Transcriptomic Analysis of Soybean Defense Mechanism Against Aphid and Nematode Co-Infestation

1. Introduction

Soybean, Glycine max (L.) Merr., is a globally important crop valued for its high-quality protein and oil [1,2]. In the Midwestern United States, the soybean aphid (SBA), Aphis glycines Matsumura (Hemiptera: Aphididae), and soybean cyst nematode (SCN), Heterodera glycines Ichinohe (Tylenchida: Heteroderidae), are the two most economically significant pests [3,4]. The soybean aphid is an aboveground phloem-feeding herbivore, whereas SCN is a belowground parasite of soybean roots. Their infestations often co-occur and can further exacerbate yield loss [5,6]. Annual economic losses in the United States are estimated at ~$4 billion for SBA and $1.3 billion for SCN [7,8,9]. Management strategies primarily rely on host plant resistance and chemical control [10,11,12]. However, challenges to sustainability have emerged. Pyrethroid resistance in SBA populations has been reported in Iowa, Minnesota, North Dakota, and South Dakota, with additional risks to non-target beneficial organisms [13,14]. Similarly, long-term deployment of SCN resistance has selected for virulent SCN populations capable of overcoming resistance genes (i.e., HG types)[15]. Although host resistance is not widely implemented for SBA, multiple virulent SBA biotypes have been identified in the U.S.; together the evolution of virulent SBA biotypes and SCN races threatens the durability of host plant resistance [15,16,17,18]. Thus, genetic insights from greenhouse-based experiments on SBA and SCN interactions are critical for identifying resistance genes and regulatory networks that can support durable pest resistance in soybean.

Despite occupying distinct niches, above- and belowground herbivores share the host through systemic tissues and can influence each other [19]. Prior studies have investigated SBA–SCN interactions in soybean using demographic approaches [5,6,20,21,22]. For example, McCarville, O’Neal, Tylka, Kanobe and MacIntosh [5] showed that SCN reproduction increased 5.24-fold in the presence of SBA and the fungus Cadophora gregata, while aphid populations decreased by 26.4% when SCN and the fungus were present. Later, McCarville et al. (2014) [6] found that SBA feeding improved soybean’s suitability as a host for SCN, although outcomes varied with cultivar and experiment duration. These demographic studies highlight the complexity of SBA–SCN interactions, but their molecular underpinnings remain unexplored.

Recent studies have substantially advanced understanding of soybean aphid biology and virulence, showing that abiotic stress such as flooding can alter fitness and transcriptomic responses of virulent versus avirulent aphid biotypes, while new chromosome-level genome assemblies of SBA now provide key resources for identifying genes underlying virulence, host adaptation, and effector functions relevant to resistance breeding [23,24]. Field-scale monitoring datasets further demonstrate that aphid population densities vary among soybean lines carrying different Rag (Resistance to Aphis glycines) genes, strengthening links between host genetics and pest pressure under real-world conditions[25]. For soybean cyst nematode (SCN), recent discoveries include novel resistance genes at the qSCN10 locus that confer broad-spectrum resistance independent of classic rhg1/Rhg4 pathways [26], alongside functional validation of defense-related genes such as GmPOD53L that enhance resistance through lignin deposition and peroxidase activity [27]. Updated surveys have also expanded the known geographic distribution of SCN, reinforcing its continued and growing agricultural significance worldwide [28]. From an applied perspective, extension surveys report more than 900 SCN-resistant soybean varieties available for the 2025 growing season, reflecting increased deployment of diverse resistance sources such as Peking in response to evolving SCN populations, while ongoing USDA-NIFA projects are developing molecular tools to detect and manage insecticide-resistant soybean aphid populations [29,30].

RNA sequencing (RNA-Seq) is a widely used approach for both qualitative and quantitative gene expression profiling, enabling detailed insights into transcript abundance and regulatory dynamics [31,32]. In this study, we analyzed RNA-Seq datasets resulted from interaction experiments using two soybean genotypes SCN-susceptible Williams 82 (PI 518671) and SCN-resistant MN1806CN—infested with SBA (biotype 1) and SCN (HG Type 0). Whole-root samples yielded over 1.1 billion reads (61.4 GB) from 47 RNA-Seq libraries. These data provide a strong foundation for dissecting the molecular basis of soybean–SBA–SCN interactions and identifying candidate genes for durable resistance.

