Transcriptome and metabolomic analyses reveal regulatory networks and key genes controlling maize stomatal development in response to blue light


 Background: Blue light is helpful in the formation of maize stomata, but the signal network is still unclear. We replaced red light with blue light in an experiment, and provided a complementary regulatory network for the stomatal development of maize by using transcriptome and metabolomic analyses. Results: The blue light led to 1296 differentially expressed genes and 419 different metabolites. Transcriptome comparisons and correlation signaling network analysis detected and identified 55 and 6 key genes, respectively, which regulate the responses of stomata to the substituted blue light in environmental adaption and signaling transduction pathways. Two genes involved in carotenoid biosynthesis and starch and sucrose metabolism participate in stomatal development as identified by transcriptome and metabolomic analyses. Conclusions: Blue light remarkably changed the network signal transduction and metabolites, and the new genes played important role in the formation of maize stomatal development network.


Background
The stomata is a specialized epidermal cell structure distributed on the surface of the leaves of all terrestrial plants. It is an important channel for the exchange of oxygen, carbon dioxide, water and heat between plants and the environment. The morphology, distribution, movement and development of stomata directly determine plant physiological activities such as photosynthesis, respiration and transpiration. Under changeable environmental condition, plant sense outside signals including gas, water, heat, and light, and collaborate with endogenous factors by nely adjusting the stomatal number to optimize the exchange of gas, water and heat to adapt to the environment. However, the unpredictability of environmental factors and the complexity of reacting internal signals bring great challenges to the research on the regulated signaling network of stomatal formation. The cell fate, structure formation, and cell cycle of the leaf stomatal lineage cells of monocotyledons, such as Gramineous maize, are quite different from those of dicotyledonous plants [1].
A large and complex regulatory signaling network controls the differentiation and distribution of stomatal cells. The main function of this network is to divide the meristem and precursor cells on time to achieve the cell cycle of stomatal lineage cell. However, the in uence of environmental factors on the stomatal formation regulatory network is di cult to explain because of the complexity of meristem development and differentiation. Light quality is one of the important factors that mediate stomatal development in plants among many environmental factors and plays an important role in the difference in the development mechanism between monocotyledons and dicotyledons [2]. As an exogenous signal, blue light is the direct initiator of signal network regulation and a key environmental factor that determines the interaction signal at other levels [3,4]. External and internal factors must be considered to explain the regulation mode of blue light on the stomatal development in maize. The regulation process of blue light on stomatal development is closely related to the constitutively photomorphogenic 1 (COP1) -elongated hypocotyl_5 (HY5) -phytochrome interacting factor (PIF) pathways regulated by cryptochrome (CRY) [5].
Blue light dominantly affects the binding of erecta (ER) -too many mouths (TMM) complex of epidermal patterning factor (EPF) transcription factor in EPF1/2 negative regulation and STOMAGEN positive regulation in maize [6]. As another signaling cascade, the mitogen-activated protein kinase (MAPK) module, which integrates internal and environmental signals, is the core of blue-mediated regulation of in the stomatal development network [2]. Blue light participates in the direct signal transduction and indirect regulation of COP1-HY5-PIF and MAPKs and thus participates in regulation networks. Many other protein signals mediated by blue light are involved in stomatal lineage differentiation. The adjustment of blue light nally affects the adjustment of terminal molecular switches. The regulation mechanism of cell cycle repeats the stomatal development process [7]. Therefore, blue light plays a particularly important role among the many environmental factors that affect the stomatal development of plants.
The known pathway is not enough to understand the whole network of stomatal development in blue light. After all, different modules regulated by environmental factors and endogenous metabolic modulators in stomatal development form a complex regulatory network. Blue light provides energy for photosynthesis, and in uences the growth and development of plants as an environmental signal.
Consequently, plants are able to fully adapt and utilize the regulated information of blue light. However, the photomorphogenesis and photochemical processes of plants are species-speci c for C4 plants in Gramineae. The mechanism of blue light signals should be further studied systematically. This project intended to reveal issues, such as intercellular signal transduction, the network and location of signaling cascades, and the regulated direction of metabolism in maize from the analyses of transcriptome and metabolome. The purpose of this experiment was to nd and identify putative functional genes that affect the stomatal development of maize from the aspects of light-dependent signal transduction and key signal cascades.

