Transcriptome Analyses and Antioxidant Activity Profiling Reveal the Role of a Lignin-Derived Biostimulant Seed Treatment in Enhancing Heat Stress Tolerance in Soybean

Soybean (Glycine max Merr.) is a worldwide important legume crop, whose growth and yield are negatively affected by heat stress at germination time. Here, we tested the role of a biostimulant based on lignin derivatives, plant-derived amino acids, and molybdenum in enhancing soybean heat stress tolerance when applied on seeds. After treatment with the biostimulant at 35 °C, the seed biometric parameters were positively influenced after 24 h, meanwhile, germination percentage was increased after 72 h (+10%). RNA-Seq analyses revealed a modulation of 879 genes (51 upregulated and 828 downregulated) in biostimulant-treated seeds as compared with the control, at 24 h after incubation at 35 °C. Surprisingly, more than 33% of upregulated genes encoded for ribosomal RNA (rRNA) methyltransferases and proteins involved in the ribosome assembly, acting in a specific protein network. Conversely, the downregulated genes were involved in stress response, hormone signaling, and primary metabolism. Finally, from a biochemical point of view, the dramatic H2O2 reduction 40%) correlated to a strong increase in non-protein thiols (+150%), suggested a lower oxidative stress level in biostimulant-treated seeds, at 24 h after incubation at 35 °C. Our results provide insights on the biostimulant mechanism of action and on its application for seed treatments to improve heat stress tolerance during germination.


Introduction
Soybean (Glycine max L. Merr.) is one of the crops most largely cultivated all around the world, with a strong production located in Brazil, Argentina, and USA [1]. Its market demand has been growing due to increased population and to different applications of this crop. Thanks to its high protein content (40-42%) as compared with other food crops, soybean represents one of the main protein sources for food and feed. Moreover, it is also largely employed for oil production, owing to the fact that it has one of the highest oil contents (18-20%) among legume crops [1]. For all these reasons, it has become essential to optimize its production, and avoid to use new lands, which would

Biostimulant Application Positively Influenced Seed Morphological Parameters and Germination Percentage under Controlled Conditions
In order to evaluate whether the biostimulant treatment had a visible influence on soybean seed morphology at the early germination phase under heat stress conditions, experiments were carried out under controlled conditions. The biometric parameters were monitored on both control and biostimulant-treated seeds, at the start of treatment (T 0 ) and at 24 h, after incubation at 35 • C, in the dark (T 1 ). Except for the weight, all the evaluated parameters showed a significantly higher percentage increase in seeds treated with the biostimulant as compared with untreated seeds ( Table 1). The biostimulant treatment prompted a significant increase in seed area (+50.34%), perimeter (+18.28%), length (+20.39%) and width (+14.74%) as compared with the control. By contrast, the weight increase was only 0.63% (not statistically significant) higher in biostimulant-treated seeds as compared to control seeds. Table 1. Biometric parameters of untreated and biostimulant-treated seeds in controlled conditions. Measurements were done before sowing and treatment (T 0 ) and at 24 h after incubation at 35 • C, in the dark (T 1 ). Values are expressed as mean ± standard deviation (SD) of five different biological replicates, each of which are composed of 45 seeds. For each row, different letters indicate statistically significant differences between untreated and biostimulant-treated samples, as calculated by t-test (p < 0.05). The last column (%∆) reports the percentage relative change between untreated and biostimulant-treated samples. The effect of the biostimulant application on soybean seed germination is shown in Table 2. The germination percentage of control and biostimulant-treated seeds was recorded at 24, 48, and 72 h after incubation at 35 • C in the dark. The biostimulant treatment prompted a significant increase in the germination percentage as compared with the control, both at 48 and 72 h. In addition, at 24 h a visible germination was not detected. Table 2. Germination percentage of untreated and biostimulant-treated soybean seeds incubated at 35 • C in the dark. Values are expressed as mean ± standard deviation (SD) of five different biological replicates, each of which are composed of 45 seeds. For each row, different letters indicate significant differences between untreated and biostimulant-treated samples, as calculated by t-test (p < 0.05). The last column (%∆) reports the percentage relative change between untreated and biostimulant-treated samples.

