The RNA-Seq data of 36 soybean root samples treated with HG Type 0 stress were obtained from our previous study (Jiang et al., 2020; Jiang et al., 2021). ‘Heinong 37’ (developed by Heilongjiang Agricultural Academy, susceptible to SCN HG type 0), ‘Dongnong L-10’, and ‘Dongnong L-204’ (provided by Northeast Agricultural University, resistant to SCN HG type 0) were planted in pots with a diameter of 13 cm in a greenhouse (day/night temperature: 27–28°C; relative humidity: 60–70%; day/night illumination time 16h/8h), collected and inoculated with SCN using the method described by Wan et al. (2015). A random complete block design (RCBD) was used with 3 replicates of 10 seedlings each. Root samples inoculated with 2000 SCN (J2) were collected and sequenced at 0, 3, 7, and 10 days after inoculation. The RNA-Seq data of these 36 samples included transcriptome data of 28,656 genes, from which 10,000 genes with an average transcriptome FPKM value greater than 2.5 were selected for co-expression network analysis.

The good “Genes Samples” function in the R packet of WGCNA was used to eliminate the genes and samples with too many missing values, and the sample-level clustering-pruning graph method was used to eliminate outlier samples—that is, to remove outlier samples that were significantly higher than other samples (Giulietti et al., 2016). The hclust function in R was used to construct the clustering tree. The Pick-Soft threshold function was used to analyze the network topology information, and the appropriate β value was selected as the soft threshold to construct the network to make the network meet the scale-free topology characteristics. The topological difference value of the network was calculated, and the genes were divided using a hierarchical clustering method to generate the gene clustering tree. The gene module with a high degree of correlation was identified by the one-step construction method; 50 was set as the minimum number of genes in the module, and the threshold of cutting height was 0.25 ( (Giulietti et al., 2016)). The gene module was identified by the cutting tree method, and modules with high similarity were combined. The gene groups were named with different colors, which is convenient for distinguishing different gene modules in subsequent gene functional recognition and visual analysis. The WGCNA package (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/) in R was used to perform the analysis, as described previously (Mason et al., 2009). The genes of MM > 0.85 and GS > 0.6 in the key modules were introduced into Cytoscape using the MMC (Maximum Margin Criterion) algorithm, and the core genes were selected according to their weight and gene network regulation position.

The soybean variety ‘Dongnong 50’ (DN50, SCN-susceptible cultivars) was used as the WT and the background plant for genetic transformation. ‘Dongnong L-204’ and ‘Heinong 37’ were used as resistant and susceptible cultivars, respectively, for gene expression pattern analysis. The seeds of soybean were germinated in vermiculite and peat soil (1:1). Plants were grown under long-day conditions (16 h light/8 h dark) at 25 ± 1°C for routine maintenance. HG Type 0 was obtained from the Soybean Research Institute of Northeast Agricultural University and was isolated and purified for many generations. A randomized complete block design (RCBD) was utilized, with 3 replicates and 10 seedlings per replicate. Each seedling was inoculated with approximately 2000 second-stage juvenile nematodes (J2s). Accordingly, mock-inoculation with distilled water was also conducted for each line as a control. All treatments and controls were watered daily to maintain soil moisture and promote uniform infection throughout the root system. Both SCN-inoculated and mock-inoculated root samples were harvested at 0, 5, 10, and 15 day post inoculation (dpi). A randomized complete block design (RCBD) was utilized with three replicates and ten seedlings per replicate. The preparation of the SCN egg suspension and its identification by acid fuchsin staining were performed following Jiang et al. (2021).

The RNA extraction methods were described previously (Jiang et al., 2021). The extracted complete RNA was reverse transcribed into cDNA. DEGs were selected and verified by qRT-PCR based on SYBR Green PCR Master Mix (TIANGEN BIOTECH, BeiJing, China) and the 7500 Fast PCR detection system. The relative mRNA level of each candidate gene was evaluated against soybean GmACTIN (GenBank Accession Number AF049106) as a reference gene. Three technical replicates were performed per gene, and the relative levels of transcript abundance were calculated using the 2^-ΔΔCT method (Jain et al., 2016). The sequences of the primer pairs were used to amplify the candidate genes ( Table S4 ). The promoter elements of the 1000 bp pro-GmHg1 sequence were analyzed using PlantCARE software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

