The increasing world population has spawned tremendous concerns about impending global food shortages. Nevertheless, the per capita increase in global food production has been greater than the increase in population since 1960, due to enormous increases in cereal crop production (Trewavas, 2001). This dramatic increase in crop yield, known as the Green Revolution, largely resulted from the development of genetically improved high-yielding crop varieties (Hedden, 2003). Among the genetic approaches used to improve yield potential, “ideal-type breeding” was used to modify plant architecture (morphology); for example, the improved rice cultivar IR24 has more tillers, reduced plant height and erect leaves. Rice breeders have harnessed natural variation in cultivars to modify rice plant architecture. For example, incorporating the recessive gene sd1 from the Chinese variety Dee-geowoo-gen, reduced plant height (Suh and Heu, 1978), and incorporating the photoperiod-insensitive gene se1 increased rice adaptability to a wide range of latitudes (Sano, 1992; Tamura et al., 1998).

Although these genetic variants played an important role in the marked increase in crop production, the demand for rice varieties with higher yield potential and greater yield stability has continued to rise due to increasing human population and changing global climate (Khush, 2005). To further enhance the yield potential of rice over that of existing high-yielding cultivars, crop scientists proposed an ideal plant architecture that includes fewer tillers with few or no unproductive tillers, more grains per panicle, and thick and sturdy stems (Khush, 1995; Li et al., 2012). Some of the genes and QTLs regulating architecture-related traits have been reported. For example, OsSPL14, which is regulated by OsmiR156, affects architecture-related traits; a point mutation in OsSPL14 perturbs OsmiR156-mediated regulation of OsSPL14 and generates an ideal-type rice plant with fewer tillers, increased lodging resistance and enhanced grain yield (Jiao et al., 2010). Marathi et al. (2012) reported 19 novel QTLs and 15 QTL hotspots affecting plant architecture, including yield-related traits, using RILs developed from Pusa1266, a semi-dwarf, high-yielding cultivar. Further work identified several other QTLs regulating plant shape, including leaf length, leaf width, tiller number and panicle length (Farooq et al., 2010; Liu et al., 2010; Wang et al., 2010; ²⁰¹¹). Modifying plant architecture for architecture-based breeding will require identification of QTLs and their causal genes in different genetic backgrounds.

QTL mapping and cloning using DNA markers are useful tools to unlock the genetic basis of complex phenotypic variation. Moreover, the application of molecular markers developed from QTL analysis enhances breeding efficiency by enabling marker-assisted selection for particular agronomic traits. However, traditional approaches for QTL analysis require substantial time and labor because of the low efficiency of polymorphic marker development for high-resolution map construction and the difficulties of candidate gene identification. Whole-genome sequencing, using next-generation technology can overcome these limitations of QTL analysis. A sequencing-based geno-typing method uses single nucleotide polymorphisms (SNPs) identified by whole-genome sequencing of a mapping population (Huang et al., 2009). This approach could substantially reduce the amount of time and effort required for genotyping of a mapping population. Recent work also proposed another approach, named ‘QTL-seq’, for rapid identification of QTLs for a given phenotype. QTL-seq works by sequencing DNAs from two populations, each composed of 20-50 individuals showing extreme opposite trait values in segregating progenies (Takagi et al., 2013). Although these approaches improve the time- and cost-effectiveness of QTL mapping, approaches to increase the efficiency of candidate gene identification have not yet been proposed.

In this study, we conducted QTL mapping for six traits related to plant architecture in rice and identified 11 main-effect and 16 minor-effect QTLs. By analysis of whole-genome sequencing data and meta-expression profiles, we identified 15 strong candidate genes for the 11 main-effect QTLs. The plant-shape QTLs and their candidate genes provide useful genetic resources for breeding of high-yielding rice cultivars.

Plant materials and growth conditions

