Previous studies have confirmed that GS3 is the major QTL that controls grain length and weight in rice, with an SNP (C→A) mutation in the second exon leading to early termination of translation and a subsequent long-grain phenotype (Fan et al., 2006; Fan et al., 2009; Mao et al., 2010). In NIL-GS3, the ORF length of GS3 was 696 bp and encoded 231 amino acids. The domain prediction identified four functional domains: the N-terminal organ size regulator (OSR) domain, a transmembrane domain (TM), a C-terminal TNFR/NGFR family cysteine-rich structure, and a VWFC structure (Figure 1A). In NIL-gs3, however, the GS3 locus mutated from C to A at 165 bp, resulting in a premature termination codon (TGC→TGA) and a sequence with only 54 amino acids (Figure 1A). The GS3 loci in NIL-GS3 and NIL-gs3 were the same as those reported in Zhenshan 97 and Minghui 63, respectively (Mao et al., 2010).

When comparing the grain length, length-to-width ratio, and thousand-grain weight, NIL-gs3 increased by 14.32, 4.69, and 16.17%, respectively, relative to NIL-GS3 (p < 0.05). Additionally, the grain width of NIL-gs3 increased by 9.17%, but we did not note any significant difference in grain thickness. Therefore, the gs3 locus contained in NIL-gs3 had a positive effect on rice grain size (Figures 1B–F; Supplementary Table S1).

To examine the cytological basis of grain size differences in the GS3 isogenic lines, we examined the morphology of the cells in the middle portion of the glumes of mature NIL-GS3 and NIL-gs3 seeds and compared them using SEM. The results showed that the glume cell length was significantly longer in NIL-gs3 than NIL-GS3 (6.50%, in the grain-length orientation), whereas the cell width decreased significantly (9.77%). Subsequently, the glume epidermal cells in NIL-gs3 appeared to be more slender, and the length-to-width ratio increased. These findings indicated that the GS3 gene regulates rice grain shape by regulating glume cell morphology (Figure 1G–J; Supplementary Table S1).

To comprehensively analyze the transcriptome changes in rice grain development following a GS3 mutation, we examined three representative stages (S:small panicle size, M:medium panicle size, and L:large panicle size) from the beginning of panicle to the beginning of the heading in the NILs. We obtained a total of 114.78 G base pairs, with an average of 6.38 G per sample. We then successfully aligned the reads to the reference genome IRGSP-1.0 using Hisat2. The results showed that the percentage of unique alignments to the genome was between 89.49 and 91.96%, with an average of 90.84% (Supplementary Table S2).

After RABT assembly, we obtained 85,411 transcripts from 44,668 genes, with 5,936 new transcripts identified in the intergenic regions (Supplementary Table S3; Figure 2). We then predicted coding potentials using CPC2 and PLEK and identified 18,343 transcripts (15,136 genes) as lncRNA candidates. Among them, 4,143 were located in the intergenic regions and 155 were located in introns. These two types of lncRNAs can be collectively called lincRNAs. Additionally, 1,255 lncRNA transcripts overlapped with gene exons, presumably because of alternative splicing (Supplementary Table S3). Quantitative analysis of the expression profiles showed that lncRNA expression, especially lincRNA expression, was lower than that of coding genes, consistent with previous studies (Lv et al., 2019).

During rice growth and development, dynamic changes in gene expression occur in different tissues and organs during different developmental stages. At the booting stage, many genes associated with cell proliferation, cell expansion, transcription factors, and hormone regulation modulate inflorescence and panicle morphogenesis (Ke et al., 2016; Wang Y et al., 2019; Zha et al., 2019). When identifying DEGs at the different developmental stages in NIL-GS3, we identified 5,974 DEGs at the S stage relative to the M stage, 2,395 at the L stage relative to M, and 9,400 at the L stage relative to the S. Among these three groups of DEGs, a total of 742 genes were commonly significantly differently expressed. We also observed this phenomenon in the three stages of spikelet development in NIL-gs3 (Supplementary Table S4). These findings clearly showed the presence of gene expression stratification at different developmental stages.

