The five soybean lines used in this study: Essex, Peking, PI 88788, PI 561310 and PI 567295 were obtained from the USDA Soybean Germplasm Collection. The soybean cyst nematode population used in this study was validated as HG type 0 (race 3) by evaluating the female index (FI) on each indicator lines at 34-DPI as described in (Niblack et al., 2002), using Essex as the susceptible control. PI 561310 and PI 567295 were additionally included to determine their resistance level to the SCN population used in this study. PI 88788 and Peking are known to possess rhg1- and Rhg4-based resistance mechanisms, while PI 561310 and PI 567295 have been previously identified as partially resistant without possessing rhg1- and Rhg4-based resistance haplotypes (Tran et al., 2019).

A greenhouse experiment was carried out at Agriculture and Agri-Food Canada, Saint-Jean-sur-Richelieu Research and Development Centre, QC, Canada. The experimental design was comprised of two treatment conditions: non-infected controls (CTR) and SCN-infected (INF) plants, for each of the five soybean lines: Essex, Peking, PI 88788, PI 561310 and PI 567295 in triplicates for a total of 60 samples. Seeds from each of the five soybean lines were germinated in moist Pro-Mix for five days in a growth cabinet held at 28˚C, with 50% relative humidity (RH) and a 12h:12h photoperiod until they reached the VE stage. Seedlings were then transferred to individual SC10U Ray Leach “Super Cell” UV cone-tainers (Stuewe and Sons Inc. Corvallis, OR, USA) filled with a moist mixture of 50:50% beach and masonry sand, then inoculated with a suspension of 2,000 SCN eggs from the SCN population of HG Type 0 (race 3). Two series of buckets were prepared to assess gene expression at 5 days post-infection and 10 days post-infection. Buckets were maintained at 28˚C in a temperature-controlled water table to a 50% RH and a 16h photoperiod. Roots were washed free of sand and flash-frozen in liquid nitrogen and stored at -80°C prior to downstream RNA sequencing.

Frozen root samples were ground in liquid nitrogen with mortar and pestle followed by a total RNA extraction using the RNeasy Plant Mini kit (Qiagen, Venlo, Netherlands) according to the manufacturer’s instructions. RNA integrity was analyzed using a 2100 Bioanalyzer (Agilent, Santa Clara, USA). cDNA libraries were prepared and paired-end (2 x 100bp) sequenced using the Illumina NovaSeq 6000 PE100-25M reads platform (Illumina, San Diego, USA) at Génome Québec in Montréal, QC, Canada.

To ensure accurate gene expression data, only samples with RIN value > 8 were kept. We used GenPipes, the main in-house framework of the Canadian Center for Computational Genomics (C3G), to perform all major processing steps (Bourgey et al., 2019). Briefly, adaptor sequences and low-quality score containing bases (Phred score < 30) were trimmed from reads using Trimmomatic (Bolger et al., 2014). The resulting reads were then aligned to the genome (Glycine max v2.1, GCF_000004515.4), using the splice-aware aligner STAR (Dobin et al., 2013). Read counts were extracted from bam files using HTSeq (Anders et al., 2015). The R package DESeq2 was used to identify differences in expression levels between the groups using negative Binomial GLM fitting and Wald statistics: nbinomWaldTest (Love et al., 2014). Similarly “ashr” was used to shrink log₂ fold changes in gene expression data (Stephens, 2017). Genes were considered differentially expressed at an adjusted p-value ≤ 0.01 and log₂fold change difference in expression ≥ 1 for this analysis.

Venn diagrams for inter-line transcriptomic response analyses were generated using: https://bioinformatics.psb.ugent.be/webtools/Venn/. The gene ontology biological process (GO: BP) data for each of the genes of interest identified in these contrasts was extracted from SoyBase (https://soybase.org) and Arabidopsis homologs were extracted from the TAIR10 database (https://www.arabidopsis.org). KEGG pathway enrichment was done using the Mapper search function from the Kyoto Encyclopedia of Genes and Genomes (KEGG) release v106.0 (https://www.kegg.jp/), with Glycine max (gmx) as the organism’s identifier specified.

