Mining of co-expression genes in response to cold stress at maize (Zea mays L.) germination and sprouting stages by weighted gene co-expression networks analysis

Experimental lines and treatments

In this experiment, different degrees of cold treatment were applied at two maize germination stages. The maize inbred lines Zhongxi 091/O2 and Chang 7-2 were subjected to RNA sequencing, which carried out by the Maize Research Institute of the Jilin Academy of Agricultural Sciences in Jilin, China. The maize inbred lines Zhongxi 091/O2 and Chang 7-2 were found to be resistant (R) and susceptible (S) (Liu et al., 2021) to cold stress at the germination and sprouting stages in previous studies, respectively. The identification data of the two lines are detailed in Fig. S1.

Germination stage: 200 plump seeds were selected and disinfected with 0.1% HgCl₂ for 10 min. The seeds were soaked at 25 °C for 12 h and then subjected to low-temperature treatment at 10 °C. The seeds were placed on a disinfected culture dish covered with a moist sponge and filter paper at a density of 20 seeds per dish. The dishes were placed separately in a smart light incubator (10 °C, 24 h dark culture, and 60% humidity). Samples were collected at five time intervals (0, 6, 12, 24, and 36 h), and several neatly germinating seeds were randomly selected. The embryos were dissected using a dissecting needle and weighed to a total weight of 0.6 g. The experiment was repeated three times to ensure accuracy.

Sprouting stage: 500 plump seeds were selected and disinfected with 0.1% HgCl₂ for 10 min. The seeds were placed on a sterilized culture dish covered with a moist sponge and filter paper at a density of 20 seeds per dish. The dishes were placed in a smart light incubator for germination (25 °C, 24 h dark culture, and 60% humidity). When the sprouts were approximately 4 cm long, 100 seedlings of the same growth size were selected and subjected to low-temperature treatment at 2 °C. The samples were collected at five intervals (0, 6, 12, 24, and 36 h) and several seeds with the same sprout lengths were randomly selected. The seedlings were cut into cubes and weighed for a total weight of 0.6 g. The experiment was repeated three times to ensure accuracy.

The total RNA of maize at the germination and sprouting stages was extracted using the polysaccharide polyphenol plant total RNA extraction kit (Exp.01032020, stored at RT) of the Shanghai Pudi Biotechnology Co., Ltd. After purification, the ends were repaired, base A was added to the 3′ end, and then the sequencing adapter sequence was added. The target size fragment after gel electrophoresis was recovered, the library required for sequencing was amplified by PCR, and the sequencing was completed using an Illumina HiSeq™ 2000 sequencer.

Clustering analysis and functional enrichment of module genes

