Illumina paired-end sequencing, assembly and annotation of maize leaf transcriptomes
Four cDNA libraries were prepared using mRNA isolated from the leaves of both water-deprived and well-watered maize plants of drought-tolerant hybrid ND476 at four developmental stages. The libraries were denoted ND1_Con (NDCa, NDCb, NDCc) and ND1_Tre (NDDa, NDDb, NDDc) (the leaves of control and treatment maize at V12 stage); ND2_Con (NDCA, NDCB, NDCC) and ND2_Tre (NDDA, NDDB, NDDC) (the leaves of control and treatment maize at VT stage); ND3_Con (NDC1, NDC2, NDC3) and ND3_Tre (NDD1, NDD2, NDD3) (the leaves of control and treatment maize at R2 stage); ND4_Con (NDC4, NDC5, NDC6) and ND4_Tre (NDD4, NDD5, NDD6) (the leaves of control and treatment maize at R4 stage). To identify genes responsive to drought stress in maize leaves at various growth stages, global gene expression profiling was performed by Illumina RNA sequencing of these libraries.
Resultantly, a total of 125.26 million raw reads were obtained. The raw sequencing data has been deposited in the NCBI Sequence Read Archive (SRA, Accession: SPR212360). After adaptors and low-quality reads were filtered out, 124.16 million clean reads were obtained, ranging from 23,029,648 to 72,636,578 for each sample. The clean reads were used for further analysis. Meanwhile, 20,432,633 (88.72%) to 64,399,057 (88.66%) clean reads were mapped onto unique positions on the maize reference genome (ZmB73_Ref-Gen_v4) (Table S1). The Q30 base percentage and GC percentages exceeded 94.46% and 54.6%, respectively (Table S1).
Subsequently, for functional annotation of the assembled transcriptome sequences, all the sequences were mapped onto public genome database with an E-value threshold of 1e-5. We annotated 179,093 (93.1%) and 137,849 (71.66%) genes in NCBI Nr database and the Swiss-Port protein database, respectively (Table S2). Based on KEGG analysis, only 82,758 genes were successfully annotated, accounting for 43.02% of the total number. In addition, 169,691 (88.21%) and 133,525 (69.41%) genes were annotated using COG and GO databases, respectively (Table S2).
Additionally, the similarities or differences of the twenty-four samples were analyzed using principal component analysis (PCA). The PCA results of all the samples showed clear separation between treatment and control samples at different stages (Figure S1). To measure the gene expression levels for the three replicates for each sample, Pearson correlation coefficients between samples were calculated. The results shown by way of heatmap revealed that each R2 between two samples was higher than 90% except for one comparison (NDC4_vs_NDDb) (Table S3). These results indicated the overall reproducibility and quality of the assay, which met the demands for further analysis.
Identification of DEGs in response to drought
In order to reveal the transcriptional responses of maize leaves to water-stressed conditions, we compared the genes identified under water-sufficient and water-deficit conditions at four different growth stages. Gene expression levels were calculated and normalized to the RPKM values. Based on this analysis, a total of 3451 DEGs were identified at four various maize growth stages. We obtained most numbers of DEGs at the V12 stage, including 1203 up-regulated and 1200 down-regulated. Meanwhile, 352 up-regulated and 298 down-regulated DEGs were identified at the VT stage. Similarly, we fished out 397 DEGs (95 up-regulated and 302 down-regulated) and 313 DEGs (112 up-regulated and 201 down-regulated) at the R2 and R4 stages, respectively (Fig. 1A).
The number of DEGs showing overlaps and specific responses under drought stress in different growth stages is visualized in Fig. 2B. A large number of DEGs were period-specific; there were 2164, 483, 307 and 198 DEGs, respectively at the four different stages. However, a limited number of common DEGs was detected. Area I represent 117 DEGs shared between V12 and VT stages after drought treatment, that is, the common DEGs identified in the two vegetative stages. Area II represents 15 DEGs shared between R2 and R4 stages after drought treatment, that is, the common DEGs identified in the two reproductive stages. There were 180 drought responsive DEGs identified in the vegetative stages (V12 or VT) and also differentially expressed at the reproductive stages (R2 or R4) after drought treatment; that is, the DEGs identified at both vegetative and reproductive stages (Area III, Fig. 1B).
To further understand the gene expressions between different stages, we performed the hierarchical clustering analysis of the identified DEGs (Fig. 1C-H). The DEGs in different stages showed different expression patterns. These results showed that there were different mechanisms of maize drought stress responses at various growth phases.
Functional annotation of DEGs using MapMan
All DEGs of four growth stages were assigned to MapMan functional categories. The DEGs were grouped into 35 BINs with putative functions (Figure S2). We found out that 720, 183, 147, and 72 DEGs of the V12, VT, R2, and R4 stages, respectively, were not assigned to any functional group (BIN 35) due to lack of annotation information (Figure S2). The DEGs of V12 stage were mainly annotated to cell wall, lipid metabolism, photosynthesis (PS), protein synthesis and degradation, abiotic stress, secondary metabolites biosynthesis and hormone metabolism (Fig. 2A). The highly enriched categories of the VT stage DEGs included lipid metabolism, amino acid metabolism and hormone metabolism (Fig. 2A). Meanwhile, the enriched categories of the R2 stage DEGs included lipid metabolism, protein degradation, and secondary metabolites biosynthesis, whereas the R4 stage DEGs related to transport, PS, hormone metabolism, secondary metabolites biosynthesis and RNA transcriptional regulation were highly enriched (Fig. 2A).
An overview analysis of DEGs was generated with MapMan tool and the drought-inducible regulated genes were classified into different regulatory processes (Table 1). The DEGs involved in PS showed up-regulated expression under drought stress conditions in all the stages, except the VT stage. A great deal of DEGs related to transport was altered in expression among the four growth stages. Moreover, a substantial number of protein kinases including serine/threonine-protein kinases, leucine-rich repeat receptor-like proteins, receptor-like kinases, phospholipases, and protein phosphatases were observed to be mostly up-regulated at the vegetative stages, whist showing down-regulation at the reproductive stages. Several plant hormones, which functions as regulatory compounds, were identified to be responsive to drought stress including as abscisic acid (ABA), auxin, brassinosteroids (BR), gibberellic acid (GA), salicyclic acid (SA), and ethylene. Majority of DEGs related to auxin and SA was down-regulated, whilst the DEGs involved in other hormones showed an increased expression under drought stress. Additionally, a great number of differentially expressed TFs were identified, such as C2H2, bHLH, HB, MYB domain, MYB-related, and WRKY domain (Fig. 2B). In the current study, more increased abundance TFs were identified in the vegetative stages compared to the reproductive stages.
We also found DEGs associated with other functional processes. For instance, the DEGs involved in maintaining redox homeostasis by a series of enzymatic compounds including thioredoxin (TRX), glutathione S-transferases (GST), and peroxidase (POD) which played a major role in protecting maize from oxidative damage. Additionally, several DEGs were annotated to stress defense. HSPs mainly showed decreased abundance, whist LEA and pathogenesis-related proteins (PRPs) had increased abundance under drought. The identification of such a great number of regulatory DEGs showed that there were multiple signaling mediators and intricate pathways in response to drought stress. Subsequently, we also obtained 10, 3, and 2 DEGs of the V12, VT, and R4 stages that were annotated to 'response to drought/salt' (BinCode: 20.2.3) (Table 1). Among them, responsive to dehydration 22 (RD22, at5g25610) mediated by ABA was annotated in the V12 and VT stages and was involved in response to desiccation. Drought-responsive family protein at3g05700 and ERD (early-responsive to dehydration stress, at4g22120) family protein may play a role in maize response to water-stress in the V12 stage. In addition, gene encoding AOC (allene oxide cyclase, at3g25780), which is involved in jasmonic acid biosynthesis, is suggested to play a functional role in maize response to drought at the V12 stage. Taken together, these differentially expressed genes were speculated to be the vital cogs in maize drought stress tolerance, and hence aroused our keen interest for further discussion.
Table 1
Classification of drought-response regulated DEGs into different categories according to MapMan annotation
Category | V12 stage | VT stage | R2 stage | R4 stage |
Up- | Down- | Up- | Down- | Up- | Down- | Up- | Down- |
Photosynthesis | | | | | | | | |
Photosystem II | 3 | 1 | 0 | 1 | 1 | 0 | 3 | 0 |
Photosystem I | 2 | 0 | 0 | 3 | 0 | 0 | 1 | 0 |
Electron carrier(ox/red) | 3 | 0 | 1 | 2 | 1 | 0 | 0 | 0 |
calvin cycle | 0 | 0 | 0 | 2 | 1 | 1 | 0 | 0 |
Transport | | | | | | | | |
Transport sugars | 8 | 3 | 1 | 2 | 0 | 1 | 0 | 3 |
Transport amino acids | 11 | 7 | 1 | 4 | 2 | 0 | 1 | 4 |
Transport ammonium | 0 | 2 | 0 | 1 | 1 | 0 | 0 | 1 |
Transport phosphate, nitrate,sulphate | 2 | 4 | 1 | 2 | 0 | 0 | 2 | 4 |
Transport metal | 3 | 2 | 3 | 0 | 0 | 3 | 0 | 0 |
Transport peptides and oligopeptides | 11 | 6 | 1 | 1 | 0 | 1 | 0 | 5 |
Transport potassium | 8 | 4 | 2 | 0 | 0 | 3 | 1 | 1 |
ABC transporters | 6 | 4 | 2 | 1 | 0 | 2 | 1 | 3 |
Transport misc | 9 | 6 | 2 | 5 | 0 | 4 | 1 | 2 |
Protein kinases and phosphatases | | | | | | | | |
Serine/threonine-protein kinase | 11 | 4 | 2 | 1 | 0 | 5 | 0 | 1 |
LRR recetor-like protein | 27 | 6 | 10 | 1 | 1 | 7 | 1 | 2 |
Receptor-like kinase | 28 | 20 | 5 | 7 | 0 | 5 | 0 | 10 |
Phospholipase | 8 | 2 | 5 | 0 | 0 | 1 | 2 | 0 |
Calmodulin | 10 | 6 | 0 | 2 | 0 | 4 | 0 | 1 |
Mitogen-activated protein kinase(kinase)(kinase) | 2 | 1 | 1 | 2 | 0 | 0 | 0 | 0 |
Protein phosphatase | 2 | 6 | 3 | 0 | 0 | 0 | 4 | 0 |
Protein phosphatase 2C | 4 | 0 | 2 | 0 | 0 | 0 | 2 | 2 |
Plant Hormones | | | | | | | | |
Abscisic acid | 14 | 2 | 8 | 1 | 1 | 3 | 9 | 3 |
Auxin | 4 | 7 | 3 | 3 | 1 | 2 | 1 | 2 |
Brassinosteroid | 11 | 4 | 4 | 0 | 0 | 0 | 1 | 1 |
Jasmonic acid | 3 | 3 | 4 | 2 | 1 | 1 | 0 | 3 |
Salicylic acid | 1 | 4 | 1 | 2 | 0 | 0 | 0 | 2 |
Gibberellic acid | 3 | 0 | 0 | 0 | 1 | 1 | 0 | 1 |
Ethylene | 6 | 4 | 2 | 2 | 2 | 2 | 0 | 0 |
Transcription factor family | | | | | | | | |
Basic Helix-Loop-Helix | 9 | 3 | 4 | 2 | 2 | 1 | 1 | 1 |
C2H2 zinc finger | 6 | 3 | 0 | 1 | 0 | 2 | 0 | 1 |
Homeobox transcription factor | 11 | 4 | 2 | 0 | 0 | 1 | 2 | 2 |
MYB domain transcription factor | 4 | 3 | 3 | 3 | 0 | 2 | 0 | 2 |
MYB-related transcription factor | 3 | 1 | 0 | 1 | 0 | 1 | 0 | 1 |
NAC domain transcription factor | 2 | 4 | 3 | 2 | 0 | 0 | 0 | 1 |
WRKY domain transcription factor | 3 | 3 | 1 | 1 | 0 | 0 | 0 | 1 |
bZIP transcription factor | 3 | 0 | 1 | 2 | 0 | 3 | 1 | 1 |
G2-like transcription factor | 2 | 0 | 1 | 3 | 0 | 1 | 0 | 2 |
Other transcription factor families | 42 | 54 | 22 | 16 | 7 | 11 | 11 | 3 |
DEGs related to detoxification | | | | | | | | |
Thioredoxin | 2 | 1 | 0 | 1 | 0 | 0 | 2 | 0 |
Glutathione S-transferases | 6 | 3 | 1 | 4 | 0 | 1 | 0 | 2 |
Peroxidase | 8 | 6 | 1 | 3 | 1 | 4 | 0 | 2 |
Ascorbate and glutathione | 5 | 5 | 1 | 0 | 1 | 0 | 1 | 1 |
Glutaredoxins | 5 | 1 | 0 | 1 | 0 | 0 | 1 | 0 |
DEGs involoved in defense | | | | | | | | |
Heat shock proteins | 4 | 20 | 8 | 7 | 0 | 1 | 0 | 0 |
Late embryogenesis proteins | 6 | 0 | 2 | 1 | 0 | 0 | 0 | 0 |
Pathogenesis-related proteins | 5 | 1 | 4 | 0 | 0 | 2 | 0 | 0 |
DEGs response to abiotic | | | | | | | | |
Response to heat | 4 | 23 | 8 | 7 | 0 | 1 | 0 | 0 |
Response to drought/salt | 6 | 4 | 3 | 0 | 0 | 0 | 2 | 0 |
Response to cold | 3 | 1 | 0 | 2 | 0 | 2 | 0 | 2 |
Secondary metabolism | | | | | | | | |
Isoprenoids metabolism | 13 | 2 | 5 | 0 | 3 | 0 | 0 | 5 |
Phenylpropanoids metabolism and biosynthesis | 14 | 7 | 6 | 3 | 0 | 2 | 0 | 12 |
Flavonoids metabolism | 9 | 15 | 1 | 2 | 7 | 1 | 2 | 4 |
Sulfur-containing metabolism | 5 | 6 | 2 | 2 | 0 | 0 | 0 | 4 |
Note: V12, VT, R1, and R4 stages imply maize crop developmental stages as vegetative 12th leaf; vegetative terminal; reproductive first stage; and reproductive fourth stage, respectively. |
Co-expression network analysis of DEGs by WGCNA
To capture crucial shifts in gene networks in maize under water-stressed conditions, we further applied WGCNA approach to perform a network-level analysis of co-expression relationships among 3451 DEGs based on their expression patterns throughout the four growth stages. After filtering, a total of 2771 DEGs were divided into 12 modules (clusters) (designated M1-M12) comprising of 32 to 1155 highly co-expressed genes (Fig. 3, Figure S3, Table S4). GO enrichment analysis of each module by BiNGO highlighted vital biological processes represented by a series of co-expressed genes. Module M1 formed the largest cluster of 1155 DEGs enriched in functions related to metabolic processes (cellular amino acid, oxoacid, organic acid and small molecule) and response to temperature stimulus (Fig. 4,Table S5). We also observed that cluster M5 had 172 DEGs enriched in functions related to metabolism process (peptide, cellular amide and protein) and biosynthetic process (peptide, amide, cellular macromolecule) (Table S5). Meanwhile, module M11, comprising a cluster of 34 DEGs, had its DEGs annotated to biosynthetic process (cinnamic acid, phenylpropanoid and carboxylic acid) and metabolism process (cinnamic acid, aromatic amino acid and benzene-containing compound) (Table S5). Modules M2, M10 and M12 showed enrichment of GO terms related to photosynthesis (Table S5). Additionally, module M10 represents DEGs that showed high-expression specifically in the VT-stage, and were enriched in GO terms associated with ion homeostasis and transport (Table S5). The black module M7 included 52 DEGs involved in response to abiotic or environmental stimulus (Table S5). Further, the brown module M3 included DEGs related to ribosome and ribonucleoprotein complex biogenesis, gene expression, and RNA processing (Table S5). Modules M4, M6, M8 and M9 did not get enriched in any GO term. By combining our gene expression pattern and GO enrichment analysis results, we concluded that the DEGs mainly participated in metabolic and biosynthetic processes, photosynthesis, ion homeostasis, transport, and response to abiotic stimulus under drought stress conditions.
Identification of hub genes within network modules
There were some genes with extremely high connectivity with other genes, and these were designated as hub genes in each network module. Owing to their central location within the network clusters, the hub genes were considered to be vital components of the networks. Selecting only the top 10% of genes that showed high connectivity degree, a total of 277 DEGs were identified as hub genes (Fig. 5A-B), including 17 TFs represented from distinct families including WRKY, MYB-related, C2H2, MYB and NAC TFs (Table S6). Seven hub genes were also identified as crucial enzymes playing a key role during maize drought stress response (Table S6). Besides TFs and enzymes, there were two HSP90 genes that were also observed to respond to water-deficit stress conditions. Six hub-genes were suggested to play crucial roles in drought stress response by taking part in photosynthesis. Interestingly, three hub-genes (Zm00001d005410, Zm00001d025920 and Zm00001d008462) have been implicated in drought stress response in Arabidopsis thaliana [43–44]. Additionally, two hub-genes (Zm00001d019363 and Zm00001d020272) have been observed vital in water-deficit stress response in maize [45–46], whereas another four hub-genes (Zm00001d047235, Zm00001d025015, Zm00001d037513 and Zm00001d013781) have been linked with drought response in Camellia sinensis L. (1), soybean (2) and wheat (1), respectively [15, 47–48] (Table S6). Further, we conducted significant KEGG pathway enrichment analysis of these hub genes by using the hypergeometric test. Resultantly, by comparing the top ten pathways that were most enriched in the hub-genes, we discovered that starch and sucrose metabolism (six genes), photosynthesis (4), linoleic acid metabolism (2), and photosynthesis - antenna proteins (2) were dominant under droughts stress conditions (Fig. 5C). Moreover, the significantly enriched GO terms related to drought response were identified, including photosynthesis (GO:0015979), response to stress (GO:0006950), and response to external stimulus (GO:0009605) (Fig. 5D) .
Validation of DEGs by quantitative real-time PCR (qRT-PCR)
To confirm the accuracy of the RNA sequencing results, we conducted a validation experiment by qRT-PCR analysis for three biological replicates. The representative DEGs were chosen based on them being highly differentiated in response to drought and reported to be related to drought resistance. The gene specific primers (Table S7) used for qRT-PCR analysis were designed using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). The results of the qRT-PCR analysis validated our findings based on RNA-seq data. Precisely, the patterns of RNA-Seq expressions on all the 12 DEGs were consistent with the qRT-PCR data, suggesting that the patterns of the RNA-seq expression on all the sampled genes were replicated by the qRT-PCR approach (Fig. 6,Table S8). A correlation coefficient (R2) (of the fold changes before and after drought treatment) of 93.01% was obtained (Figure S4), endorsing our RNA-Seq data as reliable.
Physiological responses of maize hybrid ND476 to drought treatment
To determine whether the water-limited conditions could influence the physiological activities within the maize leaf tissues, in this research, we measured Pn of the drought tolerant hybrid cultivar ND476 at different growth stages. Our analysis of the four stages showed that under well-watered conditions, at the V12 stage, Pn increased initially (from 1 to 9 days post treatment exposure), and Meanwhile, all the four growth stages generally showed significantly reduced Pn under water-deficit then decreased gradually thereafter. At the VT stage, Pn showed an increasing trend from day 1 onwards. However, under well-watered conditions, at the R2 and R4 stages, Pn exhibited a slight gradual decrease throughout the treatment exposure period, starting from day 1 (Figure S5). condition as compared to well-watered condition (Figure S5). This observation may indicate that with the increased drought exposure duration, leaf stomatal closure resulted in decreased leaf available CO2, or there was increased photo-oxidative damaged induced by an accumulation of ROS.