2. Results

2.1. Greenhouse Experiment Reveals Asymmetric Above–Belowground Interactions

We quantified SCN reproduction (egg counts at 30 dpi) and SBA population growth (5, 15, and 30 dpi) across SCN-resistant (MN1806CN) and SCN-susceptible (Williams 82; PI 518671) cultivars under four treatment conditions (control, SBA, SCN, and combined SCN:SBA). SCN reproduction differed strongly by cultivar, with substantially greater egg production on the susceptible genotype. As shown in Figure 1, SCN egg counts differed among treatments (one-way ANOVA: F = 87.44, df = 3, P < 0.001). For the SCN-resistant cultivar, the egg numbers for the SCN-only vs. SCN+SBA treatments at 30 dpi were not significantly different, whereas in the susceptible cultivar, SCN egg numbers were significantly higher when SBA was present (Figure 1a). Across both cultivars combined, SCN egg density and SBA abundance showed a strong negative relationship (linear regression: F = 143.5, P < 0.0001, R² = 0.91; Figure 1b), and SCN infestation reduced 30 dpi aphid counts by ~25% (resistant) and ~90% (susceptible; Figure 1c). Together, these greenhouse outcomes indicate that the direction and magnitude of above–belowground interactions depend on host genotype and the stage of infestation, motivating the time-resolved transcriptome analysis below.

2.2. RNA-seq Dataset Quality and Overall Structure of the Transcriptome Signal

Experimental design, library construction, sequencing, and primary read-processing (including trimming and read-mapping/quantification) were previously described and quality-validated in our Scientific Data paper [33]: root RNA-seq generated >1.1 billion reads across 48 libraries (5 and 30 dpi), with high base quality (Phred > 30), GC content ~43–45%, >99% reads retained after trimming, and strong alignment performance (mapped ~73.8–94.3%; uniquely mapped ~67.1–87.6%). After filtering, 43,122 genes were retained and variance-stabilized counts supported robust clustering and ordination; PCA indicated time as the dominant driver of expression variation (PC1 = 28%; p = 1.16 × 10⁻⁶) with treatment also contributing (PC2 = 15%; p = 2.02 × 10⁻⁸). Differential expression analysis across SBA, SCN, and SCN+SBA treatments identified 4,637 DEGs at 5 dpi and 19,032 DEGs at 30 dpi (FDR < 0.01; fold-change > 2), including 242 DEGs shared across treatments at 5 dpi and 1,535 shared at 30 dpi, indicating a large common response that shifts markedly over time.

2.3. Co-Expression Network Analysis Highlights Late Oxidative-Stress Programs

To compare early (5 dpi) and late (30 dpi) responses while reducing model complexity, we constructed weighted gene co-expression network analyses (WGCNA) separately for each time point. WGCNA identified seven co-expression modules at 5 dpi (2,000-gene network) and nine modules at 30 dpi (1,994-gene network) (Figure 2). Notably, hydrogen peroxide metabolic process and reactive oxygen species metabolic process were enriched only at 30 dpi (not at 5 dpi), consistent with stronger oxidative-stress signatures later in the interaction, pointing to stronger late-stage oxidative/ROS engagement under sustained herbivore/parasite pressure.

2.4. Large, Time-Dependent Cultivar Contrasts Reveal Shared Core Responses to SCN, SBA, and Their Combination

