Organization and Regulation of Soybean SUMOylation System under Abiotic Stress Conditions

Covalent attachment of the small ubiquitin-related modifier, SUMO, to substrate proteins plays a significant role in plants under stress conditions, which can alter target proteins’ function, location, and protein-protein interactions. Despite this importance, information about SUMOylation in the major legume crop, soybean, remains obscure. In this study, we performed a bioinformatics analysis of the entire soybean genome and identified 40 genes belonged to six families involved in a cascade of enzymatic reactions in soybean SUMOylation system. The cis-acting elements analysis revealed that promoters of SUMO pathway genes contained different combinations of stress and development-related cis-regulatory elements. RNA-seq data analysis showed that SUMO pathway components exhibited versatile tissue-specific expression patterns, indicating coordinated functioning during plant growth and development. qRT-PCR analysis of 13 SUMO pathway members indicated that majority of the SUMO pathway members were transcriptionally up-regulated by NaCl, heat and ABA stimuli during the 24 h period of treatment. Furthermore, SUMOylation dynamics in soybean roots under abiotic stress treatment were analyzed by western blot, which were characterized by regulation of SUMOylated proteins. Collectively, this study defined the organization of the soybean SUMOylation system and implied an essential function for SUMOylation in soybean abiotic stress responses.

INTRODUCTION
The conjugation of protein tags or modifiers to target proteins is one important protein post- translational modification for the maintenance of cellular homeostasis and biological activities in eukaryotes during rapid environmental changes, when physiological responses often occur both quickly and reversibly. In recent years, the small ubiquitin-related modifier (SUMO) proteins of approximately 100–115 amino acids, has emerged as an influential class of protein modifiers which has been shown to regulate various nuclear processes, including transcriptional control, subcellular trafficking, and regulation of the cell cycle (Miller et al., 2010, 2013; Park H. C. et al., 2011).The modification process of SUMOylation is mechanistically similar to the conjugation of ubiquitin, the most prominent representative of small protein modifiers with conserved structure (Rodriguez et al., 2001; Matic et al., 2010). SUMO is synthesized as a pre-protein that needs to be processed by ubiquitin-like proteases (ULP) to expose a carboxyl-terminal double glycine (di-Gly) motif. Analogous to ubiquitination, SUMO conjugation in plants is facilitated through a multi- step enzymatic pathway including E1 enzymes for SUMO activation, E2 enzymes for conjugation and E3 enzymes for ligation (Castro et al., 2012). The end result is SUMO isopeptide conjugatedvia its C-terminal diGly motif to accessible lysine (Lys) within substrate proteins. Interestingly, SUMOylation is a reversible process, and de-SUMOylation is accomplished through ULP class of SUMO proteases which release free SUMOs from protein conjugates for further conjugation cycles (Colby et al., 2006; Mukhopadhyay and Dasso, 2007).Unlike ubiquitylation, SUMOylation is not typically committing proteins to breakdown by the 26S proteasome, but can change target protein localization, repress or enhance activity of target transcription factor, regulate interactions between proteins, and mediate cross talk with other post- translational modifications (Bossis and Melchior, 2006).

SUMO conjugation was reported to be essential in Arabidopsis, and an increasing number of papers have demonstrated that mutants that were not able to attach SUMO1/2 onto substrate proteins exhibit typical phenotypes such as growth reduction, impaired salinity, drought, heat and freezing tolerance, and altered response to phosphate starvation (Roden et al., 2004; Catala et al., 2007; Miura et al., 2007b, 2011; van den Burg et al., 2010). Previous studies with rice (Oryza sativa) and poplar (Populus strichocarpa) showed that SUMOylation can be significantly enhanced upon heat, cold, high salinity, and abscisic acid exposure (Chaikam and Karlson, 2010; Reed et al., 2010; Li et al., 2013), which was always accompanied by a decrease in the pool of free SUMO and correlated with the duration and intensity of the stress. Interestingly, SUMOylation levels decrease within hours or even minutes when the heat shocked plants returned to the normal temperature, suggesting that SUMOylation acts transiently and reversibly (Kurepa et al., 2003; van den Burg et al., 2010). Recently, SUMOylation was reported to display a striking memory that dampened succeeding conjugate increases if a second stimulus was provided too close to the first one (Kurepa et al., 2003; Miller et al., 2013).Further genetic studies provided evidence that severalmembers of SUMOylation system contributed to the regulation of development and adaptation to abiotic stress conditions. SUMO1 and SUMO2 in Arabidopsis appeared to be functionally redundant, but elimination of them was embryo lethal (Park H. J. et al., 2011).

