Transcriptomic and metabolomic profiling of melatonin treated soybean (Glycine max L.) under drought stress during grain filling period through regulation of secondary metabolite biosynthesis pathways

There is a growing need to enhance the productivity of soybean (Glycine max L.) under severe drought conditions in order to improve global food security status. Melatonin, a ubiquitous hormone, could alleviate drought stress in various plants. Earlier, we demonstrated that exogenous melatonin treatment could enhance the tolerance of drought-treated soybean. However, the underlying mechanisms by which this hormone exerts drought resistance is still unclear. The present study used transcriptomic and metabolomic techniques to determine some critical genes and pathways regulating melatonin response to drought conditions. Results showed that exogenous melatonin treatment could increase relative water content and decrease electrolyte leakage in the leaves and increase seed yield under drought stress. Transcriptomic analysis showed that there were 852 core differentially expressed genes (DEGs) that were regulated by drought stress and melatonin in soybean leaves. The most enriched drought-responsive genes are mainly involved in the ‘biosynthesis of secondary metabolites’. Metabolomic profiling under drought stress showed higher accumulation levels of secondary metabolites related to drought tolerance after exogenous melatonin treatment. Also, we highlighted the vital role of the pathways including phenylpropanoid, flavonoid, isoflavonoid, and steroid biosynthesis pathways for improvement of drought tolerance in soybean by exogenous melatonin treatment. In all, findings from this study give detailed molecular basis for the application of melatonin as a drought-resistant agent in soybean cultivation.


Introduction
As a significant world cash crop, soybean (Glycine max L.) is an essential source of dietary oils and edible proteins [1]. However, soybean production levels is negatively impacted by an array of environmental stressors, the most deleterious of which is arguably drought-related as it affects every developmental stage of the crop [2]. As an abiotic factor, drought can lower agricultural productivity by 50%, and reports from climate models suggest that global warmingtriggered drought incidences will be a more regular occurrence [3]. More so, global soybean consumption is expected to increase significantly with the ongoing rise in global human population [4]. This has informed several research initiatives in developing more climate-smart agricultural protocols to ensure food sustainability. Many Agricultural scientists have recently focused on drought-tolerance mechanisms with which plants induce several defense pathways that make them resilient during drought seasons. These defense pathways could activate various molecular, biochemical, and physiological responses, depending on the type and intensity of the stress signal perception and transduction [5][6][7][8][9]. However, these may be inadequate under highly unfavorable drought stress. Therefore, the use of exogenous protocols has been proposed as an eco-friendly alternative to ensure high soybean productivity levels under drought sessions. Although irrigation is a viable option, it can have its limitations. In the United States, for example, only 9% of soybean-planted lands are irrigated. A more feasible alternative can be the modification of plant hormones to improve plant productivity and lower susceptibility to drought stress.
As a versatile hormone, melatonin (N-acetyl-5-methoxytryptamine) is a highly conserved molecule in plants and animals [10]. Its ubiquity implies that it is of high value to these organisms [11]. Melatonin suppresses the activities of plant stressors such as chemical pollutants, unfavorable temperature ranges, and pathogenic attacks. [12][13][14][15][16]. As an important plant growth regulator, it has also been linked with increasing plant drought tolerance [17][18][19][20]. Melatonin could endow plants with drought tolerance features by regulating several complex biological and molecular pathways [21][22][23]. However, studies on the use of melatonin as an exogenous inducer in improving soybean drought tolerance are scarce.
Previously, we reported that exogenous melatonin treatment could improve the tolerance of drought-treated soybean, which was possibly due to the enhanced content of osmolytes and higher antioxidant enzyme activities that reduced dehydration and lipid peroxidation [24]. However, there has been scarce reports on detailed mechanisms of melatonin for inducing drought resistance in soybean. Significant advancements in omics techniques have resolved many challenges in crop plant research [25][26][27]. In addition, transcriptomic and metabolomic analyses of drought-induced defense pathways are still poorly understood in soybeans. In this study, the transcriptomic and metabolomic profiling of melatonin for inducing resistance to drought stress in leaves of soybean was carried out, and we analyzed how melatonin response to drought stress by integrated omics strategies. It is anticipated that these findings will provide crucial insights into the theory and application for melatonin as a soybean drought-resistant agent.