Biostimulant Application Modulated the Expression of 879 Genes under Controlled Conditions
A total of 879 genes were differentially expressed in biostimulant-treated seeds with respect to untreated seeds at 24 h, after incubation at 35 • C in the dark. Among all the genes, 51 were upregulated, meanwhile 828 were downregulated by the treatment. The list of modulated genes is reported in Supplementary Table S3, whereas the Volcano plot generated by the RNA-Seq analysis is shown in Supplementary Figure S1. A hierarchical clustering tree, summarizing the correlation among significant GO terms was generated ( Figure 1A). (B) Functional interactions among the upregulated genes; (C) Functional interactions among the downregulated genes. ShinyGO ® generated the hierarchical clustering tree, using the significant GO terms listed in Table 3. GO terms with many shared genes are clustered together. Bigger dots indicate more significant p-values. The graphical representations of the functional interactions were generated by STRING software.
However, most of the genes were downregulated (Table 3), and many of them were related to response to stress (GO:0006950, 80 genes, FDR = 5.49 × e −06 ). In particular, treated seeds showed a reduced expression of genes coding for heat shock transcription factors, such as LOC100527682, LOC100786140, and LOC100789792. These genes are also part of two other GOs as follows: response to chemical (GO:0070887, 90 genes, FDR = 8.83 × e −11 ) and response to stimulus (GO:0050896, 135 genes, FDR = 2.24 × e −10 ). Moreover, these GOs grouped genes coding for enzymes involved in the redox regulation, such as glutathione-S-transferase (LOC100808374), glutaredoxins (LOC100792704), and thioredoxins (LOC100810192, LOC100800129). Among the downregulated genes, we also found some responding to chemical, such as peroxidases 3 (LOC100547872) and peroxidase 5 (LOC100811641), which are involved in removal of H 2 O 2 , oxidation of toxic reductants, biosynthesis and degradation of lignin, suberification, auxin catabolism, and response to environmental stresses. In the same category, the biostimulant treatment downregulated the expression of genes coding for a protein kinase superfamily protein (LOC100794703) which acts as a positive regulator of abiotic stress response and zinc finger transcription factors (LOC100806997). The biostimulant treatment appeared to influence the expression of genes coding for repressors of hormone signaling, such as AUX/IAA (Auxin/Indole-3-Acetic Acid) proteins (LOC100802759, LOC100791342, and LOC100799875) and DELLA (aspartic acid-glutamic acid-leucine-leucine-alanine) proteins (LOC100805968 and LOC100791952), together with protein kinases (LOC100794703) and receptor kinases (LOC102668647) acting into hormone-induced pathways. Along with a general downregulation of stress-related responses and the influence on hormone signaling, the biosynthesis of secondary metabolites (Kegg pathway 01110, 10 genes, FDR = 0.000697) was decreased in the presence of the biostimulant. Among the others, the biostimulant treatment downregulated genes coding for cytochrome P450s and enzymes such as beta glucosidases (LOC100820528, LOC100777773) and O-methyltransferases (LOC100787536) which act on phenolic compounds. The GO analysis also highlighted a global downregulation of primary metabolic processes (GO:0044238, 160 genes, FDR = 0.000933), including carbohydrate metabolic process (GO:0005975, 35 genes, FDR = 0.00117). In particular, it regulated the expression of genes which encode for enzymes involved in the promotion of cell wall organization or biogenesis (GO: 0071554, 21 genes, FDR = 0.00993) such as pectinesterases (i.e., LOC100792319, LOC100794948), glucosyl transferases (LOC100812586), and xyloglucan:xyloglucosyl transferases (LOC100778482). Moreover, the biostimulant acted on genes involved in sucrose metabolism such as sucrose synthase 4 (LOC100819730), and trehalose accumulation such as haloacid dehalogenase-like hydrolase domain-containing protein and trehalose-6-phosphate phosphatase J (LOC100803119). The glycolysis and fermentation pathways were also affected by the biostimulant treatment since both phosphofructokinase 3 (LOC100818755) and alcohol dehydrogenase 1 (LOC100800668) transcript levels were reduced. The enriched primary metabolic process GO term (GO:0044238) also grouped genes involved in protein and amino acid metabolism. In particular, our analysis showed the downregulation of glutamate decarboxylases (LOC100812201 and LOC100781791) and aspartate aminotransferases (LOC100780254), able to catalyze the biosynthesis of γ-aminobutyric acid (GABA) and aromatic amino acids, respectively [26,27]. Last, but not least, lipid metabolism was also negatively regulated by the biostimulant. For instance, the primary metabolic process GO term (GO:0044238) also contained lipid phosphate phosphatase 2 (LOC100782531), long-chain acyl-CoA synthetase 2 (LOC100806645), and lipoxygenase 4 (LOC100811820).

String Suite Analysis Revealed Two Functional Interaction Networks
In order to identify the cellular functions affected by the biostimulant treatment, and to support data analysis, all the functional interactions between the expressed proteins were investigated in both up-and downregulated datasets and STRING suite ( Figure 1). Regarding the interactions of the upregulated genes, a number of 34 nodes and 63 edges was observed (average node degree = 3.71) ( Figure 1B). Considering the value of 12 as the expected number of edges for 34 nodes, this network had significantly more interactions than expected [PPI (Protein-Protein Interaction) enrichment p-value = 1.00 × e −16 ]. A clear, but small, network was observed, related to SAM-dependent methyltransferase and SAM superfamily. Considering the interactions of the downregulated genes ( Figure 1C), 483 nodes and 1078 edges were observed (average node degree = 4.46). In addition, this network had significant enrichments (PPI enrichment p-value = 1.00 × e −16 ), 594 being the expected number of edges for 483 nodes. In this dataset, the visible highly interconnected subgroups were those linked to stimuli response and those linked to primary metabolism.

