The total DNA of 18 japonica rice was rice extracted from fresh leaves (10–20 mg) selected at the seedling stage using the Qiagen DNeasy Plant Mini Kit. The experimental procedure for DNA isolation was the same as that previously described (Stein et al., 2001). Following the recommended protocol of GBS, optimized by Poland et al., a restriction digest buffer NEB Buffer 4 and a core set of restriction enzymes, PstI (CTGCAG) and MspI (CCGG), were added to each sample to obtain a uniform distribution of cut sites across the rice genome (Poland J. A. et al., 2012). The DNA library fragment size and library concentration were analyzed using a bioanalyzer (Agilent) machine (Parson et al., 2013). To inactivate the enzymes, the samples were incubated twice at 37°C for 2 h and at 65°C for 20 min. A set of already developed adapters with 500 ng of high-quality DNA was attached to the sample ends preceding the ligation reaction. The ligated samples were pooled and selected on a 1.5% agarose gel (200–300 bp) followed by PCR-amplification in a single tube. The resulting fragments in the DNA library were sequenced using the Ion Torrent (PGM) sequencing machine with Ion Torrent kits (PGM 200 Kit v2) following the manufacturer specifications (Life Technologies, Carlsbad, CA, and U.S.A.). The sequencing data obtained in our experiment have been deposited into the NCBI short read archive (SRA) under the study accession number:SRR2758809.

