Genetic Bases of the Stomata-Related Traits Revealed by a Genome-Wide Association Analysis in Rice (Oryza sativa L.)

Stomatal density (D) and size (S) are an important adaptive mechanism for abiotic stress tolerance and photosynthesis capacity in rice. However, the genetic base of rice stomata-related traits still remains unclear. We identified quantitative trait loci (QTLs) associated with D and S on abaxial and adaxial leaf surfaces using genome-wide association analysis with 451 diverse accessions in two environments. D and S showed significant differences between indica (xian) and japonica (geng) accessions and significantly negative phenotypic correlations. A total of 64 QTLs influencing eight stomata-related traits were identified using 2,936,762 high-quality single nucleotide polymorphism markers. Twelve QTLs were consistently detected for the same traits in nine chromosomal regions in both environments. In addition, 12 QTL clusters were simultaneously detected for the same stomata-related traits on abaxial and adaxial leaf surfaces in the same environment, probably explaining the genetic bases of significant correlations of the stomata-related traits. We screened 64 candidate genes for the nine consistent QTL regions using haplotype analysis. Among them, LOC_Os01g66120 for qDada1, OsSPCH2 (LOC_Os02g15760) for qDada2.1 and qDaba2.1, LOC_Os02g34320 for qSada2.2, OsFLP (LOC_Os07g43420) or LOC_Os07g43530 for qSaba7.1, and LOC_Os07g41200 for qWada7 and qWaba7 were considered as the most likely candidate genes based on functional annotations. The results systematically dissected the genetic base of stomata-related traits and provide useful information for improving rice yield potential via increasing abiotic stress tolerance and photosynthesis capacity under stressed and non-stressed conditions through deploying the favorable alleles underlying stomata-related traits by marker-assisted selection.


INTRODUCTION
Plant stomata are microscopic pores on the epidermis of leaves formed by a pair of specialized guard cells, and some are surrounded by subsidiary cells in certain monocot species such as rice (Bergmann and Sack, 2007). Stomata act as a gateway for the efficient exchange of water vapor and CO 2 between plants and the atmosphere. Stomatal size (S) and density (D) are two major factors that improve stomatal conductance and photosynthetic capacity, which are dramatically associated with the grain yield of crops under stressed and non-stressed conditions (Xu and Zhou, 2008;Franks and Beerling, 2009). Plant stomatal density plays an important role in plant responses to environmental conditions, such as elevated CO 2 concentration , salt stress (Wei et al., 2019), drought stress (Zandalinas et al., 2018), and light intensity (Shimazaki et al., 2007). Moreover, many studies have reported that water deficit leads to an increase in D, and a decrease in S (Xu and Zhou, 2008;Buckley et al., 2020), indicating that this might be partially related to enhancing the adaptation of plants to drought.
Differences in D were observed between the abaxial (lower) and adaxial (upper) leaf surfaces. More specifically, the abaxial surface possesses more stomata compared with the adaxial surface (Ishimaru et al., 2001;Laza et al., 2010), and significant positive correlations were identified for stomatal density between both sides (Ishimaru et al., 2001). Generally, xian rice varieties tend to have higher D and smaller S than geng varieties (Laza et al., 2010). Significant negative correlations were also observed between D and S in rice (Ishimaru et al., 2001;Ohsumi et al., 2007;Laza et al., 2010). However, the genetic basis of D and S on adaxial and abaxial leaf surfaces is still poorly understood. Therefore, it is important to study the genetic basis of stomata-related traits to provide useful information for the genetic improvement of rice to enhance rice photosynthetic abilities and stress tolerance.
The use of quantitative trait locus (QTL) mapping has contributed to a better understanding of the genetic basis of many agronomically important traits. Over the last few decades, many researchers have identified QTLs for stomata-related traits, including D and S, using bi-parental populations in rice (Ishimaru et al., 2001;Laza et al., 2010). Notably, several genes controlling stomatal development in Arabidopsis, such as TOO MANY MOUTHS (TMM) (Geisler et al., 2000), STOMATAL DENSITY AND DISTRIBUTION 1 (SDD1) (Von Groll et al., 2002), MITOGEN-ACTIVATED PROTEIN KINASE3 (MPK3) and MPK6 (Wang et al., 2007), STOMAGEN , which encode putative receptors, proteases, or kinases, and act primarily by modulating the number and placement of asymmetric divisions in the stomatal cell lineage, respectively. In rice, certain key transcriptional factors (TFs), such as SPEECHLESS (OsSPCH), MUTE (OsMUTE), FAMA (OsFAMA), INDUCER OF CBF EXPRESSION1 (OsICE), FOUR LIPS (OsFLP) regulate the stomatal development stage (Wu et al., 2019). Mohammed et al. (2019) reported overexpressed rice EPIDERMAL PATTERNING FACTOR1 (OsEPF1) produced transgenic plants with reduced stomatal densities, resulting in decreased leaf stomatal conductance and enhanced water use efficiency.
A genome-wide association study (GWAS) is an effective method to dissect the genetic basis of complex traits, such as grain yield, and biotic and abiotic stress tolerances (Zhao et al., 2018;Zhai et al., 2018;Liu et al., 2019). In the present study, we reported the first GWAS to identify QTLs for stomatarelated traits in rice using a panel of germplasm resources with a broad genetic diversity. A diverse panel consisting of 451 accessions from the 3000 Rice Genomes Project (3K RGP) (3K RGP, 2014) was used to detect QTLs for stomata-related traits in rice and identify associated candidate genes by GWAS using 4.8M single nucleotide polymorphisms (SNPs) generated using high-throughput re-sequencing. The results will shed new light on the genetic basis of stomata-related traits and provide superior alleles to improve rice stomata-related traits by molecular breeding.

