Genome-Wide Association Study for Plant Height and Grain Yield in Rice under Contrasting Moisture Regimes

Drought is one of the vitally critical environmental stresses affecting both growth and yield potential in rice. Drought resistance is a complicated quantitative trait that is regulated by numerous small effect loci and hundreds of genes controlling various morphological and physiological responses to drought. For this study, 270 rice landraces and cultivars were analyzed for their drought resistance. This was done via determination of changes in plant height and grain yield under contrasting water regimes, followed by detailed identification of the underlying genetic architecture via genome-wide association study (GWAS). We controlled population structure by setting top two eigenvectors and combining kinship matrix for GWAS in this study. Eighteen, five, and six associated loci were identified for plant height, grain yield per plant, and drought resistant coefficient, respectively. Nine known functional genes were identified, including five for plant height (OsGA2ox3, OsGH3-2, sd-1, OsGNA1, and OsSAP11/OsDOG), two for grain yield per plant (OsCYP51G3 and OsRRMh) and two for drought resistant coefficient (OsPYL2 and OsGA2ox9), implying very reliable results. A previous study reported OsGNA1 to regulate root development, but this study reports additional controlling of both plant height and root length. Moreover, OsRLK5 is a new drought resistant candidate gene discovered in this study. OsRLK5 mutants showed faster water loss rates in detached leaves. This gene plays an important role in the positive regulation of yield-related traits under drought conditions. We furthermore discovered several new loci contributing to the three investigated traits (plant height, grain yield, and drought resistance). These associated loci and candidate genes significantly improve our knowledge of the genetic control of these traits in rice. In addition, many drought resistant cultivars screened in this study can be used as parental genotypes to improve drought resistance of rice by molecular breeding.


