All the 159 indica rice accessions were planted, with a randomized complete block design, in Wuhan in 2021. This experiment was replicated two times in different fields, namely environments 1 and 2. Each material was planted in one plot with 10 seedlings, row spacing was 16.7 cm × 20 cm, and one line empty between cells. At yellow ripening stage, GW (mm), GL (mm), and TGW (g) for each accession in MC and RR were measured for three times, and their averages were regarded as their trait phenotypes. GL in MC is abbreviated as GL-MC, and it is true for other traits.

The minimum, maximum, mean, standard deviation (SD), kurtosis, skewness (S_k), and coefficient of variation (CV), along with broad-sense heritability ( H B 2 ) , for all the above traits were calculated by R software lme4 v1.1.28. The H B 2 for each trait was calculated by H B 2 = σ g 2 σ g 2 + σ g e 2 / l + σ e 2 / r l × 100 % , where σ g 2 is genetic variance, σ e 2 is residual variance, σ g e 2 is the variance of genotype-by-environment interaction, l is the number of environments, and r is the number of replicates. The analysis of variance (ANOVA) for phenotypic data was conducted using the R function aov. Normal distribution test for phenotypic data was conducted using the R function shapiro.test.

The genotypic data of the 159 rice accessions used in this study consisted of two parts. The genotypic datasets of 134 accessions were derived from RiceVarMap database (http://ricevarmap.ncpgr.cn/), and the DNAs of leaves were extracted to conduct 1K Genobaits to verify their authenticity. The genotypic datasets of twenty-five modern breeding cultivars were obtained by double-terminal sequencing with coverage of approximate 10× based on illumina’s Hiseq 4000 technology sequencing platform at Novogene Technology Company. Then, extract the common SNPs from the genotype dataset of 134 public database accessions and 25 modern breeding cultivars to obtain new genotypic datasets with 2,019,008 SNPs. The software plink v1.90 was used to filter all the 2,019,008 SNPs based on minimum allele frequencies (MAFs) < 0.05 and all variants with missing call rates > 10%, where sliding window distance, step length, and R² were set as 1000 kb, 1, and 0.3, respectively. As a result, a total of 2,017,495 SNPs were used in subsequent GWAS.

All the 2,017,495 SNPs were used to conduct linkage disequilibrium (LD) analysis using popLDdecay (https://github.com/BGIShenzhen/PopLDdecay). The LD decay was determined by plotting the r² values against the genetic distance of a pair of loci (kb) for each chromosome. G-matrix and cluster analysis for all the 159 accessions were performed using the 2,017,495 SNPs by R package sommer v4.2.0 and amap v0.8.19, respectively. Principal component analysis (PCA) was analyzed using R function prcomp, and the first two principal components were plotted using the R package ggplot2 v3.3.6. ADMIXTURE v1.3.0 (http://dalexander.github.io/admixture) was used to determine population structure (Alexander et al., 2009), where the number of subgroups (K) was set from 1 to 10, and the K value corresponding to the minimum CV error is the most likely subgroup number.

A total of 2,017,495 SNPs of 159 rice accessions were used to associate with GW, GL, and TGW in two environments in MC and RR using the 3VmrMLM method and its IIIVmrMLM software (https://github.com/YuanmingZhang65/IIIVmrMLM; Li et al., 2022a; Li et al., 2022b). All parameters were set as default values. Population structure adopts the first three principal components. The K matrix was calculated using the IIIVmrMLM software. The probability threshold was set at 0.05/m = 2.48e-08 for significant QTNs and QEIs, where m was the number of markers. To reduce the loss of important candidate genes, some insignificant QTNs and QEIs with LOD score ≥ 3.0 were regarded as suggested QTNs and QEIs (Li et al., 2022a; Li et al., 2022b).

Candidate genes for grain size traits were mined based on the below steps. First, all the genes were found in the 200 kb regions of upstream and downstream around each significant QTN without previously reported gene, because the LD decay distance was 150 kb using popLDdecay. Then, RNA-seq datasets of Zhenshan 97 and Minghui 63 at endosperm 7, 14, and 21 days after pollination (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19024) were used to conduct differential expression analysis using NCBI (https://www.ncbi.nlm.nih.gov) GEO2R online tool, and the thresholds of significant difference were set as p-value < 0.05 and | Log2FC | > 1. Finally, all the differentially expressed genes (DEGs) were further analyzed by GO annotation using AgBase (https://agbase.arizona.edu), and the significant E-value was set as 10e-50. If biological process is related to the reported molecular mechanisms of grain size, the DEGs in the biological process were regarded as candidate genes.

The software plink v1.90 was used to extract all the significant SNP information (P < 0.05) after single marker genome scanning within one candidate gene and its upstream 2 kb, R v4.1.3 was used to calculate its haplotypes of the candidate gene, and the 159 rice accessions were grouped based on these haplotypes. Thus, ANOVA was performed using R function aov to test the significance of the QTN-associated trait across these haplotypes at a 5% probability level.

The averages plus standard deviations of GW-MC, GL-MC, TGW-MC, GW-RR, GL-RR, and TGW-RR in 159 rice accessions in two environments were 2.44 ± 0.33 ~ 2.47 ± 0.34 (mm), 8.41 ± 0.85 ~ 8.42 ± 0.84 (mm), 23.92 ± 3.01 ~ 24.03 ± 2.98 (g), 2.38 ± 0.29 ~ 2.43 ± 0.28 (mm), 8.00 ± 0.79 ~ 8.05 ± 0.81 (mm), 21.74 ± 2.84 ~ 22.32 ± 3.09 (g), and their coefficients of variation (CV) were 13.57 ~ 13.70, 9.91 ~ 10.12, 12.40 ~ 12.58, 11.50 ~ 12.17, 9.91 ~ 10.08, and 13.04 ~ 13.83 (%), respectively, having large phenotypic variations ( Supplementary Table S1 ). The analysis of variance was conducted and the results were listed in Supplementary Table S2 . As a result, genotypes, environments, and their interactions for all the three traits in MC and RR were significant at the 0.05 probability level ( Supplementary Table S2 ), and the H B 2 of GW, GL, and TGW ranged from 96.39% to 99.07% in MC and from 90.21% to 98.37% in RR, indicating large genetic variations ( Supplementary Table S1 ). In addition, main crop had higher trait averages than ratoon rice, especially for GL and TGW ( Figure 1 ). The phenotypes of TGW-MC and GW-RR in two environments, GL-RR in environment 2, and TGW-RR in environment 1 were found to obey normal distribution, while GW-MC and GL-MC in two environments, and GL-RR and TGW-RR in environment 1 were found to approximately obey normal distribution ( Figure 1 ; Supplementary Table S1 ).

To determine the LD decay distance, LD decay analysis was performed using all the 2,017,495 SNP markers. The r² gradually decreases with the increase of distance. When it drops to half of the maximum value, the corresponding distance is regarded as the average distance of LD decay. In this study the LD decay distance was 150 kb, when r² dropped to half of its maximum value (r² = 0.3) ( Figure 2C ).

The number of sub-populations was determined by principal component analysis (PCA), population structure analysis, and cluster analysis. The results were showed in Figure 2 . In PCA, the first two principal components separated all the 159 accessions into three subgroups: indica I, indica II, and indicia Intermediate ( Figure 2A ). In population structure analysis via the ADMIXTURE software, cross-validation (CV) error is the lowest when the number of subgroups is three ( Figures 2D, E ), which is consistent with that in cluster analysis ( Figure 2B ). Thus, the first three principal components were used in genome-wide association studies.

Known genes were searched within 200 kb upstream and downstream regions of QTNs and QEIs. Among the QTNs and QEIs, 3, 4, and 7 known genes were found in main crop to be truly associated with the above three traits, respectively, and 2, 4, and 4 known genes were found in ratoon rice to be truly associated with the above three traits, respectively. Among these known genes, 4 were simultaneously found in main crop and ratoon rice, and 3 were found across multiple traits. 3, 5, and 9 known genes were found to be around significant QTNs and QEIs for the above three traits, respectively, and 0, 1, and 2 known genes were found to be around suggested QTNs for the above three traits, respectively ( Tables 1 , 2 ).

Around the above QTNs, some known genes were simultaneously mined in single-environment analysis and two-environment joint analysis. GW5 was identified to be associated with GW-MC and GW-RR in two single-environment analyses and two-environment joint analysis, VLN2 was identified to be associated with GW-MC in two-environment joint analysis and GW-RR in the first environment analysis ( Figures 3A–C ; Supplementary Figure 1 ), GS3 was found to be associated with GL-MC and GL-RR in two single-environment analysis and two-environment joint analysis, and GW5 was found to be associated with GL-MC in the second environment analysis and two-environment joint analysis ( Supplementary Figures 2A–C, 3 ). For TGW, all known genes were separately detected in a single-environment analysis or two-environment joint analysis ( Supplementary Figures 4A–C, 5 ).

Around the above QEIs, two known genes, OsACOT and GW6a, for GL-RR were mined ( Table 2 ). Among the two known genes, OsACOT was found to be interacted with environments. In detail, its expression level under moderate soil drying treatment was higher than that under well-watered control (Teng et al., 2022) ( Table 2 ).

In the joint analysis of the MC and RR datasets, 1, 2, and 1 known genes were found to be around significant QTNs and to be truly associated with GW, GL, and TGW, respectively ( Supplementary Tables S35 ; Figures 3D, E ; Supplementary Figures 2, 4D, E ). Among these known genes, most of them were consistent with the above known genes, such as GW5, GS3, and qTGW3, but PGL2 was found only in the MC and RR joint analysis.

Around other QTNs without known genes, all the genes within 200 kb upstream and downstream regions were used to conduct differential expression analysis. All the differential expression genes (DEGs) were used to conduct gene annotation analysis. In gene annotation analysis, the significant biological processes were mainly included the below categories: cytokinin, abscisic acid and other plant hormone metabolism (e.g., Os02g0197600, Os02g0621300, Os02g0626100, and Os02g0178800), protein ubiquitination (Os07g0166800), sucrose starch metabolism (Os04g0169100, Os12g0112500, and Os08g0205900), protein phosphorylation (Os02g0126400, and Os03g0717700), and endosperm development (Os02g0538000, Os12g0277500, and Os01g0280500) ( Supplementary Tables S21, S22 ), which are highly consistent with the previously reported regulatory pathways in Zuo and Li (2014); Cai et al. (2018), and Li et al. (2019). These genes were regarded as candidate genes. As a result, there were 10, 7 and 1 candidate genes for GW, GL and TGW in main crop, respectively, and 8, 4, and 3 candidate genes for GW, GL and TGW in ratoon rice, respectively ( Supplementary Tables S21, S22 ).

For candidate genes for GW, Os02g0126400 and Os03g0717700 were predicted to be related to protein phosphorylation, Os02g0178800 and Os03g0592500 were predicted to respond to abscisic acid, Os02g0197600 was found to be related to cytokinin, Os07g0166800 was found to be related to the process of protein ubiquitination, Os02g0538000 and Os12g0277500 were found to be associated with embryonic development at the end of seed dormancy, and Os08g0205900 and Os04g0169100 were predicted to be related to sucrose metabolism and starch synthesis metabolism, respectively.

For candidate genes for GL, Os02g0197600, Os02g0621300, Os02g0626100, Os04g0514800, and Os03g0108600 were predicted to respond to cytokinin, abscisic acid, gibberellin, auxin, and ethylene, respectively, Os01g0280500 and Os02g0538000 were found to affect embryonic development, and Os04g0169100 and Os12g0112500 were predicted to be related to starch synthesis. In Wuriyanghan et al. (2009), Os04g0169100 was reported to affect the ethylene sensitivity of seeds to substantially enhance TGW of mutant. Here there is one issue pending, that is, whether the TGW increase is caused by the GL increase.

For candidate genes for TGW, Os02g0621300 was predicted to be related to response to abscisic acid, Os03g0411500 was predicted to be related to photosynthesis, Os03g0607400 was predicted to be related to positive regulation of unidimensional cell growth, and Os05g0445900 was predicted to participate in DNA demethylation, which is consistent with that in Zhou et al. (2021), in detail, Os05g0445900 encodes rice DNA glycosylase, and the mutation of DNG701 can lead to embryo retardation or abortion of part seeds (Zhou et al., 2021).

For the DEG Os03g0737000 (P=7.77E-03, log2FC=-1.29) around the QEI of chr3-30340995 for GL-MC, Os03g0737000 was predicted to be related to “response to salt stress” ( Table 3 ). In the future, new experiments are necessary to explore these novel gene-trait and GEI-trait associations.

In the joint analysis of the MC and RR datasets, there were 22 candidate genes around QTNs to be responsible for the above three traits, but only two were consistent with the above 25 candidate genes ( Supplementary Table S36 ). The significant biological processes of these candidate genes were mainly included the below categories: plant hormone metabolic pathway (e.g., Os02g0126400, Os05g0563400, and Os12g0288000), protein phosphorylation (Os03g0838100, Os08g0200500, and Os05g0514200), embryo development (Os02g0538000 and Os08g0428100), and protein ubiquitination (Os09g0294300 and Os12g0111500). In addition, we also mined two additional DEGs for TGW around QEIs, among which Os06g0154200 was predicted to be related to “positive regulation of response to water deprivation” and Os11g0600900 was predicted to be related to “response to light intensity” ( Table 3 ).

To further verify the reliability of candidate genes, we conducted haplotype analysis. As a result, 12 of the above 25 candidate genes had significant differences among the phenotypes of the traits corresponding to the haplotypes of each gene ( Figure 4 ). Among the 12 significant candidate genes, 7, 3, and 3 were found to be associated with GW, GL, and TGW, respectively, of which there are 6, 3, and 1 significant candidate genes in main crop and 2, 0, and 2 significant candidate genes in ratoon rice ( Figure 4 ). Os08g0205900 for GW was mined in both MC and RR, and Os02g0621300 was found in both GL and TGW ( Figure 4 ). It should be noted that 8 of 22 candidate genes, which were mined in the MC and RR joint analysis, were significant in haplotype analysis, and the eight genes were different from the 12 significant candidate genes in the above haplotype analysis ( Figure 4 ).

Ratoon rice is a new mode of rice planting, which can effectively save costs and increase benefits (Shen et al., 2021). Although it is generally accepted that RR has higher quality than MC, there are still many controversies in the studies on grain size related traits (Alizadeh and Habibi, 2016; Huang et al., 2020; Yuan et al., 2022b).

In this study, 64, 71, and 64 QTNs and 0, 1 and, 2 QEIs in MC, and 72, 63 and 56 QTNs and 2, 2, and 1 QEIs in RR were identified to be associated with GW, GL and TGW, respectively. Among these QTNs, 4 known genes were commonly detected in main crop and ratoon rice to be truly associated with grain size related traits, including GW5 and VLN2 for GW, and GS3 and GW7 for GL ( Table 1 ). Some known genes were found only in main crop or ratoon rice, such as qTGW3, SPL33, and OsSPL14 were detected only in main crop, and OsVPE3, GS2, ETR2, and UPA1 were found only in the ratoon rice ( Table 1 ). Among all the candidate genes, one was commonly found in main crop and ratoon rice, and 12 were detected only in main crop or ratoon rice ( Figure 4 ). No common QEIs were detected between main crop and ratoon rice.

Based on the above results, main crop can detect more known genes (14), candidate genes (10) and candidate GEIs (1) than ratoon rice (8, 4, and 0). Although some known and candidate genes can be commonly found in main crop and ratoon rice, there are still some specific candidate genes in main crop or ratoon rice. In the independent and joint analyses of the MC and RR datasets, most candidate genes and candidate GEIs were different across the two analyses. This indicated that more known and candidate genes and GEIs can be identified while the datasets in main crop and ratoon rice are simultaneously or jointly analyzed.

The candidate genes were mined by expression and GO annotation analysis, and further validated through haplotype analysis. In this study we identified five new and key candidate genes that were predicted to be closely related to the three traits, among which 3 were mined to be around QTNs in the MC and RR joint analysis ( Table 3 ), the evidence was as below.

Os02g0626100 for GL-MC, and Os02g0538000 for GW-MC were differentially expressed. In GO annotation analysis, the two genes were annotated as “response to gibberellin”, and “embryo development ending in seed dormancy”, respectively, in which these biological processes are highly consistent with metabolic pathways of important grain size traits in rice (Li et al., 2018; Li et al., 2020; Jiang et al., 2022). In haplotype analysis, significant GL/GW differences were observed across 2 and 5 haplotypes from 1 and 11 significant SNPs within the two genes and their 2 kb upstream. Thus, the two genes may be important candidate genes for GL/GW.

In the same way, Os08g0379300 for TGW, Os09g0294300 for GW, and Os12g0557800 for GL were found to be DEGs around the QTNs in the MC and RR joint analysis. In GO annotation analysis, the three gene were predicted to be related to “starch catabolic process”, “protein ubiquitination”, and “response to abscisic acid”, respectively, which have been confirmed to be important regulatory pathways of rice grain size (Li et al., 2008; Choi et al., 2018; Gao et al., 2021). In haplotype analysis, significant TGW/GW/GL difference was observed across 7, 6, and 16 haplotypes from 6, 6, and 16 significant SNPs within the genes and their 2 kb upstream. Thus, the three genes may be important candidate genes for TGW/GW/GL.

In addition, two candidate genes have been reported to be related to rice seed development. In Wuriyanghan et al. (2009), Os04g0169100, identified for GW-RR in this study, significantly increased TGW of mutants by increasing the sensitivity of seeds to ethylene. In Zhou et al. (2021), the mutant of Os05g0445900, identified for TGW-MC in this study, participated in DNA methylation process causing the endosperm of some seeds to be stunted or aborted.

Around the QEIs, OsACOT for GL-RR has been confirmed to be differentially expressed under two soil moisture treatments (Teng et al., 2022; Table 2 ).

Around QEIs in the independent analysis of MC or RR, Os03g0737000 for GL-MC was found to be differentially expressed, and its biological process in GO annotation was predicted to be related to salt stress. We speculate that Os03g0737000 may be affected by environmental factors, such as different salt treatments. Around QEIs detected in the MC and RR joint analysis, Os06g0154200 and Os11g0600900 for TGW were found to be differentially expressed, and the biological processes in their GO annotations were predicted to be related to water deprivation and light intensity, respectively. Thus, we speculate that Os06g0154200 may be affected by the moisture content of the environment and Os11g0600900 may be affected by the intensity of external light. The molecular functions of above three candidate GEIs need to be verified by subsequent molecular biology experiments.

To investigate the effect of interval length on mining known genes, two types of interval lengths were compared. One was 200 kb upstream and downstream regions of QTNs and QEIs, which was determined based on LD decay distance, while another was 1000 kb for QTNs and 1500 kb for QEIs. The results are listed in Table 1 and Supplementary Table S37 . As a result, 3, 4, 7, 2, 2, and 4 known genes around QTNs and 0, 0, 0, 0, 2, and 0 known genes around QEIs were found to be located on their corresponding 200 kb upstream and downstream regions and to be truly associated with GW-MC, GL-MC, TGW-MC, GW-RR, GL-RR, and TGW-RR, respectively, while 6, 7, 21, 7, 6, and 16 known genes for the above six traits were found to be located on 1000 kb upstream and downstream regions of QTNs, and 0, 0, 1, 2, 2, and 0 known genes for the above six traits were found to be located on 1500 kb upstream and downstream regions of QEIs. This indicates that large intervals can find more known genes. Thus, it is very important to determine a suitable interval length in mining known genes.

To investigate the effect of population structure on genome-wide association studies, we compared the results from evolutionary population structure (Liu et al., 2020), Q matrix, and PCA in this study. As a result, 323, 283, and 393 QTNs, 9, 6, 8 QEIs, and 11, 12, and 20 known genes were identified from evolutionary population, Q matrix and PCA, respectively ( Supplementary Tables S23–S28 ). Clearly, the PCA result is the best, followed by evolutionary population, and the worst is the Q matrix result in this study. Thus, population structure is an important parameter in genome-wide association study.