Metabolomic and Proteomic Analysis of Maize Embryonic Callus induced from immature embryo

The low ratio of embryonic callus (EC) induction has inhibited the rapid development of maize genetic engineering. Still, little is known to explain the genotype-dependence of EC induction. Here, we performed a large-scale, quantitative analysis of the maize EC metabolome and proteome at three typical induction stages in two inbred lines with a range of EC induction capabilities. Comparison of the metabolomes and proteomes suggests that the differential molecular responses begin at an early stage of development and continue throughout the process of EC formation. The two inbred lines show different responses under various conditions, such as metal ion binding, cell enlargement, stem cell formation, meristematic activity maintenance, somatic embryogenesis, cell wall synthesis, and hormone signal transduction. Furthermore, the differences in hormone (auxin, cytokinin, gibberellin, salicylic acid, jasmonic acid, brassinosteroid and ethylene) synthesis and transduction ability could partially explain the higher EC induction ratio in the inbred line 18-599R. During EC formation, repression of the “histone deacetylase 2 and ERF transcription factors” complex in 18-599R activated the expression of downstream genes, which further promoted EC induction. Together, our data provide new insights into the molecular regulatory mechanism responsible for efficient EC induction in maize.

. Overview of the maize callus and the metabolomic responses to induction medium in 18R and B73. (a) Photos of a single callus of 18R and B73 induced for 15 d. (b) Principal component analysis based on metabolite profiles of maize 18R and B73 inbred lines in 3 stages. "18R": 18R; "C": control; "I": stage I; "II": stage II; "III": stage III. (c) Numbers of individual mass spectrometry features that were significantly changed (fold change >1.5, Q value < 0.05) at each stage in each line. (d) Venn diagram of significantly up-and downregulated m/z in 18R and B73 at each stage. "I/C": stage I relative to control, as well as "II/C, " "III/C. " pathways. Among them, 56 metabolites were assigned to arginine and proline metabolism (Supplementary Table S2). It is well known that proline plays an important role in stress response. To some extent, immature embryos cultured on induction medium just undergo a type of stress (hormone stress). Moreover, 208 metabolites were annotated in carbohydrate metabolism including starch and sucrose metabolism, pentose and glucuronate interconversions, galactose metabolism and fructose and mannose metabolism. In addition, more than 170 metabolites were related to the lipid metabolism pathway (Supplementary Table S2). These results suggested that energy metabolism played an extremely important role in EC induction. Additionally, 50 metabolites related to signal transduction pathways were also differentially expressed between 18R and B73 (Supplementary Table S2). Biosynthesis of other secondary metabolites, metabolism of co-factors and vitamins, metabolism of terpenoids and polyketides, and xenobiotics biodegradation and metabolism pathways were enriched too (Supplementary  Table S2). These results suggested that the difference in EC induction capability between the two inbred lines was very complex and involved various metabolic pathways. To reveal how the significantly metabolites were regulated, proteomic analysis was then conducted.
Overview of the Proteomic Data Sets. Total proteins of CK, stage I, II and III of 18R and B73 were submitted to iTRAQ-based quantitative proteome analysis. A total of 5,797, 6,232 and 6,216 proteins were identified in replicates 1, 2 and 3, respectively. In this study, the 4,816 proteins identified in all 3 replicates were retained for further analysis (Supplementary Table S3). Relative to CK, the proteins with fold changes of >1.5 (or <0.67) and P values < 0.05 in any sample or in any line were classified as significantly changed proteins. In total, we identified 616 significantly changed proteins (Supplementary Table S4). In 18R, 232, 279 and 317 proteins were up-regulated in stage I, II, and III, respectively. In contrast, 168, 229 and 303 proteins were significantly down-regulated in stage I, II and III, respectively (Fig. 2a). In B73, however, 190, 230 and 235 proteins were up-regulated, and 112, 146 and 217 proteins were significantly down-regulated in stage I, II and III, respectively (Fig. 2a). At each stage, more proteins were either up-regulated or down-regulated in 18R than in B73. As shown in Fig. 2b, only 75 (22%), 128 (34%) and 118 (27%) proteins were significantly up-regulated at stage I, II and III in both 18R and B73, respectively. Similarly, only 48 (21%), 89 (31%) and 124 (31%) proteins were significantly down-regulated in three stages in both lines, respectively. These results suggested that 18R and B73 responded to induction medium by expressing different proteins, which could be responsible for the different EC induction capability.