Phenotypic Investigation
All of these genotypes were grown in Sanya (SY, 18.3 • N, 109.3 • E) during December 2016-April 2017 and Beijing (BJ, 40.2 • N, 116.2 • E) during May-October, 2017, and each accession was planted in two rows with 10 individuals in each row at a spacing of 25 cm × 17 cm for two replications. Field management followed the local standard management practices. To minimize flowering time effects on possible experiment error, batch samplings were conducted based on the heading date of each accession. At the full heading stage, eight uniform flag leaves from the main stem of eight plants (one main stem from one plant) in the middle of each plot were collected and kept in formalin-acetic-alcohol (FAA) fixative solution. Leaf impressions were made on the adaxial and abaxial leaf surfaces at the middle part of the flag leaf. The stomatal density (stomata per unit area) on the adaxial leaf surface (D ada , number/mm 2 ) and abaxial leaf surface (D aba , number/mm 2 ) was counted at 500 × magnification and repeated in five impressions from each of leaf using a scanning electron microscope (Hitachi TM3030, Tokyo, Japan). The guard cell length (L) on the adaxial leaf surface (L ada , µm) and abaxial leaf surface (L aba , µm), and guard cell width (W) on the adaxial leaf surface (W ada , µm) and abaxial leaf surface (W aba , µm), were observed and measured at 1200 × magnification and repeated for five replicates. The S values on the adaxial leaf surface (S ada , µm 2 ) and abaxial leaf surface (S aba , µm 2 ) were calculated based on the assumption that stomata are elliptical in shape, with their major axis equal to L and their minor axis equal to W, as follows Zhang et al., 2019): The average trait value of each accession was used for data analyses.

Genotyping
The 4.8M SNP genotype data of this association panel were selected from the 3K RGP (3K RGP, 2014). SNPs with a missing rate over 20% and a minor allele frequency less than 5% were removed. Heterozygous alleles were also eliminated. Finally, a total of 2,936,762 high quality SNP markers that were evenly distributed on chromosomes were used to carry out the GWAS.

QTL Mapping by GWAS
We performed a GWAS to detect SNPs that were significantly associated with all stomata-related traits using 2,936,762 SNPs and the mean trait values of the 451 accessions from each of the environments. Principal component analysis (PCA) was conducted using the efficient mixed-model association (EMMA) method in the Genome Association and Prediction Integrated Tool (GAPIT) R package (Lipka et al., 2012) to examine the population structure. The K matrix (kinship matrix) was calculated by the EMMAX software based on the Bayesian Network (BN) method using those high quality SNPs. Markertrait associations were conducted by the model of mixed linear (MLM), PCA+K, implemented in GAPIT. The critical P-value for declaring significant marker-trait association (1.0 × 10 −5 ) was calculated using GEC software based on the independent effective SNP number . We realized that the use of a single P-value threshold in GWAS could result in detection of a significant marker-trait associated SNP in one environment but not in another environment. To further examine the extent to which inconsistent associated SNP detection across the two environments (SY and BJ) actually occurred from type II errors, all significant associated SNPs identified under one condition were re-examined using the data from the other condition, with a minimum significance threshold of P < 1.0 × 10 −4.5 . The test statistics and QTL parameters associated with the QTLs were also reported as long as the QTLs reached the minimum threshold (Xu et al., 2005;Li, 2001). To estimate independently associated regions of identified QTLs, significantly trait-associated SNPs situated in one estimated LD block were defined as the same QTL. Each LD block containing the identified SNPs was evaluated with the R package "LDheatmap" (Hyung et al., 2006) according to the block definition suggested by Yano et al. (2016).

Identification of Candidate Genes
The haplotypes of all the genes located in the candidate regions for QTLs consistently detected in the two environments from the Rice Annotation Project Database 1 (Sakai et al., 2013) were conducted according to all available non-synonymous 1 https://rapdb.dna.affrc.go.jp/ SNPs located inside of these genes in the 451 rice accessions. Haplotypes containing more than 10 rice accessions were used to analyze the significant differences in phenotype. Six representative candidate genes were selected for a comprehensive analysis based on the significance of the haplotype analyses [analysis of variance (ANOVA)], their biochemically related functions, and their expression profiles. Gene expression profiles were downloaded from a rice expression profile database [RiceXPro (version 3.0)] 2 (Sato et al., 2013). Two-sided Fisher's exact tests in R were used to compare haplotype frequencies between xian and geng subpopulations.

Statistical Analysis
Differences in the mean phenotypic values among the haplotypes (containing more than 10 accessions) were evaluated using a one-way ANOVA. Duncan's multiple mean comparison test was carried out to determine the significance of any differences (5% significance level) using the agricolae package in R. Phenotypic correlation analyses of the eight stomata-traits were computed using the corrplot package in R. Variance components were evaluated using multiple-site analysis with all effects treated as random. We computed heritability across the two environments using the estimated variance components as V G /(V G + V GEI /s +V e /sr), where V G , V GEI , and V e are the variances of genotype, genotype-by-environment interaction variance, and residual error, respectively; s is the number of environments; and r is the number of replicates. All analyses were conducted with the PBTools package 3 developed by The International Rice Research Institute.

