, Phytophthora root and stem rot (PRR) caused by the soil-borne oomycete Phytophthora sojae (Kaufmann and Gerdemann) (Dorrance, 2018) is a devastating and widespread soybean disease [Glycine max (L.) Merr.]. The infection cycle begins with the introduction of motile zoospores of P. sojae being attracted to isoflavones, such as daidzein and genistein, under disease-conducive soil conditions. Upon reaching the surface of the soybean roots, these zoospores are transformed into cysts and germinate into hyphae that penetrate the plant cells within a few hours of inoculation (Morris et al., 1998). Oospores are subsequently produced in the colonized soybean roots and stems. The infected soybean plants show root and stem rot, yellowish leaves, wilting, and death, leading to substantial yield loss (Schmitthenner, 1985).
PRR has been reported to be a major disease in most soybean-growing countries, including the United States, Argentina, Brazil, China, Japan, South Korea, and Australia (Dorrance, 2018). Since its first identification in the 1950s, PRR has been a primary biotic yield-limiting factor in soybeans. Estimates indicate that PRR considerably impacted soybean yields from the late 1990s to the mid-2010s, with an average annual loss exceeding 1 million tons in North America (Allen et al., 2017; Koenning and Wrather, 2010; Wrather and Koenning, 2006, 2009; Wrather et al., 2001).
The management of PRR in soybean production primarily relies on two types of genetic resistance: R-gene-type, called Resistance to P. sojae (Rps), complete resistance and partial resistance. These genetic mechanisms govern R-gene-mediated hypersensitive responses and impede the spread of the pathogen through various plant defense mechanisms (Dorrance, 2018). Currently, through extensive screening of a large number of soybean germplasms, many resistant genotypes have been identified and used for genetic dissection of the trait aiming to identify genomic regions associated with P. sojae (Anderson and Buzzell, 1992; Athow, 1980; Athow and Laviolette, 1982; Bernard and Cremeens, 1981; Bernard et al., 1957; Burnham et al., 2003; Buzzell and Anderson, 1992; Chen et al., 2021; Cheng et al., 2017; Fan et al., 2009; Grau et al., 2004; Hartwig et al., 1968; Jiang et al., 2020; Kilen and Keeling, 1981; Kilen et al., 1974; Laviolette and Athow, 1977; Li et al., 2016, 2017; Lin et al., 2013; Mueller et al., 1978; Niu et al., 2017; Ping et al., 2016; Ploper et al., 1985; Sahoo et al., 2017, 2021; Sugimoto et al., 2011; Sun et al., 2011, 2014; Wang et al., 2021; Weng et al., 2001; Wu et al., 2011a, 2011b; Yao et al., 2010; Yu et al., 2010; Zhang et al., 2013a, 2013b; Zhong et al., 2018a, 2018b, 2019; Zhu et al., 2007). At least six Rps genes identified through these studies have been incorporated into commercial breeding programs to enhance resistance to PRR in soybeans (Dorrance, 2018).
The molecular responses of soybeans to P. sojae infection have been extensively studied. Compared with the susceptible recurrent line Williams, 10 near-isogenic lines (NILs), each containing a unique Rps gene/allele, revealed distinct transcriptomic patterns in resistance responses to P. sojae infection (Lin et al., 2014). The study identified significant differentially expressed genes (DEGs) involved in ethylene (ET), jasmonic acid (JA), and reactive oxygen species (ROS), highlighting their involvement in the defense responses of these 10 soybean NILs (Lin et al., 2014). Wong et al. (2014) reported that P. sojae-induced microRNAs and phased small interfering RNAs regulate defense-associated genes in soybeans following infection. RNA sequencing of NILs segregating at a major quantitative resistance locus (QDRL-18) predicted the upregulation of the serine-threonine kinase-coding gene as a susceptibility factor in susceptible NILs (Karhoff et al., 2022). A recent study combining high-throughput sequencing and deep learning predicted a gene regulatory network of the defense-related transcription factors, including WRKY and RAV family members, along with their binding sites in the soybean genome involved in orchestrating an immune response against P. sojae (Hale et al., 2023).
Although many Korean soybean germplasms with qualitative or quantitative resistance to P. sojae have been documented as useful genetic sources in the U.S. soybean research, the efforts mining resistant varieties and related resistance genes were only recently initiated in South Korea. Kang et al. (2019) first investigated the presence or absence of Rps-mediated resistance in Korean soybean varieties against four different P. sojae isolates. Further research has mapped a few resistance loci on chromosomes 3 and 18 using four recombinant inbred line (RIL) populations, including Daepung (DP) × Socheong2 (SC2) (Jang et al., 2020a, 2020b; You et al., 2023a, 2023b). Linkage mapping identified two significant loci on the distal positions of chromosome 18 in the DP × SC2 population (Jang et al., 2020b) that have been documented to be associated with resistance to P. sojae (Lee et al., 2013; You et al., 2023a). Linkage mapping for qualitative or quantitative traits is a conventional tool widely used to identify genomic regions associated with traits of interest. Contemporary high-throughput sequencing or chip-based single nucleotide polymorphism (SNP) genotyping generally facilitates high-resolution genetic mapping by determining more precise and narrower intervals for target loci. However, pinpointing causal genes from the identified genomic regions is often challenging (Jang et al., 2020a, 2020b). To overcome this obstacle, combined linkage mapping and transcriptome analyses have gained interest. The coupled analyses and subsequent interpretation provide more comprehensive knowledge to identify the candidate genes and attain revelatory insight into the genetic mechanisms of target traits, which have been documented in several recent studies on tolerance to drought stress, high temperature, and pathogen resistance in many crop species (Almeida-Silva and Venancio, 2021; Azam et al., 2023; Wei et al., 2021). Thus, in this study, we aimed to compare transcriptomic changes following infection with P. sojae between the parental cultivars used in the previous linkage mapping study. We pinpointed plausible candidate genes determining resistance against P. sojae by bridging the transcriptome analysis and the last linkage.
Materials and Methods
Plant materials, P. sojae inoculation, and tissue sampling
To compare global transcriptome activities in response to P. sojae isolate 2457, 10-15 seedlings of two soybean cultivars, SC2 (resistant) and DP (susceptible), were grown in 13 cm pots for 7 days. Phytophthora sojae was inoculated following a hypocotyl inoculation technique (Dorrance et al., 2004). Briefly, the mycelia of P. sojae isolate 2457 were grown on a clarified V8 medium for 10 days. The mycelia were macerated twice using a 50 ml syringe and transferred into a 10 ml syringe for inoculation. Afterward, 0.2-0.4 ml of the mycelial slurry was inoculated onto the stem tissue of the seedlings using sterile syringes. Stem tissues from the upper and lower 1 cm of the inoculation site were sampled immediately after inoculation (referred to as 0 h after inoculation; HAI). The inoculated seedlings were then incubated under dark and humid conditions (>90%), and stem tissues were sampled at 6 and 12 HAI.
RNA preparation, sequencing library construction, and sequencing
Approximately 0.1 g of P. sojae inoculated stem tissue of SC2 and DP were immediately frozen using liquid nitrogen for RNA isolation. Total RNAs of each tissue sample were isolated following the manufacturer’s protocol using the RNeasy Plant Mini kit (Qiagen, Hilden, Germany). The samples were then shipped to Theragen Bio (Seongnam, Korea) for high-throughput sequencing. Total mRNA samples were fragmented at 94°C for 8 min, which were then synthesized into cDNA using random primers. Illumina sequencing adapters were ligated into the cDNA samples using a TruSeq RNA Prep Kit v2 (Illumina, San Diego, CA, USA). Subsequently, the adapter-ligated cDNA samples were amplified using PCR. Libraries qualified using the Agilent 2100 Bioanalyzer were sequenced using the Illumina Novaseq6000 sequencing platform.
Preprocessing of raw data, transcriptome analysis, and detection of DEGs
The sequencing quality of the raw paired-end RNA sequencing (RNA-seq) reads was examined using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We then mapped the RNA-seq reads to the reference primary soybean transcriptome sequences (Glyma.Wm82.a2.v1) (Schmutz et al., 2010) and counted the number of reads mapped to the transcripts using Kallisto (Bray et al., 2016). Read counts for 17 RNA-seq libraries (one library was eliminated because of biased expression; see details in the results section) were subjected to the DESeq2 pipeline for gene expression normalization (median of ratios), correlation, principal component (PC), and DEG analyses (Love et al., 2014). For correlation and PC analyses, the read count matrix was transformed using variance-stabilizing transformation, as suggested in the DESeq2 tutorial on the Bioconductor webpage (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html). Genes that were significantly (P < 0.05) upregulated or downregulated (|log2 fold change| > 1) were selected as DEGs. The Python module pyVenn (https://github.com/tctianchi/pyvenn) was used to visualize the number of DEGs using a Venn diagram.
Functional enrichment analysis
We conducted a functional enrichment analysis for the 14 sets of upregulated or downregulated DEGs to investigate the biological implications of the selected DEGs. We first retrieved gene ontology (GO) annotations for the DEGs from the soybean reference genome annotation (Glyma.Wm82.a2.v1) from Soybase (https://soybase.org/). Next, we compared the number of DEGs to the total reference genes in a specific functional category. Statistical significance was calculated using the hypergeometric test in the Python module SciPy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html). The most significant 30 GO terms for biological processes were retrieved from each gene group and are displayed. A scatterplot in the Python module, seaborn (https://seaborn.pydata.org/index.html), was used to visualize the GO enrichment analysis results.
Candidate gene analysis
Based on the functional enrichment analysis, we identified several genes related to commonly conserved and genotype-specific disease responses. The two cultivars, SC2 and DP, used in this study were used as the parental genotypes of a RIL population in a genetic mapping study (Jang et al., 2020b). The study reported two genomic regions responsible for the resistance to P. sojae isolate 2457. We compared a list of defense-related DEGs with a list of candidate genes in the identified genomic regions to select plausible candidate genes.
Data availability
The Illumina sequencing data presented in this study can be found in the BioProject database of the National Center for Biotechnology Information (PRJNA1070321): https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1070321.
Results
Quality control and mapping statistics of RNA-seq libraries
Using RNA-seq, an average of 28,611,936 paired-end reads were produced (Supplementary Table 1). Before mapping the RNA-seq reads, the sequencing quality of all the RNA-seq libraries was evaluated. As shown in Supplementary Fig. 1, all RNA-seq libraries showed acceptable sequencing quality (Q > 20). Therefore, we mapped the raw RNA-seq reads directly to the reference transcriptome sequence of G. max (Glyma.Wm82.a2.v1) (Schmutz et al., 2010). The alignment rate of uniquely mapped RNA-seq reads ranged from 86.8 to 88.9% (Supplementary Table 1) overall examined libraries, suggesting the high quality of the RNA-seq data.
Preliminary inspection of all 18 RNA-seq libraries identified one sample from DP at 12 HAI with a biased transcriptome profile (Supplementary Fig. 2), which was excluded from the data, and the remaining 17 samples were used for RNA-seq analysis.
Global gene expression pattern of resistant and susceptible soybean genotypes in response to P. sojae
Global gene expression patterns were analyzed using correlation and PC analyses. For the 17 samples, hierarchical clustering was reasonable based on the genotype and time point (Fig. 1A). Three biological replicates of each sample (i.e., genotype by time point) were primarily grouped in all samples, and the groups were well separated, indicating that P. sojae infection experiments for tissue sampling were carried out meticulously (Fig. 1A). PC analysis showed better dispersion of samples by both timepoint and genotype. The PC1 explained the largest difference (81%) of the total variance, and global gene expression at 0 HAI was distinguishable compared with those at 6 and 12 HAI, regardless of genotype (Fig. 1B). These results indicate that the direction of transcriptional reprogramming of the R and S genotypes against P. sojae may be largely conserved despite key differences. The PC2 explained 7% of the total variance, indicating the substantial transcriptional difference between the two genotypes (Fig. 1B). These results suggested that genotype-specific resistance mechanisms enhance resistance in SC2 against P. sojae.
Analogy in transcriptomic changes between R and S genotypes following P. sojae infection
Dramatic transcriptional reprogramming occurred between 0 and 6 HAI, and the patterns of gene expression shown in PC analysis were similar in both genotypes, presumably because common molecular responses were activated within 6 HAI (Fig. 1B). To explore whether conserved defense mechanisms in the resistant and susceptible genotypes were activated, gene expression patterns of the samples at 6 or 12 HAI were compared with those at 0 HAI by genotype (grey arrows in Fig. 2). In the resistant genotype, SC2, 7,160 (Group 1) (Supplementary Table 2) and 9,228 (Group 2) (Supplementary Table 3) genes were upregulated and downregulated at 6 HAI compared with those at 0 HAI, respectively (Fig. 3A). Similarly, 6,364 (Group 3) (Supplementary Table 4) and 9,756 (Group 4) (Supplementary Table 5) genes were upregulated and downregulated, respectively, at 12 HAI compared with those at 0 HAI (Fig. 3B). In contrast, 5,848 (Group 5) (Supplementary Table 6) and 8,939 (Group 6) (Supplementary Table 7) genes were upregulated and downregulated, respectively, in the susceptible genotype, DP, at 6 HAI compared with 0 HAI (Fig. 3A). At 12 HAI in the DP, 4,370 (Group 7) (Supplementary Table 8) and 7,971 (Group 8) (Supplementary Table 9) genes were upregulated and downregulated, respectively, compared with 0 HAI (Fig. 3B). Approximately 25% of the upregulated genes were common to all four groups (SC2, 12 vs. 0 HAI; SC2, 6 vs. 0 HAI; DP, 6 vs. 0 HAI; DP, 12 vs. 0 HAI; upregulated Groups 1, 3, 5, and 7 in Fig. 3A). Likewise, approximately 35% (4,787) of downregulated genes were common among the four groups (SC2, 12 vs. 0 HAI; SC2, 6 vs. 0 HAI; DP, 6 vs. 0 HAI; 7, DP, 12 vs. 0 HAI; downregulated Groups 2, 4, 6, and 8 in Fig. 3B). These results indicated that a substantial number of DEG groups were commonly reprogrammed to exhibit defense-related responses.
Functional enrichment analyses were conducted for the eight groups of upregulated and downregulated DEGs based on GO annotation (Fig. 4, Supplementary Fig. 3) to determine the biological implications underlying DEGs. The analyses revealed that several disease response-related GO terms, including response to chitin (GO:0010200), respiratory burst involved in defense response (GO:0002679), response to wounding (GO:0009611), regulation of plant-type hypersensitive response (GO:0010363), JA-mediated signaling pathway (GO:0009867), salicylic acid (SA) mediated signaling pathway involved in systemic acquired resistance (SAR) (GO:0009862), defense response to fungus (GO:0050832), negative regulation of programmed cell death (GO:0043069) and defense response (GO:0031348), response to JA (GO:0009753) and ET (GO:0009723), JA (GO:0009695) and ET biosynthetic process (GO:0009693), and mitogen-activated protein kinase (MAPK) cascade (GO:0000165), were commonly overrepresented in Groups 1, 3, 5, and 7 irrespective of genotypes (Fig. 4, Supplementary Tables 10-13).
However, no functional category related to disease response was overrepresented in the four groups of downregulated DEG (Groups 2, 4, 6, and 8). However, some GO terms associated with photosynthesis were enriched (Supplementary Tables 14-17, Supplementary Fig. 3). These results implied that SC2 and DP underwent similar transcriptomic changes following P. sojae infection.
Most of the significantly enriched GO terms were analogous between the two genotypes; nevertheless, discrepancies in the enriched GO terms were also detected by genotype. For example, the SA-mediated signaling pathway (GO:0009863) was particularly overrepresented in DEG Groups 1 and 3 but not in DEG Groups 5 and 7 (Fig. 4). Such genotype-specific SA-mediated signaling pathways may play an essential role in conditioning resistance in the SC2.
Genotype-specific resistance mechanisms in response to P. sojae
PC analysis revealed a prominent separation in global gene expression following P. sojae infection between the R and S genotypes (Fig. 1B). To discover genotype-specific resistance mechanisms in response to P. sojae, global gene expression levels were compared between SC2 and DP at 0, 6, and 12 HAI (blue arrows in Fig. 2). As shown in Fig. 5A, 1,748, 2,416, and 2,257 genes were upregulated in SC2 at 0 (Group 9) (Supplementary Table 18), 6 (Group 10) (Supplementary Table 19), and 12 HAI (Group 11) (Supplementary Table 20). Meanwhile, 3,088 (Group 12) (Supplementary Table 21), 2,061 (Group 13) (Supplementary Table 22), and 3,629 (Group 14) (Supplementary Table 23) genes were downregulated in SC2 at the respective time points (Fig. 5B). Furthermore, approximately half of each DEG group was temporally upregulated or downregulated in SC2, except for DEG Group 13 (Fig. 5A and B). These results indicate that genotype-dependent regulated DEGs were temporally specifically controlled after inoculation with P. sojae.
Next, we performed the functional enrichment analysis of the DEGs upregulated in SC2 across the three distinct groups. The results revealed that GO terms related to plant defense mechanisms, such as response to chitin (GO:0010200), respiratory burst involved in defense response (GO:0002679), SAR (GO:0009627), defense response to fungus (GO:0050832), SA (GO:0009697), JA (GO:0009695), and ET biosynthetic processes (GO:0009693), along with response to JA (GO:0009753), were particularly prominent in Groups 10 and 11. However, these terms were absent in Group 9 (Fig. 6, Supplementary Tables 24-26).
These defense-related GO terms were also overrepresented in Group 12 genes (Fig. 6, Supplementary Table 27). These genes were characterized as DEGs downregulated in SC2 at 0 HAI. Moreover, in Groups 13 and 14, a common enrichment of GO terms associated with vital cellular processes, including plastid organization (GO:0009657), photosystem II assembly (GO:0010207), and rRNA processing (GO:0006364) (Supplementary Fig. 4) was observed. These GO terms differed from those enriched in Groups 10, 11, and 12 (Supplementary Tables 28 and 29, Supplementary Fig. 4). These observations highlight that the expression of defense-related genes was genotype-specifically low at 0 HAI but elevated within 6+ HAI only in SC2 compared with DP.
Although the enriched defense-related terms for Groups 10, 11, and 12 were broadly analogous, a considerable number of terms, such as defense response to bacterium (GO:0009816), ET-activated signaling pathways (GO:0009873), and JA-mediated signaling pathways (GO:0009867), SA-mediated signaling pathway involved in SAR (GO:0009862), regulation of plant-type hypersensitive response (GO:0010363), MAPK cascade (GO:0000165), and negative regulation of defense response (GO:0031348), were overrepresented in Group 12. Some of these GO terms were also found to be overrepresented in Groups 10 and 11 (Fig. 6). These findings indicated that the expression of genes involved in some defense responses was temporally reprogrammed in a genotype-dependent manner. Furthermore, a considerable number of GO terms enriched in DEG Groups 1, 3, 5, and 7 were also found to be enriched in DEG Groups 10, 11, and 12. This finding suggests that although common resistance mechanisms were activated in both cultivars, the enhanced resistance of SC2 was determined by genotype-specifically enhanced transcriptional regulation of cognate biological processes and relevant genes.
Candidate genes for enhanced resistance of SC2 against P. sojae isolate 2457
The R genotype SC2 exhibited improved resistance against P. sojae by genotype-specific and genotype-dependent elevated transcriptional regulation of SA-, JA-, and ET-dependent SAR within 6+ HAI in P. sojae. These biological processes represent broad-spectrum resistance rather than pathotype-specific resistance. In our previous study, we identified two genomic loci on chromosome 18 that conferred resistance to P. sojae isolate 2457 but not to isolate 40412 (Jang et al., 2020b). To further understand this isolate-specific response, we examined the overlap between DEGs and candidate genes within these regions to identify genes that responded uniquely to P. sojae isolate 2457. Among the 51 genes in the identified region spanning 452 kbp on chromosome 18 (from 2,171,155 to 2,623,988), a gene orthologous to HS1 PRO-1 2 (HSPRO2) of Arabidopsis thaliana was identified in the sets of DEGs enriched for disease response-related GO terms (GO:0000165, GO:0002679, GO:0009611, GO:0009627, GO:0009693, GO:0009695, GO:0009697, GO:0009723, GO:0009753, GO:0009816, GO:0009862, GO:0009863, GO:0009867, GO:0009873, GO:0010200, GO:0010363, GO:0031348, GO:0043069, and GO:0050832).
The transcriptional activity of GmHSPRO2 (Glyma. 18G029300) increased by more than 2-fold at 6 and 12 HAI compared with that at 0 HAI in SC2 (Fig. 7). However, in DP, the slightly increased expression of GmHSPRO2 at 6 HAI was reduced by more than 2-fold compared with that at 12 HAI. These findings are consistent with our observation of the specific upregulation of SA-mediated signaling genes in SC2 during P. sojae infection. Furthermore, the repressed expression at 0 HAI was increased at 12 HAI in SC2 but decreased in DP (Fig. 7), consistent with the genotype-dependent upregulation and downregulation of genes for SA, JA, and ET biosynthetic processes in SC2.
Discussion
Linkage mapping of a trait has been widely employed in uncovering the genetic basis of target phenotypic variations. With large sizes of mapping populations, high resolution of genetic maps, and comprehensive set of SNP, it facilitates more precise identification of genomic loci responsible for the investigated phenotypic variations than in the past. However, pinpointing the causal genes remains a considerable challenge. To address this limitation, the integration of linkage analysis with transcriptomic analyses offers a valuable strategy for identifying genetic determinants of phenotypic variations of interest (Almeida-Silva and Venancio, 2021; Azam et al., 2023; Wei et al., 2021).
Previously, a genomic region associated with resistance against P. sojae infection was identified on chromosome 18 in the DP × SC2 RIL population (Jang et al., 2020b). As a follow-up, transcriptome analyses on the two parental cultivars following inoculation with P. sojae were conducted in the current study. The present comprehensive transcriptomic analyses led us to understand transcriptional reprogramming in responses to P. sojae isolate 2457 between DP and SC2 and to pinpoint a novel candidate gene within the previously identified genomic interval in the current study.
A global comparison of the transcriptomic data of P. sojae-inoculated hypocotyl tissues of SC2 and DP at three different time points showed substantial variation in soybean transcriptome activities. Specifically, PC analysis unveiled two major factors governing variation among RNA-seq data. Considering the metadata, including cultivars and HAI of P. sojae of the samples, we concluded that the variance explained by PC1 and PC2 could be attributed to the commonly conserved and genotype-specific transcriptional responses against P. sojae, respectively. Consequently, we identified eight (Groups 1-8) and six (Groups 9-14) groups of DEGs that were explained by PC1 and PC2, respectively.
The subsequent functional enrichment analysis of each set of DEGs revealed that the defense response-related GO terms, including the JA-mediated signaling pathway, SA-mediated signaling pathway involved in SAR, response to JA, JA and ET biosynthesis, and MAPK, were commonly activated in both cultivars. However, the biological processes were differentially regulated in a genotype- and/or temporal-specific manner. These observations suggest that the improved resistance of SC2 can be attributed to genotype- and/or temporal-specific regulation of defense mechanisms. Defense-related phytohormones, such as SA, JA, and ET, play essential roles in the SAR mechanism (Clarke et al., 2000). The DEG profiling conducted by Lin et al. (2014) elucidates that the resistance against P. sojae in 10 NILs, each carrying a distinct Rps gene, was intricately linked to responses mediated by ET, JA, ROS, and MAPK. Furthermore, treatments with benzothiadiazole, a synthetic analog of SA, and 1-aminocyclopropane-1-carboxylic acid, a precursor of ET, have been demonstrated to enhance SA- and ET-signaling, respectively. This augmentation of signaling pathways contributes to conferring resistance to P. sojae in susceptible soybean cultivar, cv. Enrei (Sugano et al., 2013). In line with the findings of the two preceding reports, our GO enrichment analysis between DP and SC2 revealed a heightened expression of genes associated with SA, JA, ET, and MAPK pathways in SC2. Therefore, we believed that the enhanced resistance of SC2 against P. sojae could be attributed to the augmented JA-, ET-, and SA-mediated SAR. Given that SAR is implicated in the responses to diseases characterized by a broad spectrum of pathogenic infections and hypersensitive reactions (Bektas and Eulgem, 2015; McDowell and Dangl, 2000), the heightened disease resistance observed in SC2 may be influenced by non-specific disease resistance mechanisms.
The previous genetic mapping study identified different alleles governing the resistance in SC2 against two different P. sojae isolates (Jang et al., 2020b). The identified alleles explained 14.7% to 26.5% of phenotypic variation in the progenies of crosses between DP and SC2. This result indicated that the resistance of SC2 to P. sojae isolate 2457 is inherited by race-specific, incomplete resistance alleles, whereas typical Rps gene-mediated complete resistance is qualitatively governed (Niu et al., 2017; You et al., 2023a, 2023b; Zhong et al., 2019).
Among the DEGs associated with defense responses, the candidate gene, GmHSPRO2, was identified within the genomic interval associated with the resistance to P. sojae isolate 2457 (Jang et al., 2020b). The HSPRO2 functions were well conserved as resistance factors against various phytopathogens, including bacteria, fungi, and aphids (Aslam et al., 2022; Bose et al., 2018; Schuck et al., 2012). This gene has been recognized as a positive regulator of basal resistance in A. thaliana in response to P. syringae pv. tomato (Murray et al., 2007). Although the resistance mechanism of HSPRO2 against P. syringae pv. tomato in A. thaliana is well-established, and the contribution of GmHSPRO2 to the resistance mechanism against P. sojae in soybeans has never been reported. In Arabidopsis, the HSPRO2 gene governs the basal resistance mechanism against P. syringae by acting downstream of JA, ET, and SA signaling but does not mediate signaling to SAR. Furthermore, HSPRO2 lacks the nucleotide-binding and kinase domains, which are the major characteristics of R proteins (Murray et al., 2007). These findings suggest that HSPRO2 is not involved in broad-spectrum resistance mechanisms. Considering that the identified resistance locus harboring the GmHSPRO2 was significant against to P. sojae isolate 2457 (not to the isolate 40412) (Jang et al., 2020b), our findings suggest that the GmHSPRO2 may function as a pathotype-specific resistance gene. Consistently, the transcriptional accumulation of GmHSPRO2 mirrored the expression patterns of genes involved in the SA-mediated signaling pathway and ET, JA, and SA biosynthetic processes in the SC2. These results and the previously established molecular connection imply that GmHSPRO2 may be the causal gene in the previously identified pathotype-specific resistance allele. Further experimentation may be required to substantiate the hypothesis regarding HSPRO2-mediated resistance. To complete this, the generation of an HSPRO2 knock-out or knock-down soybean plant is necessary, but currently not available; furthermore, this process requires a considerable amount of time. Consequently, future studies will focus on confirming the function of this gene in the interaction between soybean and P. sojae isolate 2457 via the use of transgenic lines or virus-induced gene silencing.
In Karhoff et al. (2022), transcriptomic analysis of NILs identified a candidate gene (Glyma.18G026900) responsible for partial resistance to P. sojae from the previously identified major QTL on chromosome 18. The candidate gene was downregulated in inoculated resistant NILs compared with susceptible NILs at 3, 24, and 48 HAI, suggesting it may act as a susceptibility factor. In the present study, however, the expression level of Glyma.18G026900 was slightly high in DP at 0 HAI but not significantly different between DP and SC2 at 6 and 12 HAI (Supplementary Fig. 5).
Furthermore, the resistance observed in the present study may differ from that observed by Karhoff et al. (2022). The present study demonstrated that DP and SC2 exhibited susceptible and resistant reactions following hypocotyl inoculation. Consequently, the observed transcriptomic differences between DP and SC2 were attributed to comparing compatible and incompatible interactions. In the study by Karhoff et al. (2022), in contrast, the three parents OX20-8, PI 427106, and PI 427105B of the NILs exhibited susceptible reactions following hypocotyl inoculation, and the lesion lengths following root inoculation differed statistically between NILs with either resistant or susceptible alleles. Consequently, transcriptomic changes were investigated under compatible interactions between the NILs and the P. sojae isolate in Karhoff et al. (2022).
In summary, in this comprehensive RNA-seq analysis, we delved into unraveling the intricate resistance mechanisms used by SC2 against P. sojae infections. Our findings reveal a sophisticated orchestration of boosted resistance mechanisms within SC2, primarily driven by the nuanced transcriptional regulation of non-pathogen-specific SAR and the pathotype-specific resistance gene, GmHSPRO2. Further studies are needed for functional verification of GmHSPRO2 in the resistance mechanism against P. sojae is required for functional validation of GmHSPRO2. This research provides a glimpse and a profound understanding, offering invaluable insights that can be harnessed for developing soybean cultivars with robust resistance to P. sojae. Our investigation will help fortify agricultural landscapes against the pervasive threat of P. sojae, fostering sustainable and resilient soybean production in Korea and worldwide.
PDF Links
PubReader
ePub Link
Full text via DOI
Full text via PMC
Download Citation
Supplement1
Print