To validate the RNA-sequencing results, total RNA was re-extracted from root samples using the procedures mentioned earlier. Total RNA was then reverse transcribed into cDNA using the iScript synthesis kit according to the manufacturer’s protocol (BioRad Laboratories, Hercules, CA, USA). Genes found similarly expressed in both PI 561310 and PI 567295 at either 5-DPI or 10-DPI were selected and underwent RT-qPCR using the SsoAdvanced Universal SYBR Green PCR SuperMix (BioRad Laboratories, Hercules, CA, USA). The housekeeping gene Elongation factor 1b (Glyma.02G276600) was used for normalization purposes and samples were analyzed in triplicates. Fold change values were calculated using the pfaffl method (Pfaffl, 2001). Primer sequences used to amplify genes are listed in Supplementary Table 2 .

The resistance levels of the soybean lines used in this study; Peking, PI 88788, PI 561310 and PI 567295 to SCN HG type 0 were assessed as described in (Niblack et al., 2002), using Essex as the susceptible control. Essex was confirmed as susceptible to infection, with Peking and PI 88788 resistant to infection through known rhg-based resistance mechanisms; all three resistance ratings validated through our FI analysis ( Table 1 ). Two lines included in this study: PI 561310 and PI 567295 were previously found to exhibit partial resistance to SCN HG type 0 without possessing either resistant rhg alleles, both identified as moderately resistant (Tran et al., 2019). Our FI analysis ( Table 1 ) confirmed PI 561310 as moderately resistant, and classified PI 567295 as moderately susceptible. The rating differences between PI 561310 and PI 567295 are due to the cutoff for moderately resistant FI < 30; the FI values of PI 561310 and PI 567295 are quite close.

There was a total of 4,284,400,682 raw reads over all 60 RNA-sequencing datasets (five lines, two treatment conditions, two timepoints and three replicates per line). The trimmed reads were aligned to the Glycine max Williams 82 reference genome (Glycine max v2.1), with an 89 – 95% alignment efficiency of the transcriptome to the reference in the 60 samples ( Supplementary Table 1 ). The expression profiles of the resistant lines PI 88788 and Peking clustered together and were distinct from the susceptible line Essex on the two main axes of the principal component analysis ( Figure 1 ). The partially resistant lines PI 561310 and PI 567295 closely resembled PI 88788 and Peking on one principal axis, while on the other axis, PI 567295 showed a profile similar to Essex, and PI 561310 displayed a profile distinct from all other lines. Although the expression profiles of non-infected (CTR) and infected (INF) treatment replicates within individual lines clustered together, due to the significant differences between plant genotypes, they also exhibited differences, when looked individually, that became more pronounced over time.

A comparative analysis of the transcriptomic response of each soybean line to SCN infection was conducted by identifying differentially expressed genes (DEGs) in infected samples compared to non-infected samples within each line at 5-DPI and 10-DPI, respectively. Differentially expressed genes were considered significant at an adjusted p-value ≤ 0.01 and log₂ fold change difference in expression ≥ 1 for this analysis. Figure 2 shows the number of identified DEGs, with the breakdown of upregulated and downregulated gene fractions, for each line at 5-DPI and 10-DPI. Subsequent functional enrichment analyses of Gene Ontology (GO) terms considered all the genes which were differentially expressed.

Comparative analysis of DEGs across soybean lines at 5- and 10-DPI revealed distinct transcriptional responses SCN infection. The full list of inter-line DEGs is provided in Supplementary Table 4 , with key results summarized below.

In order to verify the transcriptome sequencing data, we identified seven DEGs that were observed to be co-regulated in PI 561310 and PI 567295 ( Table 4 ) with defense-related biological processes, and evaluated them across all lines by RT-qPCR. RT-qPCR expression analysis was deemed consistent with the results obtained from the differential expression data obtained from the RNA sequencing experiment ( Supplementary Figure 1 ).