We identified the Coding sequence (CDS) of the soybean GmHg1 gene in the plant database (https://phytozome.jgi.doe.gov/pz/portal.html) and amplified the full-length CDS sequence of GmHg1 from the developing seeds of ‘Dongnong L-204’ by RT-PCR. For construction of the recombinant vectors GmHg1-OX and GmHg1-GFP, 2130-bp CDS of GmHg1 were amplified from the cDNA of ‘DN50’ by overlapping PCR. Then, the PCR products were ligated to the linear vectors pCAMBIA3300 and pCAMBIA1302 using a homologous recombination system (ClonExpress II One Step Cloning Kit, Vazyme, China). For the construction of GmHg1-RNAi, the specific sequences of the 498-bp sense and antisense fragments of the GmHg1 gene were cloned into the RNA interference plant expression vector pFGC5941 through homologous recombination. According to this method (Zhao et al., 2018), the recombinant plasmid GmHg1-GFP was transformed into Arabidopsis protoplasts and analyzed under a fluorescence microscope (Leica, Germany). The recombinant plasmids GmHg1-OX and GmHg1-RNAi were transformed into Agrobacterium rhizogenes EHA105 according to a previously described method (Zhao et al., 2018), and transgenic soybean was obtained.

We used the roots of SCN-resistant cultivars (‘Dongnong L-10’, ‘Dongnong L-204’) and SCN-susceptible cultivars (‘Heinong 37’) under SCN HG type 0 stress for transcriptome sequencing, and 36 transcriptome samples were obtained. The transcriptome data of 36 samples were analyzed using WGCNA, and 10,000 genes with FPKM values greater than 2.5 were selected for soft threshold clustering. According to the expression trend, 10,000 genes were divided into 12 modules. Cluster analysis showed that the three biological repeat sequences of the samples were close to each other, and the overall effect was good ( Figure 1 ).

To ensure that it conforms to the scale-free network distribution, the WGCNA needs to choose the appropriate weighting coefficient β. The soft threshold function of Pick in the WGCNA software package was used to calculate β. Topology analysis showed that, when the threshold was β = 18, the scale-free topology fitting index (R2) was close to 0.85. This indicates that the network was close to a scale-free network ( Figure 2 ). Thus, we used β = 18 as the soft threshold for constructing a co-expression network. Using the WGCNA package in R, the genes with similar expression patterns were divided into modules; 50 was set as the minimum number of genes in the module, and the threshold of cutting height was 0.25. The gene modules were identified using the dynamic cut tree method, and the modules with high similarity were combined to obtain 21 modules ( Figure 3 ). Because of the large number of genes, the displayed network was difficult to identify in detail. The association diagram between modules was drawn, and the proximity of different color modules and the correlation of genes between modules were evaluated ( Figure S1 ). At the same time, a visual gene network was also drawn to show the degree of association of specific genes ( Figure S2 ). After generating the characteristic gene map for each module, the correlation between the characteristic genes and SCN stress was analyzed. The module most related to SCN stress was the turquoise module, with 632 DEGs, the correlation was 0.86, and the p-value was less than 0.001 ( Figure 4 ).

In the gene co-expression network, highly connected genes are called central genes, and central genes play an important role in response to stress (Koido et al., 2018). We selected genes with an MM > 0.85 and GS > 0.6 in module turquoise ( Figure S3 ), introduced them into Cytoscape, and screened them using the MMC algorithm, and selected the gene with the highest connectivity as the hub gene. The central gene GmHg1 was found in module turquoise, and the co-expression network of the hub gene GmHg1 was constructed ( Figure 5 ). The hub gene GmHg1 (Glyma.06G035000) encodes serine/threonine protein kinase. At the same time, 11 genes that were closest to the hub gene expression trend were found ( Table 1 ). They included Glyma.09G027700 (wall-associated receptor kinase-like 1), Glyma.06G109900 (C₂H₂-like zinc finger protein), Glyma.07G009700 (Hydroxyproline-rich glycoprotein family), Glyma.11G148500 (Arabinogalactan-protein 11), and Glyma.10G035400 (Basic helix-loop-helix (bHLH) DNA-binding superfamily protein). The domains of these genes have been identified as being involved in plant disease resistance or abiotic stress responses (Simon et al., 2005; Lamport et al., 2006; Meier et al., 2010; Wang et al., 2019; Hao et al., 2021). Genes with the same expression trend may have the same function. It is speculated that these genes may be involved in soybean resistance to SCN.

To understand the function of GmHg1 in soybean, we isolated the CDS fragment of GmHg1 from ‘Dongnong L-204’. The fragments comprised 2130-bp open reading frames and were predicted to encode 710 residue polypeptides. To verify the basic characteristics of GmHg1, a fusion expression vector of 35S::GmHg1-GFP was constructed and transiently expressed in Arabidopsis protoplasts. The fusion protein was expressed in the nucleus, cell membrane, and cytoplasm of protoplasts. ( Figure 6 ). At the same time, promoter element analysis of the 1000-bp upstream GmHg1 sequence was carried out ( Table S1 ). There were 16 elements in the GmHg1 promoter sequence, among which TGACG-motif and CGTCA-motif were cis-acting elements of methyl jasmonate, which responded to abiotic stress. ABRE, ARE, and TGACG-motif elements were corresponding elements for defense and response to stress. These results suggest that GmHg1 may be involved in resistance to stress.

To clarify the expression pattern of GmHg1 under SCN stress and verify the accuracy of transcriptome data, we treated the roots of the SCN-resistant cultivar (‘Dongnong L-204’) and the SCN-susceptible cultivar (‘Heinong 37’) with HG Type 0. The expression of GmHg1 was determined at 0, 5, 10, and 15dpi. Before HG Type 0 induction, GmHg1 expression in resistant varieties was higher than in susceptible varieties. After HG Type 0 induction, GmHg1 expression in SCN-resistant/-susceptible cultivars was significantly upregulated compared with that before induction ( Figure 7 ). The expression of GmHg1 reached the highest level after 10 days of HG Type 0 stress. At the same time, the upregulation multiple in SCN-resistant cultivars was significantly higher than that in SCN-susceptible cultivars. Therefore, GmHg1 may be involved in the resistance response of soybean to SCN, and the difference in GmHg1 expression may lead to the difference in SCN resistance among SCN-resistant and -susceptible cultivars.

To verify whether GmHg1 was involved in the response to HG Type 0 stress, we used ‘Dongnong 50’ as a WT and transformed the 35S::GmHg1 recombinant vector into soybean by the semi-seed transformation method (Zhao et al., 2018). We obtained T2 generation 35S::GmHg1 overexpressing soybean and identified GmHg1-OX soybean using PCR, qRT-PCR, and bar test strips ( Figure 8 ). qRT-PCR showed that the GmHg1 expression levels were significantly increased in T2 generation plants. At the same time, we constructed a pFGC5941-GmHg1-RNAi vector, interfered with the expression of the GmHg1 gene by RNAi, and created pure and stable soybean plants transformed with the GmHg1-RNAi vector in ‘Dongnong 50’. qRT-PCR showed that the expression levels of GmHg1 were significantly decreased in T2-generation plants ( Figure S4 ). The roots of transgenic soybean with GmHg1-OX and GmHg1-RNAi genes were further stressed by HG Type 0, and the phenotype was observed by acid fuchsin staining 10 days after inoculation ( Tables S2 , S3 ). Compared with ‘Dongnong 50’, the number of HG Type 0 in transgenic soybean plants with the GmHg1-OX gene decreased significantly from 5.34 to 2.9 per root. However, the number of cyst in transgenic soybean plants with the GmHg1-RNAi gene increased significantly from 5.34 to 7.74 per root. At the same time, cyst formation in soybean plants overexpressing the GmHg1 gene was severely inhibited compared with that of soybean plants that interfered with GmHg1 gene expression. The increase in GmHg1 gene expression may improve the resistance of soybean to HG Type 0. The GmHg1 gene plays a positive regulatory role in the process of soybean resistance to HG Type 0.

To identify the expression pattern of genes similar to the hub gene expression trend in transgenic materials, we detected the expression of 11 genes in soybean seedlings transformed with GmHg1-OX and GmHg1-RNAi genes. The expression of 11 genes in GmHg1-OX transgenic soybean seedlings was significantly higher than that in ‘Dongnong 50’ ( Figure 9 ). Among them, the expression of Glyma.09G027700 (wall-associated receptor kinase-like 1) increased the most, and that of Glyma.08G161700 (Seed storage 2S albumin superfamily protein) and Glyma.03G036600 (Ovate family protein 2) genes increased the least.

The expression of 11 genes in GmHg1-OX transgenic soybean seedlings was significantly lower than that in ‘Dongnong 50’ ( Figure S5 ). This indicates that these 11 genes were co-expressed with GmHg1 and may be involved in the process of soybean resistance to SCN. However, the specific resistance mechanisms of these genes and GmHg1 need to be further studied.