Gene expression pro le in light treatments
The number of gene expression and the distribution of gene expression value in samples have certain differences; thus, sample expression value (FPKM) was divided into different intervals, the number of genes expressed in the samples of different expression intervals was calculated, and the stacking histogram was drawn for display. The gene expression distribution map of the samples is shown in Fig. 1A. The correlation analysis of all transcriptome samples showed a remarkable correlation between biological replicates (Fig. 1B). The principal component analysis (PCA) of total unique genes showed that genes in each treatment were homogenously expressed at PC 1 with 59.1% variance (Fig. 1C). The clustering method was used to calculate the distance between samples to investigate the similarity between samples. Samples that belong to the same treatment are close in distance and preferentially clustered together. The cluster diagram of sequencing samples obtained according to gene expression is shown in Fig. 1D. The genes expression pro les indicate that blue light brings a large number of differentially expressed genes (DEGs).

DEGs in maize that respond to blue light
The summary of transcriptome analysis on maize leaf samples is listed in Table S2. In the RNA-seq libraries, more than 8 G data were obtained with GC > 56.23 % and Q30 > 93.92 %. approximately 92.85 % (blue) and 92.97% (red) clean reads were mapped to the reference genome of maize. In total, about 89.25 % (blue) and 88.97 % (red) of the unique genes were uniquely mapped to the reference genome. A total 1296 genes were differentially expressed, in which 890 were up regulated and 406 were down regulated when the light was changed from red to blue ( Fig. 2B and C). The DEGs were subjected to unsupervised hierarchical clustering. Genes in the same cluster may have similar biological functions for further analysis ( Fig. 2A).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the DEGs
GO analysis was carried out to analyze the functions of the DEGs that respond to the change in light condition. The summary of GO categories is depicted in Fig. 3A. According to the GO classi cation, DEGs involved in biological process occurred after light condition change. A total 308 genes were down regulated in the level 2 GO analysis. 198, 193, and 135 genes were involved in the top three GO terms under biological process, namely, "cellular process," "metabolic process," and "single-organism process," respectively; 224, 223, and 154 genes were involved in the top three GO terms under cellular component, namely, "cell," "cell part," and "organelle," respectively; and 173, 151 and 37 genes were involved in the top three GO terms under molecular function, namely, "binding," "catalytic activity," and "nucleic acid binding transcription factor activity," respectively. In total, 667 genes were up regulated in GO analysis. For biological process, 375, 367, and 328 genes involved in "cellular process," "metabolic process," and "single-organism process" were up regulated, respectively. For cellular component, 500, 500, and 375 genes related to "cell," "cell part," and "organelle" were up regulated, respectively. For molecular function, 339, 381, and 57 genes related to "binding," "catalytic activity," and "transporter activity" were enriched, respectively. The results of GO analysis were used for the functional analysis of the screened genes ( Fig. 3A).
KEGG analysis showed that DEGs were signi cantly enriched in "plant hormone signal transduction," "MAPK signaling pathway," "glycerophospholipid metabolism," "starch and sucrose metabolism," and "circadian rhythm" (Fig. 3B). The processes of "signal transduction" and "environment adaptation" of maize leaf stomata to blue light were collected. According to KEGG pathway identi cation, 53 genes were classi ed into four major enrichment classes, namely, "plant hormone signal transduction," "MAPK signaling pathway," "plant pathogen -interaction" and "circadian rhythm -plant," and six less enriched classes, namely "inositol phosphate metabolism," "phosphatidylinositol signaling system," "endocytosis," "ubiquitin mediated proteolysis," "amino sugar and nucleotide sugar metabolism," and "protein processing in endoplasmic reticulum." The enrichment term and annotation of these 53 DEGs are listed in Table 1. Gene annotations were retrieved from National Center for Biotechnology Information (NCBI) and InterPro. Four DEGs (LOC100281912, LOC107326007, LOC107403167, and LOC100192073) were not fully annotated in the KEGG analysis. Table 1 The list of 53 genes classi ed by "signal transduction" and "environmental adaptation" in KEGG. The correlation network of 53 DEGs in "signal transduction" and "environment adaptation" (according to KEGG classi cation level 2) were conducted in all transcriptome samples to fully understand the expression of the genes responding to light change in maize seedling stage (table 1). The PPI network of DEGs was constructed and classi ed by MCODE according to the correlation of DEGs ( Fig. 4A-F). Three nodes and 25 genes were obtained under the KEGG pathway, "environment adaptation" (Fig. 4A-C). 27 genes within three nodes were obtained in the "signal transduction" KEGG pathway ( Fig. 4D-F). The differences in the transcription expression of these 55 genes are shown in the heat map (Fig. 5). Three repetitive genes were deleted, and 24 genes unrelated to KEGG the pathway but interacted with the 28 annotated genes in Table 1 were nally obtained. 22 of the 24 genes can be annotated from NCBI and GO ( Table 2). The other two genes, namely, LOC100280217 and LOC100285849, were unidenti ed.  4G). LOC100285849 was identi ed as the CIPK16 (98.7 with CIPK-like 1) of the CBLinteracting protein kinase family. Its signaling pathway involved in stomatal development probably started with the catalysis of gamma-phosphoryl on protein substrates and/or calcium signals triggered by light. LOC107403167 was deduced to be a member of serine/threonine protein kinases. The protein sequence was closed, but the matching degree with the OXI1 of maize is low (query cover = 58%); thus, its gene function cannot be determined from existing results. LOC100192073 was known as the PP2C14 of the protein phosphatase 2C family. BLASTP results revealed that LOC100192073 and PP2C30 were closed, but an additional NADB_Rossmann protein domain was found in putative maize PP2C30 (Fig. 6). Therefore, LOC100192073 is not PP2C30 but a member of the PP2C family. LOC100280217 with PhrB conserved domain interacted with HY5 and was identi ed as (6 − 4) DNA photolyase (query cover = 92%), similar to UVR3_1, but was not reported in maize. Therefore, the effect of (6 − 4) DNA photolyase on the signal pathway of stomatal development in maize is not clear. LOC100281912 was close to the TIFY 11c (Per. Ident = 98.04%) of maize and other Gramineae plants but has two differences and more than two amino acids in length. Noticeably, the BLASTP results of LOC100281912 was also close to the pnFL-2 (Per. Ident = 98.03%) of maize.
2.5. Metabolite pro le and screening of differential metabolites in maize in respond to light change A total of 6272 peaks were detected by liguid chromatography-tandem mass spectrometry LC-(MS/MS) analysis, and 3204 metabolites were annotated in MS2 database. Among them, 1383 metabolites were obtained by negative ion scanning modes and 1821 metabolites were obtained by positive ion scanning modes (Fig. 7B). All metabolites were analyzed by principal component analysis and hierarchical cluster analysis determine the in uence of blue light on the change of metabolites ( Figure S1). A clear separation was found between the two sampling responses to light change. The volcanic diagram was used to screen different metabolites by visualizing the P value and fold change value from the 6272 peaks (Fig. 7A). Multi-dimensional analysis and one-dimensional analysis were used to screen different metabolites between groups after the Student's test and fold change analysis of the 3204 detected metabolites. As a result, 419 different metabolites were obtained by orthogonal partial least squares discrimination analysis (OPLS-DA) (Fig. 7B).
Hierarchical clustering and visual analysis were conducted on the expression levels of the top 50 differential metabolites according to the variable importance of projection (VIP) values to show the relationship and the expression differences of metabolites between samples more intuitively (Fig. 8A). Correlation analysis was helpful to measure the degree of correlation and the relationship between signi cantly different metabolites. Pearson correlation coe cient was used for correlation analysis, and Pearson product difference correlation coe cient was used to measure the linear correlation between two quantitative variables (Fig. 8B). The differential metabolite expression amount was visualized according to VIP value. Figure 8B shows the correlation of the top 50 differential metabolites.