To identify QTLs influencing plant shape in rice, we crossed a japonica cultivar ‘SNU-SG1’ with the tongil-type cultivar Milyang23 with a genetic make-up similar to indica (> 90%) (Yoo et al., 2007). Among 338 F₇ RILs developed from the F₂ population by single seed descent, we chose 178 F₇ RILs to conduct QTL mapping; these included 40 lines showing extreme phenotypic values for each trait. F_7:8 heterogeneous inbred family-near isogenic lines (H-NILs) derived from F₇ RILs were used to test the allele effects of candidate gene of plant shape-related QTLs (Tuinstra et al., 1997). In 2011, the RIL population and two parents were evaluated at the Seoul National University Experimental Farm in Suwon (37° N latitude), Korea. The field experiment was conducted in a randomized complete block design with two replications in two consecutive years (2010 and 2011) and field management followed normal rice agronomic practices. Seeds were sown in the seed bed of the greenhouse on April 20 and seedlings were transplanted into the paddy field on May 19 at a spacing of 0.3 m by 0.15 m with one seedling per hill. Compound fertilizer (21N-17P-17K) was applied in the same amount (24 kg/10a) and applied three times: before transplanting and at 7 and 10 weeks after transplanting.

Evaluation of agronomic traits

For the phenotypic evaluation of the selected 178 F₇ RILs, we measured these agronomic traits; plant height, tiller number, panicle diameter, panicle length, flag leaf length and flag leaf width. All of the traits were evaluated when panicles were fully emerged. Panicle diameter (0.5 cm below the panicle base) was measured by a digital caliper. Plant height was measured from the soil surface to the apex of the main stem. On the main stem, flag leaf length was measured from the junction of leaf blade and sheath to the leaf apex, and flag leaf width was measured at the widest location of the leaf. Tiller number was scored when the grains had matured fully. Of all the flowered tillers of an individual plant, the longest panicle was measured for panicle length.

Molecular markers for genotyping

Genomic DNAs were extracted from leaf tissues by a CTAB method (Murray and Thompson 1980). Of 131 DNA markers used for genotyping, 68 markers are RM-series simple sequence repeat (SSR) markers that were designed by Temnykh et al. (2001), and 35 markers are S-series sequence-tagged site (STS) markers that were designed based on the sequence differences between japonica and indica, using information from the Crop Molecular Breeding Lab, Seoul National University (Chin et al., 2007). The other 28 markers are NGS-based insertion/deletion markers designed by sequence comparison between the two parent genotypes. To develop these markers, we investigated 7 to 11 base pairs of unique insertion/deletion (InDel) regions by comparing whole-genome sequencing data of the two parents, and designed primers to get PCR product (< 150 bp) containing the InDels. The initial genetic map was constructed with SSR and STS markers, and the next-generation based InDel markers were used to further map unlinked regions in the linkage map.

Linkage map construction and data analysis

A molecular linkage map was constructed using Mapmaker 3.0 (Lander et al., 1987; Lincoln et al., 1993). Distance between markers is presented in centiMorgans (cM) using the Kosambi map function (Kosambi, 1994), and the order of markers was established using three point linkage analyses. QTL analysis was conducted with composite interval mapping implemented in Qgene 4.3.10 software (Joehanes and Nelson, 2008). Significance thresholds for QTL detection were determined with 1000 permutations, and the threshold of logarithm of odds (LOD) for each QTL ranged from 2.24 to 3.03 with 95% confidence. The QTLs explaining more than 10% of the phenotypic effect were defined as main-effect QTLs as described previously (Wang et al., 2010).

High-throughput whole genome re-sequencing

High-throughput whole-genome sequencing was performed on the two parents, SNU-SG1 and Milyang23. The Illumina Genome Analyzer II was used to generate short reads and bases were called with the Sequence Analysis Software of Pipeline version 1.4 of the Genome Analyzer (Illumina, Inc., USA). Genomic DNAs were extracted using the QIAquick PCR Purification Kit (Qiagen) from leaf tissues. Full sequencing was processed in the National Instrumentation Center for Environmental Management (NICEM) in Seoul National University, Seoul, Korea. The sequencing yielded 26X and 31X fold coverage, as well as 9 Gbp and 11 Gbp in read length, for SNU-SG1 and Milyang23, respectively (Supplementary Tables S1 and S2). Then, the 75 bp paired-end reads of the two parents were mapped to the japonica cv. Nipponbare reference genome, which was used as the reference sequence (International Rice Genome Sequencing Program, 2005). Finally, the low-quality bases (Q score in scale < 20), and sites with conflicting genotypes among reads were excluded, and only the reads aligned to unique locations in the reference genome were used for sequence construction. The total number of SNPs was 3,202,922, which include 163,933 from SNU-SG1 and 3,038,989 from Milyang23.