Biostimulant Application Reduced Oxidative Menaces
In order to gain more insight into the seed response to heat stress and H 2 O 2 production during the early phases of seed germination, the level of H 2 O 2 , the activities of ROS scavenging enzymes (superoxide dismutase (SOD), catalase (CAT) and glutathione S-transferase (GST)), and the content of non-protein thiols were analyzed after the application of the biostimulant. The obtained data were compared with those measured in non-treated samples, and results were expressed as relative value ( Figure 2). Supplementary Table S4 provides the report of the absolute values of this quantification. In general, the biochemical analysis conducted on biostimulant-treated soybean seeds at 35 • C suggested a better response to stress as compared with the control (untreated seeds). However, it seemed that the capability to contrast oxidative menaces, apparently, was not linked to changes in enzymatic activity. Indeed, after the application of the biostimulant, the enzymatic activity of SOD and GST did not show significant differences with respect to untreated seeds, and CAT activity was strongly downregulated (FC 0.28). In addition, the content of H 2 O 2 was drastically reduced (−40%) and the content of non-protein thiols strongly increased (+150%).  Supplementary Table S4. Each bar represents the mean ± standard deviation (SD). Lowercase letters indicate significant (p < 0.05) differences among the different relative measures, as calculated by one-way ANOVA followed by Tukey's test. The asterisk (*) indicates statistical differences between untreated (dotted line) and treated samples, as calculated by t-test (p < 0.05). Additional statistical information is reported in Supplementary Table S5.

Discussion
Seed treatment technologies consist of seed treatment with synthetic or natural compounds, which are aimed at enhancing the uniformity and vigor of seedlings and improving the plant tolerance to abiotic stresses [28]. The use of biostimulants to counteract the effect of abiotic stress is well-recognized and their function in promoting plant defenses against adverse environmental conditions has been studying lately [9,12]. Seed pretreatments have been revealed as a good method to contrast the ever-increasing environmental stress presence, improving the yield of productions, starting from seed germination [29]. This is a faster method as compared with conventional breeding or plant genetic modification and could be useful for seed treatment in countries, such as Brazil, where high temperature at sowing represents a limiting factor [30]. Among the crops whose cultivation occurs under this adverse environmental conditions, soybean (Glycine max) is affected in terms of poor germination, increased incidence of pathogen infection, and decreased economic value [31].
The potential effects derived from the application of the biostimulant objective of this study, was preliminary investigated, in 2017, at the Detec Experimental Station located in Brazil (Taquarituba, 23 • Table S1). During these field trials, the biostimulant applied as seed treatment on soybean, was able to positively affect plant vegetative and reproductive growth under natural heat stress conditions. The biostimulant seed treatment efficiency for improving germination and mitigating oxidative damages under heat stress conditions has been observed in our previous work on cucumber seeds [9]. In the present work, in order to evaluate the potential effects of the biostimulant in promoting heat stress tolerance also in soybean, we conducted germination trials under controlled conditions, and we investigated the potential mechanism of action using combined biochemical and transcriptomic (RNA-Seq) approaches.

The Biostimulant Promoted Soybean Germination and Growth after Heat Stress Exposure
Our preliminary field trials showed that at seven days after sowing, root length and height of plants grown from biostimulant-treated seeds were higher with respect to the control plants (Supplementary  Table S1). Moreover, the final yield observed in plants primed by the biostimulant was generally higher at harvesting time, suggesting that the pretreatment effect was prolonged at later plant growth stages. Indeed, although the biostimulant seemed not to significantly change soybean germination rate at the early stage, the final germination percentage measured under controlled conditions, appeared to be higher in biostimulant-treated seeds at 72 h after incubation at 35 • C. Likewise, the biostimulant seed treatment increased cucumber germination percentage under heat stress conditions only at 48 h after incubation at 35 • C [9]. Therefore, our results give an indication of a possible late effect of the biostimulant in promoting plant vegetative, and then reproductive growth, under heat stress conditions [32].
It has been demonstrated that in salt stress, the activation of different biochemical mechanisms related to stress response, promotes a "priming memory" in seeds which can be recruited upon a later stress exposure and provokes higher stress tolerance of germinating primed seeds [33]. Wei and co-workers observed that seven-day-old soybean seedlings grown from melatonin-treated seeds and subjected to drought stress showed an enhanced tolerance as compared with untreated seeds [34]. In accordance with our results, treatment of maize seeds with solutions of increasing concentrations of lignin carbon showed no influence on seed germination rate, but a positive biological activity on the emergence of maize plantlets [35]. Moreover, Amirkhani and co-workers showed that broccoli seeds coated with plant protein lysates enhanced seedling shoot and root growth as compared with uncoated seeds, whereas their germination was negatively affected by the treatment [36]. This negative effect could be related to a barrier for water uptake and gas exchange due to the coating during the whole germination process. Interestingly, the biostimulant also contains lignin derivatives and protein hydrolysates of plant origin. In addition, our data do not show any differences at the weight level for primed seeds at 24 h after incubation at 35 • C, thus suggesting the absence of changes in water uptake. In parallel, our RNA-Seq data analysis showed a downregulation of aquaporin gene TIP2-1 (LOC100803901). Aquaporins are proteinaceous channels known to regulate transmembrane water transport, and therefore play a key role in the imbibition process during seed germination [37]. Previous studies conducted on spinach (Spinacia oleracea L.) unprimed and osmoprimed seeds subjected to chilling and drought stress showed a downregulation of the genes coding for aquaporins [including TIP (Tonoplast Intrinsic Protein) subfamily] as compared with seeds germinated in optimal conditions [38]. There is evidence that different seed pretreatments cause a decrease in water absorption owing to the lower water potential of the solution. In this way, the treated seeds have more time to complete DNA repair processes and also reduce cellular damage, which often occurs during rapid seed rehydration in the germination process [30]. This hypothesis is also supported by RNA-Seq analysis, and in particular from the upregulation of Elongator protein 3/MiaB/NifB, SAM, and FTSJ-like methyltransferase genes (Table 3).
In addition to aquaporin downregulation, we could not exclude that the differences observed for some biometric parameters such as area, perimeter, length, and width of the seeds treated with the biostimulant could be related to a different level of gas exchange. Indeed, in the literature, seed pretreatments have been reported to induce the formation of a thick mucilage layer during imbibition, which could affect water and gas exchanges, by promoting seed swelling without increasing the seed weight [14,39]. However, this hypothesis should be confirmed in future experiments.

