The flowering date were tagged by panicles of each accession, and seeds and panicles having same flowering date (Table 1 and Supplementary Tables 1, 7) were harvested at 42 days after flowering (DAF) in the various rice accessions. The moisture content was adjusted by drying at 15°C (RH 10%) during a 7 days period, as previously described (Lee et al., 2016). The susceptibility for PHS (PHS value) was surveyed using three panicles per each accession, which were incubated at 25°C (RH 100%) for 7 days, after which the number of germinated seeds on each panicle was recorded and expressed as a percentage of the total grain number per panicle. The lower the PHS value, the lower the susceptibility to PHS and the higher the resistance. Seed germination at harvesting time (GHT) was determined by threshing the panicles and planting three replications of 50 seeds (fruits in hulls) onto moistened Whatman filter paper (10 mL of distilled water) in Petri dishes. The seeds were incubated at 25°C (RH 100%) based on International Seed Testing Association (ISTA) guidelines, and germinated seeds (radicle and coleoptile emerged from the hull) were counted daily for a period of 10 days. The cumulative number of germinated seeds was expressed as the percentage of seeds planted. Germination index (GI) was calculated as described by Basra et al. (2005).

Twenty-one samples (japonica, 14; indica/tongil, 7; Table 1), revealing similar flowering times (early August), were selected as representatives by PHS resistance (PHS resistant representatives: PHS < 20%, PHS susceptible representatives: PHS > 40%) at 42 DAF among the diverse rice genetic resources (Lee et al., 2016) in that environmental noises might reduce the resolution for the detection of PHS associated genetic loci in the field condition test. The samples were re-sequenced using the Illumina HiSeq 2500 platform to identify genome-wide variations and to detect genetic signals associated with PHS.

Genomic DNA was extracted from the 21 rice samples using a Gentra Puregene Cell Kit for Plants (Qiagen, Hilden, Germany). Library construction and sequencing data collection were conducted using Illumina’s official protocol with a 101 bp paired-end read length. Trimmomatic-0.33 (Bolger et al., 2014) was used to remove adapters and low-quality reads. The adapter-free trimmed reads were mapped against the reference genome of Oryza sativa L. (IRGSP-1.0.27) using Bowtie2 with default parameters (Langmead and Salzberg, 2012). The mapped reads were assigned to read groups and sorted using Picard version 1.138.¹ Picard was also used to remove potential PCR duplicates and to repair mismatches between each read and its mate pair. For the subsequent local realignment, base quality recalibration, and variant calling steps, Genome Analysis Toolkit (GATK) version 3.4.46 (McKenna et al., 2010) was used. Local realignment of reads was carried out to correct misalignments caused by the presence of InDels. Base quality recalibration was performed to compensate for base quality errors from empirical measurements. For variant calling, arguments including “UnifiedGenotyper” and “SelectVariants” were used. Finally, the identified variants were filtered using the “VariantFiltration” argument based on the following criteria: (1) Reads with a mapping quality of zero, MQ0 higher than 4, and MQ0/(1.0^∗DP) higher than 0.1, where DP is the unfiltered read depth, (2) FS higher than 200 to reduce false positives, and (3) Phred-scaled quality score lower than 30.

To identify specific SNPs in the PHS groups (PHS resistant and susceptible group), Fisher’s exact test implemented in SNPSift was conducted to analyze the resulting genotype count data (Cingolani et al., 2012) contrary to the linear or mixed model generally used for genome-wide association studies (GWAS) (Korte and Farlow, 2013). Information about PHS-specific groups and the genotype data assigned to the dominant and recessive models of the 2 × 2 contingency tables were used in this statistical test. To minimize the rate of false positives, multiple Bonferroni correction was conducted. SNPs were considered significant when the p-value of the total number of tests was below 0.05. Significantly enriched SNPs were annotated with SNPeff.

Thirty-nine candidate SNP loci (for genotyping diverse rice genetic resources) were selected, including SNPs in germination-related genes with hormonal action, such as ABA perception and signal transduction, as well as SNPs in genes expressed specifically in florets and grains (Supplementary Table 8). A total of 144 accessions (80 japonica accessions and 64 indica accessions) showing PHS variations (Supplementary Table 7) were used to genotype the candidate SNPs, and an additional 56 japonica accessions (mainly Korean landraces) were used to validate the regression model.

