Functional Characterization of Ubiquitination Genes in the Interaction of Soybean—Heterodera glycines

Ubiquitination is a kind of post-translational modification of proteins that plays an important role in plant response to biotic and abiotic stress. The response of soybean GmPUB genes to soybean cyst nematode (SCN, Heterodera glycines) infection is largely unknown. In this study, quantitative real-time PCR (qRT-PCR) was performed to detect the relative expression of 49 GmPUB genes in susceptible cultivar William 82 and resistant cultivar Huipizhi after SCN inoculation. The results show that GmPUB genes responded to cyst nematode infection at 1 day post-inoculation (dpi), 5 dpi, 10 dpi and 15 dpi. The expression levels of GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A were dramatically changed during SCN infection. Furthermore, functional analysis of these GmPUB genes by overexpression and RNAi showed that GmPUB20A, GmPUB33A and GmPUB24A negatively regulated soybean resistance under SCN stress. The results from our present study provide insights into the complicated molecular mechanism of the interaction between soybean and SCN.


Introduction
Ubiquitination is involved in many biological processes, including the cell cycle, hormone signaling, signal transduction, photomorphogenesis and senescence in plants. Ubiquitination modification participates in a series of reactions of the ubiquitin-activating enzyme (E1), ubiquitin-conjugating enzyme (E2) and ubiquitin-ligase enzyme (E3) [1]. Polyubiquitin-modified target proteins are generally degraded by the 26S proteasome, called the ubiquitin-26S proteasome system (UPS) [2]. The UPS is responsible for selecting, targeting and proteolyzing specific degradation substrates. Its components are particularly abundant in plants and are key hubs for regulating various plant cell processes. Hence, the UPS is essential for development and stress response in plants [1,3].
Ubiquitin ligase E3 is classified into two groups: the HECT domain family and the RING domain family. The HECT domain plays a role mainly through the sulfur ester bond necessary for the formation of the catalytic effect with ubiquitin, while the RING domain provides a residential site for E2 and substrate to enable E2 to catalyze the transfer of ubiquitin to the substrate. Recently, a new E3 family, the U-box protein family, was discovered. The characteristic of the U-box E3 ligase is that it has a conserved U-box domain of approximately 70 amino acids. Secondary structure prediction showed that the U-box was a modified RING domain-missing scaffold and conserved zinc-chelated cysteine and histidine residues in a typical RING domain. There are 64 U-box genes in Arabidopsis and 77 in rice. These genes participate in biotic and abiotic stress processes, fertility and hormone the domains of 49 GmPUBs, and these GmPUBs were divided into six categories in the present study. The C-terminus of the U-box domain of type II PUB has an ARM repeat structure, and 21 of 49 soybean PUBs belong to this type. The domain of the type III PUB contained the GKL isomerase domain, including 5 soybean PUBs. Class IV PUBs contain kinase domains, and 5 of the 49 soybean PUBs are classified as such. Class V PUBs contain only U-box domains, 16 of 49 soybean PUBs. Type VI and type VII PUBs contained only one protein. The C-terminus of the type VI PUB contains the WD40 domain, and type the VII PUB contains the TPR structure (Table S1, Figure 2).

Bioinformatics Analysis of GmPUBs
The U-box domain is highly conserved in eukaryotes. In the present study, multiple sequence alignment showed that 49 GmPUB proteins contained a complete and conserved U-box domain (Figure 1). At the same time, according to Wang's method [37], we analyzed the domains of 49 GmPUBs, and these GmPUBs were divided into six categories in the present study. The C-terminus of the U-box domain of type II PUB has an ARM repeat structure, and 21 of 49 soybean PUBs belong to this type. The domain of the type III PUB contained the GKL isomerase domain, including 5 soybean PUBs. Class IV PUBs contain kinase domains, and 5 of the 49 soybean PUBs are classified as such. Class V PUBs contain only U-box domains, 16 of 49 soybean PUBs. Type VI and type VII PUBs contained only one protein. The C-terminus of the type VI PUB contains the WD40 domain, and type the VII PUB contains the TPR structure (Table S1, Figure 2).  Then, the physicochemical properties of the amino acid sequences of 49 GmPUB family members were analyzed using the online tool Protparam of ExPASy Proteomics. The results show that the relative molecular weights of GmPUBs were concentrated in 23.65-166.32 kD (Table 1); the highest theoretical isoelectric point was Glyma12g188900 (9.11), and the lowest was GlymaU008400 (4.7). Furthermore, to compare the genetic relationships among the 49 GmPUB genes, the amino acid sequences of 49 GmPUB proteins were downloaded from NCBI and used for constructing a phylogenetic tree ( Figure 2). Furthermore, to compare the genetic relationships among the 49 GmPUB genes, the amino acid sequences of 49 GmPUB proteins were downloaded from NCBI and used for constructing a phylogenetic tree ( Figure 2).