We compared resistant vs. susceptible cultivars within each treatment at 5 dpi and 30 dpi. Using stringent thresholds (FDR < 0.01; fold change > 2), we identified 4,637 differentially expressed genes (DEGs) at 5 dpi and 19,032 DEGs at 30 dpi across SCN, SBA, and SCN: SBA contrasts. Despite this extensive divergence, a substantial subset of genes behaved consistently across all three pest contexts within each time point (Figure 3). Specifically, 242 DEGs overlapped among all treatments at 5 dpi, and 1,535 DEGs overlapped among all treatments at 30 dpi (Figure 3). At 5 dpi, the shared 242-gene set was enriched for molecular functions associated with acyl-transfer chemistry and nucleotide/ADP binding, and for pathways including circadian rhythm and flavonoid-related metabolism (Figure 4). Promoter motif enrichment suggested involvement of homeodomain and MYB-related regulatory families (Supplementary File 10). At 30 dpi, the shared 1,535-gene set was enriched for oxidation–reduction functions and broad primary/secondary metabolism, including amino sugar/nucleotide sugar metabolism and phenylpropanoid biosynthesis (Figure 5), with promoter motifs implicating AP2, bHLH, bZIP, MYB and additional plant regulatory families (Supplementary File 11). Collectively, these results define (i) an early shared signature with signaling/biosynthetic features and (ii) a later shared signature dominated by redox and metabolic remodeling.

2.5. A Subset of DEGs Coincide with Known SCN Resistance QTL Regions (QTL-Coincident Signals)

To connect expression signatures to prior genetic mapping, we intersected DEGs with a curated set of SCN QTL-associated genes. Thirteen genes from the shared core sets localized to SCN QTL regions: three from the 5 dpi shared set and ten from the 30 dpi shared set (Figure 6). Among the most notable were Glyma.18G022400 (amino acid transporter-like), Glyma.18G022500 (GmSNAP18), and Glyma.18G022700 (wound-induced protein), which were consistently upregulated (reported range ~2.53 to 5.01 log2 fold change) in the resistant cultivar relative to the susceptible cultivar across time points and treatments. These genes co-localize with the rhg1 region implicated in durable SCN resistance mechanisms, supporting biological coherence between transcriptomic and genetic evidence.

2.6. Within-Cultivar Contrasts Show Fewer DEGs in the Resistant Genotype and Identify Interaction-focused Candidates

We examined treatment effects within each cultivar to identify genes responsive to pest challenge in SCN-resistant (MN1806CN) versus SCN-susceptible (Williams 82). Across contrasts, MN1806CN generally exhibited fewer DEGs than the susceptible genotype, consistent with a more buffered transcriptional response. In the resistant cultivar under SCN+SBA, only four DEGs were unique at 5 dpi—Glyma.03G044900 (dirigent-like), Glyma.13G147600 (2OG-Fe(II) oxygenase), Glyma.16G214400 (Exo70B1), and Glyma.20G089400 (proteasome component domain)—a compact, interpretable set consistent with defense signaling/trafficking and proteostasis (Table 1). At 30 dpi, 100 DEGs were uniquely expressed in the resistant cultivar, including Glyma.04G096400 (strong induction; cystatin/cysteine-type endopeptidase inhibitor activity), multiple cytochrome P450s, and a prominent UDP-glucosyltransferase theme, consistent with hormone/JA-related turnover and specialized-metabolite modification (Table 1; Table 2).

2.7. Transcription Factor Motif Enrichment Suggests Distinct Regulatory Regimes in Early vs. Late Responses

Promoter motif enrichment analyses (300 bp upstream) of DEGs supported a time-structured regulatory architecture. In MN1806CN under combined infestation, the early response showed enrichment for motifs linked to homeodomain and WRKY-related families among induced genes, whereas late responses showed strong WRKY/TBP/bHLH signatures among repressed genes and TBP/bZIP/MYB/AT-hook/bHLH signatures among induced genes (Table 3). These patterns are consistent with layered defense regulation where early recognition and signaling cascades transition into longer-term metabolic and redox remodeling.

2.4. Pathway-Level Analysis Separates Early Lipid/Defense Signaling from Late Carbohydrate and Specialized Metabolism—

PGSEA/KEGG analyses of the most variable genes showed clear temporal partitioning of pathway activity. At 5 dpi, enriched pathways were dominated by immune signaling and wound-associated lipid processes, including plant–pathogen interaction, ubiquitin-mediated proteolysis, cutin/suberin/wax biosynthesis, α-linolenic acid metabolism, and fatty acid degradation (Figure 7, Figure 8 and Figure 9). By 30 dpi, pathway enrichment shifted toward primary metabolism (multiple carbohydrate pathways), fatty acid biosynthesis/elongation, phenylpropanoid and isoflavonoid biosynthesis, and one-carbon metabolism by folate (Figure 7, Figure 8, Figure 9 and Figure 10). This temporal progression suggests that early defense emphasizes signaling and barrier/wound chemistry, whereas late defense involves extensive resource reallocation and secondary metabolite deployment.