Phenotypic Variation of Stomata-Related Traits
The rice panel used in this study showed wide variations for all the investigated traits and most traits appeared to be normally distributed (Supplementary Table S1 and Supplementary Figure S1). The 451 accessions presented substantial variations for the tested stomata-related traits in Sanya (SY) and Beijing (BJ). Comparing the D and S values on the abaxial leaf surface with those on the adaxial leaf surface revealed that D aba (average 739.7 number/mm 2 and 683 number/mm 2 in SY and BJ, respectively) was significantly higher than D ada (average 548.2 number/mm 2 and 517.6 number/mm 2 in SY and BJ, respectively) (Supplementary Figure S2). In addition, S aba was significantly larger than S ada in BJ, whereas no significant difference was observed between the two sides in SY (Supplementary Figure S2). Based on estimates of variance components for the eight traits, most traits were controlled mainly by V G , whereas V GEI was the main source for L ada and S ada (Supplementary Table S2). The heritability of the eight traits ranged from 0.41 for L ada to 0.69 for D ada (Supplementary Table S2).
Based on the results from 3,010 rice accessions , 271 accessions were classified into the xian subpopulation, and 145 accessions were classified into the geng subpopulation for further phenotypic analysis in the two subpopulations. Xian varieties exhibited significantly higher D, but significantly smaller L, W, and S than geng varieties on the two sides in the two environments, except D ada and L ada in BJ ( Figure 1A).
In the two environments, the trends of pairwise phenotypic correlations were similar. Phenotypic correlations analysis demonstrated that significantly positive correlations of D and S could be observed between the abaxial and adaxial leaf surfaces. The D value had strong negative correlations with S for the abaxial and adaxial leaf surfaces. As expected, S showed significant positive correlations with its corresponding component traits (L and W) on the two sides of the leaves in the two environments ( Figure 1B).

Basic Statistics of Markers
A total of 2,936,762 high density SNPs with minor allele frequencies greater than 5% and missing rates of less than 20% were used in this study. The number of markers per chromosome varied from 168,424 on chromosome 9 to 359,991 on chromosome 1. The size of the chromosomes ranged from 22.9 Mb for chromosome 9 to 43.3 Mb for chromosome 1. The whole genome size was 372.9 Mb, with an average marker spacing of 127.8 bp, which ranged from 108.8 bp on chromosome 10 to 152.2 bp on chromosome 5 (Supplementary Table S3).

Identification of QTLs for Stomata-Related Traits by GWAS
Total of 723 and 570 significant SNPs associated with stomatarelated traits on the adaxial and abaxial surfaces were identified on 12 chromosomes, respectively, in SY and BJ ( Table 1 and Supplementary Figure S3). Adjacent significantly associated SNPs located in one estimated linkage disequilibrium (LD) block were defined as one QTL. Thus, a total of 64 QTLs were identified for the eight traits in the two environments ( Table 1), including 34 QTLs detected only in SY, 18 detected only in BJ, and 12 detected simultaneously in both environments, including five (qD ada 1, qD ada 2.1, qD ada 2.2, qD ada 4, and qD ada 6), one (qW ada 7), two (qS ada 2.2 and qS ada 4.2), one (qD aba 2.1), one (qW aba 7), and two (qS aba 4.2 and qS aba 7.1) for D ada , W ada , S ada , D aba , W aba , and S aba , respectively.

Coincidence of QTLs for Stomata-Related Traits on Adaxial and Abaxial Leaf Surfaces
Comparisons of 26 QTLs affecting stomata-related traits at the adaxial surface with 38 QTLs at the abaxial surface demonstrated that 12 QTL regions were simultaneously detected for the same stomata-related traits at both leaf surfaces in the same environment. These regions included two regions (8.87-9.05 Mb on chromosome 2 harboring qD ada 2.1 and qD aba 2.1, and 9.34-9.46 Mb on chromosome 9 harboring qD ada 9 and qD aba 9) for D, six regions (23.34-23.51 Mb on chromosome 1 harboring qS ada 1 and qS aba 1.2, 19.71-19.89 Mb on chromosome 2 harboring qS ada 2.1 and qS aba 2. 1, 20.28-20.49 Mb on chromosome 4 harboring qS ada 4.2 and qS aba 4.2, 4.32-4.60 Mb on chromosome 6 harboring qS ada 6 and qS aba 6. 1, 26.32-26.49 Mb on chromosome 7 harboring qS ada 7 and qS aba 7.2, and 23.56-23.74 Mb on chromosome 12 harboring qS ada 12 and qS aba 12.2) for S, two regions (21.06-21.35 Mb on chromosome 1 harboring qW ada 1 and qW aba 1, and 24.51-24.65 Mb on chromosome 7 harboring qW ada 7 and qW aba 7) for W, one region (10.74-10.85 Mb on chromosome 9 harboring qL ada 9, qL aba 9, qS aba 9.1, and qS ada 9) affecting L and S, and one region (20.35-20.61 Mb on chromosome 2 harboring qD ada 2.3, qD aba 2.2, qS ada 2.2, and qS aba 2.2) for D and S (Figure 2), which probably explain the genetic basis of the significant correlations of stomata-related traits at the two leaf surfaces.

