Identification of regulatory networks and hub genes controlling soybean seed set and size using RNA sequencing analysis

Abstract

To understand the gene expression networks controlling soybean seed set and size, transcriptome analyses were performed in three early seed developmental stages, using two genotypes with contrasting seed size. The two-dimensional data set provides a comprehensive and systems-level view on dynamic gene expression networks underpinning soybean seed set and subsequent development. Using pairwise comparisons and weighted gene coexpression network analyses, we identified modules of coexpressed genes and hub genes for each module. Of particular importance are the discoveries of specific modules for the large seed size variety and for seed developmental stages. A large number of candidate regulators for seed size, including those involved in hormonal signaling pathways and transcription factors, were transiently and specifically induced in the early developmental stages. The soybean homologs of a brassinosteroid signaling receptor kinase, a brassinosteroid-signaling kinase, were identified as hub genes operating in the seed coat network in the early seed maturation stage. Overexpression of a candidate seed size regulatory gene, GmCYP78A5, in transgenic soybean resulted in increased seed size and seed weight. Together, these analyses identified a large number of potential key regulators controlling soybean seed set, seed size, and, consequently, yield potential, thereby providing new insights into the molecular networks underlying soybean seed development.

Keywords: Seed development; seed number; seed size; soybean; transcriptomic analysis..

Fig. 1.

Stages of soybean seed development and sampling times of two soybean varieties. (A) Seed developmental stages in soybean varieties V1 and V2. Pod (a) and seed (b) at 5–7 days after fertilization (DAF) (Stage S1), seed at 10–14 DAF (Stage S2; c), seeds at 20–24 DAF (Stage S3; d), seed at maturation stage (e) and desiccation stage (f). Bars=1 mm. (B) Cross-section of the seeds of V1 and V2 at stage S1. V1 showed enlarged seed size partly due to increased cell size in the outer and inner integument compared with V2. Bars=200μm. (C) Seed lengths and widths, and widths of the outer and inner integument of V1 and V2 at stage S1. Values are presented as mean±SE. Three biological replicates were used for measurements. ** Significant difference (P<0.01), Student’s t-test. (This figure is available in colour at JXB online.)

Fig. 2.

Global gene expression profiling in early seed development stages. (A) Numbers of detected transcripts in each sample. (B) Principal component analysis of the RNA-Seq data. (C, D) Venn diagrams of differentially expressed transcripts among the four tissues of (C) V1 and (D) V2. (This figure is available in colour at JXB online.)

Fig. 3.

Numbers of DEGs in each developmental stage and tissues. (A) Numbers of DEGs between V1 and V2 at each developmental stage. (B) Numbers of DEGs in each developmental stage of V1 and V2. DEGs were filtered according to FPKM >5, FDR <0.01, log₂ fold change >1, or log₂ fold change <–1. Numbers in red or black indicate the number of up-regulated or down-regulated genes when the sample in red is compared with the sample in black, respectively. (This figure is available in colour at JXB online.)

Fig. 4.

Heatmap comparison of the DEGs at each seed developmental stage. Functional categories of significantly over-represented DEGs in S1 (A), S2 (B), S3-1 (C), and S3-2 (D). (This figure is available in colour at JXB online.)

Fig. 5.

WGCNA of genes in V1 and V2 at each seed developmental stage. (A) Hierarchical cluster tree showing coexpression modules identified by WGCNA. Each leaf in the tree represents one gene. The major tree branches constitute 12 modules, labeled with different colors. (B) Module–sample association. Each row corresponds to a module, labeled with a color as in (A). The number of genes in each module is indicated on the left. The number of transcription factors in each module is indicated by the number in parentheses. Each column corresponds to a specific tissue and replicate (R). The color of each cell at the row–--column intersection indicates the correlation coefficient between the module and the tissue type. A high degree of correlation between a specific module and the tissue type is indicated by red underline of the module name. (This figure is available in colour at JXB online.)

Fig. 6.

Coexpression network analysis of stage-specific modules. (A), (C), (E), and (G) Heatmaps showing genes in each module that were significantly over-represented in S1, S2, S3, and V1, respectively. (B), (D), (F), and (H) The correlation networks in the module corresponding to (A), (C), (E), and (G), respectively. Candidate hub genes are shown in red. (This figure is available in colour at JXB online.)

Fig. 7.

Overexpression of GmCYP78A5 (Glyma.05G019200) increased seed size in transgenic soybean. (A) RT-PCR analysis of the expression level of Glyma.05G019200 in Williams 82 and three Glyma.05G019200 overexpression transgenic lines. RNA was extracted from developing seeds at the S3 stage. (B) Seeds and (C) pods of Glyma.05G019200 overexpression transgenic lines and their non-transgenic segregants. (D) Seed weight and (E) seed dimensions of Glyma.05G019200 overexpression transgenic lines and their non-transgenic segregants. Values are reported as means±SE (n=3). * Significant difference (P<0.05, Student’s t-test). (This figure is available in colour at JXB online.)