Genotyping was performed using Fluidigm 192.24 Dynamic Array integrated fluid circuits (Fluidigm Incorporated, San Francisco, CA, United States). The IFC Controller utilizes pressure to control the valves in the chips and to the load samples and genotyping assay reagents into the reaction chambers. The EP-1 system was used as an endpoint image reader. Specific target amplification (STA) reactions were performed in a GeneAmp PCR System 9700 from Applied Biosystems.

For SNP genotyping, 5 μL of STA mix was prepared for each sample, containing 1 × Qiagen Multiplex PCR Master Mix (Qiagen, PN 206143) and 50–60 ng of genomic DNA. After diluting the STA mix, an allele-specific SNP genotyping assay was carried out with the SNPtype assays protocol in the Fluidigm 192.24 Dynamic Array IFC User’s Guide (PN 68000098 N1). After PCR amplification, the endpoint fluorescent image data were acquired on the EP-1 System. Data were analyzed using Fluidigm SNP Genotyping Analysis software to obtain genotype calls.

Since the japonica and indica ecotypes have different genetic constitutions, the accessions were separately analyzed based on ecotype (Morishima and Oka, 1981; Sun et al., 2002, 2015). Despite the scarcity of suitable molecular markers for seed germinability, some reported markers were available, such as InDel and SNP loci in qLTG3-1 controlling low-temperature germinability in rice (Fujino et al., 2008; Fujino and Sekiguchi, 2011); these markers (and modified markers) were used to elaborate the regression equation for PHS.

One-way analysis of variance (ANOVA), Duncan’s multiple range test (DMRT), correlation analysis, and regression analysis were conducted using R software (ver. 3.2.3)², and differences were considered significant when p-value < 0.05. The spatio-temporal expression patterns of genes in various tissues throughout plant growth were surveyed using the Rice Expression Profile Database (RiceXPro)³.

A total of 21 accessions with similar heading dates were sequenced using an Illumina HiSeq 2500 sequencer to identify polymorphisms between PHS resistant and susceptible accessions (Table 1). The genomes of these selected accessions were sequenced at an average of ∼25.7 × coverage, with an average of 9.59 Gb total bases (8.01–10.82) and 94.9 million reads (Supplementary Table 1). The sequences were aligned against the reference genome of Oryza sativa L. (IRGSP-1.0.27) using Bowtie2 with default parameters (Langmead and Salzberg, 2012). The alignment rate ranged from 92.8% (IT226904, Koshihikari) to 98.7% (IT 259863, KELEE), with a mean alignment rate of 96.8%. After filtering of potential PCR duplicates and correction of misalignments due to the presence of InDels, the genome-wide SNPs were detected using GATK (McKenna et al., 2010). To minimize the number of false-positive calls, we used several filtering steps before subjecting the candidate SNPs to further analyses based on the following: Phred-scaled quality score, mapping quality, quality depth, and Phred-scaled p-value. We ultimately acquired approximately 4.27 million genome-wide SNPs, with an average ratio of one SNP per 87 bp (Supplementary Table 2).

We conducted enrichment analysis using Fisher’s exact test implemented in SNPSift (Cingolani et al., 2012) to detect allelic differentiation between PHS resistant and susceptible accessions and annotated significantly enriched SNPs in each PHS grou2.p using SNPef. Among the 4.27 million SNPs detected genome-wide, approximately 21,000 and 18,000 SNP loci showed distinguished allelic distribution between PHS resistant and susceptible accessions (21 accessions, including 14 japonica and 7 indica/tongil) based on dominant and recessive models (p-value < 0.05), respectively, and the common SNPs detected in both model were 16,753 loci (Table 2 and Supplementary Tables 3, 4). Among the 14 japonica accessions, approximately 5,600 and 4,300 SNPs were detected based on dominant and recessive models, respectively, and the common SNPs detected in both models were 4,115 loci (Table 2 and Supplementary Tables 5, 6).