Candidate Gene Identification of the Important QTLs
For the eight stomata-related traits, we conducted haplotype analyses to identify candidate genes for the 12 QTLs consistently identified in the nine chromosomal regions across the two environments ( Table 1). In addition, considering that the eight stomata-related traits were significantly different between xian and geng in the two environments, except D ada and L ada in BJ (Figure 1), we conducted haplotype analysis of the targeted genes associated with the eight stomata-related traits in whole population, and the xian and geng subpopulations. A total of 315 annotated genes located in the nine chromosomal regions were selected for haplotype analysis according to the Rice Annotation Project Database, and 64 genes were finally considered as candidate genes (Supplementary Table S4).
For D ada and D aba , haplotype analysis revealed 7, 9, 11, 5, and 5 candidate genes for qD ada 1, qD ada 2.1/qD aba 2.1, qD ada 2.2, qD ada 4, and qD ada 6, respectively (Supplementary Table S4). In the region from 38.32 to 38.63 Mb on chromosome 1 harboring qD ada 1 ( Figure 3A), 45 annotated genes were selected based on the Nipponbare reference genome IRGSP 1.0 for haplotype analysis with all non-synonymous SNPs inside of the gene coding regions in the whole, xian, and geng populations. Highly significant differences in D ada were detected among different haplotypes of the seven candidate genes (LOC_Os01g66120, LOC_Os01g66130, LOC_Os01g66180, LOC_Os01g66230, LOC_Os01g66240, LOC_Os01g66250, and LOC_Os01g66260) (Supplementary Table S4 and Figure 3). Among them, LOC_Os01g66120 is identical to OsNAC6, which is a member of the NAC transcription factor gene family in rice ( Table 2; Nakashima et al., 2007). Three haplotypes were found at OsNAC6 and the frequencies of them were significantly associated with the rice subpopulations according to Fisher's exact tests (Supplementary Table S5). Haplotype (Hap) 1 (CT) exhibited significantly higher D ada than Hap3 (CC) in the whole and xian populations in the two environments ( Figures 3B,C). In addition, 87.9% of the accessions with the high-D ada Hap1 (n = 211) and 84.6% of the accessions with the low-D ada Hap3 (n = 33) belonged to the xian subpopulation ( Figure 3D). OsNAC6 is relatively highly expressed in specific organs (leaf blade, ovary, embryo, and endosperm) according   to a publicly available rice gene expression profile database [RiceXPro (version 3.0)] ( Figure 3E). In the region from 8.87 to 9.05 Mb on chromosome 2 harboring qD ada 2.1 and qD aba 2.1 (Figure 4A), highly significant differences in the mean D ada and D aba were detected among different haplotypes of the nine candidate genes (Supplementary Table S4 and Figure 4). Among them, LOC_Os02g15760 is identical to a SPCH transcriptional factor gene OsSPCH2, which regulates the initiation of stomatal lineage, as previously reported (Liu et al., 2009;Wu et al., 2019; Table 2). Haplotype analysis revealed that Hap2 was significantly associated with larger D aba and D aba than Hap3 in the xian population (Figures 4B,C). In contrast, 63.3% of the accessions with the high-D ada Hap2 and 95.1% of the accessions with the low-D ada Hap3 belonged to the xian subpopulation ( Figure 4D). The frequency distributions of Hap3 significantly differed between the xian and geng subgroups, whereas the frequency distributions of Hap2 did not significantly differ between the xian and geng subgroups (Supplementary Table S5). In addition, 11, 5, and 5 candidate genes for qD ada 2.2, qD ada 4, and qD ada 6, respectively, were identified with significant differences in the mean D ada among the different haplotypes based on haplotype analysis (Supplementary Table S4).
For W ada and W aba , a high peak of qW ada 7 on chromosome 7 was mapped together with the peak of qW aba 7. A total of 19 annotated genes were selected from the region of 24.51-24.65 Mb (140 kb) (Figure 5A), and 11 candidate genes were identified with significant differences in the mean W ada and W aba among the different haplotypes (Supplementary Table S4 and Figure 5). Among them, LOC_Os07g41200 is identical to Grain Length on Chromosome 7 (GL7), which regulates longitudinal cell elongation, contributing to an increase in grain length and improvement of grain appearance quality ; Table 2). Hap3 showed significantly larger W ada and W aba values than Hap1 and Hap2 in the whole and xian populations (Figures 5B,C). We determined that 82.3% of the accessions with high-W aba Hap3 belonged to the geng subpopulation, whereas 65.5 and 96.4% with low-W aba Hap1 and Hap2, respectively, belonged to the xian subpopulation ( Figure 5D). The frequency distributions of Hap2 and Hap3 significantly differed between the xian and geng subgroups, whereas the frequency distribution of Hap1 did not significantly differ between the two subgroups (Supplementary Table S5). GL7 is relatively highly expressed in some specific organs except leaf blade, embryo, and endosperm according to rice gene expression profile database [RiceXPro (version 3.0)] ( Figure 5E).
Regarding S ada and S aba , haplotype analysis revealed that four, eight, and four candidate genes were identified for qS ada 2.2, qS ada 4.2/qS aba 4.2, and qS aba 7.1, respectively FIGURE 2 | Distribution of quantitative trait loci (QTLs) associated with stomata-related traits on abaxial and adaxial leaf surfaces across the two environments on 12 chromosomes according to physical distance. D ada , stomatal density on the adaxial surface; D aba , stomatal density on the abaxial surface; L ada , guard cell length on the adaxial surface; L aba , guard cell length on the abaxial surface; W ada , guard cell width on the adaxial surface; W aba , guard cell width on the abaxial surface; S ada , stomatal size on the adaxial surface; S aba , stomatal size on the abaxial surface; Envir., environment; SY, Sanya; BJ, Beijing.  (Figure 6A), which contains 33 annotated genes. Four candidate genes (LOC_Os02g34260, LOC_Os02g34270, LOC_Os02g34320, and LOC_Os02g34410) were detected with significant differences in the mean S ada among the different haplotypes (Supplementary Table S4 and Figure 6). Among them, LOC_Os02g34320 encodes a basic helix-loop-helix protein ( Table 2). Five haplotypes were detected at LOC_Os02g34320 and the frequency distributions of them significantly differed between the xian and geng subgroups (Supplementary Table S5). Hap2 exhibited a significantly higher S ada value than the other four haplotypes in the whole population (Figures 6B,C). 84.5% of the accessions with the high-S ada Hap2 belonged to the geng subpopulation, 88, 87.1, and 96% of the accessions with the low-S ada Hap1, Hap4 and Hap5, respectively, all belonged to the xian subgroup, and 87.8% of the accessions with the low-S ada Hap3 belonged to the geng subpopulation ( Figure 6D). Moreover, Hap2 showed a   significantly higher S ada value than Hap3 in the xian population in both environments ( Figure 6C). A QTL cluster (qS ada 4.2 and qS aba 4.2) affecting S ada and S aba was identified in the region of 20.28-20.49 Mb on chromosome 4 in the two environments, containing 28 annotated genes. Among them, eight candidate genes (LOC_Os04g32160, LOC_Os04g32170, LOC_Os04g32200, LOC_Os04g32210, LOC_Os04g32320, LOC_Os04g32380, LOC_Os04g32480, and LOC_Os04g32570) were identified with significant differences in S ada and S aba values among the different haplotypes (Supplementary Table S4). In addition, qS aba 7.1 was fine-mapped to the region of 26.00-26.25 Mb on chromosome 7 (Figure 7A) which harbored 39 annotated genes. Highly significant differences in S aba were detected among different haplotypes for four candidate genes (LOC_Os07g43420, LOC_Os07g43510, LOC_Os07g43530, and LOC_Os07g43580) (Supplementary Table S4 and Figure 7). Among them, LOC_Os07g43530 encodes a basic helix-loop-helix protein, while LOC_Os07g43420 is identical to OsFLP, which negatively regulates guard mother symmetric cell division in rice (Wu et al., 2019) (Table 2). Haplotype analysis of LOC_Os07g43420 revealed that Hap1 was significantly associated with a larger S aba value than Hap2 and Hap3 in the whole population ( Figure 7B). The frequency distributions of these three haplotypes significantly differed between the xian and geng subgroups (Supplementary Table S5). Moreover, 79.6% of the accessions with the high-S aba Hap1 belonged to the geng subpopulation, whereas 95.4 and 91.9% of the accessions with the low-S aba Hap2 and Hap3, respectively, belonged to the xian subpopulation ( Figure 7C). OsFLP is relatively highly expressed in some specific organs except endosperm based on rice gene expression profile database [RiceXPro (version 3.0)] ( Figure 7D). For LOC_Os07g43530, Hap2 and Hap3 of LOC_Os07g43530 exhibited a significantly larger S aba value than Hap1 in the whole population ( Figure 7E). The frequency distributions of these three haplotypes differed significantly between the rice subgroups (Supplementary Table S5). In contrast, 82.3 and 78.9% of the accessions with the low-S aba Hap2 and Hap3, respectively, belonged to the xian subpopulation, whereas 97.1% of the accessions with the high-S aba Hap1 belonged to the geng subpopulation ( Figure 7F). LOC_Os07g43530 is relatively highly expressed in some specific organs except leaf blade according to rice gene expression profile database [RiceXPro (version 3.0)] ( Figure 7G).