On the other hand, over-expression of either AtSUMO1 or AtSUMO2 gene was correlated with ABA-mediated growth inhibition (Tatham et al., 2001; Aguilar-Martinez et al., 2015). Besides, null E3 ligase SIZ1 alleles (siz1-1, siz1-2, and siz1-3) in Arabidopsis displayed several phenotypes including hypersensitive to extreme temperature and drought stress, as well as enhanced tolerance to salinity stress, indicating that SIZ1-mediated SUMOylation play an important role in response to abiotic stress (Miura et al., 2007b, 2011; Cheong et al., 2009; Chen et al., 2011). Moreover, Arabidopsis mutants lacking SUMO proteases OTS1 and OTS2 showed inhibition of root growth under salt treatment (Conti et al., 2008), while transgenic rice plants over-expressing OsOTS1 gene exhibited increased salt tolerance concomitant with reduced levels of SUMOylated proteins (Srivastava et al., 2016). It indicated that manipulation of conjugation or de-conjugation of SUMO from target proteins could be an effective method in plants to cope with stress conditions.Soybean is one of the most important legume crops in the world, offering high-quality protein and vegetable oil for human and animal consumption, but the growth and development of soybean seedlings is always affected by adverse conditions (Nakashima et al., 2009; Schmutz et al., 2010). Given the potential importance of SUMOylation to stress defense, presented here is a comprehensive description of core components of SUMO system in soybean. We also investigated the responses of SUMOylation system under salt, high temperature and abscisic acid stress conditions both at transcriptional and translational levels.