Among the SNPs showing allelic differentiation between PHS resistant and susceptible accessions, the number of SNPs in genic regions was approximately 4,100 and 1,000 among the 21 total accessions (japonica, 14; indica, 6; tongil, 1) and the 14 japonica accessions, respectively (Supplementary Tables 3, 5), based on the dominant model and approximately 3,700 and 800 among the 21 and 14 accessions, respectively, based on the recessive model (Supplementary Tables 4, 6). Among the genic SNPs, intron variants were most prevalent among the 21 accessions (dominant model: 48.3% and recessive model: 48.0%), followed by exon variants (34.0 and 34.3%), 3′ UTR variants (11.4 and 11.3%), and 5′ UTR variants (6.4 and 6.4%). Among the 14 japonica accessions, exon variants were most prevalent (dominant model: 43.4% and recessive model: 44.6%), followed by intron variants (31.1 and 28.9%), 3′ UTR variants (16.5 and 16.9%), and 5′ UTR variants (9.1 and 9.5%).

We constructed a density plot of SNPs in PHS resistant versus susceptible accessions using a 1 Mb bin size along the rice genome (Figure 1). Among the SNPs detected genome-wide, an average of 55 SNPs per bin revealed allelic differentiation between PHS resistant and susceptible accessions. Chromosome-1 has the highest number of SNPs, with an average of 142 per bin, followed by chromosome-11 (78 per bin) and chromosome-12 (69 per bin). Chromosome-8 has the fewest SNPs (12 per bin). Among the japonica accessions, only an average of 13 SNPs per bin revealed allelic differentiation between the PHS resistant and susceptible groups. Chromosome-11 has the highest number of SNPs, averaging 36 per bin, followed by chromosome-12 (34 per bin) and chromosome-4 (16 per bin). The fewest SNPs were found on chromosome-10 (three per bin). We compared the locations of the SNPs with publicly reported QTLs related to PHS and dormancy along the rice chromosome, finding that some chromosomal regions contain high numbers of these SNP loci. Region-1 (Chr.1: ∼1–5 Mb, 71 genic SNPs / 184 common SNPs) has a high density of SNPs overlapping with QTL Sdr6 (Marzougui et al., 2012), qSD-1 (Miura et al., 2002), qSD1 (Gu et al., 2004, 2006) and qDEG1 (Li et al., 2011). Region-2 (Chr.1: ∼31–33 Mb, 140 genic SNPs/787 common SNPs) contains the QTL qSdn-1 (Lu et al., 2011), while Region-3 (Chr.3: ∼3–11 Mb) contains the QTLs Sdr1 (Takeuchi et al., 2003), qDT-SGC3.1 (Jiang et al., 2011) and qSD-3 (Guo et al., 2004). Region-6 (Chr.7: ∼22–25 Mb, 91 genic SNPs / 787 common SNPs) harbors QTLs overlapping with qSD-7-2 (Gu et al., 2004), qMT-SGC7.2 (Jiang et al., 2011), qPHS-7 (Dong et al., 2003) and fine-mapped Sdr4 gene (Sugimoto et al., 2010). Region-7 (Chr.12: ∼21–27 Mb, 386 genic SNPs/1,391 common SNPs) harbors the QTL qSD-12 (Gu et al., 2004) and fine mapped qSD-12 gene (Gu et al., 2010). Although several studies reported QTLs in Region-4 (Chr.5: ∼24–28 Mb, 1 genic SNPs/5 common SNPs) and Region-5 (Chr.6: ∼6–12 Mb, 0 genic SNPs/4 common SNPs), we detected few SNPs in these regions in the current study. Other chromosomal regions revealing high density of SNPs without overlapping with reported QTLs included several germination associated loci by GWAS study (Magwa et al., 2016).