3. Discussion

3.1. Host Genotype Shapes the Direction of SCN–SBA Interactions

This study is among the first to couple demographic outcomes with root transcriptomics to examine interaction effects of SCN and SBA on soybean using SCN-resistant cultivar MN1806CN and the SCN-susceptible cultivar Williams 82 (PI 518671). MN1806CN carries the Rps1k gene for resistance to Phytopthera root rot making it resistant to races 1, 4, and 17 and is susceptible to the SBA. Our greenhouse design followed the general approach of McCarville, Soh, Tylka and O’Neal [6], but differed in key parameters (cultivars, aphid density, and sampling schedule). Specifically, we used MN1806CN and Williams 82 (PI 518671), infested plants with 15 biotype 1 SBA (rather than zero, five, or ten biotype 1 SBA), and collected aphid data at 5 and 30 dpi (rather than focusing on a 30-day-after-planting). McCarville et al. (2014)[6] showed both SCN eggs and the number of females increased by 33% in SCN-resistant cultivar and reduced by 50% in the SCN-susceptible cultivar in 30 days after planting. In our greenhouse experiment, at 30 dpi, SBA feeding significantly affected the reproduction of SCN depending on the cultivars. We observed a significant difference in SCN egg counts in susceptible cultivars and no effect in resistant cultivars in the presence of SBA. This suggests the application of aphid increases the reproduction of the belowground SCN in the susceptible plants and no effects in the resistant plants. However, if we consider each treatment, the final SCN eggs counts have increased from the approximate initial counts of 2,000 eggs in both resistant and susceptible cultivars.

Previously, Heeren et al. (2012)[22] utilized resistant and susceptible lines with respect to both soybean aphid and SCN to study the interaction effect of soybean aphid and SCN in the field conditions. The effect of soybean aphid feeding on soybeans on the SCN reproduction was not observed in any of the soybean cultivars as the SCN eggs and aphid densities, less than 100 SCN eggs per 100 cc of soil and less than 10 aphids per plant for <10 days, respectively, were too low in some of the cultivars. We analyzed our data with respect to the aphid counts at 5 dpi, 15 dpi, and 30 dpi to see the trend on the SBA populations. At 5 dpi and 15 dpi, we could not observe a significant difference on the SBA counts between all types of treatments. However, at 30 dpi we observed a significant difference in the SBA counts between all types of treatments. The facilitation at lower herbivore densities and competition at higher herbivore densities might be the reason for differences in the population densities of aphids depending on the length of experiment [34,35]. Particularly, we observed 90% decrease in susceptible plants and 25% decrease in aphid counts receiving SCN in resistant plants. The decline in the SBA populations compared to the plants that did not receive SCN might be due to competition for food resources as both herbivores absorb nutrient assimilates via phloem and affect each other [34,36]. A similar pattern was also observed in a study by McCarville et al. (2014)[6] as SBA populations declined at 30 d and 60 d experiments when infested with five and ten SBA.

The negative association between SCN reproduction and SBA abundance at 30 dpi aligns with broader evidence that root herbivores/parasites can suppress aboveground phloem feeders via plant-mediated changes. For example, Pratylenchus penetrans infection in Brassica nigra reduced infestation by Pieris rapae [37]. Nematode infection reduced aphid fertility in grasses, attributed to altered phloem amino acids [38], and reduced offspring performance in Myzus persicae on Plantago lanceolata infected with P. penetrans [39]. Detrimental effects of nematodes on aphids have also been reported in Brassica systems [40] and in Ammophila arenaria with multiple nematode taxa [41]. Mechanistically, nematode-induced changes in leaf traits (e.g., cuticle waxes, leaf toughness, or water content) may reduce aphid suitability[42], and water stress can disproportionately affect phloem-feeding insects [43]. In addition, increased phenolics and related defense chemistry can contribute to reduced shoot herbivory [37,44,45]. Together, these studies support the interpretation that SCN infection can decrease SBA performance through systemic physiological and defensive shifts, while SBA effects on SCN depend strongly on host genotype and experimental context.