Rna-Seq Analysis Displayed That Heat Stress Mitigation after Biostimulant-Treatment May Be Linked to Different Molecular Pathways
Although the biostimulant treatment did not affect the early germination percentage, RNA-Seq analysis seemed to corroborate the already active effect of the biostimulant for mitigating stress in seeds germinated under high temperature conditions. First, we observed a strong downregulation following biostimulant treatment, since only 51 genes out of 879 modulated, were upregulated (Tables 2 and 3). However, a stronger downregulation, as compared with upregulation, was observed with other biostimulants. For example, the treatment with lipophilic components of the brown seaweed Ascophyllum nodosum improved freezing tolerance in Arabidopsis thaliana by inducing the upregulation of 463 genes and the downregulation of 650 genes [40]. Moreover, Contartese et al. (2016) similarly observed a higher number of downregulated genes (2415) as compared with the upregulated genes (1731) in tomato plants treated with Expando ® , a biostimulant developed to increase fruit size [41].
In this study, in the presence of heat stress, soybean seeds treated with the biostimulant showed upregulation of three proteins involved in the biogenesis of ribosomes, as well as a specific group of 12 methyltransferases acting on rRNA (Table 4, Figure 1B). DNA methylation is well known to affect plant developmental processes and to be activated in response to abiotic stress [42], strongly connected to the perception of external environmental stimuli, to changes in stress-responsive gene expression, and to pretreatment effect activation [43,44]. In contrast, RNA methylation in plants is less investigated, but accumulating reports have demonstrated that RNA methyltransferases are essential for plant growth as well as for abiotic stress responses [45]. In animals, ribosomal RNAs (rRNA), which are encoded by multiple copies of rDNA at the genomic level, display tissue-specific expression patterns, raising the chance to create combinations of rRNAs and ribosomal proteins forming ribosomes ("ribosome heterogeneity") likely able to translate subsets of transcripts [46]. Ribosomal RNA modifications (e.g., methylation) during production/maturation are at the base of the ribosome heterogeneity [46] which can regulate ribosome assembly on a group of specific transcripts [47,48]. Recently, a specific ribosome methylation event facilitating the selective translation of a cyp-29A3 gene has been reported, which by increasing, the production of eicosanoids, modulated stress resistance in C. elegans [49]. In our system, the heat stressed seeds, pretreated with the biostimulant, specifically induced rRNA methyltransferases which exhibited a high level of connectivity in the same emerging gene network ( Figure 1B). According to [49], this methylation process could, by analogy, facilitate selective translation of genes guiding the observed resistance to heat stress induced by the biostimulant treatment. This hypothesized translation-dependent layer of regulation deserves to be deeply investigated in future experiments. Moreover, few genes involved in stress response (Supplementary Table S3) are also upregulated by the application of the biostimulant. In particular, our analysis revealed the upregulation of NF-X Like 1 (LOC102669482) and pentatricopeptide repeat containing protein (LOC100788313). On the one hand, the human NF-X1 protein and homologous proteins in eukaryotes represent a class of transcription factors whose common feature is the cysteine-rich region, which possess a variable number of repeated motifs, defined as NF-X1 type zinc finger motifs. It has been demonstrated in Arabidopsis that NF-X Like 1 protein is involved in a regulatory mechanism able to improve the physiological status of plants and to support growth and survival under salt stress [50]. On the other hand, pentatricopeptide repeat (PPRs) constitutes one of the largest gene families in Arabidopsis. Several genes belonging to this family are involved in tolerance to different biotic and abiotic stresses [51]. The mitochondrial PPR protein PGN (pentatricopeptide repeat protein for germination on NaCl) was identified to positively regulate biotic and abiotic stress response. Arabidopsis plants with mutation in PGN show low resistance against necrotrophic fungi as well as towards GABA, glucose, and high salinity [52].
The potential stress mitigation induced by the biostimulant treatment could also be confirmed by the downregulation of stress response-related genes whose interconnection is high in a protein-protein network ( Figure 2). Among others, genes coding for heat shock factors (HSFs), in particular HSFB2A, B3, and B4 (LOC100789792, LOC100786140 and LOC100527682) were found to be downregulated. Currently, more than 30 genes coding for HSF divided into three classes, A, B and C, have been described in soybean [53]. Functional analyses show that some members of class B HSFs are active transcriptional repressors. The repressive activity of HSFB1 and HSFB2b was reported also in Arabidopsis and a domain, designated the B3 repression domain (BRD), was found in HSFB1 and HSFB2b at the C terminus [54]. The two HSFs, HSFB1 and HSFB2b, are able to repress the general heat shock response under non-heat-stress conditions and act via direct repression of the expression of the heat stress inducible HSFA2 [54]. The overexpression of OsHsfB2b, a member of B HSF in rice, significantly decreased drought and salt tolerance in transgenic plants, while the OsHsfB2b-RNAi plants enhanced tolerance under these stresses indicating that this protein could negatively regulate the abiotic stress tolerance. The observed downregulation of HSFs belonging to class B could be important for promoting enhancement of heat stress tolerance in the biostimulant-treated seeds germinated at 35 • C.
The RNA-Seq analysis also reveals the downregulation of different genes belonging to redox response, such as different isoforms of peroxidases and RbohB, a gene coding for a calcium-dependent NADPH oxidase that generates superoxide. These proteins are involved in responses to different stresses as revealed by reverse genetics approach and "overexpression" of Rboh genes [55]. For example, the pathogens Pseudomonas syringae and Peronospora parasitica induce an oxidative burst in Arabidopsis by enhancing AtrbohD and AtrbohF expression [56]. Interestingly, Rboh was also downregulated at 24 h in cucumber seeds pretreated with the biostimulant and incubated at 35 • C [9].