Expression Pattern Analysis of GmPUB Genes under Nematode Stress
The qRT-PCR results demonstrated that all 49 GmPUB genes changed at four time points, but the expression levels were different ( Figure 3). The expression level of six genes

Expression Pattern Analysis of GmPUB Genes under Nematode Stress
The qRT-PCR results demonstrated that all 49 GmPUB genes changed at four time points, but the expression levels were different ( Figure 3). The expression level of six genes changed significantly: Glyma01g131800, Glyma07g106000, Glyma03g088400, Glyma15g038600, Glyma02g195900 and Glyma12g188900. At the same time, according to the relationship between the proteins encoded by these genes and PUB proteins in Arabidopsis, we renamed them GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A (Figure S2). In HPZ, four genes (GmPUB16A, GmCHIPA, GmPUB23A and GmPUB24A) showed a trend of first decreasing and then increasing, and the relative expression levels were the lowest at 10 dpi. The expression level of GmPUB20A increased with inoculation time. The relative expression of the GmPUB33A gene first decreased and then increased. In Williams 82 (W82), the relative gene expression of GmPUB16A and GmPUB24A after SCN infection increased first and then decreased and peaked at 5 dpi. GmCHIPA showed a trend of first decreasing and then increasing. GmPUB23A showed an upward trend after longer induction times, while GmPUB20A and GmPUB33A showed a downward relative gene expression trend ( Figure 4). pression levels were the lowest at 10 dpi. The expression level of GmPUB20A increased with inoculation time. The relative expression of the GmPUB33A gene first decreased and then increased. In Williams 82 (W82), the relative gene expression of GmPUB16A and GmPUB24A after SCN infection increased first and then decreased and peaked at 5 dpi. GmCHIPA showed a trend of first decreasing and then increasing. GmPUB23A showed an upward trend after longer induction times, while GmPUB20A and GmPUB33A showed a downward relative gene expression trend ( Figure 4).  without soybean cyst nematode infect at 1 dpi, respectively. H1S and W1S represent the relative expressions of GmPUBs after soybean cyst nematode infection at 1 dpi, respectively. Likewise, the samples of H5C, W5C, H5S and W5S were taken at 5 dpi. The samples of H10C, W10C, H10S and W10S were taken at 10 dpi. Additionally, the samples of H15C, W15C, H15S and W15S were taken at 15 dpi. The relative expression value of GmPUBs higher than 1 and lower than 1 represent upregulation and downregulation, respectively. expressions of GmPUBs without soybean cyst nematode infect at 1 dpi, respectively. H1S and W1S represent the relative expressions of GmPUBs after soybean cyst nematode infection at 1 dpi, respectively. Likewise, the samples of H5C, W5C, H5S and W5S were taken at 5 dpi. The samples of H10C, W10C, H10S and W10S were taken at 10 dpi. Additionally, the samples of H15C, W15C, H15S and W15S were taken at 15 dpi. The relative expression value of GmPUBs higher than 1 and lower than 1 represent upregulation and downregulation, respectively. The statistical comparisons were conducted with the t-test for members of treatment and control. ns represent not significant. *, **, *** and **** indicate significant difference at the p < 0.05, p < 0.01, p < 0.001 and p < 0.0001 probability levels, respectively. The statistical comparisons were conducted with the t-test for members of treatment and control. ns represent not significant. *, **, *** and **** indicate significant difference at the p < 0.05, p < 0.01, p < 0.001 and p < 0.0001 probability levels, respectively.

Cloning of Five GmPUB Genes
Specific primers were designed according to the full-length sequence, and PCR was performed to obtain GmPUB genes. The length of the CDS region was as follows: Gm-PUB20A: 1314 bp; GmCHIPA: 837 bp; GmPUB33A: 2454 bp; GmPUB23A: 1257 bp; and GmPUB24A: 1284 bp. The electrophoresis and sequencing results show that five full-length sequences of GmPUBs were successfully obtained. The GmPUB16A was not successfully cloned in different soybean cultivars after several times ( Figure S3).

Multiple Sequence Alignment of GmPUBs
Using the amino acid sequence of the GmPUB gene as a probe, several similar proteins in Arabidopsis thaliana, Theobroma cacao, Capsicum annuum and Oryza sativa were found in the GenBank database, indicating that GmPUB proteins are ubiquitous in higher plants. DNAMAN software was used to align the sequences of the GmPUB proteins with AtPUB23 (At2g35930), AtPUB22 (At3g52450) and AtPUB24 (At3g11840). The amino acid sequences of AtPUB25 (At3g19380), AtPUB20 (At1g66160), AtPUB18 (At1g10560), AtPUB19 (At1g60190), TcPUB23 (Thecc1E022966), CaPUB1 (ABA59556) and OsPUB15 (Os08g01900) were subjected to multiple alignment ( Figure S4), and it was found that the amino acid sequences of proteins coded by the GmPUB genes were highly similar to those of PUB proteins in other species. This suggests that their function is highly conserved.

