A total of 227 O. sativa accessions were obtained from China (154), Vietnam (57), and Japan (16). The details of these accessions, including their origins and accession IDs, are summarized in Table S1. The seeds of all accessions were collected, stored and supplied by the State Key Laboratory of Crop Genetics, and Germplasm Enhancement at Nanjing Agricultural University.

Plants from the 227 accessions were grown in the paddy field of Jiangpu Experimental Station at the Nanjing Agricultural University in Nanjing (32°07′N, 118°64′E), Jiangsu Province, China, from May to October in 2013, and 2014. The field experiments conducted in each of the two consecutive years was treated as two independent environments. The field trials followed a completely randomized block design with two replicates per year. Each plot contained five rows, with eight plants in each row. The rows were spaced 20 cm apart, and each plant was spaced at a distance of 17 cm according to standard agronomic management practices.

After pistil maturation, 10 spikelets were randomly collected from five plants of each accession prior to glume opening (approximately 10:00 a.m.). For each spikelet, the STL, SBPL, and SNBPL (Figure 2A) of the fertile flower were measured with a micrometer using a stereomicroscope (10 ×, MC50, Guangdong, China), and the average values served as the measurements for the accessions.

To understand the relationship between stigma traits and grain length (GL), fully filled grains were measured. In each year, five normally developed plants from each accession were harvested from the middle of each of the plots at the grain ripening stage and dried under natural conditions for GL investigation. Ten randomly selected grains (after the awns were removed) from each accession were lined up length-wise along an electronic digital Vernier caliper (http://www.guanglu.com.cn) to measure the grain length, and the average values served as the measurements for the accessions.

Genomic DNA was extracted from young and healthy leaf blades of each accession approximately 3 months after germination using the methods described by Gross et al. (2009). According to the rice molecular map and microsatellite database published by Temnykh et al. (2000) and McCouch et al. (2002), 249 SSRs scattered on 12 chromosomes were selected. The SSR primers were synthesized by Shanghai Generay Biotech Co., Ltd. (Shanghai, China).

Each 10-μl PCR reaction contained 10 mM Tris-HCl (pH 9.0), 50 mM KCl, 0.1% Triton X-100, 1.5 mM MgCl₂, 0.5 nM dNTPs, 0.14 pM forward primer, 0.14 pM reverse primer, 0.5 U of Taq polymerase, and 20 ng of genomic DNA. DNA amplification was performed using a PTC-100™ Peltier Thermal Cycler (MJ Research™ Incorporated, USA) under the following conditions: (1) denaturation at 94°C for 5 min; (2) 34 cycles of denaturation at 94°C for 0.5 min, annealing at 55–61°C for 1 min, and extension at 72°C for 1 min; and (3) a final extension at 72°C for 10 min. The PCR products were separated on an 8% polyacrylamide gel for 1 h at 150 V and visualized using silver staining. One pair of SSR markers was used to detect one locus, and each polymorphic band at the same marker locus in the population was recorded as one allele. After screening the PAGE products, the size of each band was determined using Quantity One software (Bio-Rad Company, USA).

The genetic clusters in the 227 accessions were identified using STRUCTURE version 2.2 (Falush et al., 2007). This analysis was performed five times for each number of clusters (K) (from 2 to 10) with random starting points. We set the length of the burn-in period equal to 50,000 iterations and defined a run of 100,000 Markov Chain Monte Carlo (MCMC) replicates after burn in. A mean log-likelihood value over five runs at each K was used. If the mean log-likelihood value was positively correlated with the model parameter K, a suitable value for K could not be determined. In this situation, the optimal K value was determined through an ad hoc statistic, ΔK, based on the rate of change in [LnP(D)] between successive K values (Evanno et al., 2005). Non-admixed individuals in each genetic group were determined using a Q-matrix assignment greater than 0.9. A principal component analysis (PCA) was performed using the pcaMethods package (Stacklies et al., 2007) of R.2.11.1 (R Development Core Team, 2011) to examine the population structure. The first three principal components (PCs) were evaluated. The genetic distance was calculated based on 249 molecular markers using Nei's distance (Nei et al., 1983), and phylogenetic reconstruction was performed based on the neighbor-joining method implemented in PowerMarker version 3.25 (Liu and Muse, 2005). The tree used to visualize the phylogenetic distribution of the accessions and ancestry groups was constructed using MEGA version 4.0 (Tamura et al., 2007). The results from the three different approaches were then summarized and compared.

All basic statistical analyses were performed using the SAS package (SAS Institute Inc., Cary, NC, USA). The broad-sense heritability (HB2) was calculated based on the natural population through an analysis of variance using the formula HB2=σg2/(σg2+σe2/n), where σg2 is the genetic variance, σe2 is the error variance, and n is the number of replicates.

The polymorphic information content (PIC) was used to measure the probability that two randomly selected alleles from a population were distinct. The number of alleles per locus, the gene diversity, and the PIC values were determined using PowerMarker version 3.25 (Liu and Muse, 2005).

The average standardized individual allele size of the SSRs was determined by calculating the mean standardized size of 249 SSR loci following the method reported by Vigouroux et al. (2003). This standardization ensured that each locus contributed equally to the average individual size of the genome.

The D′ value (Farnir et al., 2000) was calculated with TASSEL 3.0 software using 1000 permutations and was used to measure the level of linkage disequilibrium (LD) between loci (Bradbury et al., 2007). Before the association analysis, rare alleles with allele frequencies of 5% or less were removed from the dataset. According to the level of LD and the genetic distance among markers with intrachromosomal combinations, the regression equation of the LD as a function of changes in the genetic distance was calculated through a regression analysis. A LD decay plot was drawn to observe the relationship between LD and genetic distance.

A Mixed Linear Model (MLM) analysis, which can significantly reduce spurious marker-trait associations (type I errors showing false positives) resulting from the population structure, was performed using TASSEL 3.0 to calculate the associations between traits and markers (Bradbury et al., 2007). The matrices Q and K were used as covariants in the MLM analysis. The Q matrix was adapted from the analysis results obtained from Structure 2.2. The K matrix (kinship matrix) was obtained from the results of the relatedness analysis using SPAGeDi (Hardy and Vekemans, 2002). According to the correction method published by Benjamini and Hochberg (1995), a false discovery rate (FDR) of 0.05 was used as a threshold for significant associations. Based on the identified association locus, the “null allele” (non-amplified allele) was used to determine the phenotypic effects of other alleles (Breseghello and Sorrells, 2006).

To detect a phylogeographic signal in O. sativa, we divided the distribution range into six regions: Vietnam (8°30′–23°22′N, 102°8′–109°24′E), Southern China (20°45′–30°08′N, 97°31′–124°34′E), Eastern China (27°12′–35°20′N, 116°18′–123°00′E), Northern China (31°23′–42°37′N, 110°21′–119°45′E), Northeastern China (38°43′–53°33′N, 118°53′–135°05′E) and Japan (30°38′–44°09′N, 130°49′–145°15′E). This delimitation was based on the latitude of each region. We performed a locus-by-locus analysis of molecular variance (AMOVA) (Weir and Cockerham, 1984) based on genetic groups delimited by the Bayesian clustering method in the program Arlequin 3.5 (Excoffier and Lischer, 2010) to statistically verify the geographical structure using SSR and standard multi-locus frequency data.

The mean value, standard deviation, skewness, and kurtosis for stigma traits measured in 227 rice varieties were calculated (Table 1). A significantly positive correlation was found between stigma traits in 2013 and 2014 (Figure S1). The materials exhibiting the minimum and maximum values are shown in Figures 2B,C. The stigmas exhibited longer SNBPLs and shorter SBPLs in some accessions (Figure 2D), and these trends were reversed in other accessions (Figure 2E). The average STL, SBPL, and SNBPL values of 227 accessions over the 2 years were 1.8 mm, 1.1 mm, and 0.7 mm, respectively. The phenotypic data for the STL, the SBPL, and the SNBPL in the studied population followed normal distributions (Table 1). Table 1 also shows that the average SBPL was greater than the average SNBPL. The broad-sense heritability values for the STL, SBPL, and SNBPL traits averaged over 2 years were 93.8, 92.8, and 88.6%, respectively. Furthermore, the three stigma traits exhibited significant positive correlations. The correlation coefficients between the STL and SBPL, the STL and SNBPL, and the SBPL and SNBPL were 0.803, 0.787, and 0.264, respectively (Table 2).

The highest average GL value over 2 years was 12.8 mm, the lowest average GL value over 2 years was 6.5 mm, and a skewed distribution was observed based on the skewness and kurtosis statistics (Table 1). The broad-sense heritability value averaged over 2 years was 98.3%. The correlation coefficients between the GL and STL, the GL and SBPL, and the GL and SNBPL were 0.710, 0.477, and 0.654, respectively (Table 2).

A marker analysis of the 227 accessions using 249 SSR molecular markers resulted in the detection of 2072 alleles. The numbers of alleles ranged from two (at locus RM7163_Chr11) to 18 (RM162_Chr6), with an average of 8.32 alleles per locus (Table S2). The genetic diversity averaged 0.699 and ranged from 0.083 (RM7163_Chr11) to 0.894 (RM162_Chr6) (Table S2). The PIC had a mean value of 0.668, ranged from 0.050 (RM7163_Chr11) to 0.885 (RM162_Chr6), and was mainly distributed between 0.509 and 0.809 (Table S2). A total of 210 markers (84%) were highly informative (PIC > 0.5), 27 (11%) were moderately informative (0.5 > PIC > 0.25), and 12 (5%) were slightly informative (PIC < 0.25).

The highest delta K values from the STRUCTURE analysis were obtained for eight clusters (Figure S2). The population structure data based on the Q matrix for each accession are summarized in Table S1, and the 227 accessions could be divided into eight subpopulations, viz. Group 1 to Group 8 and an admixed group (Figure S3A). Group 1 contained 38 accessions, which were modern improved varieties obtained mainly from Eastern China. Group 2 contained 45 accessions, which were mainly from Vietnam. Group 3 contained 20 accessions, and these were mainly obtained from Japan. Group 4 included 32 accessions, which consisted of landraces from Eastern China. Group 5 contained 21 accessions, which were mainly from Northern China. Group 6 contained 25 accessions, and these were mainly obtained from Southern China. Group 7 included 30 accessions, which were mainly from Northeastern China. Group 8 only contained four accessions obtained from North and Northeastern China. The admixed group contained the remaining 12 accessions.

Although the non-admixed individual cluster contained some admixed individuals, the neighbor-joining tree indicated that each of the eight groups corresponded to a distinct branch on the tree (Figures S3B, S4). These admixed individuals tended to be genetically similar to the group corresponding to their branch but exhibited a small amount of introgression from other groups.

The PCA results were essentially consistent with the STRUCTURE findings. The eight genetic groups were distinct in the principal component space of the three major PCs, whereas some individuals appeared mixed (Figure S3C). The three major PCs identified from the PCA explained 31.8% (11.4, 10.3, and 10.1%) of the total variance.

The following analysis did not include the accessions in the admixed group or Group 8 because Group 8 contained only four accessions.

To test whether the microsatellite size exhibited directional evolution in rice, we compared the average allele sizes in the geographically derived groups (Groups 1–7). Figure 3 shows the average individual microsatellite allele size for these seven groups. No significant differences were detected between Groups 1 and 4 (P > 0.05, two-tailed t-test), Groups 2 and 6 (P > 0.05, two-tailed t-test), or Groups 3 and 5 (P > 0.05, two-tailed t-test), whereas the rest groups exhibited significant differences, with P-values less than 0.001. The average individual SSR allele size in the subgroups from low latitudes (Groups 2 and 6) was 0.14, which is greater than the value obtained for subgroups from high latitudes (-0.04; Figure 3). These data indicate a decrease in allele size from south to north in the geographically derived groups.

A total of 18,120 of 30,876 pairs, including both interchromosomal, and intrachromosomal combinations, showed significant LD (based on D′, P < 0.05). Of these 18,120 pairs, 1530 were intrachromosomal pairs of SSR loci, which corresponds to a rate of 8.4%. The extent of LD was assessed for each of the seven groups (Table 3). For the average D′, Group 1 showed the lowest LD (0.399), whereas Group 5 exhibited the highest LD (0.468).

We used the D′ value corresponding to the intrachromosomal SSR loci and the genetic distance in each subpopulation to draw the attenuation map. Figure 4 shows that the D′ values decayed as the genetic distance (cM) increased. A regression analysis between the D′ value and the genetic distance of syntenic marker pairs revealed that the seven subpopulation genomes fitted the equation y = b lnx + c. The minimum distances of LD decay for Groups 1, 2, 3, 4, 5, 6, and 7 were 19.5 cM, 23.2 cM, 24.1 cM, 12.9 cM, 22.2 cM, 22.8 cM, and 21.7 cM, respectively. Thus, among the seven subpopulations, Group 4 exhibited the highest decay velocity with the shortest decay distance, whereas Group 3 demonstrated the lowest decay velocity.

The average F_ST among the seven subpopulations was 0.672, with the F_ST for each locus ranging from 0.197 for RM333_Chr 10 to 0.984 for RM6361_Chr 2. A results of a pairwise comparison based on the values of F_ST can be interpreted as the standardized population distances between two subpopulations. The pairwise F_ST values obtained in the present study ranged from 0.162 (between Groups 2 and 6) to 0.832 (between Groups 2 and 7), with an average value of 0.679 (Table 4). The AMOVA results indicated that 67.88% of the total genetic variation occurred between the subpopulations, whereas 32.12% occurred within the subpopulations (Table S3). These results indicate the existence of a high degree of genetic differentiation across the seven subpopulations.

A marker-trait association analysis based on an MLM revealed that six markers located on chromosomes 1, 2, 4, and 6 were associated with the STL (Table 5). The phenotypic variation explained (PVE) ranged from 2.9 to 12.5%. RM5753_Chr 6, which resides at 124.4 cM, had the highest PVE values for the STL, namely 12.5% in 2013, and 8.9% in 2014 (Table 5). Two markers distributed on chromosomes 4 and 6 were associated with the SBPL (Table 5), and of these, RM136_Chr 6 had the highest average PVE of 6.5% over the 2 years. Seven markers distributed on chromosomes 1, 2, 4, 6, and 11 were associated with the SNBPL (Table 5), and their PVE values ranged from 3.1 to 12.9%.

Thirteen markers distributed on chromosomes 1, 2, 3, 4, 6, 7, 8, 10, and 11 were associated with the GL (Table 5), and their PVE values ranged from 2.2 to 22.9%. RM6327_Chr 11, which resides at 1.7 cM, had the maximum PVE for the GL, specifically 21.8% in 2013, and 22.9% in 2014 (Table 5).

Using the entire set of accessions, we identified 28 markers associated with stigma traits and grain length, including six markers associated with the STL, two related to the SBPL, seven associated with the SNBPL, and 13 associated with the GL. Nineteen of the 28 associations were in regions where the QTL associated with the given trait had been identified by Redoña and Mackill (1998), Tan et al. (2000), Uga et al. (2003a), Aluko et al. (2004), Li et al. (2004), Agrama et al. (2007), Bai et al. (2010), Uga et al. (2010), Wang et al. (2011, 2012), Dang et al. (2015) (http://www.gramene.org/); these QTLs are listed in Table S4. Nine loci identified in this study, including one for the STL, one for the SBPL, three for the SNBPL, and four for the GL, are novel (Table 5).

The marker RM280 was co-associated with the STL, SBPL and SNBPL. Five markers, i.e., RM5389, RM450, RM7598, RM280, and RM5753, were co-associated with the STL and SNBPL. Four markers, namely RM5389, RM450, RM7598, and RM280, were co-associated with the GL and STL and with the GL and SNBPL. Two markers, RM280 and RM136, were co-associated with the GL and SBPL (Table 5).

The alleles with positive effects identified in this study were considered elite alleles for all four traits measured. A summary of the elite alleles and their typical carrier materials is shown in Table S5. The total numbers of elite alleles detected across the entire population for the STL, SBPL, and SNBPL were 11, 5, and 12, respectively. The allele RM450-135 bp showed the greatest phenotypic effect (0.235 mm) on the STL, and its typical carrier accession was Yuedao 32. The allele RM136-200 bp showed the greatest phenotypic effect (0.111 mm) on the SBPL, and its typical carrier accession was Yuexiangzhan. The allele RM450-135 bp showed its greatest phenotypic effect (0.233 mm) on the SNBPL, and its typical carrier accession was Yuedao 32.

A total of 30 elite alleles were detected for the GL. The allele RM136-200 bp exerted the greatest phenotypic effect (1.612 mm) on the GL, and its typical carrier accession was Yuexiangzhan (Table S5).

Correlations between the measured traits were observed, and the STL was significantly positively correlated with the SBPL and SNBPL. Furthermore, the GL was also significantly positively correlated with the STL, SBPL, and SNBPL. We identified one SSR marker co-associated with the STL and the SBPL, and the allele RM280-175 bp in this SSR marker locus simultaneously increased the phenotypic effect values of the STL and SBPL (Table 5 and Table S5). We identified five SSR markers co-associated with the STL and SNBPL, and the alleles RM5389-120, RM450-135, RM450-155, RM7598-115, RM280-175, RM5753-200, and RM5753-205 bp at these marker loci simultaneously increased the phenotypic effect values of the STL and SNBP (Table 5 and Table S5). We also identified four SSR markers co-associated with the STL and GL, and the alleles RM5389-120, RM450-135, RM450-155, RM7598-115, and RM280-175 bp at these marker loci simultaneously increased the phenotypic effect values of the STL and GL (Table 5 and Table S5). We further identified one SSR marker co-associated with the STL, SBPL, SNBPL, and GL, and the allele RM280-175 bp at this marker locus simultaneously increased the phenotypic effect values of the STL, SBPL, SNBPL, and GL (Table 5 and Table S5). These co-associated alleles have the correct sign with respect to trait correlations, and these data illustrate the genetic basis of trait correlations.

Based on the number of elite alleles that could be substituted into an individual plant and the expected phenotypic effects of the elite alleles that could be pyramided, the top five cross combinations for improving the STL, SBPL, SNBPL, and GL are proposed (Table S6). The elite alleles carried by the parents in excellent crosses and the corresponding phenotypic effects are listed in Table S7. Figure S5 shows the excellent parents in the superior cross for stigma traits. Certain accessions were found repeatedly in these proposed parental combinations (e.g., Yuedao 32 emerged four times in the combinations for the STL), indicating that these accessions possess unique elite alleles.

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.