3.2. Time-Resolved Transcriptomes Reveal an Early Signaling Phase and a Late Metabolic/Redox Phase

We leveraged RNA-seq to characterize root transcriptional responses to SBA, SCN, and co-infestation at two time points (5 dpi and 30 dpi). RNA-seq provides sensitive and reproducible transcript measurements across broad dynamic ranges relative to conventional approaches [46]. To reduce analytical complexity and better capture temporal dynamics [47], 5 dpi and 30 dpi datasets were analyzed separately. Across analyses, early responses (5 dpi) were dominated by defense signaling and biosynthetic activities consistent with rapid recognition and response programs, whereas late responses (30 dpi) showed large-scale transcriptional reprogramming involving redox processes, carbohydrate metabolism, and secondary metabolism. The enrichment of hydrogen peroxide and reactive oxygen species-related processes specifically at 30 dpi is consistent with a shift toward oxidative stress and redox remodeling during sustained infestation, potentially reflecting longer-term defense activation and altered root physiology under chronic pest pressure.

3.3. Candidate Mechanisms for MN1806CN Resilience Under Combined SCN: SBA Pressure

In this study, identification of DEGs in resistant and susceptible cultivars was of prime importance particularly with the treatment that both SCN and SBA received to study the genes that are differentially expressed during the interaction. We did a comparison between and within the cultivars. Upon comparison, the discordant expression of genes particularly in resistant cultivar was considered important. In total, we found 4 and 100 discordantly expressed DEGs in resistant cultivar at 5dpi and 30 dpi, respectively. At 5dpi, Dirigent-like protein, 2OG-Fe (II) oxygenase superfamily), Exo70 exocyst complex subunit, and Proteasome component domain protein were differentially expressed in the resistant cultivar that received both SCN and SBA. Dirigent (DIR) -like protein are particularly induced in different kinds of biotic such as wounding and abiotic stresses ranging from drought, cold, abscisic acid (ABA), H2O2, salinity, and osmotic stress [48,49,50]. These proteins play a crucial role in plant defenses against pathogens and lignin and lignan formation [51]. A Dirigent-like protein (Glyma.03G044900) was upregulated by 8.04 log2foldchange at 5dpi. DIR-like proteins were also upregulated during feeding of spruce (Picea spp.) stem-by boring insects (i.e., white pine weevil, Pissodes strobi) in bark tissue and defoliating insects (i.e., western spruce budworm, Choristoneura occidentalis) in green apical shoots [48].

In soybean, GmDIR22 conferred resistance to Phytophthora sojae regulating lignan biosynthesis [52]. Another gene, Glyma.16g214400 that belonged to exocyst subunit exo70 family protein B1 was upregulated by 7.50 log2fold change. The exocyst subunit Exo70B1 interacts with soluble N-ethylmaleimide–sensitive factor attachment protein receptor (SNARE) complex protein SNAP in the process of vesicular trafficking which mediates the exocytosis [53,54]. Previously, the role of α-SNAP which is one of the important genes in Rhg1 mediated SCN resistance has been uncovered in SCN resistance in many studies [55,56,57,58]. SNAP protein involves in vesicle trafficking that affects the exocytosis of food in syncytium which in turn affects the nematode physiology[55]. Another important DEG is Glyma.13G147600 (2OG-Fe (II) oxygenase superfamily) which is down regulated by 3.59 log2fold change. 2OG-Fe(II) oxygenase superfamily that constitutes flavone synthase I (FNS I), flavonol synthase (FLS), anthocyanidin synthase (ANS) and flavanone 3β-hydroxylase (FHT) which play important role in flavonoid biosynthesis [59,60]. The remaining DEG, Glyma.20G089400 (Proteasome component domain protein) was also downregulated by 1.04 log2fold change. Plant proteasomes play important role in the auxin signaling pathway, oxidative stress and hypersensitive responses which are an important component of plant defenses[61,62].