Biostimulant Treatment Mitigated the Accumulation of Reactive Oxygen Species (Ros) via a Non-Enzymatic Way
In seed physiology, reactive oxygen species (ROS) are usually considered to be toxic molecules, whose accumulation leads to cell injury and consequent problems in seed germination and development. However, there is increasing evidence that ROS, at low concentrations, can function as signaling molecules involved in a wide range of responses to various stimuli [57]. The dual function of ROS in plants mainly relies on the cellular antioxidant machinery, which involves detoxifying enzymes and antioxidant compounds. Such mechanisms can scavenge potentially toxic ROS, generally produced under stressful conditions, or rather tightly control ROS concentrations in order to regulate various signaling pathways. Among ROS, hydrogen peroxide plays a key role during germination, however too high levels of H 2 O 2 can be toxic for seeds [58]. In our study, the level of H 2 O 2 were significantly lower in the biostimulant-treated seeds, suggesting a possible role of the biostimulant for preventing its accumulation. Biostimulant-treated cucumber seeds also showed a reduction in endogenous H 2 O 2 levels [9]. The lower amount of H 2 O 2 observed in biostimulant-treated soybean seeds is correlated to CAT activity. This enzyme, along with peroxidase (PRX), is directly involved in the disruption of H 2 O 2 , and their activation indicates the effort of plants to reduce oxidative damage. A lower activity of CAT and an undetectable activity of PRX can be an indication of the biostimulant's capability for mitigating stress effects on the seed germination process. Our enzymatic data are confirmed also by gene expression results. Indeed, genes coding for peroxidases, such as PRX-3 and PRX-5, resulted downregulated along with RbohB (Table 3). Non-protein thiol compounds represent another important parameter in stress response. They are the most important source of sulfur in soybean seed, a fundamental element involved in metabolic pathways, nutritional quality, and plant productivity. Thiol metabolism plays a key role in plant growth, development, and defense against a range of environmental stresses [59]. It has been demonstrated that thiols, such as glutathione (GSH) together with its regulation in redox signaling and defense processes, are important components for the heat stress tolerance [60]. The glutathione pool has been shown to be associated with the response to heat stress of maize [61], Coleus blumei and Fagus sylvatica L. [62], Triticum aestivum [63], and Vigna radiata [64]. Soybean seeds have nearly undetectable levels of glutathione, but this compound is substituted by the tripeptide homoglutathione which functions in a similar way to GSH by showing the same redox reaction [59,65]. Glutathione-S-transferases (GST) are directly linked to thiols, since they play important roles in enzymatic thiol-dependent ROS scavenging mechanisms [66] by catalyzing the conversion of H 2 O 2 using glutathione or homoglutathione as substrates. In our study, the levels of non-protein thiols were higher in the biostimulant-treated seeds as compared with the untreated seeds, while the activity of GST did not change after biostimulant application. Differently, the biostimulant promoted an upregulation of GST activity in Cucumis sativus, at 24 h after incubation at 35 • C, while non-protein thiol level was higher at 48 h [9]. In general, the biostimulant seed treatment led to a reduction in ROS accumulation, due to a direct or indirect scavenging effect induced by the product [67,68]. On the basis of our results, we hypothesize that the biostimulant pretreatment effects on antioxidant capacity appear to be more correlated with non-protein thiol levels than with increased transcription or activity of enzymatic ROS scavengers in soybean seeds.