The FASTQ format sequence reads obtained from the sequencing machine were treated for quality control. For this, the adopter sequences and low quality reads were analyzed and trimmed. The cleaned sequences of all 18 samples were aligned onto the Nipponbare genome (IRGSP-1.0) (http://rapdb.dna.affrc.go.jp/) with bowtie 2 software (Langmead, 2010). The –M4 –very-sensitive local parameter for mapping was adopted (-D 20 -R 3 -N 0 -L 20 -i S,1,0.50), where if the alignment is greater than 4, it will output only 1 alignment at a random level. The reads that mapped successfully once onto a reference genome were considered to be uniquely mapped reads, whereas the reads mapped onto multiple locations were considered to be multiple-mapped reads. The short reads that could not be mapped onto any part of the reference were called unmapped reads. The GBS data were further used for SNP discovery using the Tassel-GBS application with default parameters (Glaubitz et al., 2014). The SNPs that evolved repeatedly were deleted from the final results.

The physical locations of the discovered SNPs were annotated by SNPEff software (Cingolani et al., 2012). The SNPs that were situated within exon, intron and other corresponding regions were classified in detail. We also grouped the identified SNPs into different classes of polymorphisms, such as transition/transversions mutations and gene-base synonymous and non-synonymous nucleotide base mutations.

In our experiment, all the plant materials belonged to Japonica rice, comprising nine CMS and nine restorer lines. These parental lines are widely used in China for the commercial breeding of three-line hybrid rice. For hybrid development, all the CMS and restorer lines were planted and manually crossed in a set of 9 × 9 combinations following the NCII mating design (North Carolina mating design II) (Shukla and Pandey, 2008). At the end of rice-cropping season, a population of 81 hybrids was obtained. The resultant hybrid seeds were grown and evaluated the next year at the Jiangpu research station, Nanjing Agricultural University, Nanjing China, with randomized complete block designs in four replications. The seeds of hybrid and their parents (restorer lines, and maintainer lines instead of CMS lines) were sown in the rice nursery in May 10, 2015. After 30 days, the seedlings of each hybrid combination were transplanted to paddy fields. Each plot in the paddy contained 4 rows with 8 plants per row. Twenty centimeter row-to-row and plant-to-plant spacing for all hybrids and their parents was maintained. Recommended field management and cultural practices were applied accordingly to the local environment. At maturity, six plants from each plot were randomly selected to calculate the phenotypic data of 12 yield-related traits. The phenotypic data of the traits were tested for significance using the statistical method of analysis of variance, along with the coefficient of variation (CV %) as described by Gomez and Gomez (1984). The GCA effect values of each CMS and restorer line for the traits were calculated by using the phenotypic data of the developed hybrids, as per the Griffing model (1956), following method-2 (Mo, 1982). The significance of the parental GCA value was tested with an LSD test at a P < 0.01.

To identify the SNP locus/loci associated with the CA of the yield-related traits, an association analysis between the identified SNPs and CA of 12 parental traits was performed using a computational software called CA screen 1.0 operated in the MATLAB language and developed by our laboratory^* (Liang et al., 2010). The script of our association method follows the principle of single-marker analysis (SMA). This method of association led to the statistically significant identification of SNPs and their effect on CA in homozygous and heterozygous associations. Moreover, we explain the principle of this association model corresponding to the research article provided, where we developed 81 F₁ combinations by crossing nine CMS with nine restorer lines. Now, at a given locus of SNP, if the parental lines of 41 F₁ hybrids possess the heterozygous SNP genotype (for example A-G, A-T), and the parental lines of the remaining 40 F₁ hybrids possess the homozygous SNP genotype (A-A or G-G). Now, after association, if the average trait value of the 41 heterozygous crosses is significantly greater or less than the average trait value of the 40 homozygous crosses, then the SNP marker associated with the CA of the trait is significantly positive or negative. If the difference in the trait value in the heterozygous association is positive and we observed a positive effect on the CA of the trait, we consider this SNP locus to be a favorable associated marker genotype of the elite CA for the trait, and vice versa.

In our study, to avoid spurious associations, the significance of the associations was assessed by using a t-test at P < 0.01 (Pradeep et al., 2007). Furthermore, coefficient of determination (R²) was calculated to determine the percentage of phenotypic variation explained by each SNP marker.

To determine the substantial number of functional SNPs, we selected associated SNPs based on their genomic positions. Only the SNP loci situated within the coding region of the genome (synonymous, non- synonymous and splice variants) were considered. Furthermore, in the NCBI, the database GENE was used to search the corresponding genes with associated SNPs for their nomenclature, chromosomal location, gene product, protein domain, and expected functions.

Using a GBS assay, approximately 16 million (16.3 × 10⁶) short reads were constituted, which comprised 4.7 GB of data. After removal of low quality reads, the sequenced data were mapped onto the Nipponbare reference genome. The mapping result indicated that 60.4% (9.86 × 10⁶) of the short reads were uniquely mapped onto the 12 chromosomes of the Nipponbare genome. Approximately 17% (2.78 × 10⁶) of the short reads were mapped more than once (Figure 1). The remaining 22.4% (3.67 × 10⁶) of the short reads could not be mapped onto any part of the reference genome.

Within the mapped sequenced reads of the nine CMS and nine restorer lines, a total of 39,001 SNPs were discovered at 11,085 sites, including 13,186 and 25,815 in the CMS and restorer lines, respectively. Despite having the same coverage in the genomes, the SNPs identified in the restorer lines were double those in the CMS lines. Of the nine CMS lines, 13,186 SNPs were identified, where the number of SNPs within the CMS lines ranged from 365 in Aizhixiang A to 3412 in Chunjiang 18A (Supplementary Table 1). The discovered SNPs in the CMS lines were distributed non-randomly on the 12 chromosomes, where chromosome 11 occupied a maximum (2740) and chromosome 5 occupied a minimum (381) of SNP numbers (Figure 2). The density distributions of the discovered SNPs in the CMS and restorer lines were explored by calculating the SNP frequency within a 200 kb genomic region. The average SNP densities between the CMS lines per 200 kb window were 7.5. Likewise, the density distribution varied across the chromosomes, where chromosome 9 possessed the highest SNP densities of 18.3, while chromosome 5 possessed the lowest SNP densities of 3.2. This phenomenon was also found to be uneven within the individual chromosome. For example, chromosome 11 in the CMS lines was found to be dense (750 SNPs) from the 25 to 28.8 MB genomic region, but scant SNPs (48 SNPs) were observed from the 11 to 16 MB genomic region. Similarly, on chromosome 6, the 10 MB region from 1.8 to 11.8 contained 542 SNPs, but on the same chromosome, the region from 16.6 to 22.4 MB had only 53 SNPs. On chromosome 9, we observed 34 high-density SNP sites with >30 SNPs per 200 kb region (Figure 3).

The discovered SNPs in nine restorer lines were counted as 25,815. These SNPs ranged from 1508 in R4179 to 5397 in LR5 (Supplementary Table 2). Due to the different chromosomal lengths, the numbers of SNPs across individual chromosomes varied. For example, for all the restorer lines, chromosome 11 was found to be rich in SNP number with 4845 SNPs, whereas a low SNP number (464) was found on chromosome 4 (Figure 4). The average densities of the SNPs among the nine restorer lines were 15 SNPs per 200 kb genomic region. Such densities across chromosomes were non-randomly distributed. The highest (35.6) SNP density was observed on chromosome 11, whereas the lowest (2.7) was observed on chromosome 4. Such densities were also found to be uneven within the individual chromosomes. In chromosome 2, the region between 5 and 22.2 MB had an unusually high (2796 SNPs) SNP density, whereas the region between 26 and 32.4 MB had a very low (24 SNPs) SNP density. The same results were observed on chromosome 11, where a high (3508 SNPs) and long SNP interval was detected between the 16 and 28.4 MB genomic region (Figure 5). A small (120 SNPs) and short interval was also found from the 0.4 to 4 MB region of the same chromosome. On chromosome 11 of the restorer lines, we identified a maximum of 57 high-density SNP sites, with >30 SNPs per 200 kb region.

The annotation of the Nipponbare rice was used to locate the distribution pattern of the identified SNPs within various genomic regions. Overall, a similar SNP distribution was observed in the CMS and restorer lines. Of the 13,186 SNPs discovered in the nine CMS lines, 7582 (57.5%) were found within intergenic spaces. The other 3212 (24.3%) and 2649 (20%) were in the intron and exon regions, respectively. Among the total SNPs in the coding regions, synonymous were counted as 1155 (8.7%), while non-synonymous were further divided into missense, 1214 (9.2%), and nonsense, 27 (0.2%) (Table 1, Supplementary Table 3). Moreover, the detected SNPs of the nine CMS lines were also categorized as transition (C/T and G/A) and transversion (G/T, T/A, A/C, C/G) nucleotide bases (Supplementary Table 5). In our results, the average ratio between the transitions and transversions was 1.4. The number of transition substitutions was significantly higher than that of transversion substitutions in the identified SNPs of all nine CMS lines (Figure 6).

Subsequently, the identified SNPs of the restorer lines were also annotated, where 14,970 (58%) were situated in the intergenic region of genome. However, 6226 (24%) and 5377 (20%) were detected within the intron and exon regions (Table 2, Supplementary Table 4). Of the SNPs in the coding regions, 2189 (8.4%) led to synonymous while 2711 (10.5%) led to non-synonymous amino acid changes with a protein-altering effect on genes. The average ratio of the transitions and transversions in the identified SNPs was also 1.4 (Supplementary Table 6). The similar result of higher transition substitutions was also observed in the identified SNPs of all nine restorer lines (Figure 7).

The GCA performance of the 18 parents of hybrid rice for 12 yield-related traits is summarized in Table 3. Among the nine CMS lines, genotype Zhe04 exhibited the highest GCA effect value for plant height, whereas A171 had the highest GCA effect for days to heading and seed width. In terms of the panicle length, the number of panicles per plants and the number of grains per panicle, CMS line 863A was a good general combiner due to the highly significant and positive GCA effect value. CMS 18A recorded the maximum GCA effect value for seed thickness and grain yield per plot, whereas 90167A revealed the highest GCA effect value in terms of seed length and 1000-seed weight. Parent AizhixiangA and Chunjiang18A had the maximum GCA effect in terms of number of panicles per plant and seed setting rate, respectively (Table 3).

The GCA assessment of the restorer lines indicated that LR27 recorded the maximum GCA effect value for plant height, whereas Shenhui254 recorded the maximum GCA effect value for days to heading. The restorer C4115 possessed higher GCA effect values for panicle length, seed thickness and seed-setting rate. In terms of seed length, seed width and 1000-seed weight, restorer Yanhui R8 recorded the highest GCA effect value. The parent LC64 exhibited a greater value for the number of spikelets per panicle, whereas LR5 exhibited a greater GCA effect for the number of filled grains per panicle and grain yield per plot. The restorer R4179 revealed the maximum GCA for the number of panicles per plant (Table 3).

The main objective of our study was to identify SNP locus/loci of CA via association analysis; the discovered 39,001 SNPs at 11,085 genomic positions were integrated with the parental CA of 12 yield-related traits at P < 0.01. We revealed a total of 362 SNP locus/loci with the CA of parental traits that caused a positive or negative effect on F₁ trait performances. The overview and detailed information of the identified SNP loci of CA for the following traits are presented in Figure 8 and Supplementary Table 7.

The exploitation of heterosis in rice has increased the global rice-grain yield. This has led to an increase in productivity (growth, size, development, fertility and yield) of the progeny over that of both homozygous parents and is a fascinating phenomenon (Hochholdinger and Hoecker, 2007). The application of heterosis has been applied not only to rice but also to other crops such as tomatoes (Williams and Gilbert, 1960), maize (Shull, 1908), wheat (Singh et al., 2004), and soybeans (Pandini et al., 2002). In plants, heterosis is considered to be a complex trait that involves the presence of potential genetic distances, combining abilities and heterotic patterns of the parents traits (Beck et al., 1990). Rice breeders emphasized that heterosis in F₁s is caused by the presence of combining ability alleles, that are present at different loci of crossing parents (Liu et al., 2002). They described how the combination of different alleles at a specific locus of crossing parents creates heterozygosity that results in the heterosis of offspring.

One part of our present study describes genome-wide SNP discovery in japonica rice. A large number of rice cultivars have been sequenced for SNP discovery to perform genetic and genomic experiments. However, only a few studies have been conducted on japonica rice (Nagasaki et al., 2010). Rice breeders have concluded that SNP discovery in japonica rice is challenging compared to that in indica rice, due to its close genetic relatedness and lower diversity (Glaszmann, 1987; Yang et al., 1994; Gao et al., 2005; Negrao et al., 2008). In our recent efforts, 18 parents of japonica rice were sequenced for SNP discovery, following optimized genotyping with a sequencing method. The use of two restriction enzymes based on a sequencing method of GBS is the most convenient, cost-effective and easy approach to the SNP discovery of large genome crop plants, including rice (Furuta et al., 2016). Nevertheless, this method also has the drawback of generating a large amount of missing data (Davey et al., 2011; Kim et al., 2016).

Here, the successful mapping of 60.4% of short reads onto the reference genome revealed 39,001 SNPs at 11,085 positions. The unique mapping of short reads plays a crucial role in the determination of SNPs (Nielsen et al., 2011). However, un-mapping reduced the chances of SNP discovery, which eventually occurs due to genomic deletions during the process of sequencing (Arai-Kichise et al., 2011). In our study, the number of SNPs among the parents varied. We observed higher polymorphisms in the nine restorer lines compared to the CMS lines. Of the 18 sequenced parents, the CMS Aizhixiang A reported a minimum number of SNPs. An in-depth analysis of the sequence genome of Aizhixiang A indicated that 70% of its sequence genome was similar to that of the Nipponbare reference. Furthermore, 25% of the missing data were also found, which further affected the chances for SNP discovery. We also found that the SNPs of CMS and the restorer lines were non-randomly distributed across 12 chromosomes and revealed several SNP-rich and SNP-poor intervals. Such extended SNP regions were similar to those reported in rice, wheat, mei (Prunus mume sieb. et zucc) and peach cultivars (Ravel et al., 2006; Subbaiyan et al., 2012; Fresnedo-Ramírez et al., 2013; Sun et al., 2013). In our study, there were also non-uniform patterns of SNPs within individual chromosomes, particularly at the terminal region of chromosome. We observed this situation at the terminal region of chromosome 11 for both the CMS and restorer lines. This may possibly arise due to the low variant rate at the centromere region of the chromosome, where recombination are typically more rare than at the chromosome terminal (Guo W. et al., 2007).

In the sequenced genome of animals and plants, SNPs do not present equally throughout the genome. Some parts of the genome revealed more SNPs than others. The discovered SNPs in our study were annotated according to their positions in the genome. The annotation of identified SNPs indicates that more than half (57%) of the SNPs are within the intergenic region, whereas 24.4% of the SNPs have been found within the intron region of the genome. Only 20.4% of the SNPs occurred inside the exon regions, with 8.6% synonymous and 9.9% non-synonymous in nature. The number of SNPs in the exon region was 37% less than that in the intergenic regions. This is possibly, because of the higher frequency of polymorphisms evolving in intergenic spaces rather than coding sequences (Guo X. et al., 2007). However, the SNPs situated inside intergenic regions are also prominent, where there may have been the possibility of the presence of some functional alleles or genes (Salvi et al., 2007). In our study, the ratio of transitions to transversions was 1.4. This ratio was in line with former sequencing experiments in flax (Kumar et al., 2012).

The concept of combining ability (GCA and SCA) was developed by Sprague and Tatum in 1942. They characterized the breeding values of inbred lines by testing and evaluating them by combining ability evaluation trails (Sprague and Tatum, 1942). In every hybrid-breeding program, a breeder's goal is to select parents with good general combining ability and crosses that have high specific combining ability in their traits. From a genetic point of view, GCA evolved under additive and additive × additive gene action, whereas SCA evolved under non-additive gene action. The GCA of inbred parents is considered more important than the SCA, in which the parents were involved in random mating. The GCA evaluation of 18 parents of japonica hybrid rice indicated that most of the CMS lines exhibited maximum and positive GCA effect values for their studied traits. Previously, the same results of good general combining ability in the CMS lines were also reported (Singh et al., 1996). In contrast, most of the restorer lines also revealed greater GCA effect for their traits. The evidence of a higher GCA of restorer lines was also found in earlier studies (Rogbell et al., 1998).

With the discovery of genome-wide SNPs, several SNP-trait associations were reported in plant species, including maize, wheat, cotton, soybean, barley and Arabidopsis thaliana (Atwell et al., 2010; Yang et al., 2010; Pasam et al., 2012; Hu X. et al., 2015; Zhang et al., 2015; Su et al., 2016). In the past, most of the association studies on rice have focused on the improvement of yield-related traits (Huang et al., 2010; Zhao et al., 2011). Very few SNP associations with CA of traits have been performed. Here, we studied the heterozygous combinations of nucleotide bases (A-T/G-C) at a locus and observed its associations with the CA of traits. Moreover, we worked directly with the heterozygous association group and revealed 362 CA loci of 12 yield-related traits of japonica rice. For those results, each of the CMS and restorer lines contained a number of favorable or unfavorable SNP locus/loci of CA in the genome that caused a positive or negative effect on the traits. Few of the associated SNP loci were found to be co-associated with the CA of more than one yield-related trait. Among the CMS lines, 18A had a maximum number of negative SNP loci of CA for panicle and grain-related traits; whereas CMS A171 had a maximum number of positive SNP loci of CA for grain thickness and day to heading. For the restorer lines, LC64 and LR27 contained maximum numbers of positive SNP loci of CA for most of the studied traits. Previously, in maize, several favorable and unfavorable marker loci of CA were detected for five yield-contributing traits (Qi et al., 2013). Of these, the positive improved, while the negative decreased, the CA of the traits. Similarly, Qu et al. identified the CA loci of agronomical traits in rice by using QTL mapping of BCRILs (Qu et al., 2012). Moreover, all of the identified and associated SNP loci in our study were spread genome wide, with some on chromosomes 2, 5, 7, 9, and 11 exhibiting a maximum positive effect for the CA of traits.

The correct identification of SNPs based on annotation results made them informative. Previous studies have explained that the SNPs within the coding region can alter the promoter activities of gene expression, transcription and translation capability (LeVan et al., 2001; Alonso-Blanco et al., 2005; Ng and Henikoff, 2006; Shastry, 2009; Batley, 2015). Such identification and characterization of gene-based SNPs could be directly utilized in any breeding program regarding crop improvement (Wang et al., 2008; Huang et al., 2010). Among the 362 SNP loci detected in this study, 32 SNP loci for CA of seven traits were situated within the reported genes (Supplementary Table 8). For CA of plant height, we found 7 SNP loci within previously identified genes. These genes encode proteins for photosynthesis, hormonal activity, lipid binding and transporter activity. Of the seven reported genes for plant height, Os03g0309200 encodes a proteins causing significant effect on photosynthesis activity, flowering time, and also have strong effects on stem and inter-node elongation in rice. Among four SNP loci for CA of grain width detected within the characterized genes, Os06g0225300 is directly involved in improving rice architecture for high yield. Overall, the biological function of most of the reported genes were involved in photosynthesis, flower and embryo development, post-embryonic development, multi-cellular organismal development, cellular processes and metabolic processes. Similarly, the molecular function of these genes revealed DNA binding, lipid binding, nucleotide binding, protein binding, DNA binding transcription factor activity, catalytic activity and kinase activity. Such SNP loci situated inside the genes could be used as a candidate marker locus for tagging genomic regions responsible for improving CA of yield traits.

In conclusion, we present the identification of 215 positive and 147 negative coding-based SNP loci of CA for heterosis breeding in japonica rice. Our study implies that the identified SNP loci of CA will increase the selection efficiency of rice breeders in the proper selection of inbred parents. These insights into rice will be productive for future work to improve the CA of parents for yield-related traits and develop superior japonica hybrids by utilizing breeding designs. We expect that, the parental lines with more positive SNP loci of CA and less negative SNP loci of CA will be useful in hybrid breeding, because they have much higher potential for CA improvement for the traits.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.