Among the approximately 20,000 SNP loci revealing distinct allelic distribution between PHS resistant and susceptible accessions, we genotyped 144 rice germplasm (80 japonica accessions and 64 indica accessions, Supplementary Table 7) using the 39 SNPs mainly located in genic or promoter region of candidate genes (Supplementary Table 8). We analyzed the differences in PHS and GI values according to genotype. Six and seven loci, including qLTG3-1 (Fujino et al., 2008; Fujino and Sekiguchi, 2011) revealed significant differences in PHS and GI values according to allele type in 80 japonica accessions (p-value < 0.05 and p-value < 0.01), respectively (Table 3). Furthermore, the existence of an awn was also significantly associated with PHS and GI values. Three significant SNP loci were associated with PHS and GI values in 64 indica accessions. Specifically, among significant SNPs associated with PHS in japonica and indica rice, V2 is a significant variant in an upstream region of the rice homolog ortholog of AtBG1, which functions in ABA conjugation in Arabidopsis. V5 is a variant in the intron region of a ROS homeostasis gene, and the S4 locus is located in the upstream region of OsVP1. S13 and S21 are missense variations in the BRCT domain-containing protein gene and the Guanine nucleotide-binding protein gene, respectively. S3 and S9 are located in a Dynamin family protein gene, representing a missense variation and an extra stop codon, respectively. S16 and S19 are missense variations in CINNAMYL ALCOHOL DEHYDROGENASE 4. Finally, S1, which is significant only in the indica group, is a missense variation in a hypothetical protein gene.

Seven new PHS- and GI-associated loci (V2, S4, V5, S13, S21, S16, and S19) in the japonica group are not included among previously reported QTLs, whereas the Dynamin family protein gene (OS03G0260000), containing the S3 and S9 variations associated with GI, overlaps with reported QTLs including Sdr1 (Takeuchi et al., 2003), qDT-SGC3.1 (Jiang et al., 2011) and qSD-3 (Guo et al., 2004). Among the three loci that are significantly associated with PHS in the indica group, S1, located on the short arm of chromosome-1, overlaps with reported QTLs including Sdr6 (Marzougui et al., 2012), qSD-1 (Miura et al., 2002), qSD1 (Gu et al., 2004, 2006), qDEG1 (Li et al., 2011). We investigated the spatio-temporal expression patterns of genes harboring significant PHS- and GI-associated SNPs by surveying a publically available expression database (RiceXPro).⁴ OsVP1 (S4 variation), BRCT domain containing protein (S13 variation), and Hypothetical protein (S1 variation) are highly expressed in embryos and endosperm compared with other tissues, whereas G-protein (S21 variation) is expressed at low levels in embryos and endosperm during seed maturation (data not shown). Dynamin family protein (S3 and S9 variation), THIOREDOXIN H-TYPE (V5 variation), and OsCAD4 (S16 and S19 variation) are more highly expressed in endosperm than in embryo tissue, while BETA-GLUCOSIDASE 36 (V2 variation) is expressed at lower levels in endosperm than in embryo tissue.

We conducted regression analysis based on the genotypes of 11 significant loci including qLTG3-1 and PHS values among 80 japonica accessions. The genotypes of significant loci and the presence/absence of awns were converted into binary data (Supplementary Table 9): the reference genotype was coded as “0,” while alternative types were coded as “1”; the absence of awns was coded as “0,” and their presence was coded as “1.” We performed regression analysis using this data set and calculated parameter estimates, p-values, and adjusted R² values accordingly. The regression model revealing the highest adjusted R² value was adopted for PHS estimation, with an adjusted R² value of 0.401 [PHS predicted value = 33.082 - 12.171(qLTG3.1) - 14.479(V2) - 11.629(V5) - 20.62(S4) + 22.544(S21) + 12.209(S3) + 10.864(S13) + 11.767(Awn)] (Table 4). The regression model is composed of eight factors, including six SNP loci (five loci – significant in PHS, one locus – significant in GI), one InDel locus of qLTG3-1, and the awn character.

To validate the estimation model for PHS, we genotyped 56 additional Korean japonica accessions and measured the PHS value at 42 DAF (Supplementary Table 10). The estimated PHS values derived from the PHS regression equation were calculated and compared to actual PHS measurements (Figure 2). The R² (coefficient of determination) value between the predicted PHS estimates and the actual measurements was 0.430 (p-values < 0.001).