Effect of light change on metabolic pathways and gene screening in the joint analysis of transcriptome and metabolome
The function and interaction of metabolites were constructed by using the KEGG database. The enrichment analysis of differential metabolites by analyzing metabolite-related metabolic pathways is helpful to understand the mechanism of metabolic pathway changes in different samples. The analysis of the differential accumulation of metabolites in stomatal development showed that 65metabolites accumulated differentially after 12 h of blue light change. Among them, 31 were up regulated and 34 were down regulated. During the light change from red to blue, the key metabolites affecting stomatal development, such as those under the "carotenoid biosynthesis" and "starch and sucrose metabolism" pathways, were analyzed (Fig. 9A).
The KEGG Markup Language (KGML) le in the KEGG database was used to map the network of genes and metabolites. This map will contribute to a more systematic study of the interactions between transcriptome and metabolomics. The results of KGML showed that "carotenoid biosynthesis" and "starch and sucrose metabolism" are the KEGG pathways of the transcriptome that affected the genes related to environmental adaptation and signal transduction (Fig. 9B). Two genes, LOC103641704 and umc1774, are involved in each of the two processes. LOC103641704 is zeaxanthin epoxidase of chloroplastic. Umc1774 is alpha-1,4-glucan phosphorylase called starch phosphorylase1.