Plant materials and experimental design
This research was carried out in 2018 at Heilongjiang Bayi Agricultural University, Northeast China (124˚19'-125˚12'E, 45˚46'-46˚55'N). Eight Suinong 26 soybean (Glycine max L.) variety seeds were sown per plastic pot (33 × 30 cm) in a 1:1 (v/v) vermiculite and perlite mixture. Watering of pots with Hoagland solution was done once daily from sowing to emergence, and after development, twice daily. Distilled water was applied every 5 d to stop salt and vermiculite/perlite accumulations. Seedlings were collected at the cotyledon stage (VC).
The growth temperature was shown in S1 Table. Planting pots during the grain filling stage (S1 Fig) were divided into well-watered and drought-stressed groups following a previously described method [28]. The former was kept at 80% relative soil water content (RSWC) while the latter was maintained at 45% RSWC. 50% of these pots were treated with 45 mL of a 100 μmol L -1 melatonin solution for 3 d [17], while the remaining were sprayed with distilled water as a control treatment. The weight of each pot was taken daily for a period of 10 d to ensure that the set RSWC limits were maintained.
There were three treatments: (1) control (CK): well-watered (80%) and sprayed with distilled water; (2) drought (D): drought-stressed (45%) and sprinkled with distilled water; (3) D-M: drought-stressed (45%) and sprayed with melatonin. A randomized complete design with three replications was adopted for this study. Full-grown leaf samples were collected from each treatment after 3 d, sprayed with melatonin solution and frozen in liquid nitrogen, and stored at -80˚C.

Determination of the relative water content and electrolyte leakage in the leaf and seed yield
The relative water content (RWC) was determined by oven drying method. Briefly, the weight of fresh soybean leaf samples was measured (W1). The leaf at the drum stage was incubated with sterile distilled water for 12 h, and the weight of the leaf was determined after water surface absorbance by filter paper (BW). Finally, it was baked to constant weight to evaluate the dry weight (W2). The following formula was used to determine the RWC: RWC = (W1-W2)/ (BW-W2).
The electrolyte leakage (E L ) was measured on the day of sampling, and the surface of the leaf was washed with distilled water. After cleaning, the leaf was incubated with sterile distilled water for 24 h. The initial conductivity (E 0 ) was measured using a conductivity meter. The tube was heated at 120˚C in a boiling water bath for 30 min, and the maximum conductivity (E max ) was determined after cooling. The E L was computed as follow: E L (%) = E 0 /E max .
The pod number per plant, number of grains per plant, 100-seed weight, and grain weight per plant at maturity were determined based on 15 representative samples in each group.

Total RNA isolation and transcriptome analysis
The EASYspin Plus kit (Aidlab, Beijing, China) was used to extract the RNA of 9 samples following the manufacturer's instructions. A total of 9 RNA-Seq libraries were prepared and sequenced as described recently [29]. The obtained raw data is available in the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/) sequence read archive (SRA) with accession number PRJNA576585. After the processing of raw reads through inhouse Perl scripts, clean reads were aligned to the soybean genome using HISTAT protocol [30]. Then the reads numbers mapped to each gene were counted using HTSeq software (0.6.1 version) [31] and gene expression levels were determined using the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) [32]. The "ggplots" R-package software was used to show the cluster interactions among samples. A significant false discovery rate-adjusted P value (FDR) < 0.05 was used as the empirical parameter to identify the differentially expressed genes (DEGs) as previously described [33]. Furthermore, Gene Ontology (GO) of the DEGs was performed using the "clusterProfiler" R-package [34]. Finally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses was carried out by Blast software [35].