Our results may suggest a role for SUMOylation in soybean abiotic stress responses, and provide valuable information for further functional studies of SUMOylation cascade components in soybean.Soybean SUMO pathway genes were identified by tblastn searches against the soybean genome at Phytozome v11.0 (https://phytozome.jgi.doe.gov/pz/portal.html), using the amino acid sequences of the known pathway genes from Arabidopsis as queries. The molecular weight and isoelectric point (pI) of each gene were calculated by ExPASy (http://www.expasy. org/tools/). Protein subcellular localization was predicted by the WoLF PSORT (http://www.genscript.com/wolf-psort.html) (Horton et al., 2007).Multiple sequence alignments of SUMO pathway components in soybean, Arabidopsis, poplar, rice, maize and tomato were carried out using clustalW with default parameters (Thompson et al., 2002). The phylogenetic tree was constructed by MEGA (V7.0.20) using the neighbor-joining (NJ) method with the following parameters: Poisson correction, pair-wise deletion and 1000 bootstrap replicates (Tamura et al., 2007). Protein domains were predicted using Pfam (http://pfam.xfam.org) and visualized using BoxShade online software (http://www.ch.embnet.org/ software/BOX_form.html) (Beitz, 2000). In addition, possible SUMOylation sites and SIM sequences were predicted using the medium settings of program GPS-SUMO v1.0.1 (Zhao et al., 2014).Cis-Acting Elements Analysis in Promoters of the Soybean SUMO Pathway GenesPlantCARE online program (http://bioinformatics.psb.ugent.be/ webtools/plantcare/html; Lescot et al., 2002) was used to predict cis-acting elements in promoters of soybean SUMO pathway genes. Sequence of 1,500 bp upstream of the start codon downloaded from Phytozome 11.0 database was used for cis- acting elements scan.The normalized transcriptome data (Reads/Kb/Millon, RPKM) of 14 soybean samples from six different tissues harvested at different growth periods was downloaded from the SoyBasewebsite (https://soybase.org/soyseq/). These samples include young leaf, flower, one cm pod, pod-shell 10-DAF (days after flowering), pod-shell 14-DAF, root, nodule and developing seeds (seed harvested at 10-DAF, 14-DAF, 21-DAF, 25-DAF, 28-DAF,35-DAF, and 42-DAF) (Severin et al., 2010).

The criteria of expressed gene was the RPKM value equal to or greater than two in the expression atlas (Belamkar et al., 2014). The RPKM normalized read count data of expressed gene was then log2- transformed and displayed in the form of heatmap using the OmicShare tools (http://www.omicshare.com/tools).Plant Growth and TreatmentsSoybean seeds (Glycine max cv Dongnong 50) were germinated◦on two layers of filter paper soaked in distilled water at 22 Cin the dark. After 2 days, germinated seedlings were transferred to vermiculite and irrigated with half-strong modified Hoagland nutrient solution (pH 5.8) in a growth chamber under normal◦conditions (25/21 C day/night temperature, relative humidityof 60–80% and 16 h light period/day at intensity of 120 µmolamplified products were added Sac I and Hind III restriction sites to clone into the pET32a expression vector. The resultant vectors were transformed into E. coli BL21 cells, which were induced for expression by adding 1 mM IPTG after achieving OD600 0.6 in◦LB broth. Cells were collected after 4 h of growth at 28 C, andcentrifuged at 4000 × g for 30 min then resuspended in PBS buffer. After sonication, the soluble protein was used directly for SDS-PAGE or purification using Ni Sepharose 6 Fast Flow resin (GE Healthcare) according to protocols in the manual.Protein Extraction and Western Blot AnalysisTotal protein samples from three biological replicates were isolated from soybean root tissues using previously described method (Cai et al., 2017). Protein concentration was determined by Bradfrod method using BSA as the standard. 15 µg of purified GmSUMO protein and each protein sample of soybean roots was separated on 12% SDS−PAGE and transferred tonitrocellulose membrane according to standard protocols. The—photons ms−1) for 2 weeks.

When the second trifoliate leavemembrane was blocked using phosphate buffered saline (PBS)unrolled, soybean seedlings were transferred to 1/2 Hoaglandsolution. Stress treatments were performed as follows: for salt treatment, the roots of soybean plants were immersed in nutrient solution containing 200 mM NaCl at room temperature; for heat◦treatment, hydroponic seedlings were transferred into a 42 Cgrowth chamber; for ABA treatment, soybean seedlings were grown in nutrient solution added 100 µM ABA. Root tissues were harvested after 0, 0.5, 1, 6, 12, and 24 h exposure to salt, heat, and ABA treatment. Three independent sets of control and stress treatment samples were collected, and each replicate represented a pooled sample of three individual plants. After collection, all the samples were immediately frozen in liquid nitrogen and stored at−80◦C until use.The total RNA samples were isolated using Ultrapure RNA Kit (CWBIO, China). cDNA was synthesized using 2 µg RNA using the HiScript II Q Select RT SuperMix for qPCR Kit (Vazyme, China). Gene specific primers for quantitative real-time RT-PCR (qRT-PCR) analysis were designed using Primer 5.0 according to soybean cDNA sequences (Table S1). The soybean actin 11 was used as internal reference gene. qRT-PCR reaction was performed using ChamQ SYBR qPCR Master Mix (Vazyme, China) and was conducted on ABI 7300 Real-time Detection System (Applied Biosystems, USA). The PCR reaction was carried out with thefollowing reaction conditions: 95 C for 20 s; followed by 40 cyclesbuffer containing 5% milk for 1 h at room temperature, then probed with AtSUMO1 primary antibody (AbcamInc., China, ab5316; a rabbit polyclonal antibody that reacts with ArabidopsisSUMO1/2) which was prepared in the PBS buffercontaining 1% BSA (1:4000) at 4 C for overnight. After removingunbound antibodies by washing with a PBST buffer (PBS containing 0.05% Tween 20), the blot was incubated with goat antirabbit IgG secondary antibody (horseradish peroxidase conjugated, Thermo FisherScientific, USA) in the PBST buffer at a dilution of 1:20000. Chemiluminescence was carried out using the ECL Plus kit according to manufacturer’s (Merck Millipore, USA) instructions.

RESULTS
Using homology searches and domain confirmation, a total of 40 genes encoding SUMO pathway members were found in soybean, including six family groups. The detailed characteristics of each family member were shown in Table 1.Our soybean list included six SUMO genes, designated as GmSUMO1 to GmSUMO6, all of which had three exons and encoded proteins harboring a SUMO-type β-grasp domain and predicted to localize in the nucleus (Table 1). Among these SUMO genes, GmSUMO2 and GmSUMO3 encodedproteins sharing 92 and 80% similarity with that fromof 95 C for 15 s, 60 C for 20 s and 72 C for 20 s. Samples forqRT-PCR were run in 3 biological replicates with 3 technical replicates and the data were represented as the mean ± SD (n= 3) for Student’s t-test analysis. The relative gene expressionwas calculated using the ∆∆Ct algorithm (Livak and Schmittgen,2001; Bustin et al., 2009).Prokaryotic Expression and Purification of Soybean SUMO ProteinsGmSMO1, GmSMO2, and GmSUMO4 genes were amplified by PCR using soybean cDNA as template. After sequencing, theAtSUMO1, while GmSUMO1 encoded a protein with 85% identity to AtSUMO2. By contrast, GmSUMO4, GmSUMO5, and GmSUMO6 encoded products with little resemblance to either AtSUMO1/2 or GmSUMO1/2/3. To better understand the evolutionary relationship between these SUMO isoforms, we searched a number of other plant genomes for related sequences and constructed a phylogenetic tree (Figure 1A). Phylogenetic analysis revealed a highly conserved SUMO group including GmSUMO1/2/3, AtSUMO1/2, ZmSUMO1a/b and conserved representatives from other species. In contrast, a “non- canonical” group that included GmSUMO4/5/6, AtSUMO3/5,and ZmSUMO2 encoded proteins shared about 30% amino acid similarities to their canonical paralogs (Figure 1A). Regardless of their dissimilarities, it is noteworthy that non-canonical members also have the C-terminal di-Gly motif necessary for conjugation (Figure 1B) which indicated their potential abilities of covalent attachment to target proteins. In support, we identified some residues important for the non-covalent binding of SUMOs to E1, E2, and SUMO-interacting motifs (SIMs, ψψXψD/S/E, D/S/EψXψψ or ψψDLT, where ψ stands for the hydrophobic amino acids and X represents any residue) conserved across the biological kingdoms (Figure 1C).

Besides these identified SUMO proteins above mentioned, two longer SUMO related proteins designated as SUMO-variants(V) were detected in soybean genome (Figure 1). GmSUMO- Vs contained the N-terminal extension of more than 200 residues upstream of the SUMO β-grasp domain (Table 1 andFigure 1B). Without an exposed diGly motif (Figures 1B,C), SUMO-v may not enter the well-studied E1/E2/E3 SUMOylation pathway, but probably have other functions to be explored in plants.The initial activation of SUMOylation is adenylation of mature SUMO on the C-terminus diGly by the heterodimeric SAE1/SAE2 E1. Soybean genome expressed two small subunit SAE1 isoforms (GmSAE1a/b) and two large subunit SAE2 isoforms (GmSAE2a/b) (Table 1). The ubiquitin-fold domain and cysteine domain important to E1 catalysis were shown in Figure S1. After being activated by E1, the bound SUMO is transferred to an active-site cysteine in the E2 by transesterification (Miura et al., 2007a). In contrast with the single SCE1 gene in Arabidopsis, the soybean genome contained four E2 genes named GmSCEa, GmSCEb, GmSCEc, and GmSCEd (Table 1). Conformed with the observationsmade earlier in the homologs AtUbc9 (Kraft et al., 2005) and OsSce1 (Nigam et al., 2008), all of these four genes were predicted to encode proteins localized predominantly to the nucleus along with SUMOs, however, whether GmSCEs function related to the nuclear pore complex is still not known and deserved further studies. Motif-scan analysis of GmSCEa-d revealed that they have highly conserved Ub-Conjugating Enzyme Catalytic (UBC) domain (spanning from 8 to 150aa residues; Figure S2A), which was common to UBC family members (Nigam et al., 2008). Phylogenetic analysis reflected that soybean SCE genes had high similarities with homologs in other species (Figure S2B) highlighted the conservation of E2 proteins. Collectively, we speculated that soybean homologs of SUMOylation system will be biochemically functional.

Although SCE itself can directly transfer SUMO to substrates, E3 ligases are often required to increase the rate of SUMOylation and increase the substrate specificity during SUMOylation process (Gareau and Lima, 2010). When searching the soybean genome using E3 genes in Arabidopsis, we got seven soybean E3 genes (Table 1), all of which contain the conserved Siz1/PIAS(SP)-RING domain which plays the role of binding to the SUMO-E2 intermediate. These E3 genes could be divided into three groups, including SIZ1, Methyl Methane sulfonate-Sensitivity protein-21/High Ploidy-2 (MMS21/HPY2), and Protein Inhibitor of Activated STAT (PIAS)-Like (PIAL) (Figure 2A). The soybean SIZ1 group contained four proteins (GmSIZ1a-d) that included the signature MIZ/SP-RING domain, as well as PINIT and SXS motifs that promoted selectiveSUMOylation (Figure 2B; Yunus and Lima, 2009), which was reported to be important for Arabidopsis SIZ1 activity (Cheong et al., 2009). Phylogenetic analysis of SIZ1 type SUMO E3 ligases from various plant species showed that GmSIZ1a-d were evolutionarily most closely to PtSIZ1a/b (Figure 2A).Thirteen SUMO protease genes in soybean were classified into four groups (Table 1 and Figure 3), all of which have the signature peptidase_C48 domain characterized by a catalytic triad of His-Asp-Cys that leads to peptide bond cleavage (Colby et al., 2006; Mukhopadhyay and Dasso, 2007). For class A, there was only one member in soybean, confirmed the notion that most plants encoded a single gene of this class (Figure S3; Novatchkova et al., 2012). For GmESD4a-d, several SIM sequences were identified that could bind SUMO (Figure S4).

For the OTS family in soybean, the predicted GmOTSa encoded protein contained intact peptidase_C48 domain, whereas two loci Glyma.08G287800 and Glyma.08G287700 in soybean genome encoded N-terminal and C-terminal of the peptidase_C48 domain respectively (Figure 3), and whether either of them could exhibit SUMO protease activity is not known and deserved further study.In addition, soybean genome contained another family of SUMO chain binding protein genes characterized by a tandem arrangement of four SIMs (Table 1 and Figure S5), which therefore have specific affinity for binding to SUMO chains. Besides, these proteins share a RING domain which was speculated to direct proteins with SUMO chains into the ubiquitin-proteasome degradation pathway (Plechanovová et al., 2011; Praefcke et al., 2012).In order to get some information of stimulus induced as well as temporal and spatial specific expression patterns of these soybean SUMO pathway genes, 1,500 bp promoter regions upstream of all of these 40 soybean genes were extracted and used for cis- acting elements searches (Table S3). As shown in Figure 4, most genes involved in SUMOylation in soybean contained several environmental stress-related elements, including cis-elements that were related to plant responses to heat (HSE), drought (MBS), low temperature (LTR), anaerobic stress (ARE), and fungus (BOX-W1). Besides, several promoters contained cis- elements related to MeJA (CGTCA/TGACG-motif), salicylic acid (TCA-element), ABA (ABRE), GA (TATC-box, GARE- motif and P-box), auxin (AUXRR-core and TGA-element), and ethylene (ERE) (Figure 4). These results suggested that these genes were likely to be transcriptionally regulated upon environmental stresses and hormone signals. Moreover, seven types of cis-acting elements were involved in plant meristem (CAT-box and CCGTCC-box), endosperm (SKn-1 motif), leaf and shoot development (HD-zip2 and as-2-box), cell cycle (MSA- like) and circadian control (circadian), which covered most plant developmental stages (Figure 4).

It suggested that SUMO pathway genes may be transcriptionally regulated during growth and development of soybean plants.Gene functions are always closely related to where and how they are transcriptionally expressed. Transcript profiles from 14 tissues/organs were collected from soybean RNA-Seq data and downloaded from the Soybase (Table S5). A total of 30 SUMO pathway genes were expressed at a variety of developmental stages in soybean (Figure 5). Among these SUMO genes, GmSUMO1 was expressed at very high levels in almost all tissues, GmSUMO2 and GmSUMO3 showed high expression levels inaerial and underground tissues, but much lower expression levels in seed developing stages. It suggested that GmSUMO1 could be functional complement to the other two orthologs. In the E1 group, GmSAE1a/b and GmSAE2a/b genes exhibited relatively higher expression levels in aerial tissues, while low expression levels during seed development. Notably, all of the E2 genes showed high expression levels in all of these 14 tissues, suggesting this family of genes may support soybean normal growth and development. Interestingly, expression levels of E3 family members in some tissues were quite different. For example, expressions of GmSIZ1c and GmPIAL1 were hardly detected in nearly all tissues, while GmSIZ1a showed relatively high expression levels in most tissues, besides, GmSIZ1b showed higher expression levels in flower and root, but GmMMS21 expressed at relatively high levels in young leaf, one cm pod and nodule (Figure 5). Furthermore, our results showed that ULP family gene GmPROa expressed at low levels in almost all developmental stages, except the young leaf.

Another interesting scenario was that GmB2d and GmESD4a/b/c were highly expressed in nodule, suggesting their potential roles in nodulation although experimental evidence is lacking for their involvement in nodule organogenesis or development. As shown in Figure 5, GmCBa and GmCBb exhibited similar expression patterns, suggesting potential redundant functions between the orthologs. Additionally, GmCBc was expressed at very low levels in almost all tissues, on the contrary, GmCBd gene showed strong expression in all of these investigated tissues. Given the results that most genes studied exhibited relatively higher expression levels in roots, the following transcriptional and translational expression profiles were performed using soybean root tissue.In order to help appreciate the potential roles and transcriptional regulations of the soybean SUMO system members under abiotic stress, we examined the expression changes of 13 SUMO pathway genes in soybean roots in response to NaCl (200 mM), heat (42◦C) and ABA (100 µM) treatments.Quantitative-PCR analysis revealed that the majority of theSUMO cascade component transcript levels were differentially regulated under salinity conditions (Figure 6A). GmSUMO2 transcripts decreased to very low levels after exposure to salt stress, whereas, GmSUMO5 transcripts showed up-regulated after treated for 6 h. The expression of GmSAE1a was slightly increased after 6 h of treatment, and was maintained nearly at initial levels after 24 h of salinity stress (Figure 6A). Among the E2 enzymes, GmSCE1a showed a transient increase, whereas, GmSCE1d transcript was significantly increased since treated for 1 h, and reached a maximum at 12–24 h which was induced about 50-fold. Among the E3 genes, GmSIZ1a began to increase after 1 h of treatment, whereas, GmMMS21 showed a transient increase and decreased when treated for 24 h, which exhibited similar expression patterns with GmOTSa gene. These increases in gene expression may contribute to the accumulation of SUMOylated proteins after exposure to salt stress. Besides,decreased significantly, thereby indicating a specific role for GmSUMO1 in ABA stress response.

All of the E2 and E3 family genes transcript levels increased after 6 h of treatment, and GmMMS21 showed a transient up-regulation followed by restoration to basal levels. Among the ULP family, GmESD4a transcript showed up-regulation 1 h after the ABA treatment, indicating that GmESD4a may be responsible for a quick de- SUMOylation process after ABA treatment.GmSUMO3, GmOTSb1, GmESD4a, GmB2d, and GmCBb showed no detectable changes in the 24 h period after salt treatment.When treated with heat stress, GmSUMO1 and GmSUMO2 were found to be induced, whereas, GmSUMO5 showed a significant reduction under heat treatment (Figure 6B). The expression patterns of GmSAE1a under heat stress were similar with that under salinity treatment. Both of the E2 genes were up- regulated under heat stress, and GmSCE1a gene exhibited a faster response than GmSCE1d gene. Among these E3 genes, GmSIZ1a showed dramatic up-regulation 6–24 h after heat treatment; however, GmMMS21 was significantly down-regulated 1 h after treatment, and returned to the initial levels after treated for 12h. In addition, similar expression patterns were observed among GmOTSa, GmESD4a, and GmB2d, which showed increased expression levels from 6 h to 24 h of treatment; while GmOTSb1 showed slightly increased from 6 h under heat treatment and after 12 h showed down-regulation. In contrast, GmCBb transcript showed increased gradually over the 24 h period of treatment. Notably, the majority of genes exhibited a transient decrease during the first 1 h of heat treatment, which indicated that SUMOylation may be depressed temporarily at this time point in soybean roots.After the ABA treatment, expression of most SUMO pathway members showed slightly up-regulation with the exception of GmSUMO2, GmSAE1a, GmOTSa, and GmOTSb1 genes which were maintained nearly basal levels during the 24 h period of treatment (Figure 6C). GmSUMO1 gene showed a transient positive response 6 h after ABA treatment and subsequently returned to initial levels.

Conversely, GmSUMO5 transcriptIn order to explore the regulation of SUMOylation in response to abiotic stress in soybean roots, western blot analysis was performed to monitor the expression changes on the translational level. To test whether anti-AtSUMO1 antibody could detect soybean SUMOs, GmSUMO1, GmSUMO2, and GmSUMO4 genes were amplified (Table S6) and cloned into pET32a expression vector for protein expression. The western blot analysis indicated that anti-AtSUMO1 antibody was able to specifically detect GmSUMO1 and GmSUMO2 proteins, but not for GmSUMO4 (Figure S6).During salt treatment, roots accumulated large amount of low molecular weight SUMO conjugates as soon as 0.5 h after salt stimulus (Figure 7A), and the accumulation of high molecular weight SUMOylated products was detected from 0.5 to 1 h of salt treatment. At the same time, the level of free SUMO was greatly reduced from 1 h of treatment, and when treated for 24 h, most free SUMOs have been exhausted in soybean roots. Consistent with the results got from poplar (Reed et al., 2010), SUMOylation by salinity treatment in soybean roots was rapidly inducible and probably has some functional implications.Rapid stimulation of SUMOylation in response to heat shock had been reported in several plant species (Chaikam and Karlson, 2010; Lopez-Torrejon et al., 2013; Augustine et al., 2016). In this study, SUMO conjugates were detected to be accumulated as soon as 0.5 h after heat stress, which was conformed with previous studies (Cai et al., 2017). Accordingly, the free SUMO pool showed decreased after 1 h of heat shock, and restored to nearly the original level after 24 h of treatment (Figure 7B).After treating with ABA, the high molecular weight SUMO conjugates accumulated significantly at 1 h after treatment, but slightly decreased at later time points, and increased from 12 to 24 h after treatment (Figure 7C), which indicated that ABA treatment effected SUMOylated protein accumulation in soybean roots. These above results further supported the hypothesis that the rapid protein SUMOylation was an essential functional component in response to abiotic stress in soybean.

DISCUSSION
Previous studies in Arabidopsis connected the protein SUMOylation to plant development and defense against environmental challenges (Castro et al., 2012). Thus, manipulation of this post-translational modification system might provide opportunities to modulate the plant stress resistance with lower inputs. However, only a few examinations have been made of SUMO pathway members in soybean (Cai et al., 2017). In this study, 40 SUMO pathway genes in soybean were identified. The soybean genome contained gene models for SUMO and its respective conjugation and de-conjugation machinery including E1, E2, E3, and SUMO protease components that were closely related to homologs in Arabidopsis. Transcriptional and translational expression profiling of soybean SUMOylation system allowed us to identify the cascade members that have crucial roles in the soybean responses to salt, heat and ABA stresses, and explain the role of SUMOylation in soybean abiotic stress responses. Compared with the eight putative SUMO copies in Arabidopsis genome, soybean genome contained six members shared a similar ubiquitin-like structural conformation characterized by a β-grasp fold. In soybean, three canonical SUMO genes were present that expressed polypeptides homologous to canonical members in other eukaryote species, which have been well recorded to conjugate to a variety of target proteins in response to environmental stimuli (Kurepa et al., 2003; Saracco et al., 2007). Among these SUMO genes in soybean, GmSUMO1 was constitutively expressed in most organs and seed developmental stages according to the RNA seq data (Figure 5), which showed similar expression patterns with the counterpart genes in Arabidopsis (van den Burg et al., 2010). On the contrary, the expression of non-canonical SUMO GmSUMO4 in soybean appeared relatively weak based on transcriptional analysis.

Moreover, GmSUMO1 transcript showed a sustained up-regulation 1 h after ABA treatment, while GmSUMO2 and -5 were not induced (Figure 6), suggesting a specific role of GmSUMO1 in ABA response. Based on the high sequence similarity and overlapping expression profiles, we speculated that GmSUMO1, -2, and -3 might be functional redundant, while whether the biological functions of GmSUMO1/2/3 were essential as their homologs in Arabidopsis still need to be identified. The roles of SUMOs under abiotic stress were inferred by the detection of increasing SUMOylation during stress (Catala et al., 2007; Conti et al., 2008; Golebiowski et al., 2009; Miura et al., 2011). In this study, rapid regulations of SUMO-conjugates in soybean roots were identified upon salinity, heat and ABA elicitation (Figure 7). Although most target proteins bear a single SUMO moiety, some are modified with a poly-SUMO chain onto a single Lys in the substrate. Previous study demonstrated that the residue Lys10 in Arabidopsis SUMO2 was required for self-SUMOylation to form the multimeric chains in vitro (Colby et al., 2006). It can be deduced that soybean SUMO proteins except for GmSUMO4 may be able to form poly-SUMO chains, given the existence of conserved Lys residue in these soybean SUMO proteins (Figure 1C). Notably, there were indeed SUMOylated products with the molecular mass of approximately 30 and 45 kDa in soybean SUMOylation profiles (Figure 7), which were speculated to be potential SUMO dimers and trimers in soybean. The poly-SUMO chain was reported to serve as degradation signal, channeling chain-modified target proteins into degradation by SUMO chain binding ubiquitin ligases (Plechanovová et al., 2011; Praefcke et al., 2012). Thus, the appearance of SUMO dimers and trimers in soybean was supposed to be an early response to stimuli that promoted reallocation of SUMO moieties to target proteins.

Besides SUMOs, soybean genome encoded all components participate in SUMOylation and de-SUMOylation cascade. The SUMO-activating enzyme E1 (SAE) enzymes facilitate conjugation of SUMO proteins through a series of reactions including adenylation, formation of E1-SUMO thioester and transfer of thioester from E1 to E2 conjugating proteins (Lois and Lima, 2005). Previous study has demonstrated that SAE2 formed a tight heterodimer with SAE1 and controlled the localization of the SUMO SAE heterodimer to the nucleus (Truong et al., 2012). In this study, nuclear locations of SAE subunits were also found in soybean genome (Table 1). In Arabidopsis, loss-of-function mutant sae2-1 exhibited developmental arrest at the early stages of embryogenesis, suggesting the essential role of SAE2 (Saracco et al., 2007). Given the comparable expression patterns between SAE1a and SAE1b, SAE2a and SAE2b (Figure 5), we suggested that there may be functionally redundancy between these paralogous SAE polypeptide. Strikingly, GmSAE1a gene was found to be enhanced by eight-fold after heat treatment for 12 h, which was presumed to contribute to the accumulation of SUMO-conjugates in responding to unfavorable environment.

GmSCEs belongs to conjugating enzyme family, characterized by a conserved UBC domain that spans most of the protein. In soybean, four sce genes showed high similarity at the amino acid levels (Figure 2), which indicated that these orthologous genes may be a result of duplication events. Putative cis-elements in 1.5 kb region upstream of GmSCE genes possessed higher abundance of stress-related elements like HSE, MBS, ARE, TC-rich repeats and CGTCA/TGACG-motif, which provided useful cues for determining the responsiveness of GmSCE genes to environmental stimuli. Indeed, our analyses showed that GmSCE transcripts were significantly up-regulated under heat, salt and ABA treatments, which was in conformity with the observations made earlier in rice (Nigam et al., 2008). Recent study provided evidence that overexpression of SaSce9 gene from a grass halophyte in plant can improve salinity and drought stress tolerance demonstrated the potential of using this gene in crop improvement (Ratna and Prasanta, 2012). Thus, biological functions of GmSCEa and GmSCEd under abiotic stress are deserved further study. Besides, soybean genome encoded three types of SUMO E3 ligases containing the signature SP-RING domain (Figure 2). The functions of E3 ligases in the regulation of environmental stress responses and developmental processes have been extensively characterized in Arabidopsis (Lois, 2010; Miura and Ohta, 2010) and rice (Li et al., 2013; Mishra et al., 2017). Recently, two identical SIZ1 homologs in soybean, GmSIZ1a and GmSIZ1b, were identified, which were reported to mediate SUMO modification of nuclear proteins and regulate vegetative growth in soybean (Cai et al., 2017).

In this study, SIZ1 gene in soybean was significantly up-regulated in response to heat and salt stress treatments (Figure 6). However, whether GmSIZ1a/b could regulate environmental stress responses in soybean is still needed to be determined. The second type of E3 ligases MMS21 was reported to function as an important regulator of cell proliferation and cytokinin in signaling during root development in Arabidopsis (Huang et al., 2009). Previous work has reported that Arabidopsis plants exposed to ABA accumulate increased levels of SUMOylated proteins dependent on MMS21 activity, because the accumulation of SUMO conjugates is significantly lower in mms21 mutant, and higher in MMS21-OE compared with WT plants (Zhang et al., 2013). In this study, we found that the expression levels of GmMMS21 gene were induced by ABA and salt treatment (Figure 6), it was presumed that GmMMS21 gene could act a pivotal part in salt stress response in soybean by regulating gene expression in the ABA-dependent manner. Moreover, there is another type of E3 ligases containing SP-RING domain and SUMO interaction motifs (SIM) which designated as PIAL1 and 2 that function as SUMO ligases specialized in SUMO chain formation. Mutant studies in Arabidopsis indicated that PIAL1 and 2 functioned in the regulation of salt and osmotic stress responses by removal of SUMOylation substrates, and may contribute to connection of SUMOylation pathway to ubiquitin proteolytic pathway (Tomanov et al., 2014).

However, the biological functions and molecular mechanisms of PIAL1 and 2 in SUMOylation cascade are still largely unknown in soybean. The up-regulation of most of E2 and E3 genes in soybean under NaCl, heat and ABA stimuli (Figure 6) may account for stress induced accumulation of SUMO conjugates which detected by western blot analysis (Figure 7).Compared with other SUMO pathway components, there are a larger number of plant SUMO proteases (Figure 3) which cleave precisely between the terminal Gly of SUMO and the substrate Lys, and play a vital regulatory role during de-SUMOylation (Hickey et al., 2012). So far only a few bona fide SUMO proteases like ESD4 and OTS1/OTS2 have been characterized using genetic, physiological, and biochemical approaches. The Arabidopsis esd4 mutant showed extreme early flowering and alterations in shoot development, suggesting that SUMOylation played roles in plant development which regulated by ESD4 (Murtas et al., 2003). Previous study suggested that the OTS class of ULP-like SUMO proteases may provide routes for boosting crop growth under stress because over-expression of OTS gene led to salt tolerance (Conti et al., 2008). However, SUMO proteases remain largely unstudied, especially in soybean and other crop plants. Similar with AtOTS1, OTS gene expression was enhanced in soybean under heat and salt stress would suggest that soybean stimulated the de-SUMOylation process in response to abiotic stresses, which had been shown in the soybean SUMOylation profiles (Figure 7).

CONCLUSIONS
In summary, a total of 40 genes involved in SUMOylation system from the soybean genome were identified. We investigated the classification, characterization and tissue-specific expression of these SUMO pathway genes, revealing SUMO pathway members broadly participate in the regulation of plant tissue development. Meanwhile, qRT-PCR analysis to 13 SUMO pathway genes suggested that these genes functioned at the transcriptional level during adaption to adverse environments. Moreover, the investigation of SUMOylation profiles in response to salt, heat and ABA treatment demonstrated that soybean responded to abiotic stress treatments by regulating the SUMO conjugated proteins. These results may contribute to the existing knowledge on the complexity of soybean post translational modifications that occur in response to Subasumstat environmental stresses. Further proteome wide identification of various targets will be paramount to defining how SUMOylation ultimately participates in stress protection.