It is well-established that Peking and PI 88788 employ distinct resistance mechanisms against SCN infection, albeit from shared loci. Peking-type resistance is conferred by low-copy number rhg1-a in an epistatic interaction with Rhg4, that leads to a rapid hypersensitive response that results in syncytial collapse and the death of nematodes at the J2 stage, as early as 2-DPI (Klink et al., 2009; Mahalingam and Skorupska, 1996). This early activation likely explains the minimal number of DEGs observed in Peking at both 5- and 10-DPI, as the key defense responses may have already been initiated and resolved by these later time-points. On the other hand, PI 88788-type resistance is attributed to the high-copy number rhg1-b, resulting in a prolonged resistance reaction in which nematodes die at the J3 and J4 stages (Acedo et al., 1984; Klink et al., 2010). This is consistent with the observed DEG profile in PI 88788, in which the numbers of DEGs increase dramatically from just eight genes at 5-DPI to over 200 genes 10-DPI, corroborating that the defense mechanisms are being triggered later in this line. Importantly, the haplotypes conferring Peking-type and PI 88788-type resistances have been found to exhibit elevated and constitutive transcript levels even in the absence of SCN infection and as corroborated in our RNA sequencing data and other independent studies, the rhg1 and Rhg4 genes are not identified as differentially expressed in Peking and PI 88788 lines between non-infected and infected treatments (Chu et al., 2022; Cook et al., 2014, 2012; Lee et al., 2015; Torabi et al., 2023; Zhang et al., 2017)

At both 5-DPI and 10-DPI, both the infected and non-infected replicates form genotype specific clusters. Furthermore, PI 561310 and PI 267295 are observed distinctly separated from one another, as well as from PI 88788, Peking and Essex, indicative of differing patterns of gene expression. PI 561310 shows a substantial increase in DEGs from 5-DPI to 10-DPI, going from 53 to 732 genes, with the majority of differential gene expression being positive at both time-points. This substantial increase in gene regulation over time implies that PI 561310 may employ defense responses that progress and escalate as the infection progresses, peaking later in the infection cycle. In contrast, PI 567295 shows a marked decrease in DEGs from 5-DPI and 10-DPI, going from 570 to 27 genes, in which the majority of differential gene expression is negative at both time-points. Previous work with PI 567295 infected with SCN had reported that J2 nematodes were able to penetrate roots at 2-DPI and initiate syncytium formation, but by 4-DPI delayed nematode development at the J2 stage and evidence of degenerating syncytium around the head of J2 SCN was found (Tran, 2022). At 8-DPI, a clear necrosis of roots were observed, indicating that the plant had initiated a defense response and produced a necrotic layer to disconnect the syncytium from the vascular cylinder starving the females and arresting further development. The transcriptomic response paired with previous cytological observations suggests that the initial response of PI 567295 to SCN infection is substantial, but stabilizes by 10-DPI; consistent with a return to a homeostatic state.

The phenylpropanoid pathway and its associated downstream pathways were found to be differentially regulated across lines in response to SCN infection to varying degrees, with the most pronounced regulation observed in PI 561310. The phenylpropanoid pathway is an important metabolism pathway in plants responsible for the biosynthesis of a variety of secondary metabolites including lignin, flavonoids and phenolic acids, from phenylalanine.