Metabolome analysis
Freeze-dried leaf samples were crushed using a mixer mill (MM 400, Retsch) with a zirconia bead for 1.5 min at 30 Hz. Afterward, 100 mg powder was extracted overnight at 4˚C with 1.0 mL 70% aqueous methanol. Following centrifugation at 10000 x g for 10 min, the extracts were absorbed (CNWBOND Carbon-GCB SPE Cartridge, 250 mg, 3 mL; ANPEL, Shanghai, China) and filtered (SCAA-104, 0.22 μm pore size; ANPEL, Shanghai, China) before LC-MS analysis. The quantification of metabolites was performed following the multiple reaction monitoring (MRM) analysis procedures using a triple quadrupole-linear ion trap mass spectrometer (QTRAP), API 4500 Q TRAP LC/MS/MS System. Analyses were carried out as described previously [36], with each group having six biological replicates.

Confirmation of the expression profiles by qRT-PCR
All samples were analyzed using the quantitative real-time PCR (qRT-PCR) technique to validate the transcriptomic results of this study. The same RNA samples for the transcriptome analysis were reverse transcribed using a PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China) following the manufacturer's instructions. qRT-PCR procedures were performed on a Real-Time PCR System (CFX96, Bio-Rad, USA) using TB Green1 Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa Co. Ltd, Dalian, China) according to the manufacturer's instructions. The primer sequences for specific genes are listed in S2 Table. Relative expression levels of selected genes were calculated using the 2 −ΔΔCt method [37], and the soybean gene GAPDH served as an internal control.

Statistical analysis
Data obtained from this research were expressed as the mean ± standard deviation (SD) and analyzed using SPSS 17.0 software (SPSS Inc., Chicago, IL, USA). The statistical significance among the various groups was determined using one-way analysis of variance (ANOVA), followed by Duncan's multiple range test. Values of P < 0.05 were considered to be statistically significant. For metabolome analysis, the relative abundances of each metabolite were converted to log values before analysis to meet normality. Principal components analysis (PCA) on normalized data was performed using the R software (version 3.1.1). Differential metabolites (variable importance in project, VIP � 1) were chosen at fold change � 2 or fold change � 0.5 statistically significant (P < 0.05) based on Orthogonality Partial Least Squares-Discriminant Analysis (OPLS-DA).

Effect of exogenous melatonin treatment on the drought tolerance
The RWC of soybean leaves on 6 and 9 d in the D or D-M group significantly decreased (P < 0.05) when compared to that on 3 d in each group. However, the decrease rate in the D-M group was lower than that in the D group ( Fig 1A). The E L of soybean leaves on 6 and 9 d in the D or D-M group significantly increased (P < 0.05) when compared to that on 3 d in each group. However, the increase rate in the D-M group was lower than that in the D group ( Fig 1B). The RWC of soybean leaves in the D group decreased when compared to the WW group ( Fig 1A). After of 3, 6, and 9 d of drought stress, it reduced by 15.83%, 32.22%, and 47.33%, respectively. This trend shows that the RWC of soybean leaves gradually decreased as the drought stress period increased. However, it significantly increased by 6.35%, 21.75% and 10.64% after 3 d, 6 d, and 9 d of exogenous melatonin treatment, respectively (P < 0.05). In Fig  1B, E L of soybean leaves in the D group was higher than that in the WW group after of 3, 6, and 9 d of drought stress, it increased by 74.15%, 99.49% and 104.60%, respectively. This observation suggests that there is a positive correlation between drought stress and leaf E L . However, it decreased significantly by 14.75%, 10.85% and 5.32% after 3 d, 6 d and 9 d of exogenous melatonin treatment, respectively (P < 0.05). Based on the results above, exogenous melatonin treatment could effectively alleviate the decrease of RWC and the increase of E L induced by drought stress.
As shown in Table 1, the pod number per plant, number of grains per plant, 100-seed weight, and grain weight per plant in the D group was significantly lower than that of the WW group (P < 0.05), indicating that drought stress could decrease seed yield. Exogenous melatonin treatment did not affect the pod number per plant of soybean under drought stress (P > 0.05) substantially. However, exogenous melatonin treatment significantly increased the  number of grains per plant, 100-seed weight and grain weight per plant, respectively (P < 0.05), which implied that exogenous melatonin treatment during grain filling stage could effectively prevent the decreases of the number of grains r per plant, 100-seed weight and grain weight per plant to improve the seed yield under drought stress.