Co-expression Networks of Significantly Changed Proteins.
To study the expression pattern, 616 significantly changed proteins were subjected to co-expression networks analysis using the WGCNA R package 14 . This unsupervised and unbiased analysis identified 14 and 12 distinct co-expression modules in 18R and B73, respectively (Fig. 2d, Supplementary Fig. S2a, and Supplementary Table S4). We then performed a modules-preservation test using repeated random sampling and modules recalculation. All the Zsummary scores of the modules of the two inbred lines were larger than 6 (Supplementary Table S5). According to previous study 15 , modules with Zsummary scores >2 were moderately preserved, which indicated that all of the modules in our study were conserved. The module preservation of 18R modules in the B73 data was calculated, which can also be regarded as a differential network analysis. The salmon, blue, green and purple modules of 18R with Zsummary scores of <2 cannot be recalculated in B73 data ( Supplementary Fig. S2b). Accordingly, 150 proteins in these 4 modules were selected as differentially expressed proteins (DEPs) between 18R and B73.
We also compared the B73 modules with the 18R modules and calculated the correlation of eigengenes value (the first principal component of the protein profiles) for each module between 18R and B73 (Fig. 2d). As expected, the proteins in most of the 18R modules overlapped with several different B73 modules, and most of the correlations resulted in very high values. Finally, 93 proteins with correlation of <0.4 were selected as DEPs, too. In total, 243 DEPs (Supplementary Table S4: red marked proteins) between 18R and B73 were chosen for further analysis.
EC Formation-Related DEP Annotation. Based on the eigengene dendrogram (Fig. 2c), the MEbrown and MEgreenyellow modules had very similar expression profiles, indicating that their genes were highly co-expressed; thus, these modules were merged and renamed as cluster 1 (C1). Similarly, MEblue, MEturquoise, MEblack and MEpink were merged and renamed as cluster 3 (C3), whereas MEred, MEyellow, MEgreen,  MEmagenta and MEtan were merged and renamed as cluster 5 (C5). MEsalmon and MEpurple were renamed as cluster 2 (C2) and cluster 4 (C4), respectively. As shown in Fig. 3a, the C1 DEPs in 18R were substantially up-regulated in stage II, and there were no significant changes in stages I and III. However, the DEPs in B73 increased at all 3 stages. These DEPs likely played an important role in the initial EC formation and were not required for stage I or III. In C2, DEPs were up-regulated at 3 stages in 18R and especially increased at stage III, whereas no significantly changes were noted in B73. The DEPs in C3 were sharply up-regulated during all stages in 18R, but not even down-regulated in B73. On the contrary, the DEPs in C5 decreased at stage I and retained a lower level at stages II and III in 18R, but did not significantly change in B73, which suggested that these DEPs are rapid and continuous response factors. The C4 DEPs did not change at stage I, but were extremely reduced at stages II and III in 18R. Still, they were sharply down-regulated at all stages in B73, which suggested that they were needed for embryo enlargement in 18R.
To understand the role of differentially expressed genes in EC, the 243 DEPs were submitted to Gene Ontology (GO) enrichment analysis. The top ten GO terms of each cluster are shown in Fig. 3b and c. These molecular functions and biological process GO terms include metal ion binding, oxidoreductase activity, protein binding, RNA binding, lyase activity, hydrolase activity, metabolic processes, response to stimulus and metal ions, oxidation-reduction processes, and organonitrogen compound metabolics (Supplementary Table S6).
Interestingly, a total of 69 DEPs in C1, C2, C3, C4 and C5 were enriched for metal ion binding or corresponding functions (P < 0.05, Table 1). Among them, four dehydrogenase family proteins were induced in all of the stages (Table 1). Likewise, three cytochrome family proteins were up-regulated, while two were down-regulated (Table 1). Six peroxidase super family proteins in C1, C2 or C3 (Table 1) were up-regulated in the process of EC formation. These DEPs may function by activating oxygen scavenging. In addition, three DEPs were glutathione S-transferase-like proteins, which are widely known to respond to stimuli or detoxification in scavenging of toxic organic hydroperoxides to protect plants from ROS damage 16 . Similarly, genes involved in energy metabolism were induced at all 3 stages (Table 1 and Fig. 3). Clearly, EC formation is a complex process that requires sufficient energy for fast growth and cell division 17 so that energy metabolism was already identified as important for maize and Arabidopsis callus induction 10 .
By GO analysis, 30 DEPs were classified into cell differentiation, division, communication and programmed cell death (Table 2). Although their P values were larger than 0.05, we selected them for further analysis because the process of EC induction does involve cell dedifferentiation, division and re-differentiation. Among them, 4 reduced and 4 increased DEPs were assigned to cell-cell junction or cell communication (Table 2). Meanwhile, 3 increased and 2 reduced DEPs were classified as cell differentiation, and 3 down-regulated DEPs as players in cell division (Table 2). Moreover, GRMZM2G408768_P01 was annotated as a cell cycle regulator, which was up-regulated at stage I in 18R but immediately reduced to normal levels at following stages (Table 2). However, the level of this protein did not increase at stage I and even decreased at stages II and III in B73.
In addition, the RNI-like protein superfamily (GRMZM2G393272_P02), including the SART-1 (GRMZM2G109284_P01) and Endochitinase PR4 proteins (GRMZM2G129189_P01) were assigned to the biological processes of meristem development, meristem initiation and somatic embryogenesis, respectively ( Table 2). These 3 DEPs of 18R were dramatically increased during all stages; however, they did not significantly change in B73, or the fold-change in B73 was smaller than in 18R (Table 2). Moreover, 4 up-regulated and 5 down-regulated DEPs in 18R (Table 2) were predicted as plant-type cell wall related genes, which were required during cell division.
KOBAS (http://kobas.cbi.pku.edu.cn/home.do) was used to perform pathway analysis. Similar to GO analysis, energy metabolism-related pathways were identified (Table 3), which suggested that EC induction required a lot of energy. As expected, amino acid synthesis and metabolism pathways were also enriched. In addition, hormone-related pathways were identified, such as plant hormone signal transduction, brassinosteroid biosynthesis, alpha-linolenic acid metabolism (which participates in jasmonic acid biosynthesis and zeatin biosynthesis), cysteine and methionine metabolism (which participates in ethylene biosynthesis) and tryptophan metabolism (which is involved in auxin biosynthesis) ( Table 3).

Plant Hormone-Related DEPs Induced by Auxin-containing Medium. Both the GO and KEGG
pathway annotation suggested that hormone signaling transduction-related pathways were activated. To identify whether the proteins are regulated by hormones, the Hormonometer program was used to compare the level of protein fold changes with indexed data sets treated by hormones 18 (http://hormonometer.weizmann.ac.il/hormonometer). In Arabidopsis (Arabidopsis thaliana), the gene expression profiles treated by methyl jasmonate, ethylene, abscisic acid, indole-3-acetic acid, cytokinin (zeatin), brassinosteroid, GA3 and salicylic acid have already been completed. To carry out this analysis, only the proteins having its Arabidopsis homologs and containing probe set identifiers (required by Hormonometer) were kept. A total of 2,750 Arabidopsis homologs of maize proteins (Supplementary Table S7) were submitted to Hormonometer analysis.
As expected, during EC formation changes of protein expression were associated with hormone signaling (Fig. 4a). Based on the dendrogram, protein expression changes at stage I were distinct from those observed at stages II and III, and changes in 18R and B73 were very comparable at each stage (Fig. 4). Moreover, there was a strong positive correlation between the proteins induced by hormone-containing medium and the Arabidopsis homologous genes induced for 0.5, 1 and 3 h with jasmonic acid and for 3 h with salicylic acid. Similarly, a positive correlation was observed between the expression changes of maize protein induced by hormone-containing medium and the Arabidopsis genes induced for 1 and 3 h by auxin. Notably, at stage II and III, the correlation of 18R at 3 h auxin treatment was enhanced compared to that of B73. A similar positive correlation was identified after cytokinin treatment for 3 h, and the correlation in 18R was enhanced compared to that of B73 for stage I. On the contrary, there was a significant negative correlation between maize proteins and Arabidopsis genes treated with ethylene for 0.5 h and 1 h, brassinosteroid for 0.5 h and 1 h, GA3 for 3 h, 6 h and 9 h. Moreover, a positive correlation, especially at stage I, was observed for maize proteins induced by induction medium and Arabidopsis genes treated with ABA for 1 and 3 h. However, a strong negative correlation was identified after ABA treatment for 0.5 h at all 3 stages in both 18R and B73. During maize leaf development, 431 expressed genes are identified to be relevant to hormone synthesis, transduction and responding 19 . In the present study, 83 expressed proteins were involved in hormone related pathways (Fig. 4b). During callus formation, 18R and B73 show different abilities in synthesis, transduction or responding of Jasmonic acid (JA), ethylene (Eth), cytokinin (CK), salicylic acid (SA), auxin, and brassinosteroid (BR). Different Hormone Biosynthesis Patterns in 18R and B73. The phenylalanine, tyrosine and tryptophan biosynthesis pathways are involved in the synthesis of chorismate, which is necessary for auxin and salicylate synthesis in plants (Fig. 5). Compared to the control, chorismate levels were significantly increased in 18R at all 3 stages, but not in B73 (Fig. 5). Furthermore, the levels of anthranilate, tryptophan, tryptamine, N-hydroxyl-tryptamine, indole-3-acetaldehyde oxime and indoleacetaldoxime dehydratase (EC: 4.99.1.6, GRMZM2G167549_P05) were all increased at least at one stage in 18R, but not in B73. These changes led to the increase of indole-3-acetic acid (auxin) in 18R, but no changes in B73. Moreover, in the auxin transduction pathway, a GH3 gene family protein, GRMZM2G410567_P01, was significantly up-regulated at 3 stages in both lines, but changed more in 18R than in B73. In the salicylic acid synthesis pathway (Fig. 5), prephenate, L-arogenate, trans-cinnamate, benzoate and tyrosine and arogenate dehydrogenase (EC: 1.3.1.78, GRMZM2G084942_P01) were all up-regulated at least in one stage in 18R, but with no significant changes in B73, which led to a significant increase of salicylate in 18R at stages II and III, but not in B73.
Different concentrations of cytokinin are usually added for in vitro callus induction. Hormone-regulated protein analysis showed that several proteins were regulated by cytokinin (Fig. 4). Therefore, we also focused on the cytokinin synthesis pathway (Fig. 7a). As expected, cytokinin trans-hydroxylase (CYP735A, GRMZM2G014395_ P01), trans-zeatin riboside, zeatin (cytokinin), and dihydrozeatin (cytokinin) were all significantly up-regulated at all stages in 18R but were unchanged in B73, which suggested that 18R and B73 have different cytokinin synthesis capabilities.
In contrast to these hormones, the brassinolide synthesis and transduction pathway were repressed (Figs 4 and 7b). The levels of 6-α-hydroxy-castasterone and brassinolide were reduced at all 3 stages and peaked at stage II in 18R, but did not change in B73. Moreover, the up-regulation of CYP734A1 (PHYB activation tagged suppressor 1, EC:1.14.-.-, GRMZM2G014395_P01) in 18R led brassinolide to be modified to 26-hydroxy-brassinolide, which further lowered its concentration, whereas there were no changes in B73. Furthermore, BRI1 (EC: 2.7.10.1, GRMZM2G066274_P01) significantly decreased in the brassinolide transduction pathway at all stages in 18R  Table 3. Pathways annotated from all DEPs.
but not in B73 imply its function in meristem development (Table 2). Additionally, GRMZM2G109284_P01, a SART-1 family protein, sharply increased at all 3 stages in 18R (Table 2); however, no changes were detected in B73. Arabidopsis MDF (belonging to SART-1 family protein), function in stem cell formation and meristematic activity 29 .    with no significantly changes in B73. However, BAK1 (GRMZM2G468756_P01, EC: 2.7.10.1) was significantly up-regulated at stage I in 18R and at all stages in B73. These results suggested that the brassinolide synthesis and transduction pathway were suppressed in 18R but not in B73. Likewise, ethylene biosynthesis was also suppressed during EC induction. As shown in Fig. 6b, except for S-adenosyl-L-methionine, which level was more reduced in B73, the level of other metabolites including L-methionine, 1-aminocyclopropane-1-carboxylate and ethylene decreased more in 18R than in B73 during EC induction. Ethylene response factors (ERF) function in ethylene-mediated signaling transduction 20 . In Arabidopsis, AtHD2A can directly interact with tomato Pti4 transcription factor (ERF family proteins) to further inhibit the expression of downstream genes by binding to the GCC-box region of promoter 21 . Similarly, Longan fruit DIERF1 also can directly interact with DIHD2 22 . HD2 is a plant-specific histone deacetylase and involved in many processes, such as seed dormancy 23 , ABA and salt stress response 24 .
Interestingly, relative to control, HDT2 (GRMZM2G100146_P01, homologous to DIHD2) and HDT3 (GRMZM2G159032_P01, homologous to DIHD2) significantly decreased both in RNA level (Fig. 8c) and protein level at all stages in 18R (Fig. 8a). Thus, we hypothesize that the ERF/HDT complex decreased because of the reduction of HDT2 or HDT3, which further led to the up-regulation of downstream genes. As expected, 21 DEPs in the blue module, with a highly negative correlation (<−0.7) to HDT2 and HDT3 at the RNA level, have at least one GCC-box in their promoter (3000 bp upstream of ATG, Fig. 8a). This suggested to us that these proteins might be directly regulated by the ERF/HDT complex. The maize homologous genes of the downstream genes of AtERF6 (homolog of DIERF1) in Arabidopsis 25,26 were also sharply up-regulated during EC formation (Fig. 8b), and half of them contained more than 1 GCC-box in their promoter region. The blocked ethylene synthesis and reduced ERF/HDT complex might lead to the up-regulation of downstream genes and additional regulation of EC induction.

Discussion
To investigate the dynamic changes of proteins and metabolites during maize EC formation cultured on hormone-containing induction medium, it is necessary to consider multiple induction stages and different inbred lines with different EC induction capabilities. Most previous studies of proteomic changes focused on embryonic callus and non-embryonic callus after induction for 4 weeks or even later in only one inbred line 9,10 . However, our metabolic and proteomic results showed that a larger range of metabolites and proteins changed significantly at early stages (Fig. 1b). We also used a genetic approach by using two maize inbreds with different EC formation capabilities and found that the two different inbred lines showed geneotype-specific responses to the induction medium. Furthermore, complementing a proteomic with a metabolomics analysis provides us with a better understanding of the maize embryonic callus induction process.
Indeed, B73 and 18R showed different cell division properties. In this study, tubulin folding cofactor A (GRMZM2G143627_P01) sharply decreased at stage I in 18R but increased at stage II in B73. Arabidopsis mutations in tubulin folding cofactor failed to form microtubules, which are essential for cell division and organized cell expansion 27 . Functional loss of these genes leads to impaired cell division with normal cell growth, which further causes formation of grossly enlarged cells. So, the sharply reduced expression of GRMZM2G143627_P01 at stage I might directly promote embryo elongation while also maintaining cell division activity.
Similarly, the abilities of meristem development, meristem initiation and somatic embryogenesis differed between 18R and B73. In Arabidopsis, despite the functional loss of ATMAX1 and ATMAX2, the first axillary shoots remain active during shoot formation 28 , which suggests that they repress axillary meristem differentiation. Thus, as an ATMAX2 homologous protein, the up-regulation of GRMZM2G393272_P02 at stages I and III in 18R Scientific RepoRts | 7: 1004 | DOI:10.1038/s41598-017-01280-8 The mdf mutant lost stem cells and lacked meristematic activity both in the root and vegetative shoot. Thus, the differences in expressions of GRMZM2G109284_P01 between 18R and B73 might be responsible for promoting callus (or stem cell) formation. Besides, Arabidopsis ATCHITIV, encoding an EP3 chitinase, was expressed in the cells surrounding the somatic embryo during somatic embryogenesis 30 . In the present study, GRMZM2G129189_P01, a gene homologous to ATCHITIV, was significantly up-regulated at all stages in 18R. However, it was increased only at stage III in B73 (Table 2). So, this gene might also play an important role in somatic embryogenesis.
The expressions of AtFLA11 and AtFLA12 were positively related to secondary cell-wall cellulose synthases in stems 31 . A homologous gene of AtFLA11, GRMZM2G022931_P01, increased sharply at stage II and stage III in both 18R and B73, whereas the fold change in 18R was higher than in B73 (Table 2), which might function in cell-wall cellulose synthases. During cotyledon callus induction, the outer cell wall was degraded. New cell walls were re-synthesized and cell wall synthesis related proteins were also reduced 16 . Cell wall loosening and/or remodeling is a necessary foundation for cell elongation and might be regulated by various hormones 32 .
It was reported that 14-3-3 proteins bind to pS/pT motifs to regulate several signal transduction pathways, such as ethylene 33 and brassinolide 34 . In Arabidopsis, 14-3-3 proteins interact with ethylene-overproducer1 (ETO1)/ETO1-like (EOLs) to rapidly degrade 1-aminocyclopropane-1-carboxylate synthase 33 to inhibit ethylene synthesis. As a 14-3-3 like protein, GRMZM2G408768_P01 was annotated as involving in cell cycle Figure 8. Expression pattern of genes potentially regulated by histone deacetylase. (a) Heatmap of proteins in the 18R blue module that show a significantly negative correlation with HDT2 (GRMZM2G100146) and HDT3 (GRMZM2G159032). No. of GCC box is the number of "GCCGCC" sequences in the promoter (3000 bp upstream of ATG). The "GCC box" was reported to be bound by ERF to regulate gene expression. ERF can directly interact with histone deacetylases to further function. (b) Heatmap of transcription changes of potential downstream genes of ERF and histone deacetylase at the 3 stages of 18R. (c) Relative expression level of HDT2 and HDT3 are shown as means ± SD (n = 3) during callus formation of 18R. *P < 0.05 (Student's t test) for differences between control and each stage. regulation and significantly increased at stage I in 18R and at all stages in B73 (Table 2), which might result in the down-regulation of ethylene (Fig. 6b) to further regulate the downstream genes (Fig. 4). Similarly, in the brassinolide transcription pathway, 14-3-3 proteins were found to bind to phosphorylated BZR1/BZR2, which efficiently retain proteins in the cytoplasm to further inhibit brassinolide levels 34 . Thus, GRMZM2G408768_P01 might act as a rapid response factor to ethylene and brassinolide synthesis at stage I to regulate the embryo expansion by controlling cell cycles.
Both proteomic and metabolomic analyses (Figs 3 and 4, Table 3, Supplementary Table S2) indicated that hormone synthesis and transduction, including auxin, cytokinin and brassinosteroids, play a role in EC induction. Combined with our previous study, the transcription levels of auxin influx carrier AUX1, auxin receptor TIR1, auxin-responsive protein IAA, ARF (auxin response factor) and GH3 (auxin responsive GH3 gene family) significantly changed during 18R EC induction 13 , which implies that the auxin transduction pathways were activated. Moreover, one GH3 protein (GRMZM2G410567_P01) was sharply upregulated at all 3 stages in both maize lines (Fig. 5b). In addition, proteins expressed in 18R showed a higher positive correlation with auxin induction than B73 (Fig. 4). Altogether, 18R showed a stronger capability of auxin synthesis and transduction when cultured on auxin-containing medium.
It is widely accepted that the balance of auxin and cytokinin content in the cell is critical for callus induction in different plants 35 . Auxin is involved in cell proliferation and elongation, whereas cytokinin promotes cell proliferation and differentiation. The cytokinin transduction pathway was also activated during EC induction of 18R 13 . The present study indicated that the production of zeatin and dihydro-zeatin (cytokinin) was increased in 18R but without a change in B73 ( Fig. 7a and b). Taken together, these findings suggest that the cytokinin synthesis in 18R is more critical for EC induction than in B73.
In response to EC induction, the synthesis of jasmonic acid was enahnced in 18R compared to B73 (Fig. 6). Furthermore, the proteins expressed during callus formation showed a significantly positive correlation with the treatment of jasmonic acid, which suggests that these proteins are regulated by jasmonic acid. Jasmonic acid is involved in biotic and abiotic stress responses, and plays especially a key role in the transduction of wound signals 36 . In addition to auxin, wound signals are generally considered to be another signal for callus induction 37 . A previous study showed that adding jasmonates to the medium could improve the callus formation ratio of Vigna mungo 38 . Combined with the findings in the present study, the synthesis and transduction of jasmonic acid in 18R seems to partially contribute to its higher EC formation ratio than B73.
Based on discussed above, we hypothesized an EC induction mechanism ( Supplementary Fig. S3). In summary, maize inbred line 18R and B73 show different stress response, energy supply, stem cell initiation, stem cell development, somatic embryogenesis, hormone synthesis, hormone transduction and hormone response abilities. All of these differences may be responsible for the high EC induction ratio of 18R.

Methods
Plant materials and tissue culture conditions. Maize inbred lines 18R and B73 were planted in a greenhouse (14/10 h light/dark, at 28 °C and 70% relative humidity). Immature embryos were harvested on the twelfth day after self-pollination. Modified N6 medium 13 was used to induce embryonic callus, and the embryos were then incubated for 15 days in darkness at 28 °C. Immature embryos induced for 0 d were used as the control (CK). During 0-15 d, 500 mg calli were sampled every 24 h in 9 biological replicates (3 for protein extraction and 6 for metabolite isolation), quickly frozen in liquid nitrogen and stored at −80 °C. Metabolome analysis. The low-molecular-weight metabolites (<1000 Da) were isolated according to Decourcelle et al. 39 . Briefly, 100 mg of embryos or calli were ground to a powder in liquid nitrogen. Then, 1 ml of methanol (HPLC grade) was added and incubated at 70 °C for 15 min, followed by centrifugation for 10 min at 14,000 g. Finally, the supernatant was harvested. The metabolites from 1-5 d were equally mixed to form sample stage I (biological replicates 6 times); likewise, stage II (6-10 d), stage III (11 to 15 d) and CK were prepared. Then, 10 μL of these samples were submitted to HPLC-MS analysis. The MS data were analyzed using XCMS and CAMERA software 40 . To obtain high quality data, robust Loess signal correction (R-LSC) 41 was used to normalize peaks to the QC samples, and the relative standard deviation (RSD) was calculated. The m/z with RSD value between 0 and 30% was subjected to principal component analysis (PCA), partial least-squares discriminant analysis (PLS-DA) and variable importance in the projection (VIP) analysis. Finally, the m/z with VIP > 1, fold change >1.5 (or <0.67) and Q value < 0.05 in any sample was retained for further metabolite identification and pathway annotation using the online HMDB database (http://www.hmdb.ca/), METLIN (http://metlin. scripps.edu/) and KEGG (www.genome.jp/kegg/) based on the differences between observed and theoretical mass smaller than 10 ppm.

Protein extraction, iTRAQ labeling and LC-ESI-MS/MS analysis.
Total protein was extracted as previously reported 42 . Then, 20 μg of the total proteins from 1 to 5 d was mixed to form samples of stage I (biological replicates 3 times); likewise, stage II (6-10 d), stage III (11 to 15 d) and CK (0 d) were prepared. Then, each mixed sample was digested with Trypsin Gold (Promega, Madison, WI, USA) at 37 °C for 16 hours and labeled according to the manufacturer's protocol for the 8-plex iTRAQ reagent (Applied Biosystems). Four milliliters of buffer A (25 mM NaH2PO4 in 25% ACN, pH 2.7) was added to reconstitute the iTRAQ labeled peptide. Then, the peptides were eluted on a LC-20AB HPLC Pump system (Shimadzu, Kyoto, Japan) with a 4.6 × 250 mm Ultremex SCX column at a flow rate of 1 mL/min with buffer A for 10 min, 5-60% buffer B (25 mM NaH2PO4, 1 M KCl in 25% ACN, pH 2.7) for 27 min, 60-100% buffer B for 1 min, 100% buffer B for 1 min and buffer A for 10 min. The absorbance of the eluent was measured at 214 nm, and 20 fractions were obtained, desalted and vacuum-dried for further MS analysis.
Each fraction was re-suspended in buffer A (2% ACN, 0.1% FA) and centrifuged at 20,000 × g for 10 min. The supernatant was loaded onto the LC-20AD nanoHPLC (Shimadzu, Japan) and then separated at 300 nL/ min starting from 2 to 35% B (98% ACN, 0.1% FA), followed by a 2 min linear gradient to 80%, 4 min at 80% B and finally returning to 5% in 1 min. After nanoelectrospray ionization, the peptides were submitted to tandem mass spectrometry (MS/MS) in the QEXACTIVE (Thermo Fisher Scientific, San Jose, CA) coupled online to the HPLC, followed by detection in the Orbitrap. High-energy collision dissociation (HCD) of the MS/MS was used to select peptides with a normalized collision energy setting of 27.0; ion fragments were detected in the Orbitrap at a resolution of 17,500. A data-dependent procedure that alternated between one MS scan followed by 15MS/ MS scans was applied for the 15 most abundant precursor ions above a threshold ion count of 20,000 in the MS survey scan with a following Dynamic Exclusion duration of 15 s. An electrospray voltage of 1.6 kV was applied. Automatic gain control (AGC) was used to optimize the spectra generated by the Orbitrap. The AGC target for full MS was 3e6 and 1e5 for MS2. For MS scans, the m/z scan range was 350 to 2,000 Da. For MS2 scans, the m/z scan range was 100-1,800.
iTRAQ protein identification and quantification. MS/MS data were searched using Mascot (2.3.02 version, Matrix Science, London, United Kingdom) against the Gramene database (Zea_maize, 105722 sequences) (ftp://ftp.gramene.org/pub/gramene/release42/data/fasta/zea_mays/pep/) to identify proteins. The search parameters were: mass tolerance of 20 ppm was permitted for intact peptide masses and 0.05 Da for fragmented ions; one missed cleavage was permitted for the trypsin digests; the potential variable modifications was Gln->pyro-Glu (N-term Q), oxidation (M) and deamidated (NQ), and fixed modifications were carbamidomethyl (C), iTRAQ8plex (N-term) and iTRAQ8plex (K); the charge states of peptides were set to +2 and +3; only peptides with Mascot probability greater than "identity" (95% confidence) were counted as identified. For protein quantitation, the proteins with at least two unique spectra in at least two replicates were kept for further analysis. The quantitative protein ratios were weighed and normalized by the median ratio in Mascot. Only the ratios with p-values < 0.05 and fold changes >1.5 in the 3 replicates were considered to be significantly changed 43 . The proteins that significantly changed in any inbred line or in any samples were selected for further analysis.
Weighted gene co-expression network analysis. All of the selected proteins from 18R and B73 were independently subjected to signed hybrid weighted correlation network analysis through the R package WGCNA 14 as previously reported 44 , using publically available scripts (http://labs.genetics.ucla.edu/horvath/ CoexpressionNetwork/Rpackages/WGCNA/). Briefly, pairwise correlation matrices based on Pearson correlations between all pairs of proteins were constructed. Next, the power β of 12 (soft-threshold of the correlation matrix) was used to construct an adjacency matrix to raise the correlation matrix and suppress noisy small correlations. Based on the results of the adjacency matrix, the topological overlap was calculated to measure the network connections of gene pairs. According to the topological overlap, hierarchical trees with average linkage were built and cut with the Dynamic Tree Cut algorithm to group genes with highly similar co-expression relationships into the same modules. The first principal component (module eigengenes) of each module was used to represent its expression level.

Module preservation analysis.
To determine whether our modules were conserved, 20 permutations of the original data were subjected to module preservation statistical analysis to calculate Zsummary scores using the R function "modulePreservation" in the WGCNA package 14 . Then, the averaged Zsummary scores combined with SD (20 replicates) were used to summarize whether a module was preserved or not. As argued in the previous study, modules with Zsummary scores of >10 are strongly preserved. On the contrary, if Zsummary scores of <2, the modules are not preserved. Zsummary scores between 2 and 10 represent modules that are moderately preserved.
Identifying differentially expressed proteins between 18R and B73. To identify the DEPs between 18R and B73, the following two methods were used. First, the significantly changed proteins of B73 were used to calculate the "module preservation" of 18R as illustrated above. The Zsummary scores here served as a measure to which each of the proteins of B73 was preserved in the 18R network. The modules with Zsummary scores below 2 were regarded as different modules; correspondingly, the proteins in these modules were identified as the DEPs between the two samples.
Secondly, we related B73 modules to 18R modules using WGCNA to obtain the correspondence of B73-specific and 18R-specific modules. Then, we used the eigengenes to represent the modules to calculate the Pearson correlations for each pair of modules between the two inbred lines. The corresponding proteins with correlations of <0.4 were also considered to be DEPs between 18R and B73.
RNA isolation and quantitative real-time PCR. Total RNA was isolated using Trizol Reagent (Life Technologies, Carlsbad, CA, the United States). RNAs from 1-5 d (sample I) were mixed with equal proportion, as well as samples from 6-10 d (sample II) and 11-15 d (sample III). The RNAs extracted from 0 d were used as control. The reverse transcription was performed as the protocol of PrimeScript RT Reagent Kit with gDNA Eraser (Takara Biotechnology Co., Ltd., Dalian, China). The qRT-PCR was performed on ABI 7500 real-time PCR System (Applied Biosystems, USA). Actin 1 (GRMZM2G126010) was used as a reference gene. The primers for HDT2 (GRMZM2G100146) amplification are F:GATTCTGATGATTCTGATGAAGGCG and R: AGAAAGAGGCGTTTTCAGAGCAT. The primers for HDT3 (GRMZM2G159032) quantitative real time PCR are F:CGGACGATTCTGATGAGGGTT and R:TAGAGGAGTTTTCAGCACGGAAC.