At 30 dpi, we found 100 DEGs that were uniquely expressed in the resistant cultivar of which only 21 genes were upregulated. Particularly, Glyma.04G096400 showed high expression (20 log2fold change) which shows cysteine-type endopeptidase inhibitor activity and possess cystatin domain. The cystatins are basically low molecular weight proteins that inhibit various exogenous proteases or digestive enzymes of invasive pests and pathogens [63,64]. It has been demonstrated that the serine protease activity of H. glycines has been inhibited by cowpea trypsin Inhibitor (CpTI)[65]. Numerous studies on the expression of both native and transgenic cystatins have shown resistance to phytonematodes in a wide range of hosts [66]. The transgenic expression of rice cystatin in eggplant improved nematode resistance in eggplant against root knot nematode, Meloidoyne incognita [66]. Three DEGs belonged to cytochrome P450s [Glyma.03G160100 (CYP94B1), Glyma.10G115900 (CYP71B34), and Glyma.12G239100 (CYP712A1)]. CYP94-genes play important role in Jasmonic Acid signaling pathway via catalyzing the sequential ω-oxidation of JA-Ile [67]. Whereas CYP71 is involved in flavonoid biosynthesis in producing isoflavone and pterocarpan derivatives such as glyceollin in soybean in pathogen-infected tissues [68]. RNA-seq analysis of two Glycine soja genotypes, PI 424093 and PI 468396B, upon infestation by HG type 2.5.7 revealed upregulation of JA, including SA, and ET pathway [69]. Upon pathway enrichment of 100 DEGs, six genes (Glyma.01G046300, Glyma.09G128300, Glyma.09G162400, Glyma.11G000500 Glyma.14G175400, and Glyma.16G158100) were enriched for glucosyltransferase activity and four genes (Glyma.06G079900, Glyma.12G089800, Glyma.03G157800, and Glyma.13G191200) were enriched for calcium binding activity. Previously, the role of glucosyltransferase has been shown in Mi-mediated nematode resistance in tomato [70] which play an important role in carbohydrate and cell wall biosynthesis [71]. The involvement of calcium/calmodulin-mediated signaling has been shown in responses to H. glycines infection in G. soja [69].

3.4. Expression Patterns at SCN QTL-Linked Genes Support Known Resistance Biology

We examined DEGs that coincided with the 251 genes that were assessed from the SCN QTLs. Remarkably, we found three genes at 5 dpi and ten genes at 30 dpi located in SCN QTLs among common genes in comparisons of resistant versus susceptible cultivars. The main purpose was to find if the resistant cultivar, MN1806CN provides rhg1 and Rhg4 mediated resistance. Three genes Glyma.18G022400 (Transmembrane amino acid transporter protein), Glyma.18G022500 [soluble N-ethylmelaimide sensitive factor (NSF) attachment protein (GmSNAP18)], and Glyma.18G022700 (Wound-induced protein WI12) that were upregulated at both 5 and 30 dpi on all treatments are important for rhg1 mediated SCN resistance [55,56,57,72] Grunwald et al., 2022; [73]. We could not find Rhg4 gene as DEG at 5 dpi and was downregulated at 30 dpi in all treatments. This informs that resistant cultivar, MN1806CN might possess rhg1 mediated SCN resistance [18].

3.5. Pathway Enrichment Supports a Working Model Integrating Defense Signaling and Resource Dynamics