Characteristics of Stomata Between
Xian and Geng and Their Responses to the Environment Asian cultivated rice (Oryza sativa L.) is classified into two subspecies, namely xian and geng. There are significant differences in many morphological and physiological traits associated with evolution of xian-geng, resulting in many distinguishing features for the two subspecies of cultivated rice (Morishima and Oka, 1981). Generally, xian varieties tend to have higher D and smaller S values than geng varieties (Laza et al., 2010). In our present study, xian cultivars exhibited significantly higher D and smaller S values than geng cultivars on both the adaxial and abaxial leaf surfaces in two environments, except for D ada in BJ (Figure 1A). In addition, the frequencies of most haplotypes at six representative candidate genes for the stomatal-related traits were significantly associated with the rice subpopulations according to Fisher's exact tests (Supplementary Table S5), suggesting D and S are related to differentiation of xian and geng and dependent of population structure. Compared with geng varieties, the higher D of xian varieties was significantly associated with high photosynthetic abilities (Weng and Chen, 1988), most probably because smaller and more numerous stomata can potentially improve stomatal conductance (Zhang et al., 2019). The leaf transpiration efficiency is regulated mainly by stomatal movement, but is also affected by stomatal size and density (Nir et al., 2014). Peng et al. (1998) demonstrated that xian cultivars exhibited 25-30% lower transpiration efficiency (ratio of photosynthesis to transpiration) compared with that of the geng cultivars. Therefore, the higher photosynthetic rate of xian cultivars might be mainly attributed to its higher D and lower gas diffusion resistance.
Stomatal density is suggested to be an adaptive mechanism in plants to environmental stresses such as light intensity, water status, and CO 2 levels (Kondamudi et al., 2018). In this study, the abaxial surface possessed more stomata than the adaxial surface (Supplementary Figure S2), which was in agreement with previous studies (Ishimaru et al., 2001;Laza et al., 2010). Increases in light intensities produced greater increases in D for the abaxial surface than for the adaxial surface in Arabidopsis (Schlüter et al., 2003). In the present study, 451 rice accessions were grown in May in BJ and December in SY. During the later growth period, the light intensity in SY was stronger than that in BJ. With the increase in light supply during growth from BJ to SY, approximately 64% of the accessions (including about 40% geng and 60% xian accessions) exhibited parallel increases in D and decreases in S on both the adaxial and abaxial leaf surfaces, where the abaxial surface increased by an average of 8.3%, which was markedly higher than the average 5.9% of the adaxial surface (Supplementary Table S1). Specifically, xian accessions, on average, exceeded 10.3 and 9% of D aba and D ada in SY compared with that in BJ, whereas geng accessions showed only 4.7 and 0.5% of D aba and D ada higher values in SY than in BJ, respectively (Supplementary Table S1). The phenomenon might mainly be attributed to different adaptations to the environments, especially the light intensity on stomata between xian and geng. Xian varieties are generally adapted to tropical lowland cultivation with higher light intensity, whereas most geng varieties are adapted to temperate climates (Khush, 1997).