Gene co-expression network analysis was performed using the WGCNA package in R software (https://cran.r-project.org/web/packages/WGCNA/index.html). To ensure compliance with scale-free network distributions, for WGCNA, the appropriate weighting coefficient β must be selected. The value of coefficient β was estimated using the pick Soft Threshold function in the WGCNA package. First, the β = 1–30 was set to calculate the corresponding correlation coefficient and mean gene connectivity. The selection criterion of β is to ensure that the square of the correlation coefficient is as close as possible to 0.8, while also preserving a sufficient level of gene connectivity. The dynamic tree cut method was used to identify the co-expression module. The Automatic Network Builder Block-wise Module was used to build the network. Here, the minimum module size was 30, merge Cut Height = 0.25 merged the modules with a similarity of 0.75, and other parameters were set by default. To obtain the biological functions and signaling pathways involved in each module, genes in each module were subjected to GO analysis (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) and KEGG analysis. The R package was then used to plot and display the results. The threshold for the GO term and the KEGG pathway was set as p < 0.05.

Weighted gene co-expression network analysis and gene network visualization

The gene co-expression network obtained from WGCNA was processed, and Cytoscape V3.6.1 was used to screen hub genes and visualize their network. Each node in the network represents a gene, and the edges represent the relationships between genes.

Real-time quantitative reverse transcription PCR (qRT-PCR) verification

Ten cold tolerance-related genes were randomly selected for qRT-PCR analysis. Total RNA was extracted from maize during the germination and sprouting stages at low temperatures using a Polysaccharide Polyphenol Plant Total RNA Extraction Kit (Shanghai Pudi Biotechnology Company, Shanghai, China). Reverse transcription was performed using the TaKaRa Reverse Transcriptase M-MLV (RNase H) reverse transcription kit (TaKaRa Bio Inc., Shiga, Japan). qRT-PCR was performed on a Quant Studio six and seven real-time fluorescence quantitative PCR instrument (Applied Biosystems, Waltham, MA, USA) using SYBR^® Green Real-time PCR Master Mix (Toyobo, Osaka, Japan), following the manufacturer’s instructions. The reaction conditions were as follows: 95 °C for 5 min, 95 °C for 15 s, 58 °C for 30 s, and 70 °C for 30 s, for a total of 35 cycles. Specific primers were designed using Primer 5.0 software and are listed in Table 1 with corn UBI as the internal reference gene. Gene expression levels were calculated using the 2^−ΔΔCT method, where ΔCT refers to the difference of the CT value between the target and the internal reference gene (Livak & Schmittgen, 2001).

RNA sequencing

A total of 230.4 and 228.3 Gb raw readings were obtained from 60 samples (2 varieties × 5 time intervals × 3 replicates) at the germination and sprouting stages, respectively. After filtering, 1,544,300,100 and 1,531,483,946 clean reads were obtained, with a Q30 percentage above 95.00%. The percentage of GC ranged from 54.04% to 57.35%, indicating that all experimental samples were highly reliable in terms of collection and sequencing results. The statistical power of this experimental design, calculated using RNA-Seq Power, was 0.84. To verify the effectiveness of cold stress treatment, the reported marker genes involved in cold stress were evaluated for their expression patterns before and after cold stress treatment, and the results showed that with increasing cold stress duration, the expression of ZmG6PDH5 (LOC100383421) increased and that of ZmDREB1A (LOC103647602) decreased, reaching the lowest value at 24 h. The results of the RNA-seq data were similar to previously reported results, indicating that the cold stress treatment carried out in this study was effective and further verified the reliability of the RNA-seq data (Fig. 1).

Principal coordinates analysis of cold stress

Principal coordinates analysis (PCA) was performed to explore the correlation between sample processing methods. In the PCA of cold stress at the two stages, the first principal coordinate primarily explained the difference in stage expression and had a large contribution, accounting for 62.96% of the total variance. This indicates that the germination and sprouting stage of stress should be analyzed separately when using WGCNA to construct a gene expression network to avoid stage differences in the stress treatment (Fig. 2A).

In the PCA of cold stress of the lines at the germination stage, the first principal coordinate primarily explained the differences in line expression, accounting for 59.96% of the total variance. In the PCA of cold stress of the lines at the sprouting stage, the first principal coordinate primarily explained the difference in line expression and had a large contribution, accounting for 58.67% of the total variance. This indicates that when using WGCNA to construct a gene expression network, the R and S lines under stress at the germination and sprouting stage should be analyzed separately to avoid line differences in the stress treatment (Figs. 2B and 2C).

Differentially expressed genes analysis

There were 546 genes that responded at the germination and sprouting stages of both lines (Fig. 3A). A total of 14,384, 7,301, 10,295 and 10,799 DEGs were observed in Resistant inbred lines at the Germination stage (GR), Susceptible inbred lines at the Germination stage (GS), Resistant inbred lines at the Sprouting stage (SR), and Susceptible inbred lines at the Sprouting stage (SS) within 36 h of low temperature stress, respectively. The number of DEGs for the R and S lines both exhibited an increasing trend with the extension of stress duration, except for SS. SS decreased after reaching its maximum value at 24 h, while the other three treatments reached their maximum at 36 h (Fig. 3B). Compared to the sensitive line Chang 7-2 (18,100), a greater number of genes were involved in Zhongxi 091/O2 (24,679). Notably, the number of up-regulated genes (14,182) was significantly higher than that of down-regulated genes (10,497) in Zhongxi 091/O2. Conversely, there was no significant difference in the number of up-regulated (9,449) and down-regulated genes (8,651) in Chang 7-2 (Figs. 3C and 3D).

WGCNA analysis

The weight values were calculated, and the appropriate soft threshold was selected to ensure that the co-expression network followed a scale-free network distribution. The results showed that the optimal soft threshold values were β = 20 at the germination stage and β = 15 at the sprouting stage. This was further used to construct a co-expression network (Fig. S2).

First, cluster analysis was carried out according to the expression of genes, and genes with a high degree of clustering were divided into a module. Then, the dynamic cleavage method was used to identify the co-expression modules. Different modules are indicated by different colors, and grey represents the genes not assigned to any module (Fig. 4). Six modules were obtained at the germination stage, among which the turquoise module contained the greatest number of genes (3,635) and the red contained the least number (43). At the sprouting stage, nine modules were formed. The turquoise module had 2,668 genes, whereas the magenta and pink had only 202 genes each (Figs. 4A and 4B).

At the germination stage, the correlation of the blue module with stress duration changed from significantly negative to significantly positive and was the most significant in Zhongxi 091/O2 at 36 h (r = 0.45, p = 0.01). A significantly positive correlation was observed between the yellow module and stress duration in Zhongxi 091/O2, which reached a maximum at 36 h (r = 0.55, p = 0.002), whereas Chang 7-2 showed no significant correlation with stress duration (Fig. 4C).

At the sprouting stage, significantly positive correlations between the turquoise and magenta modules and stress duration in Zhongxi 091/O2 were observed, which reached a maximum at 36 h (r = 0.66, p < 0.001; r = 0.68, p < 0.001). Finally, we determined that the blue and yellow modules at the germination stage and the turquoise and magenta modules at the sprouting stage were maize low-temperature-specific modules (Fig. 4D).

Enrichment analysis of specific modules under low-temperature stress

To investigate gene function within specific modules under low-temperature stress, the blue and yellow modules at the germination stage and the turquoise and magenta modules at the sprouting stage were annotated. All four modules were enriched in the response to hormone (GO:0009725) and oxygen-containing compounds (GO:1901700).

In the blue module, the genes were mainly enriched in the following: biological processes such as plant organ development (GO:0099402), developmental processes involved in reproduction (GO:0003006), and post-embryonic development (GO:0009791); molecular functions such as transcription factor activity, sequence-specific DNA binding (GO:0003700), nucleic acid binding transcription factor activity (GO:0001071), DNA binding (GO:0003677), and 2-alkenal reductase (NAD(P)⁺) activity (GO:0032440); and cellular components such as DNA packaging complex (GO:0044815), organelle membrane (GO:0031090), and plasma membrane (GO:0005886) (Fig. 5A). In the yellow module, the genes were mainly enriched in biological processes such as response to abscisic acid (GO:0009737), response to cold (GO:0009409), and hormone-mediated signaling pathways (GO:0009755) (Fig. 5B).

In the turquoise module, the genes were mainly enriched in the following: biological processes such as the abscisic acid activated signaling pathway (GO:0009738), response to jasmonic acid (GO:0009753), and response to cold (GO:0009409); molecular functions such as transcription factor activity, sequence-specific DNA binding (GO:0003700), nucleic acid-binding transcription factor activity (GO:0001071), and protein serine/threonine kinase activity (GO:0004674); and cellular components such as plasma membrane (GO:0005886), cell periphery (GO:0071944), and cell part (GO:0044464) (Fig. 5C).

In the magenta module, the genes were mainly enriched in biological processes such as regulation of cellular macromolecule biosynthetic process (GO:2000112), regulation of transcription, DNA-templated (GO:0006355) and regulation of nucleic acid-templated transcription (GO:1903506), and molecular functions such as nucleic acid binding transcription factor activity (GO:0001071), sequence-specific DNA binding (GO:0003700), zinc ion binding (GO:0008270), and transition metal ion binding (GO:0046914) (Fig. 5D).

KEGG enrichment of the four specific modules was carried out. The blue module included plant hormone signal transduction, quorum sensing, and carbon fixation in photosynthetic organisms. The yellow module included plant hormone signal transduction, mitogen-activated protein kinase (MAPK) signaling pathway-plant, and vitamin B6 metabolism. The turquoise module included the MAPK signaling pathway, plant hormone signal transduction, and phenylpropanoid biosynthesis. The magenta module included plant hormone signal transduction, plant MAPK signaling pathway, and galactose metabolism (Fig. 6).

Gene network visualization

Hub genes are typically those that have a high connectivity within a module. We selected the top 50 genes with the highest connectivity in the blue, yellow, magenta, and turquoise modules for mapping and selected the top 10 genes in terms of connectivity as hub genes. We visualized the hub genes and their associated genes and constructed a gene interaction network (Figs. S3–S5). Forty key genes including eight transcription factors tolerant to cold were predicted. Seven genes have been proven to be involved in the response to low-temperature stress (Table S1).

RT verification

qRT-PCR validation was performed for 10 genes (Fig. 7). We found that these genes have differences in their expression profiles in different lines and at different stages. The results showed that 10 candidate genes have a trend from low to high expression. At the germination stage, the expression of nine genes (except Zm00001d028401), was higher in the R line than in the S line, Zm00001d010863 and Zm00001d028401 had the highest expression at 24 h of stress, Zm00001d034099 and Zm00001d045557 had the highest expression at 12 h of stress, and the rest of the six genes all had the highest expression at 36 h. At the sprouting stage, the expression levels of Zm00001d002533, Zm00001d010863, Zm00001d022125, and Zm00001d028401 were significantly higher in the R than in the S lines. The expression of Zm00001d034099 was significantly higher in the S than in the R lines, and the remaining genes were not significantly different in terms of expression between the two lines. Compared with the germination stage, the gene expression at the sprouting stage showed no significant change across time. Overall, 10 genes exhibited similar expression patterns in qRT-PCR and RNA-seq experiments.