Identification of candidate genes for the QTLs

Of the 27 QTLs, 11 main-effect QTLs with high LOD and R² value were used for candidate gene analysis. We identified candidate genes harboring mutated sequences such as SNPs and InDels by the following process. First, a reference sequence (Nipponbare) corresponding to each main-effect QTL was transferred from the GRAMENE database (http://www.gramene.org/) into text in Microsoft Word. The genomic region associated with each QTL was selected by confidence interval (CI) that was determined by 1.5-LOD support interval (Dupuis and Siegmund, 1999). 1.5-LOD support intervals give a 95% probability of QTL location on the map and were calculated as the longest contiguous interval in which the LOD value was within 1.5-fold of its maximum. Coding sequence regions derived from SNU-SG1 and Milyang23 were interrogated to detect mutated sequences using the AutoHotkey macro program (http://www.autohotkey.com/). This process has two steps. First, the polymorphic sequences (< 15 bp) obtained from whole genome sequencing data of the two parents were copied and pasted into Microsoft Excel for sequence blast using Auto-Hotkey source codes 1 and 2 (Supplementary Fig. S1). Second, the sequences were blasted against the Nipponbare genome sequence to locate the genomic regions of polymorphic sequences in Microsoft Word using the AutoHotkey source code 3 (Supplementary Fig. S1). Among the genes harboring mutated sequences, only the genes for which the mutation caused amino acid changes were selected as candidate genes. Arabidopsis orthologous gene information for the candidate QTL genes was obtained from The Arabidopsis Information Resource (TAIR, http://www.arabidopsis.org/).

Heatmap analysis or microarray data analysis

We used the information of 1150 Affymetrix array data (platform accession number is GPL2025) from the NCBI gene expression omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) to analyze the mRNA expression pattern of candidate genes in Supplementary Table S3 (Cao et al., 2012; Jung et al., 2011). From the Affymetrix anatomical meta-expression database, we generated the log₂ average intensity data of 16 tissues/organs; callus, coleoptile, collar, leaf, flag leaf, root, shoot, stem, internode, shoot apical meristem (SAM), panicle, anther, stigma/ovary, palea/lemma, developing seed, and embryo. Among them, the data for stem and palea/lemma were only available for indica, the data for coleoptile, collar, internode, stigma/ovary and embryo were only for japonica, and the other 9 tissues/organs have the data for both. Then, the TIGR Multiple Experiment Viewer (MEV) software was used to prepare the Heatmap using the anatomical expression data (Supplementary Fig. S2; Table S3) (Eisen et al., 1998).

RNA extraction, reverse transcription (RT) and quantitative real-time PCR (qPCR)

Flag leaf tissues of SNU-SG1 and Milyang23 grown in paddy field were homogenized in liquid nitrogen. Total RNA extraction and RT-qPCR were performed as previously described (Yoo et al., 2009). Briefly, the RT products equivalent to 50 ng of total RNA and GoTaq qPCR Master Mix (Promega) were used for qPCR reactions using the Light Cycler 480 (Roche) and Roche Optical System software was to calculate threshold cycle (DCT) value. The rice Ubiquitin gene (Os01g22490) was used as an internal control. The primers used for qPCR are listed in Supplementary Table S4.

Phenotypic variation in F₇ RILs from the Milyang23/SNUSG1 cross