Comparison of QTLs Detected in This Study With Previously Reported QTLs and Cloned Genes
Of the 64 QTLs for stomata-related traits identified in this study, seven were located in the same or adjacent regions containing previously reported QTLs and cloned genes affecting stomata-related traits in rice ( Table 1). For example, qD ada 2.1 for D ada , and qD aba 2.1 for D aba in the region 8.87-9.05 Mb on chromosome 2, which was co-located with OsSPCH2 regulating the initiation of stomatal lineage in rice (Liu et al., 2009); qD ada 6 affecting D ada in the region of 26.18-26.59 Mb on chromosome 6, which harbored a previously reported QTL for D ada identified on chromosome 6 near marker R32 (Laza et al., 2010); qL aba 6 for L aba , qW aba 6 for W aba , and qS aba 6.3 for S aba , located in the region of 27.37-27.55 Mb on chromosome 6, are co-located with a QTL for S aba on chromosome 6 that overlapped marker C962 (Laza et al., 2010); qS aba 7.1 affecting S aba was mapped in the region of 26.00-26.25 Mb on chromosome 7, which harbored the previously reported qSf7b for leaf stomata frequency (Zhao et al., 2008) and OsFLP, which negatively regulates guard mother symmetric cell division in rice (Wu et al., 2019).
Many stomata-related traits are associated with photosynthesis and drought stress responsive traits (Franks and Beerling, 2009;Xu et al., 2016;Zandalinas et al., 2018). In this study, 12 and 7 QTLs for stomata-related traits were identified in the same or adjacent regions with previously reported QTLs or cloned genes affecting photosynthesis and drought stress responsive traits, respectively ( Table 1). These 12 QTLs included qD ada 2.1 for D ada and qD aba 2.1 for D aba in the region of 8.87-9.05 Mb on chromosome 2, qW aba 3 for W aba in the region of 21.25-21.45 Mb on chromosome 3, and qS aba 6.2 for S aba in the region of 23.37-23.53 Mb on chromosome 6, which mapped together with three QTLs controlling carbon isotope discrimination, which contributed to stomatal conductance (Takai et al., 2009); qS ada 4.2 for S ada and qS aba 4.2 for S aba in the region of 20.28-20.49 Mb on chromosome 4 mapped together with abaxial rolling and vein-albino leaves 4 (OsARVL4), which influences leaf morphology and the photosynthetic rate (Wang et al., 2016); qW ada 7 for W ada and qW aba 7 for W aba in the region of 24.51-24.65 Mb on chromosome 7, qL aba 3.1 for L aba in the region of 10.13-10.39 Mb on chromosome 3, and qD ada 2.2 for D ada in the region of 16.01-16.27 Mb on chromosome 2 were mapped together with qHP7b, qHP3a, and a QTL near marker RM5521 controlling leaf photosynthesis, respectively (Adachi et al., 2019); qW ada 4 for W ada in the region of 28.55-28.89 Mb on chromosome 4 mapped together with SLOW ANION CHANNEL-ASSOCIATED 1 (OsSLAC1), which is responsible for increased leaf photosynthesis caused by elevated stomatal conductance in rice (Kusumi et al., 2012); qS aba 11 for S aba in the region 18.09-18.23 Mb on chromosome 11 mapped together with a QTL cluster harboring qPn11, qGs11, and qTr11, which affect the net photosynthetic rate, stomata conductance, and the transpiration rate, respectively (Zhao et al., 2008).
Another seven QTLs included qD ada 1 for D ada in the region of 38.32-38.63 Mb on chromosome 1 and qD ada 4 for D ada in the region of 19.28-19.7 Mb on chromosome 4 mapped together with OsNAC6 (Nakashima et al., 2007) and a High-affinity potassium transporter OsHAK1 (Chen et al., 2017), which are associated with rice drought tolerance; qS ada 4.3 for S ada in the region of 34.79-34.97 Mb on chromosome 4 mapped in the adjacent region to Albino midrib 1 (AM1), which is associated with chloroplast development and drought tolerance in rice (Sheng et al., 2014); qS ada 6 for S ada , and qS aba 6.1 for S aba in the region of 4.32-4.6 Mb on chromosome 6 mapped together with a drought-responsive AP2/ERF transcription factor OsERF71, which is associated with root structure and drought resistance of rice (Lee et al., 2016); qS aba 2.3 for S aba in the region of 30.05-30.35 Mb on chromosome 2 mapped together with Grain number, plant height, and heading date2 (Ghd2), a CONSTANS-like gene that influences drought sensitivity via regulation of senescence in rice (Liu et al., 2016); and qL aba 3.2 for L aba in the region of 29.65-30.06 Mb on chromosome 3 mapped together with a Calmodulin-like gene OsCML4 for drought tolerance through scavenging of reactive oxygen species in rice (Yin et al., 2015). Allelism of the above QTLs for stomata-related traits identified in this study with previously reported QTLs or genes requires further verification via fine mapping and QTL cloning.

