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).

PlantCARE 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 SoyBase website (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 log₂-transformed and displayed in the form of heatmap using the OmicShare tools (http://www.omicshare.com/tools).

Soybean seeds (Glycine max cv Dongnong 50) were germinated on two layers of filter paper soaked in distilled water at 22°C in 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 humidity of 60–80% and 16 h light period/day at intensity of 120 μmol photons m⁻² s⁻¹) for 2 weeks. When the second trifoliate leave unrolled, soybean seedlings were transferred to 1/2 Hoagland solution. 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°C growth 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 the following reaction conditions: 95°C for 20 s; followed by 40 cycles of 95°C for 15 s, 60°C for 20 s and 72°C for 20 s. Samples for qRT-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 expression was calculated using the ΔΔCt algorithm (Livak and Schmittgen, 2001; Bustin et al., 2009).

GmSMO1, GmSMO2, and GmSUMO4 genes were amplified by PCR using soybean cDNA as template. After sequencing, the amplified 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, and centrifuged 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.

Total 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 to nitrocellulose membrane according to standard protocols. The membrane was blocked using phosphate buffered saline (PBS) buffer 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 buffer containing 1% BSA (1:4000) at 4°C for overnight. After removing unbound 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.

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 encoded proteins sharing 92 and 80% similarity with that from AtSUMO1, 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 and Figure 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 observations made 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 selective SUMOylation (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 in aerial 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 the SUMO 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, 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 12 h. 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 transcript 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.

In 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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.