For QTL analysis of the morphological traits that determine plant architecture, we used 178 F₇ RILs from the cross of Milyang23 (indica) and SNU-SG1 (japonica) (Yoo et al., 2007). We analyzed six shape-related traits that critically affect lodging and yield: plant height, tiller number, panicle diameter, panicle length, flag leaf length, and flag leaf width. We measured phenotypic variation in the two parental lines and the RILs (Fig. 1). Most of the traits showed approximately normal distribution with bidirectional transgressive segregation. We also observed a large amount of variation between the two parental cultivars for the traits, except for plant height and tiller number. Supplementary Table S5 shows descriptive statistics of the six agronomic traits. We analyzed the correlation among the traits in the F₇ RIL population (Table 1) and observed a positive correlation between flag leaf length and other traits, except for flag leaf width. We also found relatively high positive correlations between plant height and panicle length (0.64), and between panicle diameter and flag leaf width (0.53).

Linkage map construction

A coarse-scale linkage map was initially constructed over the whole rice genome using 178 F₇ RILs and 68 SSR and 35 STS markers. To construct linkage map for the unlinked regions in the initial map, we designed 28 next-generation based InDel markers using the whole-genome sequencing data. We used Mapmaker 3.0 to construct the genetic linkage map, which consisted of 12 linkage groups and covered 2,451 cM (Fig. 2). The average distance between adjacent markers was 14.5 cM, which is less than 20 cM, the minimal compulsory level for QTL mapping (Lander and Botstein, 1989).

QTL analysis for the six agronomic traits

To identify QTLs, we used a LOD threshold value and phenotypic effect (R²) (see “Materials and Methods”). As a result, we identified a total of 27 QTLs for the six traits, with 7.4 to 33.5% phenotypic effect (R²) and 2.9 to 15.8 LOD values, across all the chromosomes except chromosomes 7 and 12 (Table 2). Thus, chromosomes 7 and 12 were excluded from Fig. 2, because they did not carry any significant QTL above the threshold. A QTL with highest LOD (15.8) and phenotypic effect (33.5%) was detected for plant height on the long arm of chromosome 1 (Table 2; Fig. 2). Overall, we detected 3 to 6 QTLs for each agronomic trait. Of the 27 QTLs, 11 QTLs explained more than 10% of the phenotypic effect and were defined as the main-effect QTLs, which include Fll4 for flag leaf length, Flw5 and Flw9 for flag leaf width, Pl1 and Pl9 for panicle length, Ph1, Ph3 and Ph9 for plant height, and Tln2, Tln11-1 and Tln11-2 for tiller number (Table 2). Especially, three QTLs (Pl9, Ph1 and Tln2) explained more than 15% of the phenotypic effect, suggesting that these genetic materials provide a valuable resource for the identification of major genes affecting plant architecture.

QTL hotspots, genomic locations associated with multiple traits, may harbor critical regulators of various agronomic traits (Neto et al., 2012). Two QTL hotspots were identified on chromosomes 4 and 9 (grey blocks in Fig. 2). One, in a window of around 8 cM on the long arm of chromosome 9, harbors three main-effect QTLs (Flw9, Ph9 and Pl9 for flag leaf width, panicle height, panicle length, respectively) and two minor-effect QTLs (Pd9 and Fll9 for panicle diameter and flag leaf length) (Fig. 2). The other, in a window of around 5 cM on chromosome 4, carries three QTLs, including one main-effect QTL (Fll4 for flag leaf length) and two minor-effect QTLs (Flw4 for flag leaf width and Tln4 tiller number) (Fig. 2).

To compare the physical locations of the QTLs identified in this study and previously reported QTLs, we searched the literature for cases in which mapping populations derived from the crosses of indica and japonica cultivars were used for QTL identification. By using flanking DNA markers, we compared the physical locations between the QTLs in the current study and the previously published QTLs. Among 27 QTLs identified in this study, 13 QTLs were in accord with the chromosomal regions carrying the published QTLs (Table 3). The remaining 14 QTLs did not overlap with the published QTLs, indicating that these are potentially novel QTLs for the traits. Notably, the QTLs for panicle diameter are all newly discovered. Among the new 14 QTLs, 6 QTLs are the main-effect QTLs with above 10% of phenotypic effect, including Fll4, Flw5, Pl9, Ph3, Tln11-1 and Tln11-2 (Table 2).

Sequencing analysis for candidate gene identification

