The soybean breeding line 09-138 was developed by the Heilongjiang Academy of Agricultural Sciences (Hua et al., 2018). SCN race 4 (SCN4, HG type 1.2.3.5.6.7) and SCN race 5 (SCN5, HG type 2.5.7) were originally collected from the field and cultured from single cysts for more than 5 generations on the susceptible soybean variety Dongsheng1 in a greenhouse with 16 h of light and 8 h of darkness at 23–28°C, and then identified by race test and HG type indicator lines (Hua et al., 2018). Every other year, nematode HG types/races were reconfirmed without virulence change with time.

A nematode inoculum was prepared according to the method described by Huang et al. (2022). Plant root tissue and soil were collected 35-40 days after inoculation and put into a 2 L-beaker. The mixture was stirred vigorously with a glass rod for approximately 1 min and then precipitated for 10 s. The fluid of the supernatant was gently poured into 75/25 μm nested sieves. The mixture of cysts and root debris on the top of the sieves was rubbed with a rubber stopper to release eggs. Eggs on 25-μm sieve were rinsed with a high-pressure water faucet for 1 min and then with sterile water before collection. The collected eggs were then transferred onto six to eight layers of tissue paper supported by a metal screen on a hatching dish containing 3 mM ZnSO4 in sterile water for hatch at 28°C. J2s were then collected for inoculation after 3-4 days.

Seeds of 09-138 were sterilized by soaking in 0.5% sodium hypochlorite for 20 min and rinsed with sterile water three times. Two seeds were sown in a black plastic pot (8 cm diameter × 12 cm depth) filled with autoclaved soil and sand at a ratio of 1:1. After 4 days, two seedlings were thinned to one in each pot. An eight-day seedling was inoculated with a 1-ml suspension containing 2,000 J2 of SCN4 or SCN5. Seedlings were inoculated with 1 ml water as control. The plants were maintained in a growth chamber at a 16-/8-h day/night regime, 28°C day/22°C night, and 50% relative humidity.

Collected roots were stained 3, 6, 8, 10, and 12 days after inoculation with acid fuchsin (Byrd et al., 1983). Nematode development inside the roots was observed, and roots and nematodes were photographed under an Olympus SZX16 dissecting microscope using the Cellsens Standard image software (Olympus Corporation, Japan).

To collect roots for RNA extraction, plant roots 8 days after inoculation were washed and rinsed thoroughly with water. Three roots from each treatment were wrapped together with aluminum foil as one replication (one sample). Three replications were made for each treatment. A total of 9 samples with SCN4- and SCN5-infection and control were collected. Immediately, each prepared sample was put in liquid nitrogen to freeze it and was kept at -80°C for RNA extraction and sequencing.

Total RNA was extracted using RNAeasy Plant Mini Kit (Qiagen, United States), and RNase-free DNase (Qiagen) was used to remove DNA contamination in the total RNA. The concentration, purity, and integrity of the extracted RNA were measured with 1% agar gel (Thermo Fisher Scientific, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, United States). cDNA library construction started with 1 μg total RNA using a cDNA-PCR sequencing kit (SQK-PCS109) provided by Oxford Nanopore Technologies (ONT, Inc., United Kingdom) following the manufacturer’s instructions. Final cDNA libraries were added to FLO-MIN109 flow-cells and run on the PromethION platform at Biomarker Technology Company (Beijing, China) for sequencing. Experimental processes including sample quality testing, library building, library-quality testing, and library sequencing were performed in accordance with standard procedures provided by ONT.

Raw reads were subject to filtering with an average read quality score ≤ 7, and read length ≤ 500 bp, and ribosomal RNA mapped to the rRNA database was also discarded. Primer sequences on both ends of clean reads were searched to determine full-length, non-chimeric (FLNC) sequences. Detected FLNC transcripts were then mapped to soybean reference genome Williams 82.a2.v1¹ with minimap2 (Li, 2018) to obtain FLNC clusters, and pinfish² was applied to polish each cluster to attain consensus isoforms. All mapped reads were further collapsed using the cDNA_Cupcake package with a mini-identity of 90% and mini-coverage of 85%. When the redundant transcripts were collapsed, 5′ difference was not considered. The achieved transcripts were compared with known transcripts of reference genome Williams 82.a2.v1 utilizing gffcompare, and novel transcripts were identified in order to make supplementary genome annotation. Gene boundaries were modified, and transcripts with expressed levels ≤ 1 were filtered.

Alternative polyadenylation was identified through further analysis of FLNC by Transcriptome Analysis Pipeline from Isoform Sequencing (TAPIS) (Foissac and Sammeth, 2007). Multiple Expectation Maximization for Motif Elicitation (MEME) (Bailey et al., 2006) was used to analyze the 50-bp sequence upstream of the poly A site to detect FL motifs. The consensus sequence before remove-redundant analysis was used for fusion transcript analysis. Fusion transcript was defined under the following conditions: aligned to 2 or more sites; each site covers at least 5% of a transcript with a 1-bp minimum alignment length; the total length covers more than 95% of the total length of transcripts with at least 10k bp distance between the two sites.

Alternative splicing indicates the process of pre-mRNA treatment. Gene transcription generates pre-mRNAs with many splicing methods. Five types of alternative splicing events (3’ splice site, 5’ splice site, exon skipping, intron retention, and mutually exclusive exon) of transcripts were examined by employing the Astalavista software based on the alignment results of individual samples to the reference genome (Foissac and Sammeth, 2007). The MIcroSAtellite (MISA, an identification tool) software was used for SSR analysis, and transcripts below 500 bp were discarded.