Most Likely Candidate Genes for the Important QTLs
GWAS and haplotype analysis of candidate genes revealed 64 candidate genes governing the 12 stable QTLs in nine chromosomal regions that were consistently detected in both environments. Based on the functional annotation of candidate genes, we speculated on the most likely candidate genes of qD ada 1, qD ada 2.1/ qD aba 2.1, qS aba 7.1, qS ada 2.2, and qW ada 7/qW aba 7.
The region 38.32-38.63 Mb on chromosome 1, harboring qD ada 1, contains seven candidate genes, including OsNAC6 (LOC_Os01g66120), a member of the NAC transcription factor gene family in rice (Nakashima et al., 2007). Transgenic rice plants overexpressing OsNAC6 showed an improved tolerance to dehydration and high-salt stresses, and exhibited increased tolerance to blast disease (Nakashima et al., 2007). Many previous studies indicated that plant D plays an important role in plant responses to salt (Wei et al., 2019) and drought stress (Zandalinas et al., 2018). Moreover, water deficit leads to an increase in D and a decrease in S (Xu and Zhou, 2008), indicating that this may be partially attributed to improving the adaptation of plants to drought. In the present study, Hap1 (CT) exhibited a significantly higher D ada value than Hap3 (CC) in the whole and xian populations in the two environments ( Figure 3C). OsNAC6 is highly expressed in the leaf blade according to a publicly available rice gene expression profile database [RiceXPro (version 3.0)] ( Figure 3D). For this reason, OsNAC6 (LOC_Os01g66120) is considered the most likely candidate gene of qD ada 1.
A QTL cluster (qD ada 2.1 and qD aba 2.1) was identified in the region of 8.87-9.05 Mb on chromosome 2, with the most significant associated SNP (rs2_8970096, P = 1.0 × 10 −7 ) being detected close (approximately 87 kb downstream) to stomatarelated gene OsSPCH2. A previous study demonstrated that OsSPCH2 regulates the initiation of stomatal lineage in rice (Liu et al., 2009;Wu et al., 2019). Compared with the wild-type plant, the stomatal pattern and morphogenesis of knock-out mutants of OsSPCH2 (c-osspch2) were identical to those of the wildtype plants; however, the D value in c-osspch2 was significantly reduced (Wu et al., 2019). In this study, four haplotypes of OsSPCH2 were found, with Hap 2(G-) showing significantly higher D ada and D aba values than those of Hap 3(AG) in the xian subpopulation ( Figure 4C). We considered that OsSPCH2 (LOC_Os02g15760) is the most likely candidate gene of qD ada 2.1 and qD aba 2.1.
On chromosome 7, a high peak of qW ada 7 was mapped together with the peak of qW aba 7, containing 11 candidate genes. Among them, Grain Length on Chromosome 7 (GL7, LOC_Os07g41200), encodes a protein homologous to Arabidopsis thaliana LONGIFOLIA proteins, which regulate longitudinal cell elongation, contributing to an increase in grain length and an improvement of grain appearance quality . Compared with Hap3, Hap1, and Hap 2 both decreased W on the adaxial and abaxial leaf surfaces in the whole and xian populations in the two environments ( Figure 5C). This suggested that GL7 probably affects guard cell width in rice. Reverse genetic approaches might be used to confirm whether GL7 is the candidate gene of qW ada 7 and qW aba 7.
For qS ada 2.2, LOC_Os02g34320, encoding a basic helix_loop_helix transcription factor, is the most likely candidate gene, with significant differences in the S ada value among different haplotypes (Figure 6C). In Arabidopsis, three basic helix_loop_helix transcription factors, SPEECHLESS (SPCH), MUTE and FAMA, are essential for the initiation of stomatal lineage, termination of meristemoid fate, and the transition to guard mother cell identity, respectively (Pillitteri et al., 2007;Wu et al., 2019). Previous work indicated that OsSPCH2 (LOC_Os02g15760), encoding a helix_loop_helix transcription factor, controls the initiation of stomatal files in rice (Wu et al., 2019). In the present study, five haplotypes of LOC_Os02g34320 were detected in the whole population, with Hap 2 being associated with a significantly larger S ada value than the other four haplotypes in both environments ( Figure 6C). In addition, Hap 2 had a significantly larger S ada value than Hap3 in the geng subpopulation, suggesting that LOC_Os02g34320 is a likely candidate gene of qS ada 2.2 that probably affects S ada in rice (Figure 6).
Of the four candidate genes for qS aba 7.1, LOC_Os07g43420 (OsFLP), which encodes a functionally redundant R2R3 MYB transcription factor, negatively regulates guard mother symmetric cell division in rice (Wu et al., 2019) and is the most likely candidate gene. In addition, LOC_Os07g43530, a basic helix_loop_helix transcription factor gene, is also a likely candidate gene of qS aba 7.1. In this study, no significant difference in the two genes was detected among different haplotypes in the xian/geng subpopulations, or only one prevalent haplotype (contained in more than 10 accessions) was observed in the xian/geng subpopulations, although there were significant differences for S aba between different haplotypes in the whole population (Figures 7B,E). This indicated that xian-geng differentiation probably led to significant phenotypic differences in S aba among different haplotypes. Thus, further study is required to validate whether LOC_Os07g43530 and OsFLP are candidate genes for qS aba 7.1 or evolution-related genes of xian and geng. The above promising candidate genes associated with the eight rice stomata-related traits are valuable resources for future functional characterization and marker-assisted breeding to improve rice grain yield under non-stressed and abiotic stressed conditions after functional tests using transformation or CRISPR/Cas9 genome editing.