To understand the biological functions of the genes, PSGEA analysis was carried out. PSGEA analysis showed distinct enriched pathways at 5 dpi and 30 dpi. We observed an enrichment of Plant-pathogen Interaction, Ubiquitin mediated proteolysis, Phenylalanine, tyrosine and tryptophan biosynthesis, Cutin, suberin and wax biosynthesis, Alpha-Linolenic acid metabolism, fatty acid degradation pathways were enriched at 5 dpi. Plant-pathogen Interaction and Ubiquitin mediated proteolysis play important role in plant immunity. The interaction of plant and pathogen involves pathogen-associated molecular patterns (PAMPs) of pathogens by pattern-recognition receptors (PRRs) of host and are regulated by E3 ubiquitin ligase[62,74]. E3 ubiquitin ligase has been previously reported in involvement in phytonematodes such as Heterodera schachtii [75] and Globodera rostochiensis [76]. Phenylalanine, tyrosine, and tryptophan biosynthesis pathway are related to the shikimate pathway. It is shown that the chorismite mutase enzyme present in root-knot nematode and potato cyst nematode alter the shikimate pathway of the host plant [77]. Other pathways such as Cutin, suberin and wax biosynthesis, Alpha-Linolenic acid metabolism, and fatty acid degradation are related to pathways related to plant lipids which are important for the production of JA, cutins, and suberins in plant defense via wounding [78,79]. At 30 dpi, most of the pathways were related to carbohydrate metabolism, fatty acid metabolism including fatty acid biosynthesis, fatty acid elongation, Phenylpropanoid biosynthesis, isofalvonoid biosynthesis, and one carbon pool by folate. Phenylpropanoid biosynthesis and isoflavonoid biosynthesis pathways are particularly related to the metabolism that produces various compounds such as flavonoids, anthocynanins, lignin, suberin, salicylic acid, coumarins and furanocoumarins [60,80]. It has been shown that phloem-feeding-insect, whitefly Bemisia tabaci, when infested in Nicotiana tabacum activates the phenylpropanoid pathway [81]. We expect pathway for carbohydrate metabolism to be enriched at 30 dpi as SBA and SCN might be competing for the limited food resources of the plant [35,82,83]. The pathway one carbon pool by folate is related to folate-mediated one-carbon metabolism. It has been shown that the SHMT (GmSHMT08) gene was confirmed as the resistance gene at the Rhg4 locus that catalyzes methylene carbon of glycine to tetrahydrofolate (THF)[84,85]. GmSHMT08 changes its enzymatic properties in the resistant allele that negatively affects the folate homeostasis in the syncytium resulting in hypersensitive responses (HR) leading to programmed cell death (PCD) [26,84,86].

3.6. Practical Implications and Limitations

In summary, we presented several DEGs that were changed in the days after SCN and SBA infection during SCN susceptible and resistant soybean interactions. Many genes revealed various pathways and networks involved in the interaction effects of SCN and SBA on soybeans. Although a high number of genes were differentially expressed when compared between resistant and susceptible cultivars, however, the comparison within the cultivars exhibited the fewer DEGs conferring resistance against both SBA and SCN in resistant cultivar. One limitation of the experiment was the cultivar which was resistant to SCN was susceptible to SBA. Various GO enrichment and KEGG pathways identified several molecular mechanisms involved in SCN SBA interaction. The role of transcription factors in the SBA SCN interaction in several comparisons has been identified which can be used for future research and breeding for SCN and SBA in soybean. Further work warrants functional validation of identified DEGs.

4. Materials and Methods

4.1. Plant Material and Pest Populations

Two soybean cultivars were used: Williams 82 (PI 518671; SCN-susceptible) and MN1806CN (SCN-resistant to Heterodera glycines HG type 0/race 3). Soybean aphid (Aphis glycines) biotype 1 colonies were obtained from The Ohio State University and maintained on a susceptible soybean line (LD12-15838R). The SCN population used was HG type 0, defined by <10% reproduction on indicator lines and avirulence to major SCN resistance sources.

4.2. Greenhouse Experimental Design and Inoculation/Infestation

Experimental design, infestations, sampling, RNA extraction, library preparation, and sequencing are elaborated in the first Authors PhD dissertation work [87] and published as Data Descriptor [33]; key parameters are summarized here. A randomized complete block design (RCBD) was implemented with six blocks and factorial treatments defined by soybean genotype, soybean aphid (SBA) infestation, and soybean cyst nematode (SCN) infestation. For each genotype, treatments included control (no pests), SBA only, SCN only, and SCN+SBA. Cone-tainers were filled with a soil:sand mixture (3:1) prepared using construction sand and SCN-infested clay soil. Approximately 125 cc of mixture was added per cone-tainer (1.5″ diameter; 8.25″ depth; 164 cc volume). For SCN treatments, each cone-tainer received ~2,000 SCN eggs at planting. Cone-tainers (three seeds each) were placed into 2.0-gallon buckets filled with construction sand and maintained in a water bath to keep soil temperature between 26.7–28.9 °C to support SCN development (~30 d). Plants were grown under a 16 h light/8 h dark cycle and thinned to one plant per cone-tainer at the V2 stage. For SBA treatments, V2 plants were infested with 15 mixed-age SBA (biotype 1) using a fine-tip brush applied to the abaxial surface of the first trifoliate. Buckets were covered with mesh netting to prevent aphid movement among experimental units.