To further examine expressional changes during different developmental stages, we compared differential expression between the NIL groups at each stage (Figure 3; Supplementary Table S4). During the S stage, we identified 887 DEGs (167 lncRNAs), 434 upregulated (85 lncRNAs), and 453 downregulated (82 lncRNAs) in NIL-gs3 relative to the NIL-GS3. For the M stage, we identified 1,768 DEGs (318 lncRNAs), 852 upregulated (129 lncRNAs) and 916 downregulated (189 lncRNAs); For the L stage, we identified 1,478 DEGs (280 lncRNAs), 670 upregulated (151 lncRNAs), and 808 down-regulated (129 lncRNAs). Furthermore, the number of DEGs increased with the growth of the panicle, thus suggesting that GS3 affected the expression of more genes in the mid-late stages of panicle development (Figure 3A). When comparing the young NIL panicles at the same stage, 215 genes (66 lncRNAs) were commonly differentially expressed across the three stages. Furthermore, hierarchical clustering analysis also showed a similar expression pattern between the groups (Figures 2, 3C). The number of unique differential genes in the S, M, and L stages were 442 (90 LncRNAs), 1,091 (126 LncRNAs), and 925 (120 LncRNAs), respectively (Figure 3B). To confirm the RNA-Seq findings, we randomly selected ten representatives of differentially expressed transcripts, including two lncRNAs, and eight coding genes. The qRT-PCR results showed that three of the transcripts were significantly upregulated, and five were significantly down-regulated among S, M, and L periods. (Supplementary Figure S1). All changes were consistent with expected results from the RNA-Seq data, thus indicating a high degree of reliability.

The c-means fuzzy clustering algorithm was employed for clustering gene expression profiles in all developmental stages (Kumar and Futschik, 2007). The fuzzy clustering algorithm could assign a membership value for each gene expression profile of each cluster using different fuzziness values. As a result, nine distinct clusters of temporal patterns representing NILs and development stages were generated (Figure 4A; Supplementary Table S5). Among them, cluster 7, cluster 4, and cluster 6 were identified, which represented different development stages (S, M, and L) and significantly upregulated in the corresponding stage of NIL-gs3, whereas cluster 8, cluster 5, and cluster 1 are upregulated and showed stage-specific in three stages of NIL-GS3.

GO and KEGG enrichment analysis of these 9 clusters showed many enriched biological processes and pathways corresponding to the relevant stages (Figure 4B, Supplementary Table S6). For example, genes associated with nitrogen compound metabolic process, lipid metabolic process, carbohydrate metabolic process, and nucleobase-containing compound metabolic process are highly enriched in cluster 7, 4, and 6, corresponding to three stages (S, M, and L). Those genes associated with the organic acid metabolic process, phosphorus metabolic process, and oxoacid metabolic process were significantly enriched in cluster 7 (S stage) and cluster 4 (M stage). Three gamma subunits, GS3 (Os03g0407400), GGC2 (Os08g0456600), and DEP1 (Os09g0441900), were predominantly expressed at the early stage of panicle development and then rapidly declined (Supplementary Table S3). Interestingly, small GTPase mediated signal transduction, intracellular protein transport, and signal transmission was significantly enriched in Cluster4 (M stage). Meanwhile, starch and sucrose metabolism, Glycolysis, and phenylpropanoid biosynthesis were also enriched explicitly in the M stage. Most genes such as OsRab6A1, OsRop4, OsWD40, and OsRho associated with the G-protein signaling pathway were significantly upregulated in the M stage (Figure 5A). Besides, genes associated with plant hormone signal transduction such as OsGSK1, OsRLCK98, OsXTH9, OsGH3.3, and OsRR10 were specifically enriched and upregulated in the L stage (Figure 5B). These results suggested that G-protein and hormone signaling pathways may show a synergistic interplay in regulating panicle development in rice.

We collected a large grain Japonica variety (TD70) and a small grain Indica variety (Kasalath) (Zhang et al., 2015). We then constructed NILs (NIL-gs3 and NIL-GS3) using TD70 as a donor and Kasalath as the background by conventional six-generation backcrossing. The NILs were planted in a greenhouse at the Jiangsu Academy of Agricultural Sciences (Nanjing, China). We used the small (S, 1–3 cm), medium (M, 4–6 cm), and large (L, 8–10 cm) panicle size for specimen collection. We also collected three biological repeats for each of these periods for each of the two NILs (n = 18). Spikelets were snap-frozen in liquid nitrogen following separation and stored at −80°C. The growth and development of these spikelets were highly consistent. We obtained grain length, width, and thickness measurements after harvesting and measured the thousand-grain weight after harvesting and drying at 37°C for 1 week. We observed epidermal cell morphology within the mature grain glumes with a scanning electron microscope (SEM; EVO-LS10, Carl Zeiss Microscopy GmbH, Oberkochen, Germany).