Coding sequences (CDSs) were predicted with TransDecoder (v3.0.0; Haas et al., 2013) based on the ORF. LncRNA does not code for protein. Therefore, LncRNA in novel transcripts was predicted whether it had a coding potential by protein domain analysis including all four methods, Coding Potential Calculator (CPC) (Kong et al., 2007), Coding-Non-Coding Index (CNCI) (Sun et al., 2013), Coding Potential Assessment Tool (CPAT) (Wang et al., 2013) and Protein family (Pfam) (Finn et al., 2014). LncRNA target genes were predicted using two methods: first, depending on the location relationship between the differentially expressed lncRNA and adjacent mRNA (within 100k bp distance) expressed differentially; second, according to complementary base pairing between lncRNA and mRNA using the lncTAR tool (Li et al., 2015). TFs were detected with iTAK (Zheng et al., 2016).

Mapped full-length reads with > 5 match quality were chosen for quantification. Transcript or gene expression levels were measured in counts per million (CPM) (Zhou et al., 2014) and calculated by the following:

CPM = (read number matched the transcript)/(total read number matched referenced transcriptome) × 10⁶

The differential expression among the treatments was analyzed with DESeq2 (Anders and Huber, 2010) depending on a negative binary distribution model, consequently gaining DEGs or DETs. False discovery rate (FDR) was adjusted and controlled with the method of Benjamini and Hochberg (1995), and DEGs or DETs with log₂fold change (FC) ≥ 2 and FDR < 0.01 were chosen. A heat map for DEGs in each group was developed using the pheatmap package in R (Version 1.0.12³).

Functional annotation of the genes/transcripts was conducted by blasting with databases including NR (NCBI non-redundant protein sequences) (Deng et al., 2006), Swissprot (Apweiler et al., 2004), GO (Ashburner et al., 2000), Clusters of Orthologous Groups (COG) (Tatusov et al., 2000), euKaryotic Ortholog Groups (KOG) (Koonin et al., 2004), Pfam (Kanehisa et al., 2004), and KEGG (Mckenna et al., 2010).

Gene ontology enrichment analysis of DEGs or DETs was conducted using the GOseq R package-based Wallenius non-central hypergeometric distribution (Young et al., 2010). KEGG pathway enrichment analysis of DEGs or DETs was subject to the KEGG Orthology Based Annotation System (KOBAS) software (Mao et al., 2005). Protein–protein interactions (PPIs) for all detected DEGs were predicted using the STRING database⁴ and were visualized in Cytoscape (Shannon et al., 2003).

Samples leftover from full-length transcriptome sequencing was subject to qRT-PCR for verification of DEGs. Primers were designed by using the Primer Premier 5 software (Lalitha, 2000) and synthesized by Comate Bioscience Company Limited (Changchun, China). A total of 1 μg of treated RNA was used to synthesize the first-strand cDNA using FastKing gDNA Dispelling RT SuperMix (TIANGEN, China). PCR reactions were conducted in LightCycler^® 480 System (Roche Life Science, United States) with ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) following the manufacturer’s protocol. PCR was carried out in a 20-μl volume containing 100 ng cDNA (2 μl). PCR was performed as follows: initial denaturing for 10 min at 95°C, followed by 40 two-step cycles of 95°C for 10 s and then at 60°C for 1 min. The relative expression of tested genes was calculated with the –ΔΔCt method using ACTIN as a control. Three independent biological replicates and three technical repetitions were performed for all the experiments. The primers used are listed in Supplementary Table 1. The correlation between transcriptome sequencing and qRT-PCR was completed in Excel 2016, and independent-samples t-test was performed using SPSS 17.0.

Various development stages of the two races, SCN4 (Figure 1A) and SCN5 (Figure 1B), were observed in roots. Soybean roots were penetrated by both nematode races, and there was no obvious difference in nematode size on day 3 days (Figure 1). On day 6, most of SCN5 developed to the J3 stage but SCN4 remained in the J2 stage. Nematodes had developed from J3 to late J4 stages in roots infected with SCN5 at 8, 10, and 12 days (Figure 1B) when compared with J2, J3, or a few early J4 stages in roots infected with SCN4 (Figure 1A), confirming 09-138 is resistant to SCN4 (cyst number/plant, 13 ± SE 2.7; female index, FI = 10) and susceptible to SCN5 (cyst number/plant, 119 ± SE 8.16; FI = 40) (Huang et al., 2022). More brown spots (hypersensitive response) around some nematode feeding sites were observed in resistant roots with SCN4 than in those with SCN5 (Figure 1A).

Since an obvious difference in nematode development on day 8 was observed in 09-138, root samples from day 8 were used for full-length sequencing. An average of 6.1 Gbp of clean data for nine cDNA libraries was obtained with a range from 5.73 to 6.82 Gbp (Supplementary Table 2). The average N50 length was 1,287 bp with a mean length of 1,142 bp, an average max length of 11,276 bp (9,575-13,426 bp), and a mean quality value of 11 (Q11) (Supplementary Table 2). After rRNA was filtered out, an average of 5,296,377 clean reads (4,681,541-5,888,283) and an average of 4,258,467 FLNC read numbers (3,767,343-4,641,061) were generated with an average of 80.4% FLNC ratio (Table 1).

