Soybean gene co-expression network analysis identifies two co-regulated gene modules associated with nodule formation and development

Abstract

Gene co-expression network analysis is an efficient systems biology approach for the discovery of novel gene functions and trait-associated gene modules. To identify clusters of functionally related genes involved in soybean nodule formation and development, we performed a weighted gene co-expression network analysis. Two nodule-specific modules (NSM-1 and NSM-2, containing 304 and 203 genes, respectively) were identified. The NSM-1 gene promoters were significantly enriched in cis-binding elements for ERF, MYB, and C2H2-type zinc transcription factors, whereas NSM-2 gene promoters were enriched in cis-binding elements for TCP, bZIP, and bHLH transcription factors, suggesting a role of these regulatory factors in the transcriptional activation of nodule co-expressed genes. The co-expressed gene modules included genes with potential novel roles in nodulation, including those involved in xylem development, transmembrane transport, the ethylene signalling pathway, cytoskeleton organization, cytokinesis and regulation of the cell cycle, regulation of meristem initiation and growth, transcriptional regulation, DNA methylation, and histone modifications. Functional analysis of two co-expressed genes using TILLING mutants provided novel insight into the involvement of unsaturated fatty acid biosynthesis and folate metabolism in nodule formation and development. The identified gene co-expression modules provide valuable resources for further functional genomics studies to dissect the genetic basis of nodule formation and development in soybean.

Keywords: TILLING; gene co-expression network; nodule formation and development; soybean.

FIGURE 1

Identification of two nodule‐specific modules using weighted gene co‐expression network analysis (WGCNA). (a) Network heatmap plot showing the eigengenes of 34 modules. The gene dendrograms and module colours are displayed to the left and top of the heatmap. NSM‐1 and NSM‐2 are highlighted in dark grey and violet, respectively. (b and c) Gene networks showing co‐expression events between 304 genes in NSM‐1 (b) and 203 genes in NSM‐2 (c) The gene networks were visualized using Cytoscape v. 3.7.2. (d and e) Cis‐binding motifs enriched in the promoters of genes in NSM‐1 (d) or NSM‐2 (e).

FIGURE 2

Heatmap demonstrating the expression patterns of NSM‐1 and NSM‐2 genes across 11 tissues. RNA‐Seq count data were normalized as RPKM and used to construct gene expression heatmaps of NSM‐1 (a) and NSM‐2 (b) genes. Colour scales on the top indicate log₂ of normalized count data.

FIGURE 3

Gene ontology (GO) term enrichment analysis of nodule specific modules. (a and b) Ontology enrichment analyses of the NSM‐1 (a) and NSM‐2 (b) genes were conducted using the agriGO database with Fisher's exact test and Bonferroni multitest adjustment. Adjusted p value <0.01 was set as statistically significant. Redundant GO terms were grouped based on semantic similarity to other terms in the Uniprot database and visualized as a scatter plot using the REViGO Web server (Supek et al., 2011). The size of the bubbles reflects the abundance of the GO terms. The colour scales represent log₁₀ of the GO term adjusted p values.

FIGURE 4

TILLING mutations of Glyma.19G121600 and Glyma.14G121400 impact nodule formation and development. (a) Description of TILLING mutations of Glyma.19G121600 and Glyma.14G121400, and the corresponding amino acid substitutions. (b and c) Impact of TILLING mutations in Glyma.19G121600 and Glyma.14G121400 on nodule formation and development. Nodule number (a) and size (b) were determined 22 days after inoculation with Bradyrhizobium diazoefficiens (strain USDA110). Data were obtained from at least 10 plants from two independent experiments, which are indicated in blue and red. Dots in the graphs represent values from individual plants. Asterisks indicate statistically significant differences from the wild‐type cv. Forrest as determined by analysis of variance (p < 0.05). (d) Cross‐sections of 22‐day‐old nodules formed on wild‐type Forrest and TILLING mutants of Glyma.19G121600 and Glyma.14G121400. The pink colour in the wild‐type nodule reflects a high concentration of leghaemoglobin associated with active nitrogen fixation. In contrast, the brown colour in the nodules formed in TILLING mutants of Glyma.14G121400 and Glyma.19G121600 reflects a low concentration of leghaemoglobin associated with reduced nitrogen fixation.