DIA (Data Independent Acquisition) proteomic based study on maize filling-kernel stage drought stress-responsive proteins and metabolic pathways

Abstract

Drought is one of the most severe environmental stressors hampering maize (Zea mays L.) production worldwide. Therefore, revealing the molecular mechanisms underpinning maize drought tolerance is critical for guiding drought resistance breeding and ensuring global food security. Herein, we performed a data-independent acquisition (DIA) proteomics-based comparative analysis of the proteome responses of two contrasting maize hybrids (drought-tolerant NongDan 476 and drought-sensitive ZhongXin 978) exposed to 12 days of water-deficit (drought) treatment at the kernel filling stage. Resultantly, our DIA-based analysis approach identified a total of 826 differentially abundant proteins (DAPs), among which 7 DAPs were specifically accumulated in tolerant genotype ND 476 under drought treatment conditions. Moreover, our results showed that there were significant differences in drought stress responses between the two hybrids at the proteomic level. Predominantly, our analysis revealed that maize drought tolerance is attributed to some key proteins such as chaperones, DNA replication related enzymes, vicilin and ABA-responsive RAB17 protein, participating in several metabolic networks including starch and sucrose metabolism, terpenoid backbone biosynthesis, cyanoamino acid metabolism and carbon metabolism pathways. Our findings provide further insights into the drought stress tolerance mechanisms in maize at the grain filling period, in addition to useful reference data for maize drought resistance breeding.

Keywords:

• Maize
• drought stress
• DIA
• proteome
• DAPs
• qRT-PCR
• starch and sucrose metabolism

Introduction

Adverse environmental factors, of which drought stress represents the most severe constraint to agriculture, account for approximately 70 percent of potential yield loses across the world [1, 2]. Global warming and climate change are also predicted to affect most severely arid climates, where agricultural ecosystems are most susceptible to climatic condition variations and where marginal increases in environmental temperature are very harmful to productivity [3–5]. The development of crop cultivars with increased drought stress tolerance, by both conventional breeding and genetic engineering approaches, is central in meeting global food demands for the growing human population under such resource constrained environment [6, 7]. More precisely, understanding how plants respond to drought stress at the molecular level is useful for developing improved genotypes which would perform well under water limited conditions [8].

Maize (Zea mays L.), one of the most important food crops in the world, is very sensitive to water-deficit stress, especially during flowering, pollination and embryo development [9, 10]. Prominently, the yield loss during the grain filling period in maize is usually greater than that occurring during the vegetative phase [11]. Severe water stress during the maize grain filling period results in premature formation of the kernel black layer [12, 13]. Therefore, revealing the mechanism of maize drought response during the grain-filling stage remains important for breeding drought-tolerant maize cultivars.

In response to drought stress, a wide range of physiological responses are activated at varying levels of plant organization. Activation of these physiological responses would eventually affect the transcription factor regulation and downstream stress responsive genes expression [14, 15]. Examples of such responses to drought stress include but are not limited to the wilting, rolling and abscission of foliage; modification of the root to shoot ratio and antioxidant enzyme activities activation [16, 17].

Meanwhile, there have been many studies on the analysis of crop resistance to drought stress, but most of them have mainly focused on physiology and biochemistry [18–20] and transcriptomic analysis [21–24]. For instance, Guo et al. [25] examined the transcriptomic changes induced by drought and salt stress in 8-week-old seedlings of Picea wilsonii via RNA-seq analysis. Previously, we conducted comparative RNA-seq transcriptome studies [26–28] and identified several candidate genes and metabolic pathways for drought tolerance in maize, including secondary metabolite biosynthesis related enzymes, such as probable O-methlytransfarase 2 (OMT2; Zm00001d019613) [28]; transcription factors (TF) such as MYBs, NAC and WRKYs; detoxification and stress responsive proteins including aquaporin PIP2-4 (Zm00001d017288), HSPs and GSTs [26–28]. However, transcriptome analysis of drought tolerance mechanisms is limited because mRNA translation levels are not always related to the corresponding protein levels and proteins are the direct bearers of biological function [29, 30]. Since proteins are directly involved in plant stress responses, proteomic studies can eventually contribute to dissecting the possible relationships between protein changes and plant stress tolerance [31, 32]. Therefore, proteomic analysis of the mechanisms underlying maize drought stress tolerance becomes essential [33, 34].

Earlier, we employed the isobaric tags for absolute quantification (iTRAQ) strategy to perform proteomics analysis and identify key proteins and metabolic pathways regulating the drought stress response in maize [35–38]. We observed diverse key drought responsive proteins, including chlorophyll a-b binding proteins (Lhcb5-1 and 542320/Lhcb5-2) associated with the light-harvesting complex; lipid metabolism related non-specific lipid transfer proteins [35]; histone H2A proteins (Zm00001d006547, 100216750), carbohydrate metabolism-related protein malate dehydrogenase (MDH; Zm00001d009640) and several peroxidases [36]; histone deacetylase (HDAC), HSPs, MDH, GSTs, seed storage proteins 50kD γ-zein, Z1D alpha zein [37]; peroxidase (Zm00001d008898), stress responsive protein (100283176) and oxidoreductase family proteins [38]. Notably, these genes and proteins were mainly involved in photosynthesis antenna proteins, nitrogen metabolism, secondary metabolism related pathways, especially in the tolerant maize genotype ND476 under drought stress [28, 35, 38]. This has given us some insights into maize drought tolerance at different growth stages, although gaps in our knowledge still persist.

Besides iTRAQ proteomics, data-independent acquisition (DIA) has become a powerful tool for performing large-scale studies and comprehensive identification of drought responsive proteins in plants [39]. The DIA approach is another strategy used in mass spectrometry to determine molecular structures [40]. Unlike the data-dependent acquisition (DDA) method, which selects a fixed number of precursor ions for analysis by tandem mass spectrometry (TMS), the DIA approach alternatively uses all the ions within a selected m/z range, which undergoes fragmentation before analysis at the second phase of TMS [41, 42]. Compared with the DDA method, the DIA strategy has the advantages of high data utilization, high repeatability and high quantitative accuracy [43]. Therefore, DIA technology can obtain more accurate and rich results, suitable for protein detection of large sample size and complex systems.

