Transcriptomic insights into grain size development in naked barley (Hordeum vulgare L. var. nudum Hook. f): based on weighted gene co-expression network analysis

Abstract

Background: This study investigated the molecular mechanisms underlying grain size variation between two distinct naked barley varieties using comprehensive phenotypic and transcriptomic (RNA-Seq) analyses.

Methods: In this study, we employed a comparative transcriptomics approach to analyze two naked barley varieties: the large-grained Shenglibai and the small-grained Lalu Qingke. Our investigation focused on three critical developmental periods of grain growth (early, mid, and late grain-filling periods). By integrating longitudinal three-dimensional phenotypic data with temporal expression profiles and applying weighted gene co-expression network analysis (WGCNA), we successfully identified gene modules that co-vary with morphological expansion.

Results: Phenotypic assessments revealed that grains underwent rapid expansion during the filling period, with significant differences in grain width (GW) and thickness (GT) across all three developmental periods. In contrast, grain length (GL) remained relatively consistent by the end of the filling period. Transcriptome sequencing identified a peak in differentially expressed genes (DEGs) during the mid-filling period, indicating that the regulation of grain size development is most active in the early and mid-filling phases. WGCNA identified a blue module strongly correlated with grain size, which was significantly enriched in key metabolic pathways, including starch and sucrose metabolism. Further analysis identified seven hub genes, among which HvENO1 exhibited pronounced upregulation in large-grain varieties during the mid-to-late filling periods, closely aligning with the observed phenotypic traits. Real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR) validation confirmed the period-specific and variety-specific expression patterns of these genes, further supporting the potential of these genes as targets for improving grain size in breeding.

Keywords: Grain size; Qingke; WGCNA; qRT-PCR.

Figure 1. Grain size dynamics of two varieties during development.

(A) Dynamic changes in grain morphology of Shenglibai and Lalu Qingke at different grain-filling periods. (B) Changes in grain length of the two varieties across three developmental periods. (C) Changes in grain width of the two varieties across three developmental periods. (D) Changes in grain thickness of the two varieties across three developmental periods.

Figure 2. DEGs statistics.

The horizontal coordinates represent the different groups of differential comparisons and the vertical coordinates represent the corresponding number of up-and down-regulated genes.

Figure 3. WGCNA analysis.

(A) Scale-free topology fitting curve and mean connectivity curve. The horizontal axis represents the soft-thresholding power β, the vertical axis represents the scale-free topology model fit (R²) corresponding to the adjacency matrix transformed under the given β value. A higher R² indicates that the network better approximates a scale-free distribution (left figure); The horizontal axis represents the soft-thresholding power β; the vertical axis represents the average connectivity (k/degree) of each node in the network corresponding to the adjacency matrix transformed under the given β value (right figure). (B) Module clustering dendrogram. Branches represent individual modules. The closer two module branches are, the more similar the modules, which can serve as a basis for module merging. The vertical axis represents the clustering distance. (C) Module classification tree. Genes are divided into modules based on their expression trends. Each branch represents a gene, and each color represents a module. Genes colored in gray are those not assigned to any specific module. (D) Module-trait relationship heatmap. The heatmap illustrates the correlation between modules and specific traits. The horizontal axis represents different traits, and the vertical axis represents different modules. The numbers in the leftmost column indicate the number of genes in each module. The data on the right side of each cell shows the correlation coefficient between the module and the trait, along with the significance P-value (in parentheses). The larger the absolute value, the stronger the correlation. Orange indicates a negative correlation, while purple indicates a positive correlation.

Figure 4. GO and KEGG analysis.

(A) GO enrichment analysis. The vertical axis represents GO terms, and the horizontal axis represents the Rich factor. A larger Rich factor indicates a higher degree of enrichment. The size of the points represents the number of genes in the corresponding GO term, and the color of the points corresponds to different ranges of adjusted P-values (padjust). (B) KEGG enrichment analysis. The vertical axis represents pathway names, and the horizontal axis represents the ratio of Rich factor to the number of annotated genes (background number). A larger Rich factor indicates a higher degree of enrichment. The size of the points represents the number of genes in the corresponding pathway, while the color of the points corresponds to different ranges of adjusted P-values (padjust).

Figure 5. Co-expression network analysis of genes in the blue module.

Figure 6. Expression analysis of seven hub genes.

(A) Expression analysis of seven hub genes in transcriptomic data. S: Shenglibai; L: Lalu Qingke. 7810: HORVU.MOREX.r3.7HG0637810; 9040: HORVU.MOREX.r3.1HG0049040; 7190: HORVU.MOREX.r3.7HG0647190; 5510: HORVU.MOREX.r3.7HG0675510; 2030: HORVU.MOREX.r3.2HG0142030; 1620: HORVU.MOREX.r3.2HG0191620; 1510: HORVU.MOREX.r3.2HG0131510. (B) HORVU.MOREX.r3.1HG0049040. (C) HORVU.MOREX.r3.7HG0637810. (D) HORVU.MOREX.r3.7HG0647190. (E) HORVU.MOREX.r3.7HG0675510. (F) HORVU.MOREX.r3.2HG0142030. (G) HORVU.MOREX.r3.2HG0191620. (H) HORVU.MOREX.r3.2HG0131510. Asterisks (****) indicate a statistically significant difference (p < 0.0001).