Finally, 65,038 redundant-removed transcript sequences containing 92,079,818 bp with an N50 length of 1,665 bp, mean length of 1,415 bp, and max-length of 7,285 bp were obtained through a merging of the consistent sequences (Supplementary Figure 1). After alignment with the reference genome and filtering with known annotations of the reference genome, 1,117 novel genes and 41,096 novel transcripts were identified. The consistent sequence of each sample was used for AS analysis.

Full-length sequencing can accurately identify the structure of transcripts. APA analysis based on FLNC displayed a total of 214,760 transcripts with various numbers of poly A sites. Among the three treatments, the greatest portion of 23.8% transcripts contained > 5 poly A sites, followed by 1 poly A site (21.5%), and the least with 5 poly A sites (9.6%) (Figure 2A). The number of transcripts in each poly A site distribution showed no significant difference (P > 0.05) among the three treatments (Figure 2B). In a more dissecting way, we found that SCN4 infection made soybean generate a significantly greater (P < 0.05) mean number (840 ± SE 26) of >10 poly A sites than SCN5 infection (734 ± SE 12) and the control (698 ± SE 26), indicating APA with more poly A sites might be involved in the incompatible reaction. Enrichment in T and A were detected in the upstream and downstream of the 50-bp position, respectively (Figure 2C). Three GC-rich APA motifs (CAGGGG, GGCTGC, and GGCCGC) were identified in the 50-bp sequence upstream of the poly A sites (Figure 2D). Fusion transcript screening before the remove-redundant analysis revealed 2-12 fusion transcripts for each sample. Similarly, SCN4 infection generated a total of 21 fusion transcripts, while only 14 fusion transcripts were found for SCN5 and the control (Supplementary Table 3), suggesting that fusion transcripts might be participating in defense response.

The Astalavista analysis of novel transcripts demonstrated that all the three treatments had no significant difference in the average number of AS events, 1,925 ± SE 54.6, 2,265 ± SE 145.8; and 2,287 ± SE 130.1 for the control, SCN4, and SCN5 treatments, respectively (Figure 2E), but that the AS events varied among the treatments. Based on the average percentage of AS events among the 9 libraries, there were 40.11% intron retentions, 26.28% alternative 3’ splice sites, 18.81% exon skipping, 14.34% alternative 5’ splice sites, and 0.46% mutually exclusive exons (Figure 2F). The numbers of AS events at alternative 3’ splice sites in both nematode-infected treatments were significantly lower (P < 0.05) than those in the control, while the numbers at 5’ splice sites were significantly greater (P < 0.05) than those in the control.

A total of 22,486 SSRs were identified from 61,900 novel transcripts containing a total of 90,765,603 bp (Supplementary Table 4). Three types of SSR with mono-, di-, and tri-nucleotide accounted for 98% of total SSRs. There were 1,718 SSRs present in compound formation (hybrid microsatellite, a distance of two SSRs < 100 bp) (Supplementary Figure 2 and Supplementary Table 4).

Thus, all the tested structure variations among the control, SCN4, and SCN5 treatments demonstrated that post-transcriptional modification might be involved in nematode resistance or susceptibility in soybean.

A total number of 37,469 ORFs were obtained, including 28,759 complete ORFs. Approximately, 38 and 43% of all predicted CDS coding protein lengths were within 0-100 and 100-200 aa, respectively (Figure 3A). The corresponding complete CDS protein lengths are shown in Figure 3A.

Long Non-coding RNAs associated with plant growth, development, and stress responses have garnered widespread attention (Szcześniak et al., 2015; Urquiaga et al., 2021). In total, 288 lncRNAs were predicted by CPC, CNC, CPAT, and Pfam protein domain analyses; 24 and 90 of them were unique to the CPAT and CPC methods, respectively (Figure 3B); 47.6 (137) and 42.4% (122) of them were classified as sense_lncRNA and lincRNA (long intergenic non-coding RNA), respectively (Figure 3C). Antisense-lncRNA and intronic-lncRNA accounted for 7.3 (21) and 2.8% (8), respectively (Figure 3C). The number of cis-targeted genes (275) regulated by these lncRNAs was nearly four times greater than those (85) of trans-targeted genes. A lncRNA-Pfam protein domain analysis exhibited 164 transcripts matched to the Pfam database, including defense-responsive transcription factors WRKY, TIR, EF-hand, AUX/IAA, zf-RVT (zinc finger in the reverse transcript), zf-CCHC, BolA (bacterial stress-induced morphogen protein) and others, suggesting that lncRNAs may be involved in nematode stress response in soybean. Interestingly, 69 of 164 (42%) Pfam domains were extensin-like protein repeats from two novel transcripts, ONT.12398 (chr 9) and ONT. 12400 (chr 9), and 3 domains were extensin-like regions in ONT.14325 (chr 10). Extensins are a family of hydroxyproline-rich glycoproteins (HRGPs) in plant cell walls that can mediate resistance to viruses, bacteria, fungi, and nematodes (Deepak et al., 2010; Hirao et al., 2012). Especially, extensins identified in SCN syncytium formation (Ithal et al., 2007; Matsye et al., 2011) and nematode-related lncRNA co-expression with HRGP potentially mediating SCN infection (Khoei et al., 2021) consistently confirmed that altered cell wall composition with extensins may play roles in disease defense.