Here, in order to study the response mechanisms of maize filling- kernel drought stress at the protein level, we used a DIA-based quantitative strategy to perform proteome analysis on two contrasting maize hybrids. In addition, we also evaluated some physiological indices of these two hybrids under drought stress, and the results of this research provide further insights into the drought stress tolerance mechanisms in maize.

Materials and methods

Plant materials and treatment

Two maize hybrids with dissimilar drought tolerance (drought-tolerant NongDan 476 and drought-sensitive ZhongXin 978) were used in this study [26]. The materials were provided by the laboratory (North China Key Laboratory for Crop Germplasm Resources of Education Ministry, Hebei Agricultural University, Baoding, China). The field experiment was carried out at Qing Yuan Experimental Station (115 E, 38 N) of Hebei Agricultural University. The experiment was set up under a fully automated rain-proof shelter, and a randomized complete block design was used to plant maize hybrid lines. The size of each plot is 5 m × 5 m, and the plant spacing is 0.6 m × 0.3 m. We planted maize on 15th June 2018, and the two hybrid lines were irrigated with adequate water before grain filling. One week before the flowering period, we performed strict bagging in order to prevent maize plants from self-pollination. The maize in the experimental group was subjected to drought treatment on 23^rd August 2018, about 13 days after pollination, and continued the artificial drought for 12 days until the afternoon of 4^th September 2018, while the control group received sufficient watering. The intermediate seeds of the ears were selected for proteomic analysis. During the sample collection process, the collected samples were immediately frozen in liquid nitrogen, and then stored at −80 °C for future use.

Phenotypic characteristics measurement

Under well-watered and drought stress conditions, the phenotypic characteristics of the two hybrids at the late grain-filling stage were measured. Ten ears were randomly selected from each group and used to measure the four phenotypic data of each ear’s ear length, barren tip length, kernel rows per ear and kernels per row. Then we calculate the average and use SPSS software for significance analysis.

Protein extraction

The sample was ground to a powder in liquid nitrogen, then dissolved in 2 mL lysis buffer containing 8 mol/L urea, 2% sodium dodecyl sulphate (SDS), and 1× Protease Inhibitor Cocktail (Roche Ltd. Basel, Switzerland), followed by sonication on ice for 30 min and centrifugation at 15,000 g for 30 min at 4 °C. The supernatant was transferred to a new tube. For each sample, proteins were precipitated with ice-cold acetone at −20 °C overnight. The precipitations were cleaned with acetone three times and re-dissolved in 8 mol/L urea by sonication on ice. BCA Protein Assay Kit (Genedenovo Biotechnology Co., Ltd, Guangzhou, China) was used to determine the protein concentration of the supernatant. The protein extract quality was examined with SDS polyacrylamide gel electrophoresis (SDS-PAGE) [44].

Protein digestion and high PH reverse phase separation

For each sample, 50 μg proteins extracted from cells were suspended in 50 μL solution, reduced by adding 1 μL of 1 mol/L dithiotreitolat 55 °C for 1 h, alkylated by adding 5 μL of 20 mmol/L iodoacetamide in the dark at 37 °C for 1 h. Then the sample was precipitated using 300 μL prechilled acetone at −20 °C overnight. The precipitate was washed twice with cold acetone and then re-suspended in 50 mmol/L ammonium bicarbonate. Finally the proteins were digested with sequence-grade modified trypsin (Promega, Madison, WI) at a substrate/enzyme ratio of 50:1 (w/w) at 37 °C for 16 h.