Application in Rice Breeding for High Yield Potential Under Non-stressed and Abiotic Stressed Conditions
Under well-watered conditions, an increase in D could allow plants to improve stomatal conductance for gas exchange on the leaf surface, thus avoiding photosynthetic limitation by CO 2 supply. Positive correlations between stomatal conductance and D have been reported in previous studies (Xu and Zhou, 2008;Zhao et al., 2008). Schlüter et al. (2003) revealed that only an increase in stomatal density without altering any other internal leaf architecture could improve productivity under field conditions in Arabidopsis thaliana. Franks and Beerling (2009) demonstrated that changes toward increased D coupled with reduced S could maintain or improve the total stomatal area (caused by increased D) but could also provide a shorter diffusion path (caused by the smaller stomatal depth), potentially leading to improved maximum stomatal conductance and photosynthesis. In this study, we identified five accessions, including CX28, CX35, CX182, CX276, and CX343 that carry high-D haplotypes of LOC_Os01g66120 at qD ada 1 and OsSPCH2 at qD ada 2.1 and qD aba 2.1 (Supplementary Table S5), and three accessions CX230, CX341, and CX362 that carry small-S haplotypes of LOC_Os02g34320 at qS ada 2.2, LOC_Os07g43530 and OsFLP at qS aba 7.1 (Supplementary Table S5). Although these candidates are not causal genes, the SNPs in the gene sequences are suitable for marker-assisted selection (MAS) because of the high degree of linkage disequilibrium between them. Thus, after converting these linked SNPs into Kompetitive Allele Specific PCR (KASP) markers, MAS could be carried out to improve stomatal conductance, thus probably increasing photosynthesis by introgressing the high-D alleles (haplotypes) of LOC_Os01g66120 and OsSPCH2, or pyramiding the high-D alleles of the above two genes, and the small-S alleles of LOC_Os02g34320, LOC_Os07g43530 and OsFLP into lowstomatal density varieties.
Rice plants with fewer stomata are better able to maintain stomatal conductance and survive longer than control plants under drought and high temperature (40 • C) . Low-stomatal density rice have achieved equivalent or even increased grain yields, despite a reduced rate of photosynthesis in some stress conditions . Thus, one strategy for breeding drought tolerant varieties with fewer stomata could be achieved by introgressing the low-D alleles of LOC_Os01g66120 at qD ada 1 and OsSPCH2 at qD ada 2.1 and qD aba 2.1 into highstomatal density varieties by MAS. We identified three accessions, including IRIS 313-10211, CX88, and CX123, which carry low-D haplotypes of LOC_Os01g66120 and OsSPCH2 (Supplementary  Table S6). Therefore, these three accessions could be used as donor parents in rice breeding for drought tolerance by MAS.

CONCLUSION
The 451 accessions showed wide variations for the eight stomatarelated traits. GWAS identified 64 QTLs for the eight traits, 12 of which were consistently detected in nine chromosome regions in the two environments, and 12 QTL regions controlling the same stomata-related traits were simultaneously detected on adaxial and abaxial leaf surfaces in the same environment. A total of 64 candidate genes for the nine consistent QTL regions were identified and the most likely candidate genes for five loci (qD ada 1, qD ada 2.1/qD aba 2.1, qS aba 7.1, qS ada 2.2, and qW ada 7/qW aba 7) were inferred based on haplotype analysis and functional annotation. These results will enrich our knowledge of the genetic relationships among stomata-related traits in rice and will provide valuable information to improve rice photosynthesis and stress tolerance by deploying favorable alleles of stomatarelated traits by MAS.

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.