Location of key genes and their expression in the network of stomatal development
The basic network signal map was drawn according to the reported stomatal signaling pathway of monocotyledons (Fig. 11). The results of KEGG pathway analysis, PPI, and KGML showed that four (LOC107326007, LOC107403167, LOC100192073, LOC100281912), two (LOC100285849, LOC100280217) and two (LOC103641704, umc1774) genes were obtained from the three methods ( Fig. 10). The LOC107326007 of putative MKK7 is assumed to be located in the MAPK cascade between YDA and MPK3. LOC100285849 is close to LOC107326007 because of the interactions. LOC100280217 is located next to HY5 with the existing PPI. LOC107403167 was identi ed as a member of the OXI1 family from the MAPK pathway, and the interaction among LOC107403167, MKK4, and MPK3 were retrieved from the STRING database by PPI through the connection with PTI1. Considering the exclusion of additional domains, the protein PP2C30 of Oryza sativa was used for protein interaction selection. At last, a protein signaling pathway from PP2C30 to MPK3 was established. Therefore, we speculated that LOC100192073 is also located in the signal network near MPK3. TIFY 11c is involved in jasmonate signaling. Thus, we inferred that LOC100281912 interacts more strongly with MPK3 by COI1 [8]. The association between LOC100281912 and MPK3 was retrieved from STRING through the reported protein of Panicum virgatum. LOC103641704 interacts with the PP2C of LOC100192073, and umc1774 connects MKK4 and MPK3 by linking with PHI1. The results showed that the potential genes work on the regulatory nodes of HY5 and MPK3.
The expression of 32 genes, including 26 genes that play important roles in regulating the stomatal development of maize and six potential functional genes obtained in this experiment, were determined by quantitative polymerase chain reaction (qPCR) (Fig. 11). The qPCR results were consistent with the results of transcription expression. The 26 genes were divided into eight modules according to function location. The expression of photoreceptor genes CRY1/2, phytochrome (Phy)A1 and PhyB1/2 did not change signi cantly. The two genes in the "Ligands" module, namely, STOMAGEN and EPF2, were signi cantly up regulated. ERL1/2 in the "Receptor" module were signi cantly down regulated. MPK3 was signi cantly down regulated in "MAPK cascades". Only HY5 showed signi cant up regulation changes in the "COP1-PIF4-HY5" module in the core of the network. Among the three genes on the molecular switch, only FAMA was decreased signi cantly. The expression of POLAR and PAN1 in the "Polarity and Asymmetry" modules were signi cantly lower, whereas the expression of PAN2 was remarkably higher. The expression levels of Cyclin-A2 and CDKB1-1 in the "Cell Division" module was signi cantly up regulated. The expressions of other genes were not remarkably as shown in Fig. 11. The quantitative expression of the eight potential genes and their protein network locations are shown in the Fig. 11.
Among them, the expression of six genes increased, that of one gene decreased, and the expression of one gene was too low to be detected.

Effects of blue light on transcription expression and metabolite formation in maize leaves.
Light change has a remarkable effect on the transcription expression of maize seedling leaves. It directly leads to the change of stomatal development direction in maize. A large number of differential genes are affected by light change. The inadequate annotation of the maize genome and the number of unidenti ed proteins brought new challenges to the research on the response of maize stomatal development network to environmental change. The network can be mapped and the protein function at a speci c location can be preliminarily inferred based on other reported Gramineae plants. In this study, 1296 DEGs and 419 differential metabolites were found. Among these DEGs, 53 genes directly in uence the adaptation of plants to the environment and subsequent signal transduction. Two genes, related metabolites, and their metabolic networks were obtained by analyzing the differential metabolites of the differential genes. The differential gene indicated that maize reacted in the signal pathway after red light was converted to blue light; therefore, the variation of downstream metabolite was determined. The stomatal development of plants involved in complex network signals, including the regulation of transcription, protein and metabolic levels [9][10][11]. Our results also showed that, at the transcription level, light signal is the main initial factor that in uences downstream regulation in the reaction procession of environmental signal. At the metabolic level, the metabolic pathways of "carotenoid biosynthesis" and "starch and sucrose metabolism" may also participate in stomatal development.