Soybean cyst nematode resistance in soybean is a complex polygenic trait and typically cannot be attributed to a single molecular player or mechanism operating autonomously. Canonical resistance mechanisms function via the synergistic and concerted contributions of multiple molecular components and signaling pathways that enable the plant to mount a robust and effective defense response against pathogens. A subset of genes previously identified to confer resistance against SCN was found in PI 567295, and predominantly at 5-DPI, supporting previous cytological findings (Tran, 2022). At 5-DPI, PI 567295 was observed to have upregulated SAMT1, GmGS13A, peroxidase genes, WRKY-family genes and MYB-family genes. WRKY and MYB are two families of transcription factors that regulate signal transduction and have roles in PTI and ETI (Han et al., 2022). WRKY family members possess a highly conserved WRKY domain and act in the regulation of downstream target genes via binding of cis-acting elements in the promoter regions, whereas MYB members possess a highly conserved MYB DNA-binding domain at the N-terminus and a variable activation domain at the C-terminus, and have a wide range of regulatory functions (Du et al., 2012; Yang et al., 2017). Of the differentially expressed WRKY genes identified, Glyma.13G267600 was not only found upregulated in PI 567295 at 5-DPI, but was also found in Peking and PI 561310 at the same timepoint. Zhang et al. (2017) found that Glyma.13G267600 was the most strongly induced WRKY70 gene following infection with SCN HG type 2.5.7, increased by six-fold, in the resistant wild soybean line S54. Discordantly, most of the differentially expressed MYB genes identified in PI 567295 were downregulated, and when compared to the 31 differentially regulated MYB genes presented by Zhang et al. (2017), only Glyma.10G191000 and Glyma.06G299900 were identified and similarly regulated. In fact, when comparing the list of identified DEGs from PI 567295 with those reported in studies focused on similar objectives using different soybean lines or SCN populations, we consistently found minimal overlap. This deviation could indicate that PI 567295 employs a distinct molecular strategy in response to infection with SCN HG type 0.

Moreover, PI 567295 exhibits a distinct expression profile in response to SCN infection, both compared to the other infected lines studied, but also within the scope of the available literature. Soybean cyst nematode infection is known to induce global transcriptomic changes in soybean, with the observed trend of an increased ratio of upregulated to downregulated SCN responsive genes (Hewezi et al., 2018; Rambani et al., 2020a, 2020b). In contrast, in response to infection with SCN type 0, PI 567295 was observed to have initiated a significant transcriptional response at 5-DPI predominantly composed of downregulated genes. The only other instance of this pattern of response was found when the resistant wild soybean ecotype NRS100 was infected with SCN HG type 2.5.7 (race 5) (Kofsky et al., 2021). Although, the comparison of PI 567295 at 5-DPI with HG type 0 and NRS100 at 6-DPI with HG type 2.5.7 shows a greater quantity of expressed genes in NRS100 than in PI 567295, and unfortunately, as with other DEG analyses, gene expression similarities are limited. Ultimately, while a minor subset of DEGs in PI 567295 following SCN infection have been previously confirmed or correlated with roles in SCN resistance, the majority of the genes currently lack supporting evidence. Additionally, PI 567295 exhibited extensive differential gene regulation at 5-DPI, with 570 DEGs, that reduced drastically, to under 30 genes, by 10-DPI; a significant reversal indicative of a dynamic response.

Considering these findings, we propose that the mechanism modulating the resistant phenotype in PI 567295 may be more heavily based on epigenetic regulation. Epigenetics refers to heritable changes in gene expression that does not alter the underlying genetic sequence. Biochemical modifications, microRNAs and long non-coding RNAs modulate gene expression at the transcriptional and post-transcriptional levels on a genome-wide scale. Epigenetic regulation is a complex and dynamic mechanism that allows for rapid and reversible responses, driving phenotypic plasticity (Cadavid et al., 2023; Pikaard and Mittelsten Scheid, 2014; Yang et al., 2018). When looking at dynamic, genome-wide changes in gene expression, particularly negative regulation, microRNA (miRNA)-mediated post-transcriptional gene silencing comes to mind. Mature miRNAs are small, around ˜20 nucleotides long, non-coding RNAs that cause post-transcriptional gene silencing via cleavage of target transcripts possessing complementary sequences (Hewezi et al., 2016). Indeed, post-transcriptional gene silencing via the action of miRNAs, alongside the modulation of epigenetic marks, have been implicated in soybean-nematode interactions (Han et al., 2022, 2017; Hewezi et al., 2018; Rambani et al., 2020a, 2020b). While further research of this hypothesis falls beyond the scope of this study, the results gathered so far highlights the reaction of PI 567295 to SCN infection with HG type 0 as highly atypical and particularly intriguing. PI 567295 appears to present as a novel source of resistance unlike anything currently presented in the literature and our results provide direction for its further elucidation and characterization.