We carried out whole-genome sequencing of the two parentacultivars, SNU-SG1 and Milyang23, to obtain a cultivar-specific nucleotide polymorphism data set (see “Materials and Methods”). The polymorphic nucleotide sequences, including In-Dels, were obtained by comparing the whole genome sequences of SNU-SG1 and Milyang23 with the reference Nipponbare sequence using Pipeline version 1.4 (see “Materials and Methods”). Of the total 363,845 SNPs identified by comparing with the reference sequence, 330,317 were identified from Milyang23, and 33,528 from SNU-SG1. Sequence variation between indica- and japonica-type subspecies contributed to the large number of SNPs in Milyang23. For identification of candidate genes associated with the QTLs, we analyzed all the SNPs (including InDels) derived from SNU-SG1, but only In-Dels capable of impairing protein function were considered for Milyang23 to filter out non-causal subspecies sequence variation between indica and japonica genomes.

To perform candidate gene analysis for 11 main-effect QTLs, we analyzed the genomic region that was selected based on confidence interval for each QTL peak (Table 2). Confidence intervals were determined by 1.5-LOD support interval, which provides a QTL coverage probability of approximately 95% (Dupuis and Siegmund, 1999). The confidence intervals of the QTLs detected in the integrated map using composite interval mapping are given in Table 2. We determined the specific locations of polymorphic sequences in the individual genes by BLAST to the Nipponbare reference sequence. By these coding sequence analyses, we identified 54 genes with 656 nucleotide sequence polymorphisms in their open reading frames within the 11 main-effect QTL regions (data not shown). Among them, only the 308 SNPs that caused amino acid variation were considered as candidate genes. Thus, by analysis of the whole-genome sequencing data, we isolated 54 candidate genes carrying at least one amino acid polymorphism for the six traits (Supplementary Table S3).

Strong candidate genes for the main-effect QTLs

Among the 54 candidate genes for 11 main-effect QTLs, we further isolated strong candidate genes by analyzing their amino acid sequences and meta-expression profiles. First, by protein coding sequence analysis, we selected 11 genes as strong candidates, because they harbor potential causal SNPs leading to a premature stop codon or frameshift in their amino acid sequences (Table 4). These 11 strong candidate genes included candidates for all agronomic traits examined, except for panicle diameter, with one to seven genes for each trait. These included single genes for panicle length, plant height and flag leaf width, two genes for flag leaf length, and seven genes for tiller number (Table 4). One candidate gene (Os01g68510), encoding peptide transporter 2 (PTR2), was identified for two QTLs related to panicle length and plant height on the distal part of chromosome 1 (Table 4). The largest effect QTL Ph1 for plant height was also detected with 33.5% R² value in this location on chromosome 1 (Tables 2 and 3). Meanwhile, QTL effect can be also caused by differences in gene expression between alleles due to mutations in regulatory regions including promoters and introns (Meuwissen, 2010). Thus, we next analyzed meta-expression profiles of the 54 candidate genes using Affymetrix array data (Supplementary Fig. S2), which provide the indica or japonica-specific gene expression (Cao et al., 2012; Jung et al., 2011). As we used the indica (Milyang23) and japonica (SNU-SG1) cultivars to create the QTL mapping population, the genes expressed differentially in the indica and japonica rice cultivars would be considered as strong candidates for the corresponding traits. By Heatmap analysis, we found that seven genes are differentially expressed in at least one tissue type with more than two-fold differences between indica and japonica (Table 5). Interestingly, all the probes corresponding to the seven genes were expressed preferentially in the flag leaves of japonica rice. To test whether these genes are differentially expressed in SNU-SG1 and Milyang23, we performed RT-qPCR analysis using flag leaves. Of the seven genes, five genes were preferentially expressed in SNU-SG1 with more than two-fold difference compared to Milyang23 (Supplementary Fig. S3), in accordance with meta-expression profile. However, we failed to amplify the other two genes (Os02g40190 and Os02g40240), possibly due to very low transcript abundance. This result strongly supports the idea that at least five genes identified by meta-expression profile analysis are strong candidate genes for the corresponding QTLs. Notably, three genes (Os02g40190, Os02g40200 and Os02g40240) for Tln2 were identified as strong candidates by both coding sequence and meta-expression analyses, and encode known receptor-like protein kinase precursors (Tables 4 and 5). Thus, taken together our data provide a total of 15 strong candidate genes by two different analyses.