Phylogenetic Analysis of GmPUBs
The amino acid sequences of PUB proteins from different plants were downloaded from NCBI, and the phylogenetic analysis of GmPUBs ( Figure 5) showed that six GmPUB proteins (GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A) showed relatively high similarities with four Arabidopsis PUB proteins (AtPUB18, AtPUB19, AtPUB22 and AtPUB23). Among them, GmPUB16A, GmPUB20A, GmPUB23A and Gm-PUB24A were clustered in the same large branch, whereas other two GmPUB proteins (GmCHIPA and GmPUB33A) were clustered in the other two branches. On the branches of GmPUB16A, GmPUB20A, GmPUB23A and GmPUB24A, four GmPUB proteins were closely related to AtPUB16, AtPUB20, TcPUB23 and AtPUB24, clustered in the same branch, and evolved in four different directions. In GmCHIPA and GmPUB33A, the former has high homology with AtCHIP, while the latter has high homology with JrPUB33.
performed to obtain GmPUB genes. The length of the CDS region was as follows: GmPUB20A: 1314 bp; GmCHIPA: 837 bp; GmPUB33A: 2454 bp; GmPUB23A: 1257 bp; and GmPUB24A: 1284 bp. The electrophoresis and sequencing results show that five fulllength sequences of GmPUBs were successfully obtained. The GmPUB16A was not successfully cloned in different soybean cultivars after several times ( Figure S3).

Multiple Sequence Alignment of GmPUBs
Using the amino acid sequence of the GmPUB gene as a probe, several similar proteins in Arabidopsis thaliana, Theobroma cacao, Capsicum annuum and Oryza sativa were found in the GenBank database, indicating that GmPUB proteins are ubiquitous in higher plants. DNAMAN software was used to align the sequences of the GmPUB proteins with AtPUB23 (At2g35930), AtPUB22 (At3g52450) and AtPUB24 (At3g11840). The amino acid sequences of AtPUB25 (At3g19380), AtPUB20 (At1g66160), AtPUB18 (At1g10560), AtPUB19 (At1g60190), TcPUB23 (Thecc1E022966), CaPUB1 (ABA59556) and OsPUB15 (Os08g01900) were subjected to multiple alignment ( Figure S4), and it was found that the amino acid sequences of proteins coded by the GmPUB genes were highly similar to those of PUB proteins in other species. This suggests that their function is highly conserved.

Phylogenetic Analysis of GmPUBs
The amino acid sequences of PUB proteins from different plants were downloaded from NCBI, and the phylogenetic analysis of GmPUBs ( Figure 5) showed that six GmPUB proteins (GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A) showed relatively high similarities with four Arabidopsis PUB proteins (AtPUB18, AtPUB19, AtPUB22 and AtPUB23). Among them, GmPUB16A, GmPUB20A, GmPUB23A and GmPUB24A were clustered in the same large branch, whereas other two GmPUB proteins (GmCHIPA and GmPUB33A) were clustered in the other two branches. On the branches of GmPUB16A, GmPUB20A, GmPUB23A and GmPUB24A, four GmPUB proteins were closely related to AtPUB16, AtPUB20, TcPUB23 and AtPUB24, clustered in the same branch, and evolved in four different directions. In GmCHIPA and GmPUB33A, the former has high homology with AtCHIP, while the latter has high homology with JrPUB33.

Construction of Vectors for Overexpressing and Silencing GmPUBs
Specific primers were used for PCR amplification of the recombinant overexpressing vectors, and 1% agarose gel electrophoresis was performed to check the size of PCR products. The sizes of PCR products were as follows: pOE-GmPUB20A: 1593 bp; pOE-GmCHIPA: 1116 bp; pOE-GmPUB33A: 2733 bp; pOE-GmPUB23A: 1536 bp; and pOE-GmPUB24A: 1563 bp ( Figure S5). This indicated that five plant recombinant vectors overexpressing the GmPUB were successfully constructed.