There were 5,337 TFs predicted, including 176 families of TFs (transcription factor, transcription regulator, and protein kinase) with a number ranging from 1 to 296 (Supplementary Table 5). The three treatments, the control, SCN4 and SCN5, contained 1,804, 1,787, and 1,746 TFs, respectively. The most abundant TFs were WRKY (296), AP2/ERF-ERF (APETALA2/Ethylene Responsive Factors, 276), NAC (NAM, ATAF, and CUC, 230), and GRAS (GAI, RGA, and SCA, 217) (Figure 3D).

Function annotation of transcripts obtained from the alternative splicing analysis exhibited 128,648 known and 40,090 novel transcript isoforms (Supplementary Table 6) and 56,865 known and 859 novel genes (Supplementary Table 7). The number of annotated transcripts in all and new isoforms is displayed in Figure 3E. Among them, NR annotation of species distribution indicated that 86.8, 8.9, 0.95, and 0.4% of the transcripts were aligned to G. max, G. soja, Phaseolus vulgaris, and Medicago truncatula, respectively. The GO-enrichment analysis showed that 13,483 out of 100,994 (13.5%) G. max isoforms and 4,350 out of 31,102 (14.0%) novel isoforms were associated with response to stimulus (Supplementary Figures 3A,B).

The CPM density distribution of transcript expressions of samples is displayed in Supplementary Figure 4A, and the dispersion degree of expression level distribution in a sample is shown with a boxplot in Supplementary Figure 4B. The evaluation of the Pearson correlation coefficient among the biological replicates indicated that the r range was from 0.682 to 0.963 (Supplementary Figure 4C). A clear separation of the control treatment from the SCN4- and SCN5- infected treatments is shown in principal component analysis (PCA) (Supplementary Figure 4D).

In total, 2,255 DEGs and 4,167 DETs were identified for all three comparisons; the upregulated and downregulated DEGs or DETs are listed in Table 2, and volcano plots for the DEGs and DETs are shown in Supplementary Figure 5. There were no overlapping DEGs (Figure 3F) or DETs (Figure 3G) found between CK-SCN4-up and CK-SCN5-down and vice versa. The annotated numbers of DEGs and DETs with functional annotation database are listed in Supplementary Table 8. Only 7 DEGs were found between the SCN4- and SCN5- treatments, including 3 DEGs identified in CK-SCN5, Glyma.04G061500 (protein serine/threonine activity, GO:0004674), Glyma.20G205700 (serine-type endopeptidase inhibitor activity, GO:0004867), and Glyma.U029900 (function unknown, plasma membrane, GO:0005886), and 4 DEGs were identified in CK-SCN4, Glyma.03G120700 (calmodulin binding, GO:0005516), Glyma.10G093900, Glyma.11G099300 (a structural constituent of ribosome, GO:0003735), and Glyma.19G125300 (calmodulin binding, GO:0005516).

By GO annotation analysis, even though only 7 DEGs were found between the SCN4- and SCN5- infected soybean roots, the SCN4-treated roots displayed 986 DEGs, while the SCN5-treated roots showed 913 DEGs when each was compared with the control. The top GO classifications for both CK-SCN4 and CK-SCN5 treatments in BP, MF, and CC are displayed in Supplementary Figures 6A,B, respectively. The top three enriched BPs for both treatments were a response to growth hormone, auxin-activated signaling pathway, and multidimensional cell growth. The top three CCs were plant-type cell wall and intracellular membrane-bound organelle for CK-SCN4, cell part, intracellular part, and intracellular membrane-bound organelle for CK-SCN5. The top two enriched MFs for both CK-SCN4 and CK-SCN5 were quercetin 3-O-glucosyltransferase (GTF, EC:2.4) activity and coniferyl-alcohol GTF activity, and the third MF type was pectinesterase (PE, EC:3.1.1.11) activity for CK-SCN4 and quercetin 4’-O-GTF activity for CK-SCN5. Various numbers of enriched DEGs for BP, MF, and CC were observed between CK-SCN4 and CK-SCN5 (Figure 4A). Especially, the numbers of SREs (arrowhead pointing in Figure 4A) in CK-SCN4 are greater than those in CK-SCN5, e.g., response to stimulus (SCN4/SCN5, 613/565), signaling (225/181), immune system process (121/84), detoxification (41/32), and binding (582/493) (Figure 4A).

In further analysis with BP annotation of SREs, a greater number of upregulated DEGs than that of downregulated ones was found in both CK-SCN4 and CK-SCN5 (Figure 4B and Supplementary Table 9). For example, the total number of the up-/down-regulated DEGs added together for each detected SRE was 8,780/6,662 for CK-SCN4 and 7,113/5,932 for CK-SCN5, indicating that more SRE genes were induced in the incompatible reaction than in the compatible reaction (Figure 4 and Supplementary Table 9). The BP response to stimulus included all kinds of biotic and abiotic stresses or stimuli (Supplementary Table 9), e.g., pathogens or pests, organic or inorganic chemicals, hormones, temperature, endogenous or external stimulus, defense response, cell death, hypersensitive response, immune response, gene silencing, and MAP kinases (Figure 4B).