Complement network of stomatal development in maize by PPI and KGML.
PPI plays an important role in the transcription and signal transmission between cells. Determining the interactions between proteins can provide information about the functional structure and relationship of cells [12]. Key network genes (26 genes) and the genes (27 genes) obtained from the KEGG signaling pathways "environmental adaptation" and "signal transduction" were used as templates to search for PPI.
Environmental light mediates the signal connection between activator and deactivator. 55 differential proteins interact to form six clusters with their respective nodes (Fig. 4). 12 out of 28 genes unannotated in the KEGG pathway represent the protease, such as MAPK [13], serine/threonine phosphatase [14], deoxyhypusine synthase [15], and 4-coumarate-CoA ligase [16], which can change the signaling response to light. The key node starts from the photoreceptor of stomatal development signal and is connected with the downstream of Suppressor of phyA (SPA) as a partner of COP1 [17]. Calcium signal regulates the signaling activity of the protease (Table 2), which in turn affects the cell cycle [18]. In addition, SLAH is expressed in guard cells [19]. Six uncharacterized genes that interact with HY5 and the MAPK cascade were obtained after clustering. These results indicate that these 28 genes were directly or indirectly involved in cell differentiation or stomatal separation induced by light regulation.
Two annotated genes were obtained by KGML analysis, which links the metabolites with the corresponding transcription level changes. In most plants, light promotes carotenoid biosynthesis and storage. Carotenoid accumulation is positively correlated with photomorphogenic development [20]. Cryptochromes transduce blue light signal by regulating the stability of transcription factor HY5. HY5 is a direct activator of the PSY gene, which accumulates carotenoids [21]. Meanwhile, strigolactone shares biosynthetic precursors with ABA, both of which are derived from carotenoids [22,23]. They can interact with COP1-HY5 and affect stomatal development [24]. The MAPK signaling pathway is closely related to glucose [25,26]. Therefore, as precursors, the synthesis and transportation of starch and sucrose affect the MAPK pathway, and subsequently affect the stomatal development.
3.3 Regulation of key nodes of HY5 and MAPKs by blue light was the determining factor for the differences in stomatal development.
Transcription factor HY5 is a hierarchical regulator of transcription cascades in photomorphogenesis that acts on the downstream of many photoreceptors and serves as a signal bridge between light and metabolite [27]. First, HY5, as a target of COP1/SPA complex, is regulated by light and metabolites at the protein and transcription levels [28][29][30]. Second, HY5 bridges light and metabolites signals and promotes light morphogenesis at the transcription level by activating the gene transcription of metabolite signaling pathways [31]. The reported work shows that photolyase genes are controlled by the transcription of HY5/HYH and induced by light exposure [32]. (6 − 4) photolyases reverses the formation of photodimer and uses blue light energy to restore the original base to avoid DNA damage [33]. Therefore, the increased expression of LOC100280217 may be related to blue light induction and increased HY5 expression. The detection of DNA damage is an important process to resist and tolerate environmental factors that cause damage. DNA photolyases take part in the network of proteins that is employed to recognize damage and initiate a signaling cascade that inhibits the progression of the cell cycle [34]. Although the expression level of (6 − 4) DNA photolyase was increased under blue light, it did not form cell cycle arrest by inhibiting CDKB1-1 [34]. This result indicates that HY5 may play a dominant role in the signal pathway of the HY5-(6 − 4) DNA photolyases cascade.
MPK3 in the MAPK cascade is a key signal regulator that controls organ development and responds to abiotic stress [35,36]. In the stomatal development regulated by blue light, MPK3 is the center of ligand signal transmission to molecular switches. Several potential functional genes identi ed in this study were simultaneously associated with MPK3 at PPI levels. PP2C interacts with ABA receptors and downstream SnRK2s. MPK2, which is the substrate of SnRK2s, is activated by MKK3 simultaneously in response to ABA [37,38]. Therefore, LOC100192073 is involved in the ABA signaling pathway and the MAPK cascade reaction at the same time. LOC107326007 was identi ed as the homologous gene of MKK7, which is key node in the MAPK cascade that controls stomatal development and transmit upstream signals to MPK3/4 [39]. The CBL-CIPKs of LOC100285849 form a signal cascade with MKK7/9 as induced by Ca 2+ [40]. Therefore, we inferred that LOC100285849 competes with MKKs on calcium channels and participates in stomatal development. LOC107403167 was identi ed as a member of the STPK family, but it only has one Pkinase domain whereas OXI1 has 2 Pkinase domain as determined through BLASTP. The results showed that the function may be similar, but the differential domains suggest that they are not the same protein. However, interactions with MPK3 and known signal relationships between the STPK family and cyclins indicate that this gene may be deeply involved in stomatal development [41]. OXI1 affects cyclins and cyclin-dependent kinases (CDKs) by the H 2 O 2 pathway [42]. LOC100281912 has a tify domain and a CCT2 domain, and its protein sequence is similar to that of TIFY 11c of the TIFY JAZ subfamily. Therefore, it may be a homologous allele of the TIFY family of maize. JAZ of the TIFY transcription factor subfamily protein may act as a repressor of jasmonate responses and mediate the interaction between JAZ proteins, a transcription repressor, and COI1 protein, a receptor [43]. Ca 2+ , MAPKs, and other factors constituted the JA signal cascade [44]. This composition may be the reason why LOC100281912 interacts with MPK3 protein and participates in stomatal development.
LOC103641704 and umc1774 were characterized in maize genome sequence. 103641704 is zeaxanthin epoxidase, which is linked to 100192073 through interaction with SnRK1. umc1774, which co-occurs with MKK4, MPK3, and phi1 in genomes, is alpha-1,4 glucan phosphorylase and an important allosteric enzyme in carbohydrate metabolism. Transcription and metabolism analyses showed that these two genes in uence stomatal development through their respective signal networks and form bridges between protein signals and sugar networks.
3.4 Speculation of the signaling pathway of blue light affecting the stomatal development of maize.
Key nodes in the known stomatal development network of monocotyledon plants were selected to construct the signal transmission network. These nodule genes and DEGs obtained from transcriptome analysis were used for PPI relationship construction and clustering, but no DEGs were obtained for relationship mapping except HY5. Therefore, the signal network in maize transcriptome concentrated on MAPKs and related signal cascades when light turned blue. The depicted signal networks show an interpretable signal process. The increase in ligand protein activity will promote ligand-receptor pair binding in the next step [45]. As the head of the MAPK cascade, YDA, which has no signi cant differences, suggests that competitive binding between ligands and receptors may still be conservative.
The remarkable increase of HY5 expression in the COP1-HY5-PIF4 pathway may be the reason for the signi cant decrease in MPK3 expression at the end of the MAPK cascade [2]. Existing research on monocotyledons plants could not prove the regulatory relationship between MPK3 and the molecular switch FAMA. However, FAMA has a negative regulation on cyclins and CDKs [7], which may be how blue light promotes the increase in the stomatal density and stomatal index of maize [6,46].