We further investigated the orthologs of these candidate genes in Arabidopsis thaliana, using Gramene (http://www.gramene.org/) and TAIR (http://www.arabidopsis.org/) databases. We found six Arabidopsis orthologs in which four genes have been functionally characterized (Table 4). In Pl1/Ph1 QTL, At1g69870 encodes a low affinity nitrate transporter (NRT1.7), which functions in nitrate remobilization in leaf (Fan et al., 2009). In Tln2 QTL, At1g69440 encodes ARGONAUTE 7, a member of the ARGONAUTE family, which functions in the regulation of developmental timing (Adenot et al., 2006). In Fll4 QTL, At2g 18800 encodes Xyloglucan endotransglycosylase/hydrolases (XTHs) which functions in carbohydrate metabolism (Paul and Janet, 1999). In Tln1-1 QTL, At2g07050 encodes cycloartenol synthase 11 (Corey et al., 1993) which functions in the biosynthesis of brassinosteroids (Asami and Yoshida, 1999). This identification of promising candidate genes provides useful genetic resources for the further isolation of the causal gene associated with each QTL.

Validation of QTL regions using H-NILs

To further validate the QTL regions associated with 15 strong candidate genes, we developed H-NILs that segregate for the phenotypes of the corresponding traits (Tuinstra et al., 1997). The genotypic and phenotypic analyses of the progeny H-NILs allow identification of the candidate genes underlying QTLs (Routaboul et al., 2012). In 178 F₇ RILs, we selected heterogeneous lines for at least one of the flanking marker loci in the QTL regions (Fig. 3). By phenotypic analysis using 40 F_7:8 individual plants derived from each heterogeneous line, we found several H-NILs for six main-effect QTLs that harbor strong candidate genes (Supplementary Table S6). Most of the H-NILs, including F₇ RIL #72-derived ones which showed highest variation in tiller number, displayed transgressive segregation (Fig. 3, Supplementary S4). These support that the target genes are localized within the corresponding QTL regions, enhancing the fidelity of strong candidate genes. Further linkage analysis of the markers and traits in the H-NIL mapping populations may lead to isolation of the genes underlying the QTLs.

Next-generation sequencing data enhance the efficiency of QTL analysis

Agronomically important crop traits are generally controlled by multiple genes (QTLs) that produce continuous variation in phenotype. Conventional QTL analysis remains time-consuming and labor-intensive, mainly because it requires development of polymorphic markers for linkage analysis. To overcome this limitation, several research groups have recently introduced next-generations sequencing technology as a new strategy for genetic mapping. For example, high-throughput genotyping of recombinant populations can be performed by whole-genome sequencing of individual RILs (Huang et al., 2009). This sequence-based genotyping has proven faster for genetic map construction and more accurate for determination of recombination breakpoints than marker-based genotyping methods (Wang et al., 2011). Another group performed whole genome re-sequencing of 50 rice cultivars and investigated genome-wide variation in high-quality candidate SNPs obtained from re-sequencing data for each cultivar (Xu et al., 2012). In this study, we modified conventional mapping methods by whole-genome sequencing of the two parental cultivars to improve QTL mapping efficiency. To develop polymorphic DNA markers for linkage analysis, we generated InDel markers using whole-genome sequencing data; these markers were designed based on nucleotide sequence polymorphisms between the two parents. The InDel markers showed much higher efficiency (74.5%) in polymorphic marker development compared to conventional SSR (18.1%) and STS (21.1%) markers, which were designed from sequence polymorphisms between japonica and indica genomes (Fig. 4). These results indicate that application of a whole-genome sequencing approach to the two parental cultivars substantially enhanced the efficiency of QTL mapping with relatively lower cost for whole-genome sequencing compared to high-throughput approaches.

High-resolution mapping to identify candidate genes for the detected QTLs is another time-consuming and laborious step in conventional QTL analysis. Thus, instead of performing fine mapping to narrow down the QTL region to less than 100 kb, we directly analyzed the genomic regions for the QTLs that were selected based on a 1.5-LOD confidence interval (Table 2), and identified candidate genes for each QTL of interest. Our sequence analysis found an average of five candidate genes harboring 23 causal SNPs in the confidence interval-based QTL regions by comparing genomic sequences of SNU-SG1 and Milyang23 with the reference Nipponbare. In conventional QTL analysis, fine mapping has to be conducted to narrow down the QTL region to 50 to 100 kb through development of large numbers of polymorphic DNA markers in each QTL region as well as development of mapping populations such as heterogeneous inbred family-near isogenic lines (Tisne et al., 2008) and backcross inbred lines (Li et al., 2011). Whole-genome sequencing-based QTL analysis enabled us to isolate candidate genes for all the QTLs identified from initial mapping population at the same time, while conventional QTL approaches need to develop a high resolution-mapping population for each QTL to find strong candidate gene(s). Furthermore, the combination of whole-genome sequencing-based QTL and meta-expression profile analyses enable us to isolate candidate genes for the previously reported QTLs by sequencing of the two parents that were used to construct mapping populations. Therefore, we propose that this advanced QTL analysis approach, sequencing two parent cultivars, is an efficient method for candidate gene identification, as it skips the fine-mapping procedure.

Comparison of the identified QTLs with previously reported QTLs for agronomic traits

Some QTLs were identified here and in previous studies of plant architecture that examined traits such as tiller number and panicle diameter. For example, the four QTLs, Tln2, Tln4, Tln6 and Tln9, which influence tiller number, shared a similar genomic location with QTLs identified in previous studies (Jiang et al., 2012; Liu et al., 2012; Qu et al., 2012; Rahman et al., 2008; Yan et al., 1998) (Table 3). This consistency in the QTL regions associated with the tiller number among independent studies indicates that the locations of alleles for the trait are conserved across different genetic and environmental backgrounds. By contrast, the other two QTLs, Tln11-1 and Tln11-2, have not been previously identified. All the QTLs we identified for panicle diameter (Pd3, Pd6, Pd8, Pd9 and Pd11) are novel, with no correspondence to previously reported QTLs for this trait (Table 2), indicating that these QTLs may be a potentially new set of alleles specific to this genetic background.

Analysing candidate genes and inferring their functions from Arabidopsis homologs

By the whole-genome sequencing data analysis, we found 54 candidate genes in all the detected QTL regions related to six agronomic traits (Supplementary Table S3). By further sequence and meta-expression profile analyses, We designated 15 out of 54 genes as strong candidates for 11 QTLs. Of the 15 strong candidate genes, 11 genes identified by amino acid sequence analysis carry polymorphic SNPs that lead to a frameshift or premature stop codon in the coding sequence (Table 4). For example, in the Pl1 and Ph1 QTL region (Fig. 5A), five genes containing polymorphic SNPs in their coding sequences were initially found by sequence analysis (Supplementary Table S3). Of these SNPs, a 1-bp deletion in the coding region of PTR2 occurred in the Milyang23 allele, leading to a premature stop codon (Fig. 5C), whereas simple SNPs, which do not cause premature stop codon, were found for the other four genes. These suggest that PTR2 is the strongest candidate gene for the Pl1 and Ph1 QTL (Supplementary Table S3). In the same way, we identified 10 strong candidate genes for the other five agronomic traits (Figs. 4B and 4D; Supplementary Table S3).

In addition to mutations in coding sequences, nucleotide mutations in the promoter, 5′-UTR, 3′-UTR and intron of many eukaryotic genes affect morphological traits by regulating gene expression at transcriptional or post-transcriptional levels (Lu et al., 2008; Rose et al., 2008; Yu et al., 2007). For example, TAC1, which affects tiller angle, is highly expressed in the indica rice line ‘IR24’ which displays a relatively spread-out plant architecture, but the expression of TAC1 is low in the introgressed line, IL55, derived from the japonica rice ‘Asominori’ which displays extremely erect tillers (Yu et al., 2007). This differential expression results from a nucleotide change in an intron, which affects 3′-UTR splicing. Recently, Jung et al. (2013) analyzed publicly available Affymetrix array data and identified 490 genes preferentially expressed in japonica rice and 104 genes in indica rice. Some of the genes that are preferentially expressed in each subspecies involve development of morphological traits such as plant height, panicle, culm and root (Jung et al., 2013). Thus, using the Affymetrix array data, we analyzed the expression patterns of 54 candidate genes in the indica and japonica rice samples of various tissue/organs and found seven genes that are differentially expressed between japonica and indica; this japonica-preferred expression was confirmed in SNU-SG1 by RT-qPCR analysis (Supplementary Fig. S3). Notably, three of the differentially expressed genes (Os02g40190, Os02g40200 and Os02g40240) were identified by both the sequence and expression analyses (Tables 5 and Supplementary Fig. S3). Together with 11 candidate genes that were identified by sequence analysis, an additional four genes from the meta-expression analysis between japonica and indica cultivars would be strong candidates for the causal genes for the corresponding QTLs.

The functions of the uncharacterized genes in rice can generally be predicted from the functions of orthologous genes in Arabidopsis, as proteins in these species generally have conserved molecular and biochemical functions. For example, the functions of CO and FT, flowering promoters in photoperiod pathway, are well conserved between Arabidopsis and rice and most likely in other species as well (Hayama and Coupland, 2004). Among the 15 strong candidate genes proposed in this study, the functions of four Arabidopsis orthologs have been reported (Corey et al., 1993; Fan et al., 2009; Maris et al., 2009) (Table 4). The Arabidopsis nitrate transporter NRT1.7 (Fan et al., 2009), an ortholog of PTR2, plays an important role in source-to-sink remobilization of nitrate. In addition, a rice gene, Short Panicle 1 (SP1) encoding a putative peptide transporter (PTR) has been reported to regulate inflorescence branch elongation; the knock-down sp1 mutant is defective in rice panicle elongation, producing short-panicle phenotypes (Li et al., 2009). These also indicate that PTR2 is a promising candidate gene related to the two QTLs detected for panicle length and plant height in this study. Another candidate gene Os11g18194, located in the region of Tln11-1 (Fig. 5B), is a rice ortholog of Arabidopsis CAS1, encoding a cycloartenol synthase, which is involved in the biosynthesis of brassinosteroids (Corey et al., 1993). Mutations in the genes involved in the BR biosynthesis pathway modify tiller number (Fujii et al., 1991), leaf size (Morinaka et al., 2006) and yield per plant (Choe et al., 2001). Moreover, increased levels of brassinosteroids cause a significant effect, increasing tiller number and seed weight (Wu et al., 2008). These suggest that OsCAS1 is a promising candidate gene for Tln11-1 related to tiller number. Lastly, another candidate gene Os04g51450, located in the region of Fll4, has an orthologous gene At2g18800 (XTH21), encoding a member of the xyloglucan endotransglucosylase/hydrolases (XTH), which is involved in cell loosing process at the preliminary stage of cell elongation in shoot and root meristems (Liu et al., 2007). Genes related to cell elongation, such as ANGUSTIFOLIA in Arabidopsis, play a role in leaf longitudinal development (Tomohiko et al., 1996). This supports the possible involvement of Os04g51450 in Fll4 QTL regulating flag leaf length.

In this study, we developed an improved approach for DNA marker development and candidate gene identification by whole-genome sequencing of the two parent cultivars of our RILs. This approach substantially increased the efficiency of linkage analysis by enabling marker development by simple comparison of sequence data from the two parents. In addition, the application of whole-genome sequencing to candidate gene analysis reduced the cost and time for candidate gene identification by skipping the high-resolution mapping step. By QTL analysis, we detected a total of 27 QTLs for six agronomic traits. The 15 novel QTLs, two QTL hotspots and 15 strong candidate genes identified in this study provide valuable resources for molecular breeding of new high-yielding plant varieties in rice.