Since the response to growth hormone was the first top enriched GO-BP for both treatments mentioned above, the response to each hormone was inspected and compared. Different responses to eight types of hormones were identified in both resistant and susceptible reactions compared with the control (Figure 4C and Supplementary Table 9). Although the hormones responded to infection with both nematode races, a greater number of upregulated DEGs than down-regulated ones were found in the resistant reaction to SCN4 in SA- (up/down, 31/17), JA- (28/12), abscisic acid- (ABA, 23/17), ET(13/12), gibberellic acid- (GA, 12/9), brassinosteroid- (BR, 7/6), and cytokinin- (CTK, 4/2) mediated signaling pathways and SA- (16/9) and JA- (6/1) induced systemic resistance except for auxin-activated signaling pathway (IAA, 8/10) (Figure 4C). On the contrary, in the susceptible CK-SCN5 reaction, a lower number of upregulated than downregulated DEGs were found in SA- (16/22), ABA- (13/22), and ET- (5/7) meditated signaling pathways, and SA- (8/14) or JA- (0/1) induced systemic resistance. These results suggested that more SREs with more upregulation were activated by the incompatible reaction with the SCN4 infection when compared with the SCN5 infection.

The KEGG annotation showed that 273 and 292 unigenes were involved in 96 pathways in CK-SCN4 and 91 pathways in CK-SCN5, respectively (Supplementary Table 10). Two treatments shared 83 pathways; 13 and 8 unique pathways were detected in CK-SCN4 and CK-SCN5, respectively, e.g., SNARE interactions in vesicular transport (3 DEGs) and oxidative phosphorylation (4 DEGs) only in CK-SCN4 and N-Glycan biosynthesis (4 DEGs) and fatty acid degradation (1 DEG) only in CK-SCN5.

Phenylpropanoid biosynthesis was the top enriched KEGG pathway for both CK-SCN4 (Figure 6A) and CK-SCN5 (Figure 6B). The following top 5 enriched KEGG pathways were photosynthesis-antenna proteins, caffeine metabolism, cutin, suberin and wax biosynthesis, and AGE-RAGE signaling pathway in diabetic complications for CK-SCN4 (Figure 6A), and starch and sucrose metabolism, photosynthesis-antenna proteins, nitrogen metabolism, and AGE-RAGE signaling pathway in diabetic complications for CK-SCN5 (Figure 6B). Among the top 20 enriched KEGG pathways, glycine, serine and threonine metabolism, monoterpenoid biosynthesis, other glycan degradation, phagosome, plant hormone signal transduction, pyrimidine metabolism, and ubiquitin-mediated proteolysis were unique for CK-SCN4; five metabolisms (nitrogen, galactose, cysteine and methionine, ether lipid, and inositol phosphate), and valine, leucine, and isoleucine degradation were unique for CK-SCN5 (Figure 6B). Interestingly, the KEGG pathway circadian rhythm-plant was enriched in the top 8th for CK-SCN4 and in the top 14th for CK-SCN5 (Figures 6A,B).

The PPI analysis among all the 2,255 DEGs indicated a ratio of 7,063/2,632 (2.7 ×) of protein-protein networks for CK-SCN4/CK-SCN5, including 866/194 activation (4.5 ×), 2715/1381 (2 ×) binding, 781/280 (2.8 ×) catalysis, 401/104 (3.9 ×) expression, 462/96 (4.8 ×) inhibition, 874/166 (5.3 ×) ptmod (post translational modification), and 964/411 (2.3 ×) reaction, indicating soybean resistant response to the SCN4 infection activated more protein–protein interactions than a susceptible response to the SCN5 infection, especially in post translational modification, inhibition, and activation.

Two genes, Glyma.03G120700 (GmVQ5) and Glyma.19G125300 (GmVQ70), were identified as significant DEGs with a 2.7-3.2-fold increase in CK-SCN4 than in CK-SCN5, encode calmodulin binding protein (CaMBP) (GO:0005516) carrying conserved VQ (Valine-Glutamine) motif (FxxhVQxhTG), which is capable of interacting with the transcription factor WRKY DNA-binding domain to play vital roles in stress response (Wang et al., 2014). Thus, the interactions between VQ and WRKY proteins encoded by DEGs were explored based on the Pfam domain; 6 (GmVQ5, GmVQ8, GmVQ24, GmVQ43, GmVQ44, and GmVQ70) out of 26 VQs and 19 out of 74 WRKYs annotated were detected in either CK-SCN4, CK-SCN5, or SCN5-SCN4 (Figure 8A). Interestingly, GmWRKY15, encoded by the upregulated Glyma.02G232600, homologous with Arabidopsis AtWRKY33 (WRKYGQK/WRKYGQK) in the plant–pathogen interaction pathway, was identified in CK-SCN4, which could interact with GmVQ24 (AtVQ21, Glyma.06G124400), and both were linked to GmVQ44 (AtVQ25, calcium-binding protein, Glyma.09G111800) with downregulation, which could bind with the two homologous CaM-binding proteins GmVQ5 and GmVQ70 (Figure 8B). VQ and WRKY numbers were designed according to Yu et al. (2016), Wang et al. (2019), respectively.