Total RNA was extracted from the spikelet tissue using a TRIzol® RNA Purification Kit (Invitrogen, Carlsbad, CA, United States). We determined RNA purity, concentration, and integrity using a NanoDrop 2000 Spectrophotometer (ThermoFisher Scientific, Waltham, MA, United States), Qubit® 4.0 Spectrophotometer (Life Technologies, Carlsbad, CA, United States), and Agilent 2,100 Bioanalyzer® (Agilent Technologies, Palo Alto, CA, United States), respectively. The mRNA was then enriched using Oligo (dT) magnetic beads, and double-stranded cDNA (ds-cDNA) synthesis was performed. The ds-cDNAs were then purified using Agencourt AMPure XP beads (Beckmann Coulter, Pasadena, CA, United States), followed by end repair, polyadenylation, and sequencing connector. Finally, we used AMPure XP beads to screen for fragment size and polymerase chain reaction (PCR) amplification before constructing the sequencing library. The library was sequenced on an Illumina NovaSeq 6,000 sequencer in the PE150 mode (paired-end 150 bp).

We cleaned the raw sequences by removing adapter sequences, and low-quality reads with Q20 (a Phred quality score of 20 represents an error rate of 1 in 100) using the Fastp program (Chen et al., 2018), and eliminated any sequences ≤50 bp. The cleaned reads were then mapped to Nipponbare reference genome (https://rapdb.dna.affrc.go.jp, IRGSP-1.0 version) with the Hisat2 program (Kawahara et al., 2013; Pertea et al., 2016). We extracted unique reads (MAQ ≥20) for analysis using SAM tools v1.6 (Li et al., 2009). We updated the genome annotation with a reference annotation-based transcript (RABT) method, based on IRGSP-1.0 GTF annotation file (https://rapdb.dna.affrc.go.jp), using the StringTie program (Pertea et al., 2016). We extracted updated transcripts using the Gffreads program and performed a potential coding prediction using the CPC2 (Kang et al., 2017) and PLEK programs (Li et al., 2014). We considered transcripts with a length greater than 200 bp and consistently identified them as non-coding using both algorithms to be lncRNAs.

We determined the expression abundance of each transcript based on the fragments per kilobase of transcript per million fragments mapped (FPKM) value and identified unexpressed or false positive transcripts using the largest FPKM value < 1 in multiple samples (Trapnell et al., 2010). We identified differentially expressed genes (DEGs) at different stages of panicle development using DESeq2 software (Love et al., 2014) with a Log2 (fold-change) ≥1 and q-value <0.05 required. We then visualized DEG commonalities between samples using a Venn plot.

Expressed genes from three development stages in NILs were grouped into different clusters using the Mfuzz package in R with fuzzy c-means (Kumar and Futschik, 2007). The different clusters represented the different expression patterns during rice panicle development. The gene expression values were standardized to have a mean value of zero and one standard deviation for each gene profile. The transformed expression values were then clustered using the fuzzy c-means clustering algorithm by the Mfuzz package.

The ClusterProfile package (Yu et al., 2012) was used for enrichment analysis and visualization using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, with an FDR <0.05 required for significance. Heat maps of the gene expression profiles were generated using R packages.

To validate the reliability of the RNA-Seq analysis, we performed quantitative real-time PCR (qRT-PCR) using ten randomly selected genes. The primers were designed using Primer Premier v.5.0 software (Premier Biosoft International, Palo Alto, CA, United States). We synthesized the first-strand cDNA using a TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransGen Biotech, Beijing, China). We conducted the qRT-PCR on an ABI StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, United States) using SYBR Premix Ex Taq (TaKaRa, Dalian, China). We used Actin1 (Os03g0718100) as an internal standard to calculate the relative expression level and determined expression values based on the relative quantitative method (Livak and Schmittgen, 2001).