Transcriptomic characteristics of all samples and differential expression analysis
Transcriptomic tools were used to study drought-induced changes. Here, soybean leaves were collected from the three groups, sequenced using RNA-Seq protocol and subjected to bioinformatics analyses. All treatments showed high correlation between biological replicates (R 2 > 0.96), thus validating its reliability (

DEGs affected by drought stress and exogenous melatonin treatment
The    (Fig 2A).
The top 20 enrichments of biological pathways in the D/D-M comparison include 'biosynthesis of secondary metabolites', 'isoflavonoid biosynthesis', 'phenylpropanoid biosynthesis', 'metabolism of xenobiotics by cytochrome P450', 'metabolic pathways', 'drug metabolismcytochrome P450', 'ovarian steroidogenesis', 'flavonoid biosynthesis', 'steroid hormone biosynthesis', 'drug metabolism-other enzymes', 'alanine, aspartate and glutamate metabolism', 'ubiquinone and other terpenoid-quinone biosynthesis', 'phenylalanine metabolism', 'glutathione metabolism', 'cortisol synthesis and secretion', 'isoquinoline alkaloid biosynthesis', 'tryptophan metabolism', 'monoterpenoid biosynthesis', 'plant-pathogen interaction' and 'cyanoamino acid metabolism' (Fig 2B). The genes that differentially expressed in the D/D-M comparison but not in the WW/D comparison were enriched in these metabolic pathways, besides 'biosynthesis of secondary metabolites' pathway, indicating that they may play important roles in the rapid adaptive response to drought stress by melatonin in soybean. Interestingly, 'biosynthesis of secondary metabolites' was common to both groups, suggesting that exogenous melatonin treatment could partially treat these changes triggered by drought stress. These pathways were then further analyzed in this report.

Metabolic analysis of soybean in response to drought
Metabolites were detected in soybean leaves using metabolome techniques to study droughtinduced changes. PCA analysis results show that all the groups were distinct but the distance between the WW and D-M groups were closer than between the WW and D groups (S6 Fig). A total of 706 compounds were detected, including amino acids, lipids, organic acids, sugars, alkaloids, amines, flavonoids, and terpenoids. In the D group, a total of 204 differentially-accu-

PLOS ONE
differentially accumulated in the D/D-M comparison but not in the WW/D comparison were enriched in these metabolic pathways, except for 'biosynthesis of secondary metabolites' pathway, indicating that they may play important roles in the rapid adaptive response to drought stress by melatonin in soybean. In addtion, it was observed the 'biosynthesis of secondary metabolites' pathway were both regulated by drought stress and exogenous melatonin treatment. It may suggest that melatonin treatment could partially reverse the 'biosynthesis of secondary metabolites' pathway, which was induced by drought.

Association of transcriptomic and metabolomic changes involved in crucial biological pathways
In this section, we focused on the connection between drought-responsive gene expression and metabolite changes to gain more understanding of the physiological changes in tolerance to drought stress in soybean. Most genes encoding key enzymes and crucial metabolites had positive correlation (quadrant 3 and 7), other genes encoding some essential proteins and metabolites had negative correlation (quadrant 1, 2, 4, 6, 8 and 9), implying that these enzymecoding genes may play critical roles in the formation or breakdown of important metabolites that lower drought stress (S8 Fig). Further studies are warranted to validate this hypothesis.

Effect of melatonin on the biosynthesis of secondary metabolites in soybean during drought stress
The expressions of many structural enzyme-coding genes of the phenylpropanoid pathway were significantly affected by drought stress (Fig 4). These genes were, however, not involved in phenolic acid biosynthesis. Therefore, there were no differentially accumulated metabolites related to phenolic acids such as, p-coumaric, ferulic, and caffeic acids. Exogenous application of melatonin decreased p-coumaric accumulation, but increased the expression of the two genes encoded for phenylalanine ammonia-lyase [EC:4.3.1.24] and a gene annotated as transcinnamate 4-monooxygenase [EC:1.14.14.91].
The detailed expression patterns of the genes involved in flavonoid biosynthesis pathways and their metabolites are presented in Fig 5. Compared with the WW group, the concentrations of isoliquiritigenin, butein, butin, biquiritigenin, garbanzol, naringenin chalcone, naringenin, eriodictyol, dihydrokaempferol and quercetin were lowered in the D group. Interestingly, these flavonoids' levels were elevated after exogenous melatonin treatment. In addition, the expression of most structural flavonoid genes, such as chalcone synthase Isoflavonoid biosynthesis was also affected by drought stress and exogenous melatonin treatment (Fig 6). In the steroid biosynthesis pathway (Fig 7), among the DEGs, genes related to β-sitosterol biosynthesis, such as farnesyl-diphosphate farnesyltransferase [EC:2.5.1.21], sterol 24-C-  -] were upregulated under drought stress, but the D group have a lower content of β-sitosterol than that of the WW group. Exposure to melatonin treatments significantly increased the density of βsitosterol.

Validation of uniqueness expression using qRT-PCR
The transcriptomic analysis results of this study were further validated using the qRT-PCR technique. Here, 6 DEGs, including lysosomal acid lipase, sterol 22-desaturase, isoflavone-7-O-methyltransferase, CYP81E1_7, vestitone reductase and UGT72E, were selected. The expression patterns of these DEGs were in consonance with the earlier observed transcriptomic data (Fig 8).

Confirm the metabolomic profile using HPLC
The metabolomic changes were further confirmed using the HPLC technique. Here, 4 differentially accumulated metabolites, including quercetin, genistein, glycitein and β-sitosterol were selected, because they were differentially accumulated metabolites in the flavonoid biosynthesis, isoflavonoid biosynthesis and steroid biosynthesis pathways from KEGG analysis, which play vital roles for improvement of drought tolerance in soybean by exogenous melatonin treatment. The representative chromatograms of HPLC were shown in S9-S11 Figs. The concentrations of quercetin, genistein, glycitein and β-sitosterol in the D group were lower than those in the WW group (Fig 9). However, melatonin treatment significantly increased their concentrations in the leaves when compared with the D group (P < 0.05). These results were consistent with the metabolomics profile earlier reported that quercetin, genistein, glycitein and β-sitosterol in the WW/D comparison were up-regulated metabolites, while were down-regulated metabolites in the D/D-M comparison.

Discussions
As an important world legume, soybean (Glycine max L.) is a major source of oil and protein for millions of people and livestock. However, drought stress triggered by climate change factors can result in significant losses as well as the reduction of seed quality [38]. Therefore, it is necessary to explore effective measures to alleviate the effects of drought stress on soybean. In 1995, melatonin, a plant growth regulator (PGR), was isolated from higher plants and functions in growth stimulation and boosting plant resistance [39]. Notably, the negative effects of drought stress were reversed by exogenous melatonin treatment in some plants, including apple [40], grape [41], cucumber [42], tomato [43], wheat [21] and maize [44]. In our previous study, it was demonstrated that exogenous melatonin treatment could improve the tolerance of drought-treated soybean [24]. However, the specific mechanisms by which melatonin alleviates drought stress is still unclear. In the present study, we reported the role of melatonin in resistance to drought stress using transcriptomic and metabolomic analysis, which could identify candidate genes and key metabolic pathways involved in drought response in soybean.
Our results showed that the exogenous melatonin treatment could increase the relative water content and decrease the electrolyte leakage in leaves and increase seed yield under drought stress. Data analyses using transcriptomic and metabolomic techniques can give more information into complex interactions between and within several factors. Thus, we analyzed how melatonin response to drought stress by integrating transcriptomic and metabolomic approaches. A combination of transcriptomic and metabolomic analysis could highlight the vital role of exogenous melatonin for the improvement of drought tolerance in soybean. The KEGG enrichment analysis was carried out to show the biological pathways of the DEGs and differentially-accumulated metabolites. It is noteworthy that pathways involved in 'biosynthesis of secondary metabolites' including phenylpropanoid, flavonoid, isoflavonoid, and steroid biosynthesis pathways were significantly enriched in both the WW/D and D/D-M comparisons through transcriptome and metabolome. Secondary metabolites abound in different levels in different plants, and their presence is most noticeable under stress conditions. External application of phenolic compounds can confer some drought tolerance characteristics, although these vary from one plant species to another [45]. In this study, two genes encoding phenylalanine ammonia-lyase (PAL) expression were down-regulated under drought stress, but increased after exogenous melatonin treatment. PAL expression has been linked to drought stress responses previously. In addition, PAL can convert phenylalanine to fuel the phenylpropanoid metabolism pathway [46]. However, drought stress did not change the phenolic acid contents of samples. At the level of phenolic acids under abiotic stresses, previous data indicated some possible reasons. Bettaieb et al. demonstrated that the production of phenolic acids could be in response to the level of water deficit [47], and this may explain why there was no change in phenolic acid levels under drought stress as observed in this study.
Aside from playing integral roles in plant development and reproduction, flavonoids are also important for plant stress defense [48]. This feature has been shown previously under induced stress conditions in Ligustrum vulgare [49], Scutellaria baicalensis [50], and rice [51]. Our study showed that melatonin treatment improved flavonoids accumulation. Flavinoids may also shield plants from conditions imposed by water deficiency [52]. Naghizadeh et al. recently reported that exogenous melatonin application mitigates the adverse effects of drought stress on morpho-physiological traits and secondary metabolites including flavonoid in Moldavian balm (Dracocephalum moldavica) [53]. Furthermore, Liang et al. reported that exogenous melatonin administration could delay the senescence of Kiwifruit Leaves by regulating the antioxidant capacity and biosynthesis of flavonoids [54]. The possible dual roles of flavinoids have also been posited, suggesting that these compounds maintain high antioxidant properties by blocking the production of reactive oxygen species (ROS) [55] and acting as ROS scavengers once ROS are produced [52]. These findings indicated that melatonin might have a role in the promotion of the antioxidant capacity of soybeans through secondary metabolism regulation.
Many plant species contain an active phytosterol called β-sitosterol, which functions in maintaining the lipid bilayers membrane integrity. Previous studies have also implicated this compound in conferring plants with stress tolerance features. Recently, the treatment of tomato seeds with β-sitosterol conferred heightened resistance during extreme temperature conditions, implying that this phytosterol can have drought-resistant characteristics [56]. Its external application significantly lowered stress caused by high salt levels in pepper (Capsicum annum) and sunflower (Helianthus annuus), as shown by growth and physiological analyses [57,58]. Metabolomics and physiological analyses revealed that β-sitosterol could improve white clover growth processes and resistance to water deficits [59]. In the steroid biosynthesis pathway, exposure to melatonin treatment significantly increased the density of β-sitosterol, indicating that melatonin provided beneficial effects in alleviating drought stress by regulating β-sitosterol production.

Conclusion
In summary, findings from this study showed that exogenous melatonin treatment could increase the relative water content and decrease the electrolyte leakage in leaves and increase seed yield under drought stress. Key metabolic pathways and metabolic products induced by exogenous melatonin treatment under drought conditions were also reported. Exogenous melatonin application could mitigate the adverse effects of drought stress by modulating 'biosynthesis of secondary metabolites' pathways including phenylpropanoid, flavonoid, isoflavonoid, and steroid biosynthesis, which regulate soybean leaf response to water deficit. Furthermore, this study gave insights into the molecular mechanisms regulating drought stress reduction by melatonin application in soybean. It also provides useful evidence for developing melatoninbased agents to protect soybean crops from drought conditions. Supporting information S1