AtWRKY33 can interact with a VQ protein called mitogen-activated protein (MAP) kinase substrate1 (MSK1, AtVQ21), a substrate of MPK4 (MAP kinase 4), which can be activated in the presence of pathogen attack or flagellin, and then AtWRKY33 is released to induce the expression of phytoalexin deficient 3 (PAD3) in the nucleus, and thereby to increase defense response (Andreasson et al., 2005; Cheng et al., 2012). In addition, MPK3 and MPK6 can also interact with AtWRKY33 together, leading to an increase in phytoalexin-related gene expression (Alves et al., 2014). The overexpression of AtWRKY33 resulting in decreased susceptibility to beet cyst nematode (Heterodera schachtii) (Ali et al., 2013) supports that AtWRKY33 may be involved in soybean cyst nematode resistance. To find the kinase of the substrate GmVQ24 (MSK1) in the annotation data, all expressed MAP kinases were searched and no homologous AtMPK4 was found; interestingly, only two homologous MPK3s encoded by Glyma.U021800 and Glyma.12G073000, each with a 1-fold increase in expression level in CK-SCN4, were identified to interact with GmWRKY15 and GmVQ24. Further searching of other DEGs encoding MAPKs exhibited that two genes encoding MAP kinase kinase 2 (MPKK2) had a 1.5-to 1.8-fold reduction in the expression level only in CK-SCN4 (Figure 8B). Two more MPKKKs encoded by Glma.17G173000 and Glyma.17G177900, which could bind to MPK3 and MPKK2, were found with a 1.2-fold increase of gene expression in CK-SCN4 (data not shown). All these results demonstrated that MPKKK/MPKK/MPK/WRKY-VQ-CaMBPVQ might work together in soybean to defend against SCN4 infection. Based on these data, a model with a VQ–WRKY interaction in response to SCN4 infection was established with two possibilities, one is GmWRKY33 release by MPK3 activation or any other kinase to induce phytoalexin production, which activates plant defense as described in Arabidopsis (Alves et al., 2014); the other is GmVQ44 binding to GmVQ24 (MSK1) to suppress GmVQ44 expression, which results in high expression of CaMBP GmVQ5/70 to induce SCN resistance (Figure 8C).

Another interesting protein–protein interaction was detected in the plant hormone signal branch cytokinin pathway with two upregulated CK-SCN4 DEG genes (Figure 8D), Glyma.06G187000 (ARR6, Log₂FC = 2.2) and Glyma.17G093900 (ARR5, Log₂FC = 1.5), which are classified as two-component response regulator ARR-A (type A Arabidopsis response regulator) family composed of an inner membrane-spanning histidine kinase and a cytoplasmic response regulator to allow organisms to sense and respond to changes in environmental stimuli (Stock et al., 2000). In CK-SCN4, the ARR5-like and ARR6-like proteins were capable of interacting with a mitochondrial chaperone, monothiol glutaredoxin-S15 (GRXS15), and BolA2. GRXS15 encoded by the downregulated DEG Glyma.08G209900 (Log₂FC = -1.8), and BolA2 encoded by a downregulated DEG Glyma.05G188300 (Log₂FC = -1.6) could connect with the transcription factor GATA zinc finger (Glyma.06G086400, Log₂FC = 1.7) (Figure 8C). Mitochondria GRAX15 was reported to play key roles in iron-sulfur protein maturation in the plant (Moseler et al., 2015) and a transcription regulator BolA2 containing a helix-turn-helix (HTH) motif for nucleic acid binding (Kasai et al., 2004). AtBolA3–GRXS17 interaction plays a key role in suppressing abiotic tolerance (Cheng et al., 2011; Qin et al., 2015). In a similar manner, inhibition of the gene expression of both BolA2 and GRAX15 with SCN4 infection resulting in resistance denoted that BolA2–GRAX15 interaction may negatively regulate soybean resistance to SCN4. The function of BolA2–GRAX15 interaction has not been reported yet, but it is known that BolA2 is nucleon-cytoplasmic and interacts with GRAX15 (Couturier et al., 2014). Furthermore, BolA2 was linked to the lncRNA detected above. Therefore, the interactions of ARR5/6-GRAX15-BolA2-GATA might indicate a new complex interaction partaking in plant defense in response to SCN4 infection.

The expression levels of 23 DEGs obtained from full-length transcriptome sequencing were compared with those obtained from qRT-PCR. The expression correlation (R²) between full-length-seq and qRT-PCR was up to 0.6958 (Figure 9A), and the comparison of relative expression levels between SCN4 and SCN5 demonstrated that the qRT-PCR data almost matched with the full-length-seq data (Figure 9B). The relative expression levels of four WRKY and three VQ genes are also displayed in Figure 9C.

In this study, full-length transcriptome analysis on soybean was first conducted to compare responses of the same soybean genotype incompatible and compatible to different H. glycines races. Obviously, the advantages of full-length transcriptome sequencing include high-cost performance, high throughput, no GC specificity, and base bias, long sequencing reads, accurate quantification at the transcriptome level, accurate identification of structure characteristics (e.g., AS, fusion gene, APA) without breaking sequences or gene structures, and DEG/DET analysis at once. Compared to NGS, ONT needs less read numbers to cover the same amount of transcripts. For example, an average of 3-5 × greater number of clean reads was obtained by NGS (Zhang et al., 2017; Neupane et al., 2019; Miraeiz et al., 2020) when compared with the average number of clear reads in this study. In addition, the identified novel genes can generate new AS events.

The change in numbers of AS events at 5’ splicing sites or at 3’ splicing sites after both nematode infections and various transcripts among the three treatments indicated that the nematode infections did trigger differential AS events to generate multiple transcripts that might result in plant defense or susceptibility. A few cases have been reported that AS regulates plant stress responses and/or adaptations (Bedre et al., 2019; John et al., 2021; Martín et al., 2021). For instance, the tobacco gene N resistant to tobacco mosaic virus (TMV) can be alternatively spliced to produce two transcripts that are required to contribute to complete resistance to TMV (Dinesh-Kumar and Baker, 2000). AS regulates ABA and light signaling pathways to coordinate plant growth and stress response (Thiruppathi, 2020).