Conclusion
The research on the leaf stomatal development of monocotyledons is too complex to clarify. The in uence of speci c signals, such as blue light, on signal transduction is unpredictable and variable because of the complexity of the signal cascade. Using a combination of transcriptome and metabolomic analyses can help us sort out the context in a complex network. In this study, we obtained eight genes, corresponding signal transduction, and two related metabolite networks. The identi cation of these genes and networks will help complete the existing research on the stomatal development of maize.

Plant material and growth conditions
Maize seed Xianyu 335 is a production of DuPont Company (US). The maternal part is PH6WC and the paternal part is PH4CV. Maize seedlings were cultivated in a climate-controlled cultivation chamber, which Illumina sequencing platform HiSeqTM 2500 and 125 bp/150 bp paired-end reads were generated. Raw reads of sequencing were processed into clean reads by ltering low quality reads. Then clean reads were mapped to the reference genome of maize B73_RefGen_v4. DESeq was used to analyze differential expression genes (DEGs). P value < 0.05 and foldChange > 2 or foldChange < 0.5 was set as the threshold for signi cantly differential expression. The DEGs were blastx to the related species to establish a predicted protein interaction network. The Cytoscape v 3.7.1 app MCODE was used to develop the considerable modules in the network [47]. Advanced options were 2° cutoff, 0.2 Node Score Cutoff, and 5 K-Core.