The peptide mixture was re-dissolved in buffer A (buffer A: 20 mmol/L ammonium formate in water, pH 10.0, adjusted with ammonium hydroxide), and then fractionated by high pH separation using Ultimate 3000 system (ThermoFisher scientific, MA, USA) connected to a reverse phase column (XBridge C18 column, 4.6 mm × 250 mm, 5 μm, (Waters Corporation, MA, USA). High pH separation was performed using a linear gradient, starting from 5% B to 45% B in 40 min (B: 20 mmol/L ammonium formate in 80% ACN, pH 10.0, adjusted with ammonium hydroxide). The column was re-equilibrated at the initial condition for 15 min. The column flow rate was maintained at 1 mL/min and the column temperature was maintained at 30 °C. Ten fractions were collected; each fraction was dried in a vacuum concentrator for the next step.

DIA: nano-HPLC-MS/MS analysis

The peptides were re-dissolved in 30 μL solvent A (A: 0.1% formic acid in water) and analyzed by on-line nanospray liquid chromatography with tandem mass spectrometry (LC-MS-MS) on an Orbitrap Fusion Lumos coupled to an EASY-nLC 1200 system (Thermo Fisher Scientific, MA, USA). Briefly, 3 μL peptide sample was loaded onto the analytical column (Acclaim PepMapC18, 75 μm × 25 cm) with a 120-min gradient, from 5% to 35% B (B: 0.1% formic acid in ACN). The column flow rate was maintained at 200 nL/min, and the column temperature at 40 °C. The electrospray voltage of 2 kV versus the inlet of the mass spectrometer was used.

The mass spectrometer was run under data independent acquisition mode, and automatically switched between MS and MS/MS mode. The parameters were: (1) MS: scan range (m/z) = 350–1200; resolution = 120,000; AGC target = 1e6; maximum injection time = 50 ms; (2) HCD-MS/MS: resolution = 30,000; AGC target = 1e6; collision energy = 32；stepped CE = 5%. (3) DIA was performed with a variable isolation window; each window overlapped 1 m/z, and the window number was 60.

Protein identification and quantification

Raw data of DIA were processed and analyzed by Spectronaut X (Biognosys AG, Switzerland) with default parameters. Retention time prediction type was set to dynamic iRT. Data extraction was determined by Spectronaut X based on the extensive mass calibration. Spectronaut Pulsar X determines the ideal extraction window dynamically depending on iRT calibration and gradient stability. False discovery rate (FDR) cutoff on precursor and protein level was applied at 1%. Decoy generation was set to mutated, which is similar to scrambled but will only apply a random number of AA position swamps (min = 2, max = length/2). All selected precursors passing the filters were used for quantification. The average top 3 filtered peptides which passed the 1% FDR cutoff were used to calculate the major group quantities. After Student’s t-test, differentially expressed proteins were filtered if their Q value was Q < 0.05 and the absolute AVG log2 ratio > 0.58.

Functional annotation and enrichment analysis of DAPs

We performed gene ontology (GO) annotation to assign GO terms to the DAPs using the Blast2GO web-based program (<seurld>https://www.blast2go.com/</seurld>). Further, the DAPs were searched against GO (Gene Ontology, http://www.geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomics, http://www.genome.jp/kegg/) databases to obtain their biological functions. The GO analysis [45] functionally categorized the identified DAPs into their BP, MF and CC involvement in response to drought stress. Moreover, the DAPs were assigned to various metabolic pathways using the KEGG pathway analysis. Further, significant GO functions and KEGG pathways enrichment analysis were examined within DAPs with the hypergeometric test, Q (Bonferroni-corrected p-value) < 0.05 set as the statistically significant threshold.

RNA extraction, cDNA synthesis and qRT-PCR analysis

We used the M5 Quickspin (Mei5 Biotechnology, Beijing, China) to extract total RNA. Total RNA was extracted from non-stressed and stressed kernels of the two hybrids (ND 476 and ZX 978) and prepared for quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis. To generate cDNA templates, 2 µg of total RNA was reverse-transcribed in a total volume of 20 µL, using the PrimerScript^TM RT reagent Kit (TAKARA, Dalian, China) according to the manufacturer’s instructions. We randomly selected 15 genes and designed gene-specific primers (Supplemental Table S1) for qRT-PCR using Primer Premier 5 Designer software. We used LightCycler® 96 for qRT-PCR, using TB Green ® Premix Ex Taq^TM II (Tli RNaseH Plus) (TAKARA, Dalian, China). Each total 20 µL qRT-PCR reaction mixture comprised 1 µL of template cDNA, 0.8 µL of each primer, and 10 µL of 2 × TB Green Premix Ex Taq II (TAKARA, Dalian, China). The amplification program was as follows: 95 °C for 30 s followed by 40 cycles of 95 °C for 5 s and 60 °C for 20 s, 65 °C for 15 s. The maize gene GAPDH (accession no. X07156), with forward (GAPDH-F: 5′-ACTGTGGATGTCTCGGTTGTTG-3′), and reverse (GAPDH-R: 5′-CCTCGGAAGCAGCCTTAATAGC-3′) primers, was used as an internal control. We performed three biological replicates on each sample. The relative mRNA abundance was estimated using Livak and Schmittgen’s 2^−ΔΔCt method [46].

Table 1. Drought-responsive DAPs observed specifically in ND 476.

No.	Gene ID^1	Description^2	Log2 fc^3	Expr.^4	p Value^5	Pathways^6
1	Zm00001d017494_T001	Predicted: beta-expansin 4 isoform X1	0.8782	Up	2.28E-04	–
2	Zm00001d042025_T001	Ubiquitin-fold modifier 1	0.6784	Up	2.06E-04	–
3	Zm00001d034074_T004	Starch phosphorylase1	−0.7863	Down	1.06E-06	Starch and sucrose metabolism
4	Zm00001d034074_T008	Starch phosphorylase1	−0.7863	Down	1.06E-06	Starch and sucrose metabolism
5	Zm00001d034074_T009	Starch phosphorylase1	−0.7863	Down	1.06E-06	Starch and sucrose metabolism
6	Zm00001d034074_T013	Starch phosphorylase1	−0.7863	Down	1.06E-06	Starch and sucrose metabolism
7	Zm00001d034074_T018	Starch phosphorylase1	−0.7863	Down	1.06E-06	Starch and sucrose metabolism

¹ Gene ID, the name of the transcript corresponding to the protein.

² Description, annotated biological functions.

³ log2fc, is expressed as the logarithmic value (log 2) of intensities of proteins.

⁴ Expr., Up, up-regulated; Down , down-regulated.

⁵ p value, statistical level (using Student’s t-test) below < 0.05, at which protein differential expression was accepted as significant.

⁶ Pathways, metabolic KEGG pathways.

Results and discussion

Phenotypic differences between two hybrids

To validate the previous observations that ZX 978 is drought-sensitive and ND 476 is drought-tolerant, the two maize hybrids were treated with or without moisture deficit stress in a fully automated rain-proof shelter. We then recorded the phenotypic data of sample ears from two treatment groups of ND 476 and ZX 978 (Figure 1a,b). From the phenotypic perspective, before and after drought treatment, the four phenotypic indicators including ear length, barren tip length, kernel rows per ear and kernels per row did not change significantly in ND 476, which indicated that ND 476 is a comparatively drought-tolerant hybrid line (Figure 1c–f). Although the kernel rows per ear of ZX 978 did not change, the ear length was significantly reduced, the barren tip was significantly increased, and the kernels per row were significantly reduced, which indicated that ZX 978 is a drought-sensitive hybrid line (Figure 1c–f).

Figure 1. Phenotypic characterization of the two maize hybrid lines (drought-tolerant ND 476, T; and drought-sensitive ZX 978, S) ears’ responses to drought stress. Observations/measurements were made after 12 days of water-sufficient (control, C) or water-deficit (drought, D) conditions. (a-b) Ear phenotypes; (c) ear length; (d) ear bare tip length; (e) kernel rows per ear; (f) number of kernels per row. Data are presented as mean values ± standard errors (n = 3). Different letters on error bars indicate significant differences at p < 0.05; For c-f, each replication is an average for the measurement of 10 ears.

Phenotypic charactersitics could not show obvious changes in ND476 genotype before and after drought treatment, suggesting that the perception of the stress by ND476 could have been mild, which could not warrant significant phenotypic modifications. This is consistent with the expected results of ND 476 as a drought-tolerant hybrid line and ZX 978 as a drought-sensitive hybrid line. Therefore, we can infer that under drought conditions, water deficit will lead to short ear length and slow growth of maize; increased barren tip length, insufficient grain filling and reduced number of kernels per row.

Inventory of maize proteins identified by DIA-based analysis approach

The 830,311 spectra obtained by DIA proteomic strategy were matched to known spectra using Spectronaut software, and 50,587 precursors and 36,717 peptides were identified. Consequently, 5276 protein groups and 13,659 proteins were identified from the 36,717 peptides. Further analysis showed that of those 36,717 peptides, 2941 (8.01%) were < 1000 Da, 23,814 (64.86%) were 1000–2000 Da, 8049 (21.92%) were 2000–3000 Da, and 1913 (5.21%) were > 3000 Da in molecular weight (Supplemental Figure S1a). The distribution of the peptide lengths indicated that over 80% of the peptides had lengths between 7 and 20 amino acids (Supplemental Figure S1b). In terms of number of peptides defining each protein, our analysis showed that among the 5276 protein groups, 83.95% (4429 protein groups) had proteins composed of at least two peptides, whereas the remaining 16.05% (847 protein groups) were composed of only one identified peptide (Supplemental Figure S1c). Additionally, we used the local normalization method in Pulsar software to normalize the peak intensity of the whole sample. From the distribution of quantized values of normalized peptides, the signal intensity of most samples reached the same response intensity (Supplemental Figure S1d). Therefore, the DIA sequencing results qualified for use in the subsequent analyses.

Analysis of DAPs observed in different comparisons

Similar to the approach used previously [38], comparative proteomic analysis was used to investigate the changes of maize filling-kernel protein profiles in ND 476 (drought-tolerant, T) hybrid and ZX 978 (drought-sensitive, S) hybrid under drought stress condition. A pairwise comparison of before and after treatments (drought, D, and control, C) was performed in ND 476 (TD_TC) and ZX 978 (SD_SC) individually. In addition, a comparative study on the drought stress proteome was performed between the tolerant hybrid and sensitive hybrid, under drought (SD_TD) and under water-sufficient (control) (SC_TC) conditions, giving four comparison groups (Figure 2a). In the drought-tolerant hybrid line, 36 proteins showed differential abundance before and after drought treatment (TD_TC), 8 of these DAPs were up-regulated, 28 were down-regulated (Figure 2a). In the drought-sensitive hybrid line, we identified 130 DAPs before and after drought treatment (SD_SC), 28 of these DAPs were up-regulated whilst 102 were down-regulated (Figure 2a). In total, 826 DAPs were found among the four comparison groups (Figure 2a,b).

Figure 2. Analysis of differentially abundant proteins (DAPs) identified in four experimental comparisons. (a) Total number of DAPs identified in each experimental comparison group. Red (up-regulated) means DAPs with increased differential abundance. Green (down-regulated) means DAPs with decreased differential abundance; (b) Venn diagram analysis of DAPs identified in the four comparisons. DAPs uniquely expressed in TD_TC (I), SD_TD (II), SD_SC (III) are indicated with arrows. Area IV shows 18 overlapping DAPs within line. Overlapping regions of the Venns indicate DAPs shared between/among corresponding groups. (c) Clustering analysis of DAPs in SD_TD comparison. Each row represents a significantly abundantly expressed protein. ZXD1-3 (drought-sensitive ZX 978, S; drought, D) refer to the biological replicate number for ZX 978, whilst NDD1-3 (drought-tolerant ND 476, T; drought, D) refer to the replicate number for ND 476. The DAPs were clustered based on the differentially expressed levels. The scale bar indicates the logarithmic value (log2 fold change) expression of the DAPs. Red (log2 fold change > 0) indicates up-regulated DAPs and blue (log2 fold change ＜ 0) indicates down-regulated DAPs.

According to Figure 2b, it can be seen that area I encompasses the 7 specific DAPs of drought-tolerant hybrid ND 476 (Table 1); area II represents 148 specific DAPs of SD_TD and analysis of DAPs revealed that under drought stress, the number of DAPs in ZX 978 was significantly higher than that in the drought-tolerant hybrid ND 476 (Figure 2c and Supplemental Figure S2); area III is a specific DAPs of the sensitive hybrid line ZX 978, with 45 (Table 2). Area IV represents the 18 DAPs shared by TD_TC and SD_SC, that is, the overlapping drought responsive DAPs within the two hybrid lines. Of these 18 common drought responsive DAPs, all were down-regulated in tolerant hybrid ND 476; whereas 2 were up-regulated and 16 down-regulated in sensitive hybrid ZX 978 (Table 3).

Table 2. Drought-responsive DAPs observed specifically in ZX 978.

No.	Gene ID^1	Description^2	Log2 fc^3	Expr.^4	p Value^5	Pathways^6
1	Zm00001d029631_T001	Cystatin 8 precursor	1.7236	Up	1.77E-03	–
2	Zm00001d009888_T001	DNA primase	1.3119	Up	1.36E-03	Pyrimidine metabolism; Purine metabolism; DNA replication
3	Zm00001d003869_T001	Predicted: adenylate isopentenyltransferase 5, chloroplastic-like	1.2548	Down	4.37E-04	Zeatin biosynthesis
4	Zm00001d047548_T001	17.0 kDa class II heat shock protein	1.2546	Up	8.31E-06	Protein processing in endoplasmic reticulum
5	Zm00001d047553_T001	17.0 kDa class II heat shock protein	1.2546	Up	8.31E-06	Protein processing in endoplasmic reticulum
6	Zm00001d018107_T001	Phi-1-like phosphate-induced protein precursor	1.02	Up	6.52E-05	–
7	Zm00001d047203_T001	Uncharacterized protein LOC100276254	1.0061	Up	8.31E-04	–
8	Zm00001d020364_T001	Early flowering 4 protein	0.8773	Up	1.21E-03	–
9	Zm00001d029720_T001	ABC-type Co2+ transport system, permease component	0.8375	Up	1.26E-03	–
10	Zm00001d041822_T001	Betl-9 protein precursor	0.8	Up	2.99E-08	–
11	Zm00001d049683_T001	Predicted: L-gulonolactone oxidase-like	0.7541	Up	5.51E-07	–
12	Zm00001d019277_T001	Predicted: defensin-like protein P322[Nicotiana tabacum]	0.7357	Up		–
13	Zm00001d003311_T001	2-Oxoglutarate-dependent dioxygenase DAO [Dichanthelium oligosanthes]	0.6171	Up	6.89E-06	–
14	Zm00001d031128_T004	5-Methyltetrahydropteroyltriglutamate– homocysteine methyltransferase 1	−0.5876	Down	2.02E-04	Biosynthesis of amino acids; Cysteine and methionine metabolism; Selenocompound metabolism
15	Zm00001d048176_T003	Vicilin-like seed storage protein	−0.5956	Down	3.41E-14	–
16	Zm00001d048176_T004	Vicilin-like seed storage protein	−0.5956	Down	3.41E-14	–
17	Zm00001d048176_T007	Vicilin-like seed storage protein	−0.5956	Down	3.41E-14	–
18	Zm00001d048176_T009	Vicilin-like seed storage protein	−0.5956	Down	3.41E-14	–
19	Zm00001d048176_T014	Vicilin-like seed storage protein	−0.5956	Down	3.41E-14	–
20	Zm00001d052056_T005	T-complex protein 1 subunit alpha	−0.6016	Down	3.48E-04	–
21	Zm00001d026632_T001	Stem-specific protein TSJT1	−0.6029	Down	2.12E-08	–
22	Zm00001d028755_T002	Ran-binding protein 1 homolog a	−0.6406	Down	1.81E-03	–
23	Zm00001d014325_T001	Uncharacterized LOC100839725 [Brachypodium distachyon]	−0.6475	Down	1.91E-03	–
24	Zm00001d033623_T002	Lipoxygenase	−0.6877	Down	9.49E-06	Linoleic acid metabolism
25	Zm00001d047567_T001	Embryonic cell protein 63	−0.6877	Down	4.75E-04	–
26	Zm00001d046700_T001	NAD(P)-binding Rossmann-fold superfamily protein	−0.7527	Down	1.35E-03	–
27	Zm00001d046700_T002	NAD(P)-binding Rossmann-fold superfamily protein	−0.7527	Down	1.35E-03	–
28	Zm00001d046700_T003	NAD(P)-binding Rossmann-fold superfamily protein	−0.7527	Down	1.35E-03	–
29	Zm00001d033683_T001	Phospholipase A2	−0.8001	Down	6.61E-04	Glycerophospholipid metabolism; Ether lipid metabolism; alpha-Linolenic acid metabolism; Linoleic acid metabolism; Arachidonic acid metabolism
30	Zm00001d033683_T003	Phospholipase A3	−0.8001	Down	6.61E-04	Glycerophospholipid metabolism; Ether lipid metabolism; alpha-Linolenic acid metabolism; Linoleic acid metabolism; Arachidonic acid metabolism
31	Zm00001d034801_T001	11-Beta-hydroxysteroid dehydrogenase-like 5	−0.8327	Down	1.48E-03	–
32	Zm00001d006725_T001	Eukaryotic peptide chain release factor subunit 1-1	−0.944	Down	1.50E-03	mRNA surveillance pathway
33	Zm00001d005283_T001	Uncharacterized protein Loc100278520	−0.9502	Down	1.92E-03	–
34	Zm00001d005284_T001	Uncharacterized protein Loc100278520	−0.9502	Down	1.92E-03	–
35	Zm00001d032519_T002	Predicted: apomucin-like isoform X2	−0.9787	Down	4.90E-08	–
36	Zm00001d015406_T001	Succinate-semialdehyde dehydrogenase mitochondrial	−1.071	Down	1.84E-04	Alanine, aspartate and glutamate metabolism; Butanoate metabolism
37	Zm00001d015405_T001	Succinate semialdehyde dehydrogenase	−1.1428	Down	8.62E-05	Alanine, aspartate and glutamate metabolism; Butanoate metabolism
38	Zm00001d015405_T002	Succinate semialdehyde dehydrogenase	−1.1428	Down	8.62E-05	Alanine, aspartate and glutamate metabolism; Butanoate metabolism
39	Zm00001d040026_T001	Bowman-Birk type trypsin inhibitor	−1.3178	Down	3.01E-05	–
40	Zm00001d040026_T002	Bowman-Birk type trypsin inhibitor	−1.3178	Down	3.01E-05	–
41	Zm00001d040026_T003	Bowman-Birk type trypsin inhibitor	−1.3178	Down	3.01E-05	–
42	Zm00001d040026_T012	Bowman-Birk type trypsin inhibitor	−1.3178	Down	3.01E-05	–
43	Zm00001d029218_T001	Protein SRG1	−1.5563	Down	1.04E-03	–
44	Zm00001d003938_T001	Steroleosin	−1.8606	Down	2.56E-04	–
45	Zm00001d017022_T001	Predicted: late embryogenesis abundant protein D-34-like [Setaria italica]	−2.0245	Down	8.16E-05	–

For full description of the column items, refer to Table 1 captions above.

Table 3. Overlapping drought-responsive DAPs between ZX 978 and ND 476.

No.	Gene ID	Description	ZX 978			ND 476			Pathways
			log2 fc	expr.	p Value	log2 fc	expr.	p Value
1	Zm00001d033447_T005	Globulin-1 S allele	−1.4572	Down	2.71E-16	−0.7086	Down	5.49E-08	–
2	Zm00001d023516_T001	Salt stress-induced protein	1.6153	Up	3.16E-08	−0.7246	Down	2.18E-07	–
3	Zm00001d036242_T001	Flowering locus T protein	1.2334	Up	5.95E-07	−0.7246	Down	2.18E-07	Circadian rhythm - plant
4	Zm00001d037894_T001	RAB17 protein	−1.7073	Down	3.94E-05	−0.7246	Down	2.18E-07	–
5	Zm00001d045042_T001	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
6	Zm00001d045042_T002	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
7	Zm00001d045042_T003	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
8	Zm00001d045042_T005	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
9	Zm00001d045042_T006	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
10	Zm00001d045042_T007	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
11	Zm00001d045042_T010	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
12	Zm00001d045042_T011	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
13	Zm00001d045042_T013	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7246	Down	2.18E-07	Starch and sucrose metabolism
14	Zm00001d045042_T014	Sucrose synthase 1	−0.9582	Down	4.49E-05	−1.2293	Down	1.93E-05	Starch and sucrose metabolism
15	Zm00001d045042_T015	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.8872	Down	6.06E-05	Starch and sucrose metabolism
16	Zm00001d045042_T016	Sucrose synthase 1	−0.9582	Down	4.49E-05	−0.7822	Down	1.20E-04	Starch and sucrose metabolism
17	Zm00001d003801_T001	Vicilin-like seed storage protein	−1.6559	Down	1.90E-04	−3.0691	Down	5.00E-04	–
18	Zm00001d009431_T003	Farnesyl pyrophosphate synthase1	−1.1879	Down	1.13E-03	−0.7454	Down	6.90E-04	Terpenoid backbone biosynthesis

For full description of the column items, please refer to Table 1 captions above.

Expression of ubiquitin associated proteins in ND 476

Among these drought-responsive proteins whose abundance was significantly altered under drought stress in tolerant hybrid ND 476, we observed that ubiquitin-fold modifier 1 (Ufm1) was up-regulated. Ufm1 is a ubiquitin-like modifier and conjugates to target proteins in cells [47, 48]. In addition to being a recognition signal, ubiquitin regulates cellular processes directly by altering the biophysical properties of the substrate, suggesting that ubiquitination facilitates the unfolding process and improves the substrate degradation efficiency [49], hence Ufm1 is up-regulated under drought stress, which enhances the rate of protein degradation. Moreover, Xia et al. [50] found that the ubiquitin ligase was up-regulated in maize response to drought stress.

The response of maize to drought stress is a complex process involving DNA repair process, hydrolysis process, signal transduction, intracellular environment stability, and transcription and translation, etc. [51]. In particular, both genetic analyses and the rapid accumulation of small ubiquitin-like modifiers (SUMO) conjugates in response to various adverse environmental conditions suggest that SUMOylation plays a key role in the stress responses [52]. SUMO-conjugating enzyme (SCE) mediates addition of SUMO group to various cell proteins, through process referred to as SUMOylation [53]. SUMO, in conjunction with intracellular proteins, actively participate in several biological functions such as altering signal transduction in the cytoplasm, transcription regulation, DNA repair, as well as regulation of nuclear material importation and subnuclear compartmentalization [52, 54]. Moreover, the ubiquitin-dependent proteolytic pathway degrades most proteins and is the major proteolytic mechanism in eukaryotic cells [55]. Thus, here, the up-regulation of SCE1 implies that cells promote the proteins and enzymes involved in protein ubiquitination in order to protect themselves under drought stress [35].

Other proteins for drought tolerance in ND 476

Among the dominating up-regulated DAPs in ND 476, we observed 17.5 kDa class II heat shock proteins (see TC_TD comparison, Supplemental Table S2). Heat shock proteins (HSPs) are generally activated in plants’ response to external stimuli emanating from a wide array of abiotic stressors such as water deficit, nutrient deficiency, oxidative stress, among others [56–59]. Thus, HSPs are key members of the plant stress proteins involved in stress responses [60]. Particularly, HSPs provide chaperone functions to stabilize newly synthesized proteins by maintaining correct folding, or by assisting the refolding of damaged proteins [61]. Molecular chaperones facilitate the stabilization of other macromolecular structures under stress conditions [62]. By stabilizing partially unfolded proteins, HSPs help the easy transportation of other proteins and cellular materials across membranes [63, 64]. Additionally, HSPs help preserve other proteins in their functional confirmations, thereby offering plant protection against environmental stresses [65]. A previous study has demonstrated that HSPs are greatly accumulated in alfalfa (Medicago sativa L.) leaves in response to salinity stress [66].

In our study, we found that starch phosphorylase 1 (STP-1) was down-regulated. STP-1 is largely known for the phosphorolytic degradation of starch [67]. Moreover, drought stress caused a marked reduction in the starch content of the kernel [68], shortened the grain filling period, and increased the grain filling rate [69]. Therefore, STP-1 was down-regulated to resist the decrease of starch content in kernel caused by drought stress, and reduced material consumption to protect the kernels from drought damage. In addition, Prasch et al. [70] suggested that impaired starch breakdown was accompanied by decreased stomatal opening in Arabidopsis thaliana, resulting in improved drought tolerance of Arabidopsis plants.

Enhancement of DNA replication in ZX 978

In the present study, adenine phosphoribosyltransferase 2 (APRTase 2) and DNA primase were up-regulated in response to drought stress (see SC_SD comparison, Supplemental Table S3). APRTase is an enzyme encoded by the APRT gene [71], and being part of the type I PRT enzyme family, it is involved in the nucleotide salvage pathway [72]. APRTase provides a mechanism for the production of adenine [73]. It was found that the increase of adenine in wheat significantly increased the levels of free proline, amino acid and soluble sugar, which increased the photosynthesis, water status and chlorophyll content under drought stress [74]. Therefore, the up-regulation of APRTase 2 in ZX 978 is used to combat drought stress by regulating osmotic potential and enhancing photosynthesis.

DNA primase is an enzyme involved in the replication of DNA and is a type of RNA polymerase. DNA primase plays an essential role in controlling cell cycle progression in response to DNA damage [75]. Faria et al. [76], studying Medicago truncatula, suggested that desiccation-tolerant seedlings may apply distinct strategies to survive in drought, either by avoidance or further repair of the DNA damage. Hence, DNA primase is an up-regulated enzyme that mitigates the damage caused by drought by repairing damage in DNA replication.

Overlapping drought responsive proteins between ND 476 and ZX 978

Venn diagram (Figure 2b, Area IV) analysis showed that 18 DAPs were common between TD_TC and SD_SC. All the 18 DAPs (Table 3) were down-regulated in tolerant hybrid ND 476 in response to drought treatment. Conversely, among these 18 common DAPs, two (salt stress-induced protein and flowering locus T protein) were up-regulated, whilst 16 (globulin-1 S allele, RAB17 protein, sucrose synthase 1, farnesyl pyrophosphate synthase1 and vicilin) were down-regulated in sensitive hybrid ZX 978 in response to drought treatment (Table 3).

Flowering locus T protein (FT), a member of the phosphatidylethanolamine-binding protein family, functions as a florigen or a component of florigen [77]. Long-distance signal called florigen is produced in the leaf under inductive day length conditions, and is transported to the shoot apex where it triggers floral morphogenesis [78]. Regulation of the transition process from vegetative to anthesis is vital in plants’ drought stress acclimation [79, 80]. The drought-escape response allows plants to adaptively shorten their life cycle to make seeds before severe stress leads to death [80]. Therefore, FT is up-regulated in ZX 978 to promote seeds generation to shorten the life cycle to achieve drought escape. Similarly, Su et al. [81] found that the FT was up-regulated in Arabidopsis in response to drought stress. In addition, FT is down-regulated in ND 476 maybe to slow the long-distance transportation of the signal, thereby saving energy to resist drought.

Vicilin is used as a storage protein and commonly found in plants like pea or lentil [82]. It has been suggested to be an allergen in allergic reactions to peas [83]. Vicilin is synthesised as preprovicilin with subsequent removal of a signal peptide and a C-terminal peptide as well as post translational endo-proteolytic cleavage [84]. Additionally, the down-regulation of vicilin in the two hybrid lines was consistent with the inhibition of starch biosynthesis under drought stress in plants. Similarly, globulin-1 S allele, a type of polymorphic Globulin-1 gene, encodes a storage globulin [85] that is down-regulated in two hybrid lines.

In the present study, sucrose synthase 1 (SS-1) was significantly down-regulated in response to drought stress. It has been noted that premature formation of black layer following cold temperatures, drought or defoliation might all be explained on the basis of reduced sucrose availability [13]. It was suggested that SS plays a key role in the regulation of carbon metabolism and, therefore, of nitrogen fixation under drought stress conditions [86, 87]. Moreover, SS, which reversibly catalyzes sucrose synthesis and cleavage, plays pivotal roles in regulating sucrose flux [88].

The maize abscisic acid-responsive RAB17 protein localizes to the nucleus and cytoplasm in maize cells, involved in plant responses to stress [89]. Based on the similar domain arrangements of Rab17 and a nuclear localization signal (NLS) binding phosphoprotein, interaction of Rab17 with NLS peptides was found [90]. Thus, Rab17 may play a role in nuclear protein transport [91]. Moreover, binding of the NLS peptide to Rab17 is dependent upon phosphorylation by protein kinase CK2 [92]. Under drought conditions, the reduction in RAB17 protein may be due to reduced protein transport to save energy against drought stress, similar to our previous findings where the tolerant maize genotype YE8112 plants had to reduce the synthesis of redundant proteins to help themselves save energy for drought stress endurance [35]. Figueras et al. [93] observed a protective effect of RABl7 protein in vegetative tissues under osmotic stress conditions.

GO annotation and functional classification of the drought responsive DAPs

We performed GO annotation to assign GO terms to the DAPs. Further, GO functional classification of the GO-term-assigned-DAPs into biological processes (BP), molecular functions (MF), and cellular component (CC) categories was carried out. For the tolerant hybrid line ND 476 specific DAPs, the most significantly enriched GO terms in the MF category were GO: 0008194 (UDP-glycosyltransferase activity), GO: 0035251 (UDP-glucosyltransferase activity), GO: 0016758 (transferase activity, transferring hexosyl groups), GO: 0046527 (glucosyltransferase activity) and GO: 0016757 (transferase activity, transferring glycosyl groups) (Supplemental Table S4).

The significantly enriched GO terms in the sensitive hybrid line ZX 978 (SC_SD) were GO: 0010466 (negative regulation of peptidase activity), GO: 0030162 (regulation of proteolysis), GO: 0043086 (negative regulation of catalytic activity), GO: 0045861 (negative regulation of proteolysis), GO: 0050790 (regulation of catalytic activity), GO: 0051336 (regulation of hydrolase activity), GO: 0051346 (negative regulation of hydrolase activity) and GO: 0052547 (regulation of peptidase activity) in the BP category (Supplemental Table S5); GO: 0004866 (endopeptidase inhibitor activity), GO: 0035251 (UDP-glucosyltransferase activity), GO: 0061134 (peptidase regulator activity), GO: 0061135 (endopeptidase regulator activity) and GO: 0030414 (peptidase inhibitor activity) in the MF functions (Supplemental Table S5).

Additionally, the amount of DAPs in the two comparison combinations (TC_TD, SC_SD) annotated into each biological function is shown in Figure 3. Among the tolerant hybrid line ND 476 (TC_TD) -specific DAPs, metabolic process was the most enriched BP function; and catalytic activity was most prominent in the MF category (Figure 3a). Among the sensitive hybrid line ZX 978 (SD_SC)-specific DAPs, metabolic process was the most common biological processes; whereas catalytic activity was prominent in the MF category (Figure 3b).

Figure 3. Gene ontology (GO) functional classification of differentially abundant proteins (DAPs) under drought conditions. Most significantly enriched GO terms in tolerant hybrid line ND 476 (a) and sensitive hybrid line ZX 978 (b). Note: The number above each bar graph shows the enrichment factor of each GO term. Red indicates up-regulated and blue indicates down-regulated.

KEGG pathway enrichment analysis of DAPs

To further analyze the functional consequences of the drought-responsive DAPs, we mapped them to the Kyoto Encyclopedia of Genes and Genomes (KEGG, <seurld>https://www.genome.jp/kegg/</seurld>) database and the DAPs were assigned to various biological pathways. The top three enriched pathways in ND 476 were starch and sucrose metabolism (22 proteins), ubiquitin mediated proteolysis (2) and RNA transport (2) (Figure 4a). However, the composition of the enriched KEGG pathways in ZX 978 differed significantly, with cyanoamino acid metabolism (14 proteins), nitrogen metabolism (12), starch and sucrose metabolism (15), glyoxylate and dicarboxylate metabolism (12) and carbon metabolism (14) being the most enriched pathways (Figure 4b).

Figure 4. KEGG pathway enrichment analysis of the DAPs; TD_TC comparison (a); SD_SC comparison (b). Note: The bubble size represents the number of DAPs enriched in the corresponding pathway. The vertical coordinate shows KEGG pathway. The horizontal coordinate shows the enrichment factor (rich factor ≤ 1), and the enrichment factor indicates the number of differentially abundant proteins (DAPs) participating in a KEGG pathway as a proportion of proteins involved in the pathway in all identified proteins. The significance of the enrichment of the KEGG pathway is based on Student’s t-test, q < 0.05. The color gradient represents the size of the q value; the color is from blue to red, and the nearer red represents the smaller q value, the higher the significant level of enrichment of the corresponding KEGG pathway.

Most significantly enriched metabolic pathways in the two genotypes’ response to drought stress

Using the hypergeometric test, KEGG pathways that had a q value < 0.05 were considered to be significantly affected by drought stress. Consequently, we observed that ‘starch and sucrose metabolism’, ‘terpenoid backbone biosynthesis’ and ‘ubiquitin mediated proteolysis’ were the most significantly enriched pathways in drought-tolerant hybrid ND 476 (Figure 4a). For those pathways, sucrose synthase1 and starch phosphorylase1 were down-regulated (Supplemental Table S2). These results resonate well with our GO analysis results where sucrose synthase1 were most significantly enriched in the prominent GO terms (catalytic activity and metabolic process), suggesting their role as the major drought resistant signature in ND 476 hybrid. Moreover, the sucrose content increased and the starch content decreased concomitantly, indicating that drought facilitated the partitioning of photosynthates into sucrose [94].

On the other hand, ‘cyanoamino acid metabolism’, ‘carbon metabolism’ and ‘nitrogen metabolism’ were the significantly enriched metabolic pathways in the sensitive hybrid ZX 978 (Figure 4b). Under drought conditions, the grain filling rate of maize is accelerated, but the time is short, which will lead to a decline in yield [95]. Because of this phenomenon, maize under drought stress needs more complex molecular mechanisms to maintain its own stability.

Quantitative real-time PCR (qRT-PCR) analysis

To confirm the accuracy of the DIA-seq data, we performed a verification experiment using quantitative real-time PCR (qRT-PCR). Fifteen genes were selected from different groups of DAPs (Supplemental Table S1). We found that the expression levels (log₂ FC) (FC, fold change) of 15 genes measured by qRT-PCR were not significantly different from those measured by DIA-seq (Figure 5; Supplemental Table S6). In addition, we conducted a correlation analysis of qRT-PCR date and DIA-seq date of 15 genes, showing that the correlation coefficient (R²) was 83.56% (Supplemental Figure S3). Thus, this confirmed the accuracy and reliability of our DIA-seq results.

Figure 5. Confirmation of DIA-seq results by quantitative real-time PCR (qRT-PCR). Quantitative RT-PCR analysis of the expression patterns of the maize kernel genes encoded by differentially abundant proteins (DAPs) from different comparisons: DAPs specific to TD_TC (a) and SD_SC (b). Note: The y-axis represents qRT-PCR (green) relative expression levels (log₂-fold change) and log₂-fold change of the DIA-seq data (orange). Negative values of expression level mean that the genes were down-regulated in response to drought stress. Maize gene GAPDH (accession no. X07156) was used as the internal reference. Error bars represent the SE (n = 3).

Comparative analysis of DIA and iTRAQ proteomics derived results of maize hybrid cultivars’ drought stress response at the grain-filling stage

More recently [38], we performed comparative iTRAQ proteomics analysis of the same maize hybrids (ND476 and ZX978) responses to drought stress at the same (kernel filling) growth stage. Our comparative analysis of the current (DIA) study and our most recent proteomics study [38] revealed some similarities and variations with regards to the detection capacity of the two proteomic analysis approaches and drought stress response evaluation in maize. Firstly, in terms of the detection capacity of the two strategies, the number of spectra (815,688) that were matched to known spectra using Mascot software (for iTRAQ proteomics) was close to the number of spectra (830,311) that was matched by Spectronaut software (for DIA proteomics). This may suggest that both tools can be effectively used for proteomics analysis with comparably similar success. However, a higher number of unique peptides (87,018) was identified by iTRAQ [38] as compared to significantly fewer peptides (36,717) identified by DIA in this paper. Additionally, the iTRAQ approach detected a higher number of DAPs (1655), although it had fewer (9895) total proteins [38]. Conversely, the DIA approach detected significantly fewer DAPs (826) from a significantly higher number (13,659) of total proteins (this study). Subsequently, the number of DAPs identified to be specifically accumulated in the tolerant genotype ND476 under drought stress conditions were significantly fewer (seven) by the DIA approach (Table 1) as compared to twenty by the iTRAQ method (see also Table 1 of Dong et al. [38]). This observation may suggest DIA’s high data utilization and quantitative accuracy [43]. However, both strategies can be employed in proteomics analysis with effectiveness. Further detailed comparisons of this nature will need to be performed to be able to draw up an authoritative conclusion.

With regards to the key drought responsive protein and metabolic networks identified from the two studies, our iTRAQ proteomics analysis results revealed that cellular redox homeostasis maintenance associated- and stress defense proteins were the predominant contributors of maize drought tolerance through their participation in purine metabolism, proteasome and plant hormone signal transduction metabolic pathways [38]. However, our current DIA analysis results attribute maize drought tolerance to stress responsive chaperones, DNA replication related enzymes and seed storage (vicilin) proteins being involved in starch and sucrose metabolism, terpenoid backbone biosynthesis and carbon metabolism metabolic pathways. Although both studies similarly identify redox homeostasis maintenance and stress defense proteins (particularly peroxidases, oxidoreductases and chaperones) critical for maize drought tolerance, they showed varied metabolic pathways in which these proteins were involved. This variation could possibly be explained by the differences in the number of DAPs that were significantly enriched in different pathways. Considering that the DIA approach had fewer DAPs identified, as a result of its rigorous screening, it is possible that some proteins that were excluded from analysis could have obviously caused differences in the significant pathways enrichment of DAPs. Overall, both iTRAQ and DIA approaches are effective and vital tools for performing proteomics analysis of drought stress responses that can be used with confidence either singly or complementarily. Meanwhile, our findings derived from both previous [38] and current studies add to the body of knowledge related to, and offer useful reference data for, maize drought stress tolerance breeding.

Conclusions

The DIA-based analysis of the maize filling-kernel proteomic responses of two contrasting hybrid cultivars exposed to a 12-day drought treatment at the grain filling stage identified a total of 826 DAPs. Among them, we screened out 7 DAPs specifically accumulated in the tolerant genotype ND 476 and 45 DAPs uniquely expressed in the sensitive genotype ZX 978. Our analysis revealed that the maize drought tolerance is related to some key proteins such as chaperones, DNA replication related enzymes, flowering locus T protein and vicilin, ABA-responsive RAB17 protein. These proteins chiefly participate in starch and sucrose metabolism, terpenoid backbone biosynthesis, cyanoamino acid metabolism and carbon metabolism pathways. Our findings provide useful reference data for maize drought stress tolerance during the grain filling stage, as well as offer some theoretical basis for maize drought resistance breeding.

Conflicts of interest

The authors declared no conflict of interest. Furthermore, the funders played no part in designing, data collection, analysis and interpretation of results, writing, editing or decision making process related to the publication of this study.

Additional information

Funding

This research was supported by the National Key Research and Development Project of China (Selection and Efficient Combination Model of Wheat and Maize Water Saving, High Yield and High-Quality Varieties) (Grant No. 2017YFD0300901).