Interestingly, we found that more abundant APA events per gene were identified with > 5 poly A sites (23.8%) than those in 1-5 poly A sites in soybean, while 2 poly A sites were found as the most abundant APA event distribution in sorghum (Chakrabarti et al., 2020), and 1 poly A site for the tree species Liriodendron chinense in the magnolia family (Tu et al., 2021). The poly A site distribution in plants can be shifted by abiotic stresses (Yan et al., 2021), and novel stress-specific cis-elements in intronic poly A sites may contribute to abiotic stresses (Chakrabarti et al., 2020). In this study, nematode infection increased the > 5-poly A site distribution than control, especially, the > 10 poly A sites in the incompatible reaction further confirmed that APA may be involved in nematode stress response. Generally, stress-responsive APA plays a key role in regulating abiotic stress and biotic stress through crucial stress-responsive genes and pathways (Ye et al., 2019; Chakrabarti et al., 2020; Yan et al., 2021). For instance, Ye et al. (2019) found that APA in rice responds to heat stress tolerance as a negative regulator, to Cd stress by regulating DNA repair and cell wall formation, and to disease stress (bacterial blight, rice stripe virus, and rice blast) by regulating chlorophyll metabolism. Comparisons of APA sites linked with annotated genes or pathways between control and nematode-infected treatments will uncover the APA function in soybean associated with nematode resistance or susceptibility. Furthermore, the difference in the fusion transcripts and lncRNAs among the treatments demonstrates that full-length transcriptome sequencing is a powerful tool to analyze post-transcriptional modification.

Both the GO and KEGG enrichment analyses consistently confirmed that stress response elements and associated pathways (plant hormone signal transduction and plant–pathogen interaction) contribute to plant innate immune response to SCN. Surprisingly, when all the DEGs in these groups/pathways identified in this study were compared with those listed by Zhang et al. (2017), Miraeiz et al. (2020), only 1-3 overlapping genes were found; the only DEG, Glyma.11G207000, encoding leucine-rich repeat-containing protein was found in all three studies. The difference might be caused by various SCN–soybean interaction systems and be affected by inoculation time, genotypes, nematode types, and possible nematode inoculation density. For example, a transcriptome analysis was conducted as early as 8 h post inoculation with HG type 0 on soybean Peking, G. soja PI 468916, Fayette, and Williams 82 by Miraeiz et al. (2020), and later sedentary phase at 3, 5 and 8 days (pooled samples) with HG type 2.5.7 on two G. soja genotypes by Zhang et al. (2017), and 8 days with HG Type 2.5.7 and HG Type 1.2.3.5.6.7 on the same genotype 09-138 in this study. Most transcriptome profiling studies are examined within 2-10 days since SCN feeding site establishment generally required 48 h after inoculation, and syncytium formation and syncytium collapse happens within 2-10 days after inoculation in resistance response (Klink et al., 2007; Kandoth et al., 2011).

Although uniquely differential expression genes to SCN4 or SCN5 infection were identified on 09-138, most DEGs were expressed in both resistant and susceptible reactions with only a small difference (Figure 7), indicating qualitative and quantitative traits. For example, in the plant–pathogen interaction pathway, upstream PPR-CERK1 and downstream Ca²⁺-dependent signal components CDPK, Rboh (triggering reactive oxygen burst, ROS), CaM/CML, HSP, and PR1 were all identified in both compatible and incompatible roots, and downstream WRKYs and Pti6 with upregulation were only detected in the incompatible response (Figures 7E,F), while these components are involved in PTI and/or ETI activity. In the plant-hormone pathway, all phytohormones except for brassinosteroid were activated with 1-4 key DEGs, while the phytohormone networks of the JA, ET, and SA signaling pathways are required for PTI and ETI as well (Cui et al., 2015). Naveed et al. (2020) found that both resistant and susceptible plants showed almost identical transcriptome responses, but that the resistant plants achieved high-amplitude transcriptional reprogramming several hours earlier than the susceptible plants through a defense phytohormone signaling network. Here, we only detected one time point in the later infection stage, and a time series examination will reveal more about how the pathways work together to defend against nematode attacks. The PTI-ETI continuum concept with crosstalk between PTI-ETI signal components effectively activating plant immune responses has the increasing attention of biologists (Mine et al., 2018; Naveed et al., 2020; Yuan et al., 2021), and this concept may also be applicable for SCN–soybean interaction based on the identified signal components in these pathways. In addition, surprisingly, even though 09-138 contained the Peking-rhg1a locus, there were not any DEGs found in that region, suggesting an alternative resistance mechanism in 09-138.

Transcription factors as transcriptional regulators function by binding to the promoter region of target genes and regulating plant response to environmental stress, e.g., altering the expression of cascades of defense genes (Chen et al., 2002; Alves et al., 2014). The novel transcripts identified exhibited more than 5,000 TFs, with the top 5 being WRKY, AP2/ERF-ERF, NAC, GRAS, and bHLH. WRKY TFs, as the largest family of transcriptional regulators, were consistently identified in almost all transcriptome analyses of SCN infection (Ithal et al., 2007; Kandoth et al., 2011; Mazarei et al., 2011; Wan et al., 2015; Zhang et al., 2017; Song et al., 2019; Jiang et al., 2020; Miraeiz et al., 2020). The WRKY family, associated with stress responses in soybean, has been examined, e.g., in response to soybean rust disease (Phakopsora pachyrhizi) (Bencke-Malato et al., 2014), salt stress (Yu et al., 2016), dehydration and salt stress (Song et al., 2016), and soybean cyst nematode (Yang et al., 2017). However, 19 WRKY-DEGs to SCN4 or SCN5 identified in this study were not found in previous studies for SCN infection (Yang et al., 2017; Zhang et al., 2017; Miraeiz et al., 2020), indicating that these identified WRKY factors may be specific to the line 09-138 or the nematode races used in this study.