Soybean Seed Primary Metabolism and Hormone Signaling Are Affected by the Biostimulant Treatment
Surprisingly, the biostimulant treatment also affected the primary metabolism of soybean seeds by inducing a dramatic global downregulation of primary metabolism process (GO:0044238, 160 genes, (Figure 1), whose activation was directly related to the promotion of germination [69]. Concerning carbohydrate metabolism, sucrose synthase is one of the most important enzymes directly involved in the germination process, since it catalyzes the reversible cleavage of sucrose to glucose and fructose. Its downregulation induced by the biostimulant seems to confirm the reduction in the isolation of sugar-nucleotide precursors for structural and storage polysaccharide biosynthesis [70], which indeed appears to be negatively affected by the biostimulant treatment. Moreover, some of the carbohydrate metabolism regulated genes control trehalose accumulation, which has a role in improving the abiotic stress tolerance [71]. Interestingly, the glycolysis pathway is also affected by the biostimulant treatment, i.e., a decrease in the transcript level of phosphofructokinase 3, encoding for the enzyme which catalyzes the phosphorylation of D-fructose 6-phosphate to fructose 1,6-bisphosphate by ATP, the first committing step of glycolysis. Moreover, alcohol dehydrogenase 1 downregulation suggests a reduced alcohol fermentation in seeds treated with the biostimulant. Together with the carbohydrate metabolism, protein and amino acid metabolism is also affected by the biostimulant treatment. Seed storage proteins are critical to provide amino acids for protein synthesis and energy production, particularly glutamate and aspartate, which are the most abundant amino acids in soybean seeds. In particular, our analysis showed the downregulation of glutamate decarboxylases and aspartate aminotransferases, both encoding for enzymes whose activity is promoted by imbibition and whose role is essential along seed germination [69]. While glutamate decarboxylase (LOC10081220) promotes the accumulation of GABA [27], aspartate aminotransferase (LOC100780254) directly promotes the biosynthesis of aromatic amino acids [72]. Since imbibition is a critical step in seed germination, the downregulation of the primary metabolism genes observed with the biostimulant-treated soybean seeds might be attributable to the lower water absorption linked to the seed treatment process. However, this slowing down in germination has also been demonstrated to be important for DNA repair and cellular damage prevention [73].
The hormone mediated signaling pathway (GO:0009755) has a key role in modulating seed germination, development, and response to stress [74]. In particular, during soybean germination, auxins and ABA (Abscisic Acid) appear to repress germination [75], while gibberellin and ethylene promote germination [69]. Our transcriptomic analysis highlighted a complex influence of the biostimulant treatment on seed hormone cascade. Concerning auxin signaling, the biostimulant treatment appeared to significantly reduce the expression of different AUX/IAA proteins (IAA29, IAA19, and ATAUX2), which are known to act as repressor of early auxin response genes under low auxin concentration [74]. On the contrary, the biostimulant seemed to positively affect the gibberellin signaling pathway by negatively modulating the expression of genes coding for DELLA proteins (LOC100805968, GAI and LOC100791952, RGL3), which act as repressor of gibberellin signaling. Last but not least, ABA signaling was also regulated by the biostimulant treatment, in terms of decreasing the expression of genes encoding for protein kinases (LOC100794703) and receptor kinases (RK1) acting along its signaling pathway.

Experiments in Controlled Conditions
Soybean seeds were treated by following the protocol provided by Embrapa in 2005 [76], described in details in Campobenedetto et al. [9], and currently used in Brazil for the seed treatment with different products, including phytochemicals and biostimulants [77]. Briefly, 20 mL of the biostimulant solution was diluted in distilled water in order to reach the final volume of 80 mL. Then, the diluted biostimulant solution was added drop by drop to 250 g of dried seeds kept in continuous shaking, until the complete and visible distribution of the product on the seed surface was obtained. After the treatment, the seeds were dried at room temperature, and then employed for the experimentations. In order to perform germination tests, 45 of them were placed in glass Petri dishes (20 cm of internal diameters) containing two filter papers saturated with 15 mL of distilled water. Finally, the Petri dishes were incubated for 24 h at 35 • C, in the dark. After the incubation time, the seeds were collected, grinded with mortar and pestle using liquid nitrogen, and then stored at −80 • C for further analyses. The experiment was repeated four more times, in order to have five different biological replicates. Moreover, a parallel experiment, where distilled water was used instead of the biostimulant, was conducted and used as control.

Evaluation of Biometric Parameters under Controlled Conditions
Morphological analysis on seeds germinated in controlled conditions was conducted in order to evaluate differences between biostimulant-treated and untreated seeds. Each seed was weighed and area, perimeter, length, width, and weight were measured by using the SmartGrain software [78], and then statistically evaluated comparing treated and control seeds by a Student's t-test at p < 0.05 using SYSTAT 10 software. All measurements were taken immediately after treatment (T 0 ) and after 24 h of incubation at 35 • C, in the dark (T 1 ). The changes were calculated on values obtained from the difference between T 1 and T 0 . In addition, germination percentages were calculated for control and biostimulant-treated seeds at 24, 48, and 72 h of incubation at 35 • C, in the dark.

Total Rna Isolation
Total RNA was isolated from powdered seeds, as previously described, by using TRIzol ® reagent (Thermo Fisher Scientific, Waltham, MA, USA) [79]. Moreover, in order to increase RNA yield, the manufacturer's protocol was modified according to the suggestion reported by Wang [80]. After isolation, total RNA was purified by using the RNeasy ® Mini Kit (Qiagen, Hilden, Germany), following the RNA clean-up protocol. RNA concentration was evaluated by spectrophotometry (Ultrospec 3000, Amersham Pharmacia Biotech, Little Chalfont, UK). Sample quality was checked by using the RNA 6000 Nano kit and the Agilent 2100 Bioanalyzer (Agilent Technologies, San Diego, CA, USA), according to manufacturer's instructions.

RNA-Seq Analysis
Sequencing libraries were generated using a NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) at Novogene (Hong Kong) and index codes were added to attribute sequences to each sample. The library were sequenced on an Illumina HiseqX platform and 150 bp paired-end reads were generated. Three biological replicates were used for RNA-Seq analysis. Raw data (raw reads) of fastq format were processed and clean reads were obtained by removing reads containing adapter and low-quality reads from raw data, running Sickle (https: //github.com/najoshi/sickle), and trimmed for the presence of residual adapter sequences through Scythe (https://github.com/vsbuffalo/scythe). All the downstream analyses were based on the clean data with high quality. Reference genome and gene model annotation files were directly downloaded from Genbank (https://www.ncbi.nlm.nih.gov/genome/?term=glycine+max). The index of the reference genome was built using Bowtie v2.2.3 and paired-end clean reads were aligned to the reference genome using TopHat v2.0.12. We selected TopHat as the mapping tool because it could generate a database of splice junctions based on the gene model annotation file, and thus a better mapping result than other non-splice mapping tools. HTSeq v0.6.1 was used to count the reads numbers mapped to each gene. A similarity matrix of the control and treatment (+biostimulant) data was built up calculating 1 cosine distance (http://amp.pharm.mssm.edu/clustergrammer). Differential expression analysis of two conditions (three biological replicates per condition) was performed using the DESeq R package (1.18.0). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p < 0.05 found by DESeq were assigned as differentially expressed. A gene ontology (GO) enrichment analysis of differentially expressed genes was implemented using the enrichment term engine (GO terms, KEGG pathways/INTERPRO domain) implemented in STRING (https://string-db.org) and GO terms were visualized with a hierarchical clustering tree function summarizing the correlation among significant GO terms, using ShinyGO (http://bioinformatics.sdstate.edu/go). STRING database collects, scores, and integrate all publicly available sources of protein and protein-protein interaction information. Moreover, it implements many basic services such as GO, KEGG, and INTERPRO enrichment functions together with more complex analyses (interactome). To characterize the function of genes differentially expressed in treated plants, the soybean genes were organized in lists of upregulated and downregulated genes and compared to Arabidopsis proteome (TAIR version 10), using reciprocal blast best hits, to find species orthologs. The list of orthologs was submitted to STRING and enrichments were recorded when terms were more enriched in the set of query proteins than the background, considering a FDR value less than 0.05 [81]. An interactome map was built using DEG (up-/downregulated) based on the STRING ® database, with known and predicted protein-protein interactions (PPI) and the networks were built accordingly. The number of observed edges was compared with the expected number of edges, which gives how many edges is to be expected if the nodes are to be selected at random, and PPI enrichment p-values were calculated. A small PPI value indicates that the nodes are not random and that the observed number of edges is significant.

Evaluation of Antioxidant Enzyme Activity of Seeds Incubated under Controlled Conditions
Total proteins were extracted, according to a previously reported protocol [41]. Briefly, 500 mg of powdered seeds were used for the extraction, and all the steps were carried out at 4 • C. The extraction buffer included 50 mM sodium phosphate (pH 7.5), 250 mM sucrose, 1.0 mM EDTA, 10 mM KCl, 1 mM MgCl 2 , 0.5 mM phenylmethylsulfonyl fluoride (PMSF), 0.1 mM dithiothreitol (DTT), and 1% (w/v) polyvinylpolypyrrolidone (PVPP) in a 1:10 (w/v) proportion. The homogenate was mixed by pipetting, and then centrifuged 20 min at 25,000 g and at 4 • C. The supernatant, containing total proteins, was used for enzymatic assays. Soluble protein content was evaluated by the method of Bradford [82] using bovine serum albumin (BSA) as a standard. 4.5.1. Superoxide Dismutase Activity (Sod, EC 1.15.1.1).
SOD activity was evaluated, as previously described [83]. Briefly, 1 mL containing 50 mM sodium phosphate buffer (pH 7.8), 13 mM methionine, 75 µM nitro blue tetrazolium (NBT), 2 µM riboflavin, and 0.1 mM EDTA was incubated with an appropriate dilution of enzyme extract. The samples were placed 30 cm under a light source (4000 lux) for 15 min. In parallel, two different blanks were prepared, i.e., one without the enzyme extract that was placed under the light and that totally developed the reaction; the second one, containing the enzyme extract but placed in the dark, with the aim to avoid the reaction, and therefore was used as control. The absorbance of both samples and blanks was detected at 560 nm. 4.5.2. Catalase Activity (CAT, EC 1.11.1.6).
Cat Activity was evaluated, as previously described [84]. Briefly, the absorbance at 240 nm was monitored for 120 sec with the aim to evaluate the change due to the decreased absorption of The enzyme activity was evaluated by monitoring the absorbance variation at 340 nm for 15 min of the 1-chloro-2,4-dinitrobenzene (CDNB) used as substrate [65]. The assay was performed in 1 mL of reaction solution containing 100 mM potassium phosphate buffer (pH 7.0), 1 mM reduced glutathione (GSH), 1 mM 1-chloro-2,4-dinitro-benzene (CDNB) (10 mM CDNB dissolved in 50% acetone stock solution), and enzyme extract. The reaction was started by adding CDNB.

Evaluation of Non-Protein Thiol Content of Seeds Incubated under Controlled Conditions
The assay was carried out by mixing 500 µL of crude extract, obtained as described before, to 100 µL of 25% (w/v) trichloroacetic acid (TCA). The samples were centrifuged at 12,000× g for 20 min, at 4 • C. Then, 300 µL of supernatant were added to 2.7 mL of 0.6 mM 5,5'-dithiobis (2-nitrobenzoic acid) (DTNB) prepared in 0.1 M sodium phosphate buffer (pH 8.0). The spectrophotometer was used to quantify the sulfhydryl content. The absorbance was detected at 412 nm [65].

Evaluation of Hydrogen Peroxide Levels of Seeds Incubated under Controlled Conditions
The hydrogen peroxide level was detected, according to the protocol described for the first time by Velikova and colleagues [85]. Briefly, powdered seeds (0.5 g) were homogenized with 5 mL of 0.1% (w/v) TCA. The samples were centrifuged at 12,000 g for 15 min and 0.5 mL of supernatant was added to 0.5 mL 10 mM potassium phosphate buffer (pH 7.0) and 1 mL 1 M KI. The absorbance was read at 390 nm and the H 2 O 2 content was determined based on a standard curve.

Statistical Analysis
For biometric and antioxidant system analyses, five biological replicates were used for the statistical treatment of the data. Results are expressed as mean values ± standard deviation (SD). Significance of differences observed in datasets was tested by one-way ANOVA followed by Tukey's test, or by Student's t-test at p < 0.05 using SYSTAT 10 software.

Conclusions
The results of this study displayed how seed pretreatment with the biostimulant had a positive effect on soybean grown under heat stress. In particular, under controlled conditions, this biostimulant was able to enhance the germination rate of soybean seeds incubated at 35 • C, although the effects were visible only at later stages (72 h after imbibition). Furthermore, RNA-Seq analysis carried out on seeds incubated 24 h at 35 • C showed the modulation of about 900 genes (51 upregulated and 828 downregulated). Interestingly, more than one third of the upregulated genes encoded for different methyltransferases, proteins involved in abiotic stress response and in various processes, including DNA repair, a mechanism crucial for a successful seed germination. All of them seemed to act on rRNA raising the hypothesis of a translation-dependent layer of regulation, which deserves to be deeply investigated. Some metabolic pathways, such as carbohydrate metabolism, response to stress, and hormone signaling were negatively regulated by the treatment. This unexpected downregulation might be related to the seed pretreatment activity of the biostimulant to decrease the water uptake in favor of cellular damage prevention. In addition, the antioxidant system activity (ROS scavenging enzymes and non-protein thiol content) evaluated on the same seeds employed for RNA sequencing, revealed a lower amount of hydrogen peroxide in biostimulant-treated seeds, correlated to lower activities and lower expression levels of the corresponding detoxification enzymes, thus confirming the protective action of this biostimulant.