Metabolites extraction and LC-MS analysis
80 mg accurately weighed sample was transferred to a 1.5 mL Eppendorf tube. Two small steel balls were added to the tube. 20 µL of 2-chloro-l-phenylalanine (0.3 mg/mL) dissolved in methanol as internal standard and 1 mL mixture of methanol and water (

Data processing and differentially accumulated metabolites identi cation
The acquired LC-MS raw data were analyzed by the progenesis QI software (Waters Corporation,Milford, USA) using the following parameters. Precusor tolerance was set 5 ppm, fragment tolerance was set 10 ppm, and retention time (RT) tolerance was set 0.02 min. Internal standard detection parameters were deselected for peak RT alignment, isotopic peaks were excluded for analysis, and noise elimination level was set at 10.00, minimum intensity was set to 15 % of base peak intensity. The Excel le was obtained with three dimensions data sets including m/z, peak RT and peak intensities, and RT-m/z pairs were used as the identi er for each ion. The resulting matrix was further reduced by removing any peaks with missing value (ion intensity = 0) in more than 50 % samples. The internal standard was used for data QC (reproducibility).
Metabolites were identi ed by progenesis QI (Waters Corporation, Milford, USA) Data Processing Software, based on public databases such as http://www.hmdb.ca/; http://www.lipidmaps.org/ and selfbuilt databases. The positive and negative data were combined to get a combine data which was imported into R tools package. Principle component analysis (PCA) and (orthogonal) partial leastsquares-discriminant analysis (O)PLS-DA were carried out to visualize the metabolic alterations among experimental groups, after mean centering (Ctr) and Pareto variance (Par) scaling, respectively. The Hotelling's T2 region, shown as an ellipse in score plots of the models, de nes the 95% con dence interval of the modeled variation. Variable importance in the projection (VIP) ranks the overall contribution of each variable to the OPLS-DA model, and those variables with VIP > 1 are considered relevant for group discrimination. In this study, the default 7-round cross-validation was applied with 1/seventh of the samples being excluded from the mathematical model in each round, in order to guard against over tting. The differential metabolites were selected on the basis of the combination of a statistically signi cant threshold of variable in uence on projection (VIP) values obtained from the OPLS-DA model and p values from a two-tailed Student's t test on the normalized peak areas, where metabolites with VIP values larger than 1.0 and p values less than 0.05 were considered as differential metabolites.

Gene excavation and location and RT-PCR expression analysis
48 genes were selected to perform qPCR according to reported regulating stomatal development key genes. Total RNA was isolated using the RNAqueous® Total RNA Isolation Kit AM1912 (Life Technologies Corp., Grand Island, New York, USA). The RNA yield was determined using a NanoDrop 2000 Spectrophotometer (Thermo Scienti c, Waltham, USA), and the integrity was evaluated using agarose gel electrophoresis with ethidium bromide stain. Quanti ed reactions were performed in a GeneAmp® PCR System 9700 (Applied Biosystems, Waltham, USA). RT-PCR was performed using LightCycler® 480 II Real-time PCR Instrument (Roche, Basel, Switzerland) with QuantiFast® SYBR® Green qPCR Master Mix (Qiagen, Düsseldorf, Germany). Each sample was prepared in three copies for analysis. The primer sequences were designed in the laboratory and synthesized by Generay Biotech (Generay, Shanghai, China) based on the mRNA sequences obtained from the NCBI database. GAPDH was used as a reference gene. The sequence for primers used in qPCR were listed in table S1. Genes were annotated by Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KO, and GO databases, and then applied into GO enrichment analysis and KEGG pathway analysis.

Drawing and Statistical analysis
R toolkit is used to draw the pictures of transcriptome analysis and metabonomic analysis. Heatmap were drawn with the Tbtools [48]. The analysis of phylogenetic tree was carried out by MAGA X [49], and depicted by EvolView online tools [50]. Signal network were depicted by Adobe Illustrator CC 2019 (Adobe Systems Incorporated, San Jose, CA, USA). Interaction networks of genes, proteins and metabolites were drawn by Cytoscape 3.5.1 [47]. ANOVA was used to analyze the signi cant differences between the measured data by comparing their means. The signi cance level was set at 0.05 (α).

Competing interests
The authors declare that they have no competing interests.

Funding
Funding number: JAT190164. Funding name: study on the in uence of light environment on the yield and quality of amaranth. Funding sourse: research project of young teachers of Fujian education department.
Authors' contributions TL contributed to the conception, writing and part of data analyses of the study; XZ performed the experiment and part of data analyses of the study.
All authors have read and approved the manuscript.