INTRODUCTION
Rice is one of the most important staple foods and plays an important role to ensure food safety. However, rice production consumes copious amounts of fresh water. In China, 49% of all fresh water resources are used for rice production (Zhang, 2007). Due to environmental damage and anomalous climatic variation, drought has become one of the main limiting factors for global grain production (Pennisi, 2008). Throughout Asia, 20% of all rice production areas are affected by drought per year (Gowda et al., 2011). This strengthens the disparity between supply and demand of rice production. Therefore, enhancing the drought resistance of rice is vital to reduce the effect of drought on rice production.
Drought, which is a soil water deficit, can result in insufficient moisture for a plant to adequately grow and complete its life cycle. Water shortage will seriously affect rice growth and lead to a series of physiological and biochemical changes. Many studies focused on plant drought resistance. These studies reported that drought usually leads to a delayed flowering time, restrained plant growth, and ultimately loss of yield in rice (Liu et al., 2005;Yue et al., 2005;Zou et al., 2007). The reduction in yield is due to a reduced spikelet number or reduced fertility of spikelet caused by drought stress (Ekanayake et al., 1989;Liu et al., 2010). In addition, drought can cause the change of antioxidant system and osmotic regulatory system in rice, accumulating antioxidants, and osmotic regulation substances (Ouyang et al., 2010).
With the development of molecular genetic technologies, more and more studies focused on the discovery of drought resistant genes and their function. These results are useful for breeders to develop drought resistant cultivars. QTL (quantitative trait loci) mapping used to be the main method to discover drought resistance loci. For example, many drought resistance-related QTLs were identified by using a recombinant inbred line (RIL) population derived from a cross of lowland indica rice Zhenshan 97 and upland japonica rice IRAT109 in previous studies (Yue et al., 2006;Liu et al., 2010). Based on the above drought resistance-related QTLs, 17 near-isogenic lines (NILs) were constructed and phenotypic variations of these NILs were investigated under drought and normal conditions, among them, qFSR4 was fine mapped for spikelet number, flag leaf width and root volume, further analysis showed that NARROW LEAF 1 (NAL1) regulating leaf width was located in the qFSR4 interval (Qi et al., 2008;Ding et al., 2011). Some studies have shown that NARROW LEAF 1 can affect root development and drought resistance in rice (Fujita et al., 2013;Cho et al., 2014). In addition, root traits are very important traits for drought resistance in plant. Dro1, a major QTL involved in deep rooting of rice was cloned by fine-mapping, it can increase rice yield under drought conditions (Uga et al., 2011(Uga et al., , 2013. And qRT9, a QTL controlling root thickness and root length was mapped an 11.5 kb interval in upland rice . Although linkage analysis for drought resistance has made some achievements in rice, many QTLs/genes of drought resistance remain hidden. Due to a rapid development of sequencing techniques, genome-wide association studies (GWAS) became a new method for gene mining of target traits. Compared to QTL analysis, GWAS is based on natural populations, and can detect multiple alleles at the same site (Flint-Garcia et al., 2003). GWAS have widely been used in human genetic studies as well as in plants, such as Arabidopsis, maize, sorghum, and rice (Atwell et al., 2010;Huang et al., 2011;Tian et al., 2011;Riedelsheimer et al., 2012;Morris et al., 2013;Chen et al., 2014;Wen et al., 2014;Yang et al., 2014Yang et al., , 2015Wang et al., 2015;Yano et al., 2016). Identification of the allelic variation underpinning the phenotypic diversity in rice will result in enormous practical implications for rice breeding.
Genome-wide association study (GWAS) has been widely used for the genetic analysis of rice. Huang et al. (2010) constructed the first high-density haplotype map in rice via resequencing of 517 landraces and cultivars, and they identified 37 strong association loci for 14 agronomic traits via GWAS. Their study was considered as the coming of age for GWAS in rice (Clark, 2010). The genetic architecture of rice metabolism was dissected via GWAS, and five candidate genes were identified or annotated . Recently, a high-throughput phenotyping platform was developed for GWAS, achieving good results and further promoting GWAS for the dissection of genetic mechanisms for other traits (Yang et al., 2014).
In previous study, GWAS has been used to dissect drought resistance in maize (Lu et al., 2010;Wang et al., 2016;Zhang et al., 2016). However, a systematic examination of drought resistance in rice via GWAS was still called for. In this study, we conduct GWAS on a natural rice population under well watered and drought stress conditions, exploring the drought resistance of rice. The collection contains a wide range of genetic variation , highlighting its suitability for mapping studies aimed at detecting genetic variation and segregation in a diverse set of rice varieties. Using this strategy, 18, 5, and 6 loci could be associated with plant height, grain yield per plant, and drought resistant coefficient. Nine known genes were identified for three traits. The function of these genes is either directly or indirectly related to the three target traits. We also discovered new loci that contribute to these three traits and that were missed by previous studies. Among these, OsRLK5 was preliminarily confirmed to positively regulate yield-related traits under drought condition.

Plant Material and Field Experiment
The plant material consisted of 270 accessions of rice landraces and cultivars, collected from Asia, Africa, and America. This population has previously been used in GWAS for mesocotyl elongation and ratio of deep rooting (Lou et al., 2015;Wu et al., 2015). All accessions were tested under two water regimes: well watered and drought stress, in a drought resistance screening facility (Luo, 2010) at the Baihe Experimental Station of the Shanghai Agrobiological Gene Center (31 • 15 ′ N, 121 • 10 ′ E, 4 m altitude) in 2011 and 2012. The plant material was arranged by single factor randomized block design. Seeds derived from a single plant, from which the genomic DNA was extracted for sequencing, was used for field trials. Staged sowing was used according to growth durations. There were 22 hills per row with a space of 18 cm between rows and 16 cm between hills. Fertilizer application and pest control were identical to normal field management. The method for drought treatment was the same as reported in our previous study (Liu et al., 2005). Drought stress was implemented at the early booting stage and lasted for a total of 35 days. During the treatment, drip irrigation was provided to keep the plants growing well in the well-water treatment regime every day and stop watering in drought stress regime.

Measurement of Soil Water Content, Plant Traits, and Statistical Analysis
Soil moisture content was measured every 3 days for both drought stress regime and irrigational regime during drought stress. Plant height (PH) was investigated before harvest and grain yield per plant (GY) was investigated after harvest. The drought resistant coefficient (DRC) was calculated as the ratio of the grain yield per plant under drought stress regime to the grain yield per plant under water regime. Statistic analysis was conducted using SPSS software (version 19.0). Mixed model was used for ANOVA of phenotypic data. Genotype and treatment were treated as fixed factor while year was treated as a random factor. Type sums of squares were used for ANOVA for unbalanced data in Table 2 (Langsrud, 2003).

Genotype
Genomic DNA (gDNA) was extracted from a single plant and used for sequencing. A total of 270 accessions were genotyped via re-sequencing and using an Illumina HiSeq2000. Among these, genotypic data of 101 accessions were generated by Shanghai Agrobiological Gene Center, and the genotypic data of the remaining 169 accessions were provided by the Huazhong Agricultural University Yang et al., 2014). Paired-end sequence reads were mapped to a rice reference genome sequence of japonica cv. Nipponbare (MSU v6.1) using the software BWA, then used for SNP identification, following the procedures described by Wu et al. (2015).

GWAS Analysis
The genome-wide association mapping was conducted via the efficient mixed-model association (EMMA) method using the GAPIT software package in R (Lipka et al., 2012). For this study, a total of 1019,883 SNP markers were used for GWAS. 144,995 SNPs, with less than 10% missing data in this natural population, were used for kinship calculation among individuals and principal component analysis (PCA) to adjust the population structure. The genome-wide threshold was set at p = 9.81E-07, calculated via the formula: 1/total number of SNPs, which was widely used in plant GWAS studies (Wen et al., 2014;Wang et al., 2015;Bai et al., 2016). We furthermore evaluated the extent of local LD (linkage disequilibrium) for each significant SNP. The extended region, where LD between nearly SNPs and lead SNP (with the lowest p-value) decayed to r 2 = 0.6, was defined as the local LD-based QTL interval (Yano et al., 2016). To test whether the significant associated loci respond to drought stress, we performed a two-way ANOVA for each significant locus to test the significance of the interactive effect between locus and water status (normal water and drought stress). The set of loci with significant interactions with the water status were defined as QTL responsive to drought stress .

The Deletion Mutants
The CRISPR/Cas9 method was used to generate deletion mutants for the verification of candidate gene function. The mutant of OsRLK5 was obtained via the CRISPR/Cas9 method as described by Zhang et al. (2014), using the target sequence GAAAGATCCCGAAGTGGATATGG. We obtained homozygous individuals with 1 bp deletion within target sequences (GAAAGATCCCGAAGT-GATATGG) at the CDS region.
The mutant of OsGNA1 was gained via the CRISPR/Cas9 method as described by Ma et al. (2015), with the target sequence of GGGGCACGTCGAGGACGTCGTGG. Several individual homozygous or heterozygous plants were obtained. The homozygous mutants had the target sequence GGGGCACGTCGAGGAC-TCGTGG, and those homozygous and heterozygous mutants with 1 bp deletion at the CDS region were used in further experiments. The primers used for identification of mutants and haplotype analysis of OsGNA1 were listed in Table S1.

Sequence Analysis of OsGNA1
We sequenced the promoter and coding region of OsGNA1 in the natural population. Sequence alignment was conducted with ClustalX 1.83. Cluster analysis was conducted via MEGA 6.0 (Tamura et al., 2013). LD analysis was conducted using Haploview 4.2 (http://www.broad.mit.edu/mpg/haploview/; Barrett et al., 2005).

The Change of Soil Water Content during Drought Stress in 2011 and 2012
Drought treatment lasted for a total of 35 days. In this period, the absolute water content of the soil declined from 23.03 to 7.78%. Under the drip irrigation regime in 2011, values were always higher than 18.88%. The absolute water content of the soil declined from 22.98 to 8.07% under drought stress regime, while always remaining higher than 19.01% under the drip irrigation regime in 2012 ( Figure S1). The change of soil water content showed a similar trend in both years.

Phenotype Statistics and Analysis of Variance of Traits
The descriptive statistics for three traits are shown in Table 1. The coefficient of variation among different genotypes ranged from 25.70 to 28.98, 30.51 to 84.49, and 69.04 to 85.09% for PH, GY, and DRC, respectively. The ANOVA results revealed significant differences of all three traits for genotypes, treatments, and years ( Table 2). These results indicate large phenotypic variation in this natural population, possibly increasing the chance to discover candidate genes for related traits via genetic analysis.

Genome-Wide Association Study
Principal component analysis (PCA) analysis revealed two divergent groups belonging to two subspecies of cultivated rice ( Figure S2), the indica subspecies and the japonica subspecies, respectively. This is consistent with known information of the rice germplasm .
We controlled population structure by setting top two eigenvectors and combining kinship matrix for GWAS in this study. A total of 18, 5, and 6 loci were detected for plant  height, grain yield per plant, and drought resistant coefficient, respectively ( Table 3).
Genome-wide association study (GWAS) for plant height: 13 associated loci were detected under normal water conditions and 13 associated loci under drought stress while eight loci were detected under both conditions (Figure 1, Table 3). These QTLs explained a range from 3.12 to 31.58% of phenotypic variation. The associated SNP sf0123707238 on chromosome 1 provided the largest contribution to PH variation. qPH1.3 and qPH1.4 were detected under both normal water and drought stress conditions in both years, explaining the range from 8.85 to 18.96% of the phenotypic variation. Additionally, two-way ANOVA revealed that all loci that were detected for plant height showed no significant interaction between locus and water status, indicating that these QTLs of plant height were not responsive to drought stress. The peak signal of the GWAS loci often appeared in the region near but not within the known genes . We found the similar cases of five known genes involved in plant height: OsGA2ox3 (Lo et al., 2008), OsGH3-2 (Du et al., 2012), sd-1 (Spielmeyer et al., 2002), OsDOG/OsSAP11 (Giri et al., 2011;Liu et al., 2011), and OsGNA1 (Jiang et al., 2005). These genes were located within or nearby four QTLs intervals and were identified for PH via GWAS (Figure 1, Table 4). Additionally, according to genome annotation information (MSU 6.1) and KEGG pathway database, LOC_Os09g08130 and LOC_Os10g06710 related to the indole alkaloid biosynthesis pathway, located near qPH9.1 and qPH10, respectively ( Table 3).
Genome-wide association study (GWAS) for grain yield: a total of five QTLs were detected for grain yield per plant. Three QTLs (qGY8, qGY9, and qGY10) were detected under drought stress and two QTLs (qGY1 and qGY5) were detected under normal water conditions (Figure 2, Table 3). These QTLs explained the range from 15.64 to 45.61% of phenotypic variation. The associated SNP sf0506803017 on chromosome 5 contributed the highest to GY variation ( Table 3). qGY8 and qGY10 were detected for GY under drought stress in both 2011 and 2012. In addition, the two-way ANOVA revealed all loci detected for grain yield showing a significant interaction between locus and water status, indicating that these QTLs of grain yield were likely responsive to drought stress.
OsCYP51G3 (LOC_Os05g12040) was found to be located near the lead SNP sf0506803017 (Figure 2A, Table 4) and involved in steroid biosynthesis, regulating plant height and seed setting rate (Xia et al., 2015). OsRRMh (LOC_Os09g34070) was 157 kb away from the significant association signal (sf0920270490) ( Figure 2B, Table 4) and has been reported to regulate fertility rate and number of spikelet per panicle (Liu and Cai, 2013).
Genome-wide association study (GWAS) for drought resistant coefficient: six QTLs were detected for DRC, explaining a range of the phenotypic variation from 9.38 to 23.56% (Figure 3, Table 3). Since the DRC reflected relative yield between drought stress and normal water condition, these six QTLs for DRC should be involved in the response to drought stress. Two known genes are close to the association intervals (Figure 3), i.e., OsPYL2 (Tian et al., 2015) near the lead SNP sf0207287189 and OsGA2ox9 (Lo et al., 2008) near the lead SNP sf0225263251.

Candidate Genes Involved in Drought-Response Loci
The purpose of rice drought resistance breeding is to improve grain yield under drought conditions. DRC reflects the yield changes (mostly decreases) of plants suffering from drought stress. A total of 162 candidate genes were discovered based on 10 GWAS loci (QTL interval, r 2 of LD > 0.6), involved in response to drought stress as well as RNA-seq data under contrasting water status (unpublished data) (Table S2). Among them, 63 genes were up-regulated more than 2-fold and 99 genes were down regulated to less than half their expression. The expression levels of most candidate genes under drought stress conditions tended to return to normal levels after re-watering.
OsGNA1 plays an important role for plant height and root length. It is noticeable that the sf0918825447 at the promoter of OsGNA1 (LOC_Os09g31310) was detected for plant height under drought stress conditions with a p-value of 4.18E-07 in 2012. It was 5 kb away from the lead SNP (sf0918830330) for plant height under drought stress. It remains unclear whether OsGNA1 regulates plant height, although it has been   shown to regulate root development (Jiang et al., 2005). To confirm whether OsGNA1 regulates plant height, we generated a gRNA construct and introduced it into the Nipponbare to knock out the OsGNA1 gene via a CRISPR/Cas 9 strategy (Ma et al., 2015). Several homozygous or heterozygous mutant plants were obtained. The seedling height and root length of both heterozygous and homozygous mutants were significant shorter compared to non-mutant transgenic plants (Figure 4). These results demonstrate that OsGNA1 positively regulates plant height and root length. Elevated expression levels were beneficial to plant height and root length. We sequenced OsGNA1 from −891 to 748 bp in 184 accessions of the natural population and identified 38 polymorphic loci (36 SNPs and 2 indels) (Table S3). Among them, one indel (+48) and two SNPs (+77, +427) at the CDS region caused insertion/deletion of two amino acids, substitution of Pro (CCG) to Leu (CTG), and Pro (CCG) to Ser (TCG), respectively. Cluster analysis divided 184 accessions into five groups (Figure 5). LD analysis of 38 polymorphism loci is depicted in Figure 6. We performed haplotype analyses based on 38 polymorphism loci of OsGNA1 and classified the genotype into five haplotypes (Hap) ( Table S3). Hap1 and Hap5 mainly contained the indica subpopulation, while Hap2 mainly contained the japonica subpopulation. At adult stage under drought stress Hap2 had the lowest PH, and Hap5 had the highest PH and highest DRC. At seedling stage, Hap5 had the longest root length (Table S3). The one amino acid mutation caused by +427 SNP only existed in Hap5 (Table S3).
Mutation of OsRLK5 reduced the tolerance of plants to drought stress. OsRLK5 (LOC_Os02g13430) was a differentially expressed gene in the qDRC2.1 interval (Table S2), we generated a gRNA construct that was introduced into the Nipponbare to knock out OsRLK5 via CRISPR/Cas 9 technology . Rate of water-loss of detached leaves and grain yield related traits were identified in both homozygous mutants and wild-type plants. As depicted in Figure 7, the mutant plants lost water faster than the wild-type plants. Compared to wildtype plants, plant height, panicle length, panicle neck length, seed-setting rate, spikelet numbers per panicle, and grain yield per plant were significantly reduced in the mutant plants under drought stress conditions (Figure 8). These results demonstrate that OsRLK5 is responsible for the GY under drought stress.
Additionally, using genome-wide association study for eigenvectors (EigenGWAS; Chen et al., 2016), five SNPs in the loci of OsRLK5 were significantly detected for the third eigenvector, implying this gene was under selection within this population. None or few SNPs were detected for all three largest eigenvectors for OsGNA1 and eigenvector 1 and 2 for OsRLK5 (Table S4). Roughly, eigenvector 1 is responding to the population stratification of indica and japonica subspecies. Further analysis is needed to understand the selective effect on OsRLK5 gene among subgroups determined by eigenvector 3.

DISCUSSION
Drought is one of the main abiotic stresses affecting the production of most field crops. Drought stress can occur at any crop growth stage and can affect productivity to variable degrees depending on the onset time, duration, and intensity of drought. Drought stress typically suppresses plant growth at the vegetable stage and has been reported to directly reduce production at the reproductive stage (Yue et al., 2005(Yue et al., , 2006. Agronomists and breeders have made significant progress in improving the drought resistance of crops, however, the genetic and molecular basis for drought resistance in crops remains largely unclear. In this study, drought resistance of the rice germplasm was evaluated using contrasting soil moisture regimes. Combining high density SNPs and phenotypes of this germplasm, we conducted GWAS for plant height, grain yield per plant, and drought resistant coefficient. Five, two, and two known genes for the same three traits were found within or close to the associated loci, respectively. These associated loci and candidate genes could enhance the knowledge for understanding the mechanisms underlying rice drought resistance and could be used in improvement of drought resistance of rice via molecular breeding. In addition, many drought resistant germplasm accessions were obtained in this study that should be valuable in rice breeding program facing water scarcity. Compared to linkage analysis, GWAS greatly accelerated both speed and accuracy of locating QTL or candidate genes. For example, OsSPL13 has been quickly identified based on genomewide association study for the grain size and expression profiling in panicles between small-grain and large-grain varieties. Further analysis indicated that OsSPL13 positively regulates cell size in the grain hull, resulting in enhanced grain length and yield in rice (Si et al., 2016). In this study, nine reported functional genes were found to be near the significant GWAS loci ( Table 4). The function of these genes is consistent with their associated traits. The discovery of the novel drought resistant gene OsRLK5 and the new functions of OsGAN1 in this study further support GWAS as a viable approach to quickly identify new genes that are contributing to the complex traits.
Numerous factors affect plant height of rice. Gibberellin (GA) is one of the most important determinants (Monna et al., 2002;Sasaki et al., 2002;Spielmeyer et al., 2002). In this study, several candidate genes (SD-1, OsGA2ox3, and OsDOG), controlling rice plant height were involved in synthesis or signal transduction of GA. For example, SD-1 encodes a gibberellin 20-oxidase, a key enzyme in the biosynthesis of gibberellin catalyzing (GA53 → GA44 → GA19 → GA20) (Monna et al., 2002;Sasaki et al., 2002;Spielmeyer et al., 2002). Interestingly, the GWAS approach also found that SD-1 is important for plant height under drought stress condition (Figure 1), suggesting that SD-1 may be involved in drought resistance. OsGA2ox3 was also Values represent mean ± SD. Non-m, non-mutant (n = 13), He-m, heterozygous mutants (n = 4), Ho-m, homozygous mutants (n = 6). "**" represent significances at p < 0.01. involved in the GA metabolic pathway, and it has been reported that OsGA2ox3 regulates plant height and tiller number (Lo et al., 2008). Another candidate gene for rice height, OsDOG, has been reported to act as a new dynamic equilibrium regulator for the gibberellin metabolism. It has furthermore been reported that OsDOG negatively regulates cell elongation and plant height in rice (Liu et al., 2011). OsSAP11 was an A20/AN1 zinc-finger gene and contained stress-associated proteins. Over-expression of OsSAP11 can enhance the tolerance to water deficit and salt stress in transgenic Arabidopsis plants (Giri et al., 2011). OsDOG/OsSAP11 was identified within the important locus for plant height, especially as a significant signal with lower p-value under drought stress conditions.
Auxin plays an important role in plant height. In this study, OsGH3-2 was identified at the qPH1.2 interval, encoding an enzyme catalyzing the IAA conjunction to amino acids. OsGH3-2 has been reported to be involved in the regulation of plant height, root and stomatal development, as well as in the modulation of ABA level and drought resistance (Du et al., 2012). Two more candidate genes (LOC_Os09g08130 and LOC_Os10g06710) that were detected for plant height are also related to the auxin metabolism. LOC_Os09g08130 was annotated as indole-3-glycerol phosphate synthase and was located upstream of the indole alkaloid biosynthesis pathway. LOC_Os10g06710 was annotated as amidase, involved in the metabolic pathway of indole acetic acid. Indole acetic acid is an important phytohormone and plays a crucial role in cell division, differentiation, and elongation, in root development and plant height regulation (Petersson et al., 2009;Du et al., 2012;Lu et al., 2015).
In particular, a significant SNP (sf0918825447) with p-value of 4.18E-07 for plant height is located at the promoter region of OsGNA1. The associated peak signal was detected with a lower p-value under drought stress condition compared to normal water condition, implying its contribution to drought resistance. A previous study demonstrated OsGNA1 to play an important role in root development. Compared to wild type rice plants, mutant gna1 plants exhibited reduced root elongation, shorter taproot, lateral root and root hairs (Jiang et al., 2005). Similarly, the mutant of OsGNA1 exhibited shorter root length and seedling height compared to wild type plants in this study. These results suggest that OsGNA1 can regulate plant height and drought resistance. A haplotype analysis revealed that Hap5 was likely to improve drought resistance in rice because Hap5 had the highest PH at the adult stage under drought stress and had the longest maximum root  Values represent mean ± SD, n = 5. "*" and "**" represent significances at p < 0.05, p < 0.01, respectively. length at seedling stage. Moreover, Hap5 had the highest DRC (Table S3).
In this study, the candidate genes for drought responsive associated loci were discovered together with a set of RNA-Seq data for rice plants under contrasting moisture regimes (unpublished). A total of 162 genes showed significantly different expression levels in leaves between drought treated and normal water treated plants. Two functional genes, OsPYL2 and OsGA2ox9 were discovered in this study. OsGA2ox9 was involved in GA metabolic pathways, and a previous study reported it to control plant height, tiller number, and root length in rice (Lo et al., 2008). Roots play an important role for absorbing water and nutrients. Plants can benefit from longer roots for water uptake from deeper soil layers as drought stress occurs. These above results imply that OsGA2ox9 was likely affecting drought resistance via regulation of root growth in rice. OsPYL2 was reported to be an ABA receptor (Tian et al., 2015). ABA is involved in the adaptation of plants to abiotic stress. The combination of ABA and the ABA receptor can activate the ABA signal transduction and trigger the corresponding physiological reaction, including stomatal closure (Cutler et al., 2010;Kim et al., 2010). Thus, OsPYL2 is likely to regulate stomatal behavior to reduce water loss, then to enhance the drought resistance. Furthermore, our preliminary research verified the function of OsRLK5 to be positively related to the leaf water content and grain yield under drought stress (Figures 7, 8).
Two known genes (OsCYP51G3 and OsRRMh) were mapped for grain yield in this study, and the functions of these genes correlate well with the grain yield trait. OsCYP51G3 has been reported to mediate the biosynthesis of phytosterols and brassinosteroids, and to regulate seed setting rate (Xia et al., 2015). OsRRMh plays an important role during the transition from vegetable to reproductive phase in rice. Compared to the wild-type, the OsRRMh RNAi lines exhibit enlarged panicles while over-expression of OsRRMh lines lead to lower fertility rates and less number of spikelets per panicle (Liu and Cai, 2013). Values represent mean ± SD, n = 5. "*" and "**" represent significances at p < 0.05, p < 0.01, respectively.
Genome-wide association study (GWAS) based on largescale re-sequencing provides a powerful tool to discover genetic variants that can be used for crop improvement, including improvement of drought resistance. GWAS has been successfully used at a high resolution to uncover associations that involve complex traits from inbred lines of field-grown maize and rice landraces and cultivars with genetic variants Lu et al., 2010;Chen et al., 2014;Wen et al., 2014;Wu et al., 2015;Yano et al., 2016). Using high density SNPs at a wholegenome level, GWAS provides high-resolution genetic mapping that can narrow down the associated regions to candidate genes. Our study demonstrates that known genes as well as novel genes and loci can be dissected via GWAS analysis. Moreover, candidate genes discovered in this study provide useful information for further studies in drought resistance of crops. Functional identification of these candidate genes will improve our understanding of the mechanisms of crop responses to drought stress.

AUTHOR CONTRIBUTIONS
LL, HL, HM, and XM designed this study. XM, KX, SC, TL, and XL performed the experiments. XM, FF, and HW analyzed the data. XM, HL, and LL drafted the manuscript.

ACKNOWLEDGMENTS
We thank the National Key Laboratory of Crop Genetic Improvement at the Huazhong Agricultural University for providing the seeds and genome sequence data of the rice germplasm collection.