WRKY factors in plants are divided into five groups (I, IIa + IIb, IIc, IId + IIe, and III), and the interaction between the WRKY domain and its partner domain (e.g., VQ) is involved in signaling, transcription, and other important biological processes (Chi et al., 2013; Alves et al., 2014). The WRKY–VQ interaction has been well-studied in Arabidopsis. As mentioned above, AtWRKY33 (group I)/MSK1 interaction activates plant defense gene expression (Andreasson et al., 2005; Qiu et al., 2008; Alves et al., 2014). AtWRKY25, a homolog of AtWRKY33, is also able to interact with MSK1 and MPK4 in the absence of a pathogen (Cheng et al., 2012). Additionally, AtWRKY33 is able to interact with two other VQ proteins, AtSIB1 (sigma factor binding protein 1, AtVQ23) and AtSIB2 (AtVQ16), in the nucleus to activate resistance to a necrotrophic pathogen, Botrytis cinerea (Lai et al., 2011). In this study, the lack of homologous MPK4 indicates a possible alternate interaction mechanism participating in nematode defense response. In the soybean–pest interaction system, only one GmVQ58 was identified as a negative regulator for soybean resistance to the common cutworm (Spodoptera litura Fabricius), and GmVQ58 could interact with GmWRKY32 (Li X. et al., 2020).

The CaM family is composed of ubiquitous Ca²⁺-binding proteins or calcium sensor proteins, which can play a key role in cellular signaling cascades coupling various environmental stimuli (Zeng et al., 2015). For instance, AtCaMBP25/AtVQ15 was found as a negative effector regulating osmotic stress tolerance during seed germination and seedling growth (Perruc et al., 2004). CaM-mediating signaling can regulate ROS homeostasis directly and indirectly (Zeng et al., 2015). Here, we found DEGs CaM/CML in the plant-pathogen pathway and two CaMBPGmVQ5/VQ70s, which directly or indirectly interact with other VQs and WRKYs to form complex GmVQ5/70-GmVQ44-GmWRKY15-GmVQ24 (MSK1) (Figure 8B). Only the incompatible reaction could activate GmVQ5/70, GmWRKY15, and GmVQ24, and suppress the connector GmVQ44, indicating that VQ–WRKY interactions partake in plant defense response to SCN4. Additionally, the high expression level of GmWRKY15 increasing SCN resistance is consistent with that the homologous AtWRKY33 overexpression decreasing H. schachtii susceptibility (Ali et al., 2013).

The phenylpropanoid pathway has been considered a ubiquitous defense response against pathogens including nematodes, and this pathway is always enriched in SCN–soybean interaction (Edens et al., 1995; Dixon et al., 2002; Zhang et al., 2017; Li et al., 2018; Singh et al., 2019; Miraeiz et al., 2020). Approximately 75% of DEGs in this pathway were peroxidase, and more than 80% of the enriched peroxidase genes could be induced after nematode infection either in susceptible or in resistant response, which is matched with previous reports (Miraeiz et al., 2020), suggesting that peroxidase plays a central role in response to nematode attack. Peroxidase contributes to plant defense by modifying the cell wall composed of lignin, suberin, feruloylated polysaccharides, and HPRG (extensins), enhancing ROS production and phytoalexin productions (Pandey et al., 2017). The key enzymes linked to lignin biosynthesis and uniquely expressed in the incompatible (e.g., HST, Glyma.15G002600) or compatible (e.g., beta-glucosidase, Glyma.15G031300) reaction will be potentially targeted genes for further studies to elucidate plant defense or pathogenesis.

SCN5 infection-induced more DEGs with more upregulated genes in cell wall modification and carbohydrate metabolic process than SCN4 infection but with some quantitatively overlapping DEGs, indicating unique and accumulative DEGs in carbohydrate biological process including metabolism, biosynthesis, catalyst and transport work together to make SCN susceptibility or resistance. Miraeiz et al. (2020) also reported that susceptible W82 showed less responsive genes, more downregulated genes in carbohydrate metabolism and transport proteins, and some overlapping genes with the resistant genotype 8 h post inoculation. Nevertheless, again, there are less common DEGs found between 8 h post inoculation (Miraeiz et al., 2020) and 8 days (this study), a beta-glucosidase was downregulated at 8 h but upregulated at 8 days in this study, suggesting spatially and temporally differential expression. In the starch and sucrose metabolism pathway, sucrose synthase, hexokinase, beta-amylase, and polygalacturonase were suppressed more, and beta-glucosidase, beta-fructofuranosidase, and galacturonosyltransferase were induced more in the compatible interaction than in the incompatible interaction, proving the complexity for nematode pathogenicity. Hofmann et al. (2007) demonstrated that sucrose supply to H. schachtii-induced syncytia depends on the apoplasmic pathway in the early stage during syncytium formation and on the symplasmic pathway in the later stage when syncytia are linked to the phloem. The change of a series of enzymes during a metabolic process may explain various enzyme expression patterns in different SCN–soybean interaction systems both spatially and temporally.