4.3. Aphid Counts and SCN Egg Quantification

Soybean aphid populations were counted on each plant at 5, 15, and 30 days post-infestation (dpi). SCN eggs were quantified at 30 dpi. For egg extraction, soil and roots were washed and the suspension was passed through an 850 µm sieve and retained on a 250 µm sieve. Females/cysts collected on the 250 µm sieve were macerated using a motorized rubber stopper to release eggs; eggs were recovered on a 25 µm sieve after passing through a 75 µm sieve. Eggs were suspended in 50 mL water and counted under a compound microscope using a 1 mL aliquot as a representative sub-sample.

4.4. Root Sampling, RNA Extraction, Library Preparation, and Sequencing

Whole-root tissue was collected at 5 and 30 dpi by snap-freezing in liquid nitrogen and storing at −80 °C. Frozen roots were ground to a fine powder in liquid nitrogen and total RNA extracted using the PureLink RNA Mini Kit (Invitrogen). RNA was DNase-treated (TURBO DNase, Invitrogen), integrity assessed by 1% agarose gel electrophoresis, and concentration measured using a NanoDrop 2000 (Thermo Fisher Scientific). cDNA libraries were prepared with the NEBNext Ultra II RNA Library Prep Kit (96 single index) and sequenced on an Illumina HiSeq 3000 platform (single-end, 100 bp) at the Iowa State University Sequencing Facility.

4.5. RNA-Seq Preprocessing, Mapping, Quantification, and Normalization

Raw read quality was assessed using FastQC (v0.11.3) and summarized with MultiQC (v1.3). Low-quality bases were removed with Btrim64 (v0.2.0) (quality threshold >20; 5-bp sliding window). High-quality reads were quantified by mapping to soybean primary coding sequences (Wm82.a2.v1; Phytozome) using Salmon (v0.9.1). Downstream analyses were conducted in iDEP (v0.81). Genes were filtered at ≥0.5 counts per million (CPM) in at least one sample, normalized using DESeq2, and transformed using DESeq2’s regularized log (rlog) for clustering/ordination.

4.6. Differential Expression and Candidate-Gene Analyses

Differentially expressed genes (DEGs) were identified using DESeq2 with |fold change| ≥ 2 and FDR ≤ 0.01. Effects of cultivar, treatment, and their interaction were tested using the model: expression ~ cultivar + treatment + cultivar: treatment. DEG annotations were obtained from SoyBase.

4.7. Co-Expression Networks, Functional Enrichment, and Motif Analysis

To reduce analytical complexity, count matrices were analyzed by time point for network analyses. Weighted gene co-expression networks were constructed in WGCNA using the top 2,000 most variable genes with soft-threshold power = 5 and minimum module size = 20. Functional enrichment of DEGs (GO terms and KEGG pathways) was performed using ShinyGO, REVIGO, and iDEP. Promoter motif enrichment was assessed by scanning 300 bp upstream regions of DEGs for transcription factor binding motifs (ShinyGO and iDEP). Pathway-level patterns were visualized using MapMan, with functional binning assigned via Mercator. PGSEA was used for gene set enrichment analyses with an FDR-based significance cutoff.

4.8. Statistical Analyses for Phenotypes

Aphid and SCN counts were analyzed in GraphPad Prism (v8.0.2). SCN egg counts (30 dpi) were analyzed by one-way ANOVA followed by Tukey’s multiple comparisons test. Aphid counts (5, 15, and 30 dpi) were analyzed using two-way ANOVA with the Geisser–Greenhouse correction, followed by Tukey’s multiple comparisons. Linear regression was used to evaluate relationships between SBA density and SCN egg counts at 30 dpi.