Gene Expression and Biomass Analysis of Hairy Roots
The positive soybean hairy roots were detected and screened with a handheld lamp 3415RG (Luyor, Shanghai, China) to visualize GFP expression. Hairy roots with strong GFP fluorescence signals were reserved and used for further testing, and the non-GFP hairy roots and the main roots were removed ( Figure S7).
qRT-PCR was used to verify the expression level of the GmPUB gene in the hairy roots overexpressing and RNA interfering with the GmPUB gene. The outcomes displayed that the copy number of the GmPUB gene in the hairy roots overexpressed the GmPUB gene more than the control plant (EV). The pOE-GmPUB20A and pRNAi-GmCHIPA were upregulated by 472.50 and 150.44 times, respectively. Accordingly, the copy number of the GmPUB gene in hairy roots that interfered with the expression of the gene was lower than that in the control plant (EV), and the downregulation multiple was between 0.14-0.67. Among them, the most obvious downregulation was observed for GmPUB20A and GmPUB23A. In addition, compared with transgenic soybean hairy roots harboring empty vectors, the number of lateral roots of hairy roots overexpressing the GmPUB gene increased significantly, while the number of lateral roots of hairy roots RNA interfering with the expression of the GmPUB gene decreased significantly (Figures 6 and 7).  One-way ANOVA was applied in SPSS to analyze the difference of lateral roots in EV and transgenic hairy roots. Multiple t-tests were applied to detect the significant differences of the relative expression level of all GmPUBs in EV and other transgenic hairy roots on Graphpad prism 9. ns represent not significant. * and **** indicate significant difference at the p < 0.05 and p < 0.0001 probability level, respectively. Different characters mean significant differences found at p < 0.05 level. One-way ANOVA was applied in SPSS to analyze the difference of lateral roots in EV and transgenic hairy roots. Multiple t-tests were applied to detect the significant differences of the relative expression level of all GmPUBs in EV and other transgenic hairy roots on Graphpad prism 9. ns represent not significant. * and **** indicate significant difference at the p < 0.05 and p < 0.0001 probability level, respectively. Different characters mean significant differences found at p < 0.05 level. One-way ANOVA was applied in SPSS to analyse the difference the number of lateral roots in EV and transgenic hairy roots. Multiple t-tests were applied to detect the significant differences of the relative expression level of all GmPUBs in EV and other transgenic hairy roots on Graphpad prism 9. ns represent not significant. **** indicate significant difference at the p < 0.0001 probability level. Different characters mean significant differences found at p < 0.05 level.

GmPUBs Regulate Soybean Resistance to SCN in Hairy Roots
After the hairy roots were inoculated with SCN J2, the samples were collected at 15 dpi, and the nematodes in soybean roots were stained by the sodium hypochlorite-acid fuchsin method. The nematodes in each stage were examined and counted under stereomicroscope. As shown in Figure 8, it was found that the total number of nematodes in soybean hairy roots overexpressing GmPUB20A, GmPUB33A and GmPUB24A was higher than EV, while the total number of nematodes in soybean hairy roots overexpressing GmCHIPA and GmPUB23A was lower than that in EV. The results show that overexpression of GmPUB20A, GmPUB33A and GmPUB24A could weaken the resistance of soybean plants under nematode stress, while overexpression of GmCHIPA and GmPUB23A could enhance the resistance of soybean to SCN. Multiple t-tests were applied to detect the significant differences of the relative expression level of all GmPUBs in EV and other transgenic hairy roots on Graphpad prism 9. ns represent not significant. **** indicate significant difference at the p < 0.0001 probability level. Different characters mean significant differences found at p < 0.05 level.

GmPUBs Regulate Soybean Resistance to SCN in Hairy Roots
After the hairy roots were inoculated with SCN J2, the samples were collected at 15 dpi, and the nematodes in soybean roots were stained by the sodium hypochloriteacid fuchsin method. The nematodes in each stage were examined and counted under stereomicroscope. As shown in Figure 8, it was found that the total number of nematodes in soybean hairy roots overexpressing GmPUB20A, GmPUB33A and GmPUB24A was higher than EV, while the total number of nematodes in soybean hairy roots overexpressing GmCHIPA and GmPUB23A was lower than that in EV. The results show that overexpression of GmPUB20A, GmPUB33A and GmPUB24A could weaken the resistance of soybean plants under nematode stress, while overexpression of GmCHIPA and GmPUB23A could enhance the resistance of soybean to SCN.
A further study of the resistance of soybean hairy roots with GmPUB RNA interfered with SCN, as shown in Figure 9, compared with EV; the total number of nematodes in soybean hairy roots after RNA interference with GmPUB genes decreased significantly. Among them, the total number of nematodes in soybean hairy roots after interference to GmCHIPA and GmPUB23A decreased significantly. The above results showed that soybean plants that interfered with the expression of the GmPUB gene had strong resistance.  A further study of the resistance of soybean hairy roots with GmPUB RNA interfered with SCN, as shown in Figure 9, compared with EV; the total number of nematodes in soybean hairy roots after RNA interference with GmPUB genes decreased significantly. Among them, the total number of nematodes in soybean hairy roots after interference to GmCHIPA and GmPUB23A decreased significantly. The above results showed that soybean plants that interfered with the expression of the GmPUB gene had strong resistance. One-way ANOVAs were applied in SPSS to analyze the difference in the number of nematodes in EV and other transgenic hairy roots, and the Holm-Sidak method was used to determine the statistics. Different characters mean significant differences found at p < 0.05 level.

Discussion
Ubiquitin ligase E3 is an important factor that determines substrate specificity in the ubiquitin protein degradation pathway. The plant U-box (PUB) protein family contains the U-box domain, which is a class of ubiquitin ligase E3 with important functions [38][39][40] and plays a wide regulatory role in plant growth and development, biological and abiotic stress, fertility and hormone signaling pathways [4]. There are few reports about soybean GmPUB responses to nematode stress. In the process of soybean and SCN interaction, the relationship between GmPUB resistance to nematode stress and its dependent signaling pathways is unclear.
Recently, 125 U-box genes were identified in the soybean genome [37]. Gene expression analyses can provide critical information about the potential functions of GmPUB genes [41]. Researchers analyzed the expression profile of GmPUB genes under SCN stress using publicly available RNA-seq datasets [42]. Based on the analysis of the RNA-seq datasets, 36 soybean GmPUB genes were found to have significantly changed after incubation with SCN. These results support the notion that GmPUB likely plays a vital role in soybean immunity against SCN [41]. Accordingly, this study was based on publicly available RNA-seq data and previous transcriptome data from our research group, and six soybean GmPUB genes (GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A) were screened from 49 GmPUB genes by qRT-PCR. The differential gene expression of these six genes at different time points after SCN infection indicates that they may be related to the interaction of soybean and SCN. Meanwhile, the public RNA-seq data investigated the response of soybean roots to virulent and avirulent SCN inoculation for 6 and 8 days, while this study examined the response of genes in susceptible and resistant soybeans at four time points (1, 5, 10 and 15 dpi), which expanded our understanding of the functions of ubiquitination-related genes in soybean immunity. Most soybean cultivars resistant to SCN are derived from limited resistance sources, and SCN races have begun evolving to overcome SCN resistance [41,43]. Therefore, studying the function of these six genes in the process is necessary.
In plants, U-box is often coupled with ARM (armadillo repeat) to form a U-box/ARM structure, which interacts independently or jointly with other proteins. The U-box/ARM protein is a special protein in higher plants [44]. There were 41 U-box/ARM proteins in Arabidopsis and 43 U-box/ARM proteins in rice, but there was no U-box/ARM protein in Chlamydomonas vulgaris [45]. Some proteins with U-box/ARM functional domains have been identified in other higher plants. These proteins were the ARC1 protein in rape [46], ACRE276 protein in tobacco [22] and PUB4 protein [47], PHOR1 protein in potato [48], CMPG1 protein in celery (Petroselinum crispum) [49], BG55 protein in mangrove plant Bruguiera gymnorrhiza [50] and SPL11 protein in rice [16]. Some scholars believe that the U-box/ARM protein plays an important role in the plant hormone response, disease resistance and other aspects [51,52]. Among these proteins, U-box proteins such as SPL11 in rice [18,19,53], ACER276 in tobacco and AtPUB17 in Arabidopsis are related to plant disease resistance [22,27]. SPL11 is involved in the primary defense of rice against pathogen infection and negatively regulates cell death [16,19]. Contrary to the negative regulation of SPL11 on disease resistance, CMPG1 in tomato, tobacco and celery can positively regulate plant disease resistance [22]. Meanwhile, Amador et al. reported that PHOR1 (photoperiod-responsive1), a U-box gene found in Solanum tuberosum, is a positive regulator of the GA signal transduction pathway [48]. Kim et al. found that NtPUB4 in tobacco is involved in plant development and cytokinin signal transduction [47]. However, studies on the disease resistance of the U-box/ARM protein in soybean have not been reported. In the present study, 49 soybean GmPUB genes differentially expressed under nematode stress were analyzed by bioinformatics, and 49 U-box proteins were divided into six categories. The U-box/Arm structure is the second type, and 21 of the 49 soybean PUBs belong to this type of protein. The GmPUB16A (Glyma01g131800) gene was highly expressed at different time points after inoculation with SCN. The silencing of GmPUB16A in hairy roots decreased the number of SCN. It indicates that the U-box/Arm structure proteins may play a regulatory role in plant disease resistance. These results provide a theoretical basis and reference for further studying the mechanism of U-box proteins in soybean and cyst nematode interaction.
The CHIP protein located in the cytoplasm, with auxiliary chaperone molecules and E3 activities. Its E3 activity depends on the C-terminal U-box domain [53]. Through its E3 activity and molecular chaperones, such as Hsc70/Hsp70 and Hsp90, CHIP synergistically changes the Hsc70/Hsp70-and Hsp90-mediated regulation of signaling pathways in the protein folding and degradation balance and the quality control of proteins in cells [21,54]. The Arabidopsis AtCHIP has three peptide repeats and one U-box domain. Its expression is induced by various environmental stresses, such as the increase in expression under the stimulation of high temperature, low temperature or strong light and the overlapping repair or degradation of these abnormal proteins by regulating Hsp70 and Hsp90 to reduce the damage of environmental stress on plants [44,55]). Other studies have found that the AtCHIP protein can promote the degradation of abnormal proteins in cells under strong light stress using the subunit ClpP4 of chloroplast protease [56]. The AtCHIP protein is also related to the signal transduction pathway of abscisic acid (ABA) stress, and plants overexpressing the AtCHIP gene are very sensitive to ABA stress [57]. The AtCHIP protein has become a typical example for studying the physiological functions and mechanisms of the AtPUB protein. In this study, GmCHIPA (encoded by Glyma03g088400) was highly homologous to AtCHIP. In addition, the GmCHIPA gene was expressed at different time points after inoculation with SCN J2. A further study of its function showed that overexpression of GmCHIPA reduced the number of nematodes. Silencing this gene also causes a decrease in the number of nematodes. However, this research has not yet clarified the specific mechanism by which GmCHIPA regulates soybean cyst nematode resistance. We hypothesized that the overexpression and silencing of GmCHIPA are involved in different regulatory pathways between soybean and cyst nematode interaction. Therefore, understanding the specific substrates of GmCHIPA and the upstream and downstream signaling pathways involved needs to be further investigated.
As positive or negative regulators of the plant immune response to various pathogens, U-box proteins play an important role in the plant defense response and programmed cell death [39]. Many types of U-box proteins play a negative regulatory role in plant disease resistance. The U-box E3 gene Spl11 was isolated from the indica rice variety IR68, and its mutant spl11 could increase the nonspecific resistance of rice to Magnaporthe grisea and Xanthomonas oryzae pv. Oryzae [16]. Similarly, Liu et al. reported that the U-box E3 ligase SPL11 and its Arabidopsis homologous gene PUB13 in rice played a negative regulatory role in PCD and defense responses [19]. Arabidopsis U-box E3 ubiquitin ligases AtPUB22, AtPUB23 and AtPUB24 negatively regulate the PTI response induced by PAMPs such as elf18 and chitin. Studies have shown that the PTI response was significantly enhanced in the three mutants of pub22/pub23/pub24, and the three mutants significantly improved the resistance to live vegetative pathogens. At the same time, Pseudomonas syringae pv. tomato inoculated with pub22/pub23/pub24 mutants showed enhanced resistance to pathogens, increased ROS content and upregulated expression of RbohD genes involved in ROS formation [13,14,58]. In this paper, GmPUB23A (encoded by Glyma02g195900) and GmPUB24A (encoded by Glyma12g188900) are homologous to AtPUB23 and AtPUB24. Functional studies have found that the overexpression of GmPUB24A weakened the resistance of soybean under nematode stress, and the overexpression of GmPUB23A enhances the resistance of soybean to SCN. Furthermore, this study interfered with the GmPUB24A and GmPUB23A in soybean hairy roots. The results show that soybean roots interfered with the expression of the GmPUB24A and GmPUB23A genes, conferring strong resistance to SCN. Similarly, in GmCHIPA, soybean's resistance to SCN was enhanced in both GmPUB23A overexpression and RNA interference hair roots. We speculate that the reason for this situation is that plants under two treatments may also be involved in different regulatory pathways of resistance to nematodes. In addition, Wang et al. showed that the soybean GmPUB8 gene, which was highly homologous to GmPUB23A, played an important role in abiotic stress. Its expression was induced by exogenous ABA and NaCl, regulating flowering time through the photoperiod pathway and playing a negative regulatory role under plant drought stress [37]. Similarly, the tolerance of Arabidopsis pub22 and pub23 mutants to salt stress was enhanced, while pub22pub23 double mutants showed an additive effect [59].
Unlike the above U-box proteins that play a negative regulatory role, some plant U-box proteins play a positive regulatory role in plant disease resistance. AtPUB17, AtPUB20 (AtCMPG1) and MAC are positive regulatory factors in plant disease resistance, and these three proteins are involved in the regulation of ETI in the innate immune response. Among them, the expression of AtPUB20 was induced by PAMP and had autoubiquitination activity in vitro. The substrate of AtPUB20 is the G-protein β subunit AGB1 [23,60]. To regulate host immunity, soybean GmPUB1-1 interacts with many other RxLR effectors, including Avr1d [61]. Lin et al. confirmed that GmPUB13 is involved in the virulence mechanism of Phytophthora effectors. The virulence mechanism of GmPUB13 is that Avr1d inhibits the activity of the GmPUB13 E3 ligase by competing with E2 to stabilize the sensitive factor GmPUB13 and promote Phytophthora infection [11]. In this study, GmPUB20A (encoded by Glyma07g106000) was homologous to AtPUB20 (AtCMPG1), and its function was further studied. The overexpression of GmPUB20A under nematode stress weakened the resistance of plants. In contrast, interfering with the expression of GmPUB20A more significantly enhanced the resistance of soybean to cyst nematodes. These findings highlight the involvement of U-box proteins in plant immune responses. In addition, the transmission of defense signals and the participation of UPS (ubiquitin-26S proteasome system) components in the process from pathogen perception to pathogen immune response and effector immune response also proved the multiple roles of UPS in biological stress response.
Although many ubiquitin ligases have been identified, their function is still unclear. The mining and identifying of more upstream regulatory factors and downstream substrates of ubiquitin ligases are essential for understanding their functions. Currently, studies on the function of U-box genes mainly use yeast two-hybrid screening of downstream proteins interacting with them [62] and reverse genetics based on RNA silencing or gene knockout [7]. In addition, the same U-box ubiquitin ligase may involve various biological processes. For example, the U-box ubiquitin ligases AtPUB20 and AtPUB23 are involved in plant disease resistance, growth and development. The AtCHIP protein is related to both temperature stress and abscisic acid stress. This action mode of ubiquitin ligase also provides a new basis for studying some cross-signaling pathways. Although research on PUB is a great challenge for plant science, it also helps us to understand the proteinprotein interaction in plants and many physiological activities in cells, especially the function of such proteins in plant disease resistance and stress resistance. GmPUB20A, GmPUB33A and GmPUB24A played a vital role in soybean resistance to cyst nematodes in the present study, and several problems need to be solved immediately, such as finding the substrates of GmPUB genes and clarifying the process of ubiquitin proteins' recognizing substrate proteins. The results of these future studies will be more helpful in explaining the mechanism of plants against pathogen infection through immune defense and provide new ideas for the integrated management of nematode disease.

Nematode Population, Plant Materials and Growth Condition
The race 3 (HG 0) population of Heterodera glycines (soybean cyst nematode, SCN) was propagated on the susceptible soybean cultivar Liaodou 15. Soybean (Glycine max), Williams 82 (W82) and "Huipizhi" (HPZ) were selected as susceptible and resistant cultivars to SCN race 3, respectively. Soybean seeds were surface sterilized with 1% NaClO for 10 min and were washed at least five times with sterilized water. Then, the seeds were germinated in plastic pots (length × width × height = 22 × 15 × 7 cm) containing wetting vermiculite with Hoagland nutrient solution. Four days later, uniform seedlings were selected and transferred into PVC tubes (height × diameter = 10 × 3 cm) containing equal ratios of sterilized sand and soil in a climatic chamber (16 h light/8 h dark, temperature 23-26 • C, humidity 50%). The plants were watered once every three days with Hoagland nutrient solution.

Bioinformatics Analysis of GmPUBs
The accession number and coding sequence (CDS) of soybean GmPUB genes were obtained from SoyBase (https://www.soybase.org, accessed on 1 October 2020), and the molecular weight and isoelectric point (PI) of the gene product were analyzed online by Protparam (http://web.expasy.org/protparam/, accessed on 4 April 2021). The protein sequences of the GmPUB genes were downloaded from the NCBI website. DNAMAN software was used for multisequence alignment, and the NJ algorithm of MEGA 7.0 was used to generate the phylogenetic tree.

Expression Analysis of GmPUBs
Quantitative real-time PCR (qRT-PCR) was performed to detect the relative expression of the 49 GmPUB genes of the susceptible cultivar William 82, and the resistant variety Huipizhi after soybean cyst nematode race 3 infection in this study (approximately 2000 second-stage juveniles (J2) were inoculated at the root of each soybean seedling, and uninoculated plants were used as controls). SCN-infected soybean roots were collected at 1 dpi, 5 dpi, 10 dpi and 15 dpi. RNA was extracted from soybean roots using the Reagent RNA extraction kit (CWbio, Beijing, China) according to the manufacturer's instructions. Then, RNA was used for first-strand cDNA synthesis by using the PrimerScript TM RT Reagent Kit with gDNA Eraser (TaKaRa Bio Inc., Dalian, China). qRT-PCR was performed by using a TaKaRa SYBR Green PCR Master Mix kit, and the qRT-PCR system and procedures were performed according to the manufacturer's instructions. All qRT-PCR experiments were performed on a CFX Connect Real-Time Quantitative PCR system (Bio-Rad, Hercules, CA, USA) with the manufacturer's instructions. The soybean ubiquitin 3 gene GmUBI-3 (GenBank accession D28123.1) was used as an internal reference gene, with three parallel and three biological repeats per sample. The results were analyzed by the 2 -∆∆CT method, and all the specific primers used for qRT-PCR are shown in Supplementary Materials S9.

Cloning of 6 GmPUB Genes
cDNA of soybean Williams 82 was used for gene cloning. First, RNA was extracted from soybean roots by using the Reagent RNA extraction kit (CWbio, Beijing, China) according to the manufacturer's instructions. Then, RNA was used for first-strand cDNA synthesis (unified concentration of 1000 ng µL −1 ) using PrimerScriptTM RT Reagent Kit with gDNA Eraser (TaKaRa Bio Inc., Dalian, China). The full-length primers CDS-F and CDS-R were designed according to the 6 GmPUBs (Supplementary Materials S9). The fulllength sequences of these GmPUBs were obtained by amplifying leaf cDNA as a template.

Construction of Vectors for Overexpressing and Silencing GmPUBs
The plant expression vector pNINC2RNAi with a GFP tag was constructed using pCAMBIA3301 as the backbone [63]. Primers were designed according to the cDNA sequence of soybean GmPUBs and the sequence of AscI and AvrII restriction sites in pNINC2RNAi and homologous recombination requirements (Supplementary Materials S9). Then, the full-length CDS of 5 GmPUBs was cloned into pNINC2RNAi via PCR technology.
For GmPUB16A, GmPUB20A, GmCHIPA, GmPUB33A, GmPUB23A and GmPUB24A RNAi, we synthesized the RNAi sequences (including sense sequence, loop and antisense sequence) in BioRun Biosciences Co., Ltd. (BioRun, Wuhan, China). Then, the synthetic fragment was inserted into the vector pNINC2RNAi by restriction digest cloning to obtain the 6-plant interference recombinant vector pNINC2RNAi-RNAiGmPUBs of the GmPUB gene.
After screening positive clones, the freeze-thaw method transformed the recombinant vectors' overexpressing and silencing GmPUBs into Agrobacterium rhizogenes K599.

Transformation Process of Soybean Hairy Roots
The seeds of soybean Williams 82 were disinfected and germinated in sterilized vermiculite watered regularly with Hoagland nutrient solution. Six days later, young seedlings with unexpanded cotyledons were used for infection with A. rhizogenes K599.
The soybean cotyledonary nodes were inoculated with K599 strains carrying cognate plasmids, pNINC2RNAi-OEGmPUBs, pNINC2RNAi (empty vector, EV) and pNINC2RNAi-RNAiGmPUBs, for overexpression or silencing by syringe [64]. After hairy roots were grown, the following parts of the hairy roots were subtracted ( Figure S1). The plants with hairy roots were transferred into a new plastic culture pot containing sterile vermiculite for further culture, and nutrient solution and sterile water were regularly irrigated. Approximately ten days later, the hairy roots were screened with a handheld lamp (Luyor, Shanghai, China) to visualize GFP expression. Hairy roots carrying strong GFP signals were reserved and transferred into plastic culture pots for further culture.

Gene Expression and Biomass Analysis of Hairy Roots
The samples were collected from hairy roots cultured for 15 d, and the effects of overexpression and silencing were confirmed by detecting the expression level of the GmPUB gene by qRT-PCR (the primers used for qRT-PCR are shown in the Supplementary Materials S9). The soybean ubiquitin 3 gene GmUBI-3 (GenBank accession D28123.1) was used as an internal reference gene, with three parallel and three biological repeats per sample, and the results were analyzed by the 2 -∆∆CT method. At the same time, the biomass of hairy roots cultured for 15 d was indicated by the number of lateral roots of soybean-transformed hairy roots. Each treatment was repeated five times.

Investigation of SCN Development in GmPUB-Overexpressing and -Silenced Soybean Hairy Roots
When the hairy root length was approximately 5-10 cm, it was transferred into PVC tubes (height × diameter = 10 × 3 cm) containing equal ratios of sterilized sand and soil in a climatic chamber (16 h light/8 h dark, temperature 23-26 • C, humidity 50%). The plants were watered once every three days with Hoagland nutrient solution. After the hairy roots were restored to health, the second-stage juveniles (J2) of SCN were inoculated into the hairy roots induced by pNINC2RNAi-OEGmPUBs, pNINC2RNAi and pNINC2RNAi-RNAiGmPUBs. Each plant was inoculated with 500 J2s, and hairy roots without nematodes were used as a control. The nematodes in soybean roots were stained by the sodium hypochlorite-acid fuchsin method at 15 dpi with SCN. The nematodes at different ages were examined and counted under an anatomical microscope and microscope. Each treatment was repeated five times.

Statistical Analyses
Statistical analyses were conducted using IBM SPSS STATISTIC v.22 (Armonk, NY, USA) and GraphPad Prism 9 software (GraphPad Inc., San Diego, CA, USA). Independent t-tests were applied in SPSS to analyze the difference in the number of nematodes and lateral roots. Multiple t-tests were applied to detect the significant differences in the relative expression levels of all GmPUBs in the pairs HC vs. HS, WC vs. WS, EV vs. OE and EV vs. RNAi on GraphPad Prism 9, and the Holm-Sidak method was used to determine the statistical significance with alpha = 0.05.

Conclusions
This study's results clearly indicate the characteristics of 49 GmPUB genes in soybean and their expression patterns under SCN stress. Six soybean GmPUB genes with significant responses to SCN infection were obtained. Additional functional studies of the six genes showed that GmPUB20A, GmPUB33A and GmPUB24A negatively regulated soybean resistance. These data help to further reveal the regulatory mechanism of soybean GmPUB participating in soybean-SCN resistance.