RNA-Seq and Comparative Transcriptomic Analyses of Asian Soybean Rust Resistant and Susceptible Soybean Genotypes Provide Insights into Identifying Disease Resistance Genes

Asian soybean rust (ASR), caused by Phakopsora pachyrhizi, is one of the most destructive foliar diseases that affect soybeans. Developing resistant cultivars is the most cost-effective, environmentally friendly, and easy strategy for controlling the disease. However, the current understanding of the mechanisms underlying soybean resistance to P. pachyrhizi remains limited, which poses a significant challenge in devising effective control strategies. In this study, comparative transcriptomic profiling using one resistant genotype and one susceptible genotype was performed under infected and control conditions to understand the regulatory network operating between soybean and P. pachyrhizi. RNA-Seq analysis identified a total of 6540 differentially expressed genes (DEGs), which were shared by all four genotypes. The DEGs are involved in defense responses, stress responses, stimulus responses, flavonoid metabolism, and biosynthesis after infection with P. pachyrhizi. A total of 25,377 genes were divided into 33 modules using weighted gene co-expression network analysis (WGCNA). Two modules were significantly associated with pathogen defense. The DEGs were mainly enriched in RNA processing, plant-type hypersensitive response, negative regulation of cell growth, and a programmed cell death process. In conclusion, these results will provide an important resource for mining resistant genes to P. pachyrhizi infection and valuable resources to potentially pyramid quantitative resistance loci for improving soybean germplasm.


Introduction
The challenge of maintaining and increasing the global food supply is a significant concern for human civilization in the 21st century.Maximizing agricultural crop yields has always been a formidable task due to the presence of pathogens.Soybean (Glycine max), being one of the most important crops globally, is extensively cultivated and serves as the primary legume crop that supplies 70% of the world's protein meal as well as food oil and is a sustainable fuel source (www.soystats.com,accessed on 1 May 2023).The occurrence of Asian soybean rust (ASR), caused by the obligate biotrophic fungus Phakopsora pachyrhizi, leads to significant yield losses of up to 90% in all soybean-producing countries [1].ASR exhibits a wide host range and has the ability to infect over 150 plants belonging to 53 different genera, and soybean is the best host [2].The primary pathogenic forms of rust bacteria are summer spores.The dissemination of soybean rust is facilitated by wind and characterized by strong fulmination, rendering control efforts challenging.This disease manifests severely in regions with high temperatures and humidity and has emerged as a primary affliction of soybean crops in tropical and subtropical locales [3].In South America, soybean rust is estimated to affect approximately 80% of the total cultivated area of soybean crops, with its prevalence being ubiquitous in Latin American regions, where an estimated 210 million metric tons of soybeans are expected to be produced by 2022/23, as per incomplete statistics (https://apps.fas.usda.gov/psdonline/app/index.html, accessed on 14 March 2023), and on average, represent a gross production value of USD 115 billion per season (https://www.ers.usda.gov/dataproducts/season-average-price-forecasts.aspx,accessed on 14 March 2023).Chemical control is the only method by which to control the disease due to the limited resistance available in soybean germplasm.In Brazil, the estimated cost for controlling P. pachyrhizi is approximately USD 2.8 billion per season (http://news.agropages.com/News/NewsDetail-33194.htm,accessed on 14 March 2023) [4].To address this issue, the most cost-effective and efficient approach is to select and breed superior varieties with a high yield and stable and persistent resistance, which can effectively control ASR in soybeans and other crops [5].
Over 30,000 soybean resources have been identified for their resistance to rust [6][7][8][9][10][11].However, only eight genetic loci associated with resistance have been identified thus far, and a limited number of 20 resistance resources with distinct genetic characteristics have been discovered.Furthermore, it has been observed that existing ASR isolates have the ability to overcome the resistance provided by six loci [12,13].Eight major resistance loci that confer resistance to Phakopsora pachyrhizi are Rpp1 [14], Rpp2 [15], Rpp3 [16], Rpp4 [17], Rpp5 [18], Rpp6 [19], Rpp7 [20] and Rpp6907 from the resistant material SX6907 in our lab [21].A limited number of ASR-resistant genes have been successfully cloned.Although certain cultivars can acquire rust resistance through backcross transfer, conventional breeding techniques heavily depend on the accessibility of disease-resistant resources and are hindered by pathogen evolution, posing significant challenges to substantial advancements.Given the limited genetic diversity and absence of resistant germplasm, addressing this lethal disease has emerged as a formidable predicament for scientific researchers worldwide.
Plants demonstrate their immune response through a series of signaling pathways that encompass recognition, signal transduction, and downstream defense reactions, such as the production of antibacterial compounds [22].RNA-seq technologies are instrumental in investigating the intricate biological processes at a transcriptional level [23], as they have significantly enhanced gene discovery efficiency while reducing time, labor, and cost.Recently, RNA-Seq was employed to elucidate the transcriptional resistance mechanisms to plant pathogens (e.g., [24][25][26]; however, up to now, the response of soybean to P. pachyrhizi still has only a few reports.Panthee et al. [27] discovered notable variations in gene expression during different growth stages of susceptible soybean crops in response to P. pachyrhizi by utilizing whole-genome Affymetrix microarrays.Tremblay et al. [28] employed mRNA-Seq to scrutinize the pathogen's transcript abundance at specific time points, namely 15 s, 7 h, 48 h, and 10 d after the inoculation of susceptible soybean leaves with P. pachyrhizi, with the objective of obtaining fresh insights into the transcript abundance in soybean and the pathogen, including the uredinial stage, during the infection.Hossain et al. [1] have provided an elucidation of the variation in Rpp3 gene-mediated resistance to P. pachyrhizi by focusing on the major genes of a phenylpropanoid pathway through RNA-seq.Therefore, it is crucial to concentrate on the primary infection stages for the examination of the defense responses.A comparative analysis of transcriptional profiling of soybean varieties with contrasting disease resistance will offer a comprehensive understanding of the response of resistant varieties to P. pachyrhizi infection and enhance the breeding of resistant varieties for more effective control of P. pachyrhizi.
In a prior investigation, numerous Chinese soybean germplasm resources were assessed to ascertain their level of resistance to P. pachyrhizi.Of these resources, only SX6907 exhibited immunity to local mixed P. pachyrhizi isolates [7].The current study endeavors to enhance comprehension of the molecular mechanism underlying the resistance conferred by SX6907 against P. pachyrhizi.To achieve this objective, we are conducting a comparative analysis of the transcriptome of resistant and susceptible P. pachyrhizi accessions to elucidate the molecular foundation of resistance and susceptibility.The primary objective of this study was to identify genes that are potentially implicated in conferring resistance against P. pachyrhizi.The findings of this investigation are expected to furnish the molecular underpinnings for the development of resistant cultivars and the formulation of disease management tactics while also augmenting the current knowledge on mechanisms of disease resistance.

The Phenotypic Response of Soybean Accessions to P. pachyrhizi
Leaves of SX6907 and tianlong1 seedlings were inoculated with P. pachyrhizi isolate SS4.At 10 dpi (days post-inoculation) of P. pachyrhizi, it was observed that SX6907 displayed high levels of resistance and no visible lesions, while numerous lesions spread on those of tianlong1 (Figure 1A).According to the results of the seedling incubation assay, SX6907 showed a strong immune response against isolate SS4, whereas tianlong1 was completely susceptible.
in conferring resistance against P. pachyrhizi.The findings of this investigation are expected to furnish the molecular underpinnings for the development of resistant cultivars and the formulation of disease management tactics while also augmenting the current knowledge on mechanisms of disease resistance.

The Phenotypic Response of Soybean Accessions to P. pachyrhizi
Leaves of SX6907 and tianlong1 seedlings were inoculated with P. pachyrhizi isolate SS4.At 10 dpi (days post-inoculation) of P. pachyrhizi, it was observed that SX6907 displayed high levels of resistance and no visible lesions, while numerous lesions spread on those of tianlong1 (Figure 1A).According to the results of the seedling incubation assay, SX6907 showed a strong immune response against isolate SS4, whereas tianlong1 was completely susceptible.
An examination of P. pachyrhizi growth on SX6907, a resistant cultivar, and tianlong1, a susceptible cultivar, demonstrated a nearly identical developmental process at 24 h postinoculation.The formation of P. pachyrhizi structures within leaf tissues, including germ tube (Gt), appressorium (Ap), primary invasive hypha (pih), penetration hypha (ph), invasive hyphae (ih), and sporogenous hyphae (sph), occurred at this time point in both SX6907 and tianlong1.However, significant differences in P. pachyrhizi development were observed between SX6907 and tianlong1 from 3 to 10 days post-inoculation.At 10 days post-inoculation, SX6907 exhibited substantially larger fungal colonies and a greater abundance of pathogen-feeding structures (Figure 1B).In contrast, the fungal development in SX6907 was very limited.An examination of P. pachyrhizi growth on SX6907, a resistant cultivar, and tianlong1, a susceptible cultivar, demonstrated a nearly identical developmental process at 24 h post-inoculation.The formation of P. pachyrhizi structures within leaf tissues, including germ tube (Gt), appressorium (Ap), primary invasive hypha (pih), penetration hypha (ph), invasive hyphae (ih), and sporogenous hyphae (sph), occurred at this time point in both SX6907 and tianlong1.However, significant differences in P. pachyrhizi development were observed between SX6907 and tianlong1 from 3 to 10 days post-inoculation.At 10 days postinoculation, SX6907 exhibited substantially larger fungal colonies and a greater abundance of pathogen-feeding structures (Figure 1B).In contrast, the fungal development in SX6907 was very limited.

Analysis of RNA-Seq Data
In order to comprehend the molecular mechanisms of resistance, we collected three replicates of P. pachyrhizi-treated and controlled Pp-R and Pp-S samples at 6 hpi, 24 hpi, and 10 d to capture the transcriptional changes during the biotrophic and transition phases.
After adaptor removal and quality filtering, the results of the RNA-seq analysis generated a total of 147.2934Gbp clean bases of data (983.9058M reads) from the eight transcriptome samples, including resistant control (R-0 hpi), resistant treatment 6 h (R-6 hpi), resistant treatment 24 h (R-24 hpi), resistant treatment 10 d (R-10 dpi), susceptible control (S-0 hpi), susceptible treatment 6 h (S-6 hpi), susceptible treatment 24 h (S-24 hpi), and susceptible treatment 10 d (S-10 dpi).A rigorous standard was employed in the data cleaning process, whereby all sequences with a Phred score of 20 were eliminated, resulting in a set of clean reads ranging from 37 to 45 million across all samples (Table 1).The clean reads from all 24 libraries were subsequently aligned to the soybean genome Wm82.a4.v1, with an overall mapping coverage of 70-93% observed across all samples.The percentage of uniquely mapped reads ranged from 69% to 90% for all samples.Subsequent analysis was conducted solely on the uniquely mapped reads, with the comprehensive sequencing and mapping statistics presented in Table 1.The transcriptome data of all samples qualified for further analysis to identify the differentially expressed genes (DEGs).Genes with reads per kilobase million (RPKM) values exceeding 1 in three replicates for each group were considered to be expressed.The treated and control groups of both genotypes at three time points collectively expressed more than 30,000 genes, with a varying number of specifically expressed genes in each group (Supplementary Table S1).In total, 258,567 expressed genes were identified across all eight groups, with 132,190 genes expressed in the Pp-R samples and 126,377 genes expressed in the Pp-S samples.

Validation of Selected Differentially Expressed Genes (DEGs) Using qRT-PCR
In order to verify the differential expression analysis based on RNA-seq data, the transcriptional levels of the selected DEGs were evaluated via qRT-PCR analysis.These DEGs were deemed to be potentially involved in defense against P. pachyrhizi and exhibited distinct expression profiles across both the resistant and susceptible samples treated with P. pachyrhizi and the control samples.A total of 10 DEGs were selected for qRT-PCR analysis (Figure 2).The results of the qRT-PCR analysis were in good agreement with the RNA-seq data (Figure 2), thereby validating the reliability of the gene expression levels obtained through the DGE analysis.

Validation of Selected Differentially Expressed Genes (DEGs) Using qRT-PCR
In order to verify the differential expression analysis based on RNA-seq data, the transcriptional levels of the selected DEGs were evaluated via qRT-PCR analysis.These DEGs were deemed to be potentially involved in defense against P. pachyrhizi and exhibited distinct expression profiles across both the resistant and susceptible samples treated with P. pachyrhizi and the control samples.A total of 10 DEGs were selected for qRT-PCR analysis (Figure 2).The results of the qRT-PCR analysis were in good agreement with the RNA-seq data (Figure 2), thereby validating the reliability of the gene expression levels obtained through the DGE analysis.

Comparative Analysis of tianlong1 and SX6907 Transcriptomes
This present study compared the whole genomic gene expression of tianlong1 and SX6907 at four distinct time points.The analysis revealed that the average expressional intensity (measured in fragments per kilo bases per million reads (FPKM) > 1) of whole genomic genes was highest at 24 hpi for SX6907 and at 6 hpi for tianlong1, in comparison to the other time points (Figure 3).
The quantification of the differentially expressed genes (DEGs) between SX6907 and tianlong1, used as the control, was conducted at each time point.The expression patterns of the DEGs were classified into upregulated and downregulated groups.In the 0 h

Comparative Analysis of Tianlong1 and SX6907 Transcriptomes
This present study compared the whole genomic gene expression of tianlong1 and SX6907 at four distinct time points.The analysis revealed that the average expressional intensity (measured in fragments per kilo bases per million reads (FPKM) > 1) of whole genomic genes was highest at 24 hpi for SX6907 and at 6 hpi for tianlong1, in comparison to the other time points (Figure 3).libraries, the count of upregulated genes was lower compared to that of the downregulated genes.However, after 6 h of inoculation, the count of the upregulated genes exceeded that of the downregulated genes.The highest count of DEGs was observed 10 days post-inoculation (Figure 4A).The present study utilized Venn diagrams to compare the differential expression of genes between SX6907 and tianlong1 (used as control) after inoculation of rust fungi.The results indicated that at 10 dpi, 1041 upregulated and 1930 downregulated DEGs were specifically detected, while at 24 hpi, only 70 upregulated and 54 downregulated DEGs were specifically detected.Similarly, at 6 hpi, 806 upregulated and 619 downregulated DEGs were specifically detected.Notably, the least number of differentially expressed genes was detected at 24 hpi.Overall, a total of 74 upregulated and 102 downregulated DEGs were detected at all four time points (Figure 4B).The results indicate that some The quantification of the differentially expressed genes (DEGs) between SX6907 and tianlong1, used as the control, was conducted at each time point.The expression patterns of the DEGs were classified into upregulated and downregulated groups.In the 0 h libraries, the count of upregulated genes was lower compared to that of the downregulated genes.However, after 6 h of inoculation, the count of the upregulated genes exceeded that of the downregulated genes.The highest count of DEGs was observed 10 days postinoculation (Figure 4A).libraries, the count of upregulated genes was lower compared to that of the downregulated genes.However, after 6 h of inoculation, the count of the upregulated genes exceeded that of the downregulated genes.The highest count of DEGs was observed 10 days post-inoculation (Figure 4A).The present study utilized Venn diagrams to compare the differential expression of genes between SX6907 and tianlong1 (used as control) after inoculation of rust fungi.The results indicated that at 10 dpi, 1041 upregulated and 1930 downregulated DEGs were specifically detected, while at 24 hpi, only 70 upregulated and 54 downregulated DEGs were specifically detected.Similarly, at 6 hpi, 806 upregulated and 619 downregulated DEGs were specifically detected.Notably, the least number of differentially expressed genes was detected at 24 hpi.Overall, a total of 74 upregulated and 102 downregulated DEGs were detected at all four time points (Figure 4B).The results indicate that some The present study utilized Venn diagrams to compare the differential expression of genes between SX6907 and tianlong1 (used as control) after inoculation of rust fungi.The results indicated that at 10 dpi, 1041 upregulated and 1930 downregulated DEGs were specifically detected, while at 24 hpi, only 70 upregulated and 54 downregulated DEGs were specifically detected.Similarly, at 6 hpi, 806 upregulated and 619 downregulated DEGs were specifically detected.Notably, the least number of differentially expressed genes was detected at 24 hpi.Overall, a total of 74 upregulated and 102 downregulated DEGs were detected at all four time points (Figure 4B).The results indicate that some DEGs are dependently regulated on the infection time, and some are independently regulated.
A heat map was generated through hierarchical cluster analysis of the expressed genes across all samples, which unveiled distinct expression patterns of certain diseaseresistant genes between the Pp-R and Pp-S accessions.Among the 74 upregulated DEGs across all samples, 56 clustered together to reveal higher expression levels in all four R samples compared to the very low expression in all S samples (Figure 5).Among the rest of the clusters, 18 genes were seen as highly induced in susceptible accessions upon P. pachyrhizi treatment.

FOR PEER REVIEW 7 of 24
DEGs are dependently regulated on the infection time, and some are independently regulated.
A heat map was generated through hierarchical cluster analysis of the expressed genes across all samples, which unveiled distinct expression patterns of certain diseaseresistant genes between the Pp-R and Pp-S accessions.Among the 74 upregulated DEGs across all samples, 56 clustered together to reveal higher expression levels in all four R samples compared to the very low expression in all S samples (Figure 5).Among the rest of the clusters, 18 genes were seen as highly induced in susceptible accessions upon P. pachyrhizi treatment.In the P. pachyrhizi-treated samples, 104 genes were found to be the strongest, upregulated in Pp-R, and comparatively less induced in Pp-S at all 6 hpi, 24 hpi, and 10 In the P. pachyrhizi-treated samples, 104 genes were found to be the strongest, upregulated in Pp-R, and comparatively less induced in Pp-S at all 6 hpi, 24 hpi, and 10 dpi.Some genes have an annotation correlated with disease resistance being listed (Table 2).Fifteen were LRR and NB-ARC domain-containing disease-resistant proteins; two genes were the rust resistance kinase Lr10 isoform.A total of 133 genes were found to be highly downregulated in Pp-R and comparatively less induced in Pp-S at all 6 hpi, 24 hpi, and 10 dpi.A total of 17 were LRR and NB−ARC domain-containing disease-resistant proteins; 7 were the LRR receptor protein kinase; 2 were rust resistance kinase lr10-related protein; 2 were programmed cell death protein; and 2 were the TMV resistance protein N isoform (Table 2).
Table 2. Summary statistics of important differentially expressed genes at all 6 hpi, 24 hpi, and 10 dpi.

Gene Ontology (GO)-Based Analysis of the DEGs
Gene set enrichment analysis was employed to identify the functional gene categories (GO terms) and determine the significant biological processes in SX6907 compared to tianlong1 under P. pachyrhizi infection.The results revealed that defense response was the functional category that was significantly over-represented at all four time points (Table 3, Supplementary Table S2).Additionally, a heat map was generated to display the 30 expression profiles of the core defense genes at all four time points (Figure 6).Notably, 11 NB−ARC domain-containing R genes were expressed at high levels solely in Pp−R in both the control and treated samples.Fourteen R genes containing NB-ARC domains were expressed at a low level, exclusively in Pp-R, in both the control and treated samples.Three receptor-like protein kinases rich in cysteine were found to have inducible expression.The categories of response to stress and response to stimulus were over-represented at 0 hpi, 6 hpi, and 24 hpi.In contrast, the categories of secondary metabolic process, protein phosphorylation, oxidation-reduction process, and response to biotic stimulus were over-represented among the transcripts at 10 dpi.Significant differential induction of these core defense components in Pp−R compared to Pp−S makes them the potential candidates responsible for resistance against P. pachyrhizi in Pp−R accession.

KEGG Pathway Analysis of the DEGs
The Kyoto Encyclopedia for Genes and Genomes (KEGG) enrichment analysis was conducted to identify the main metabolic pathways associated with plant defense in the SX6907 infected by P. pachyrhizi.The DEGs were mapped to canonical pathways in the KEGG database.Based on the number of enriched DEGs, we identified the top 20 pathways (Figure 7).A total of 379 DEGs were mapped with the metabolic pathways related to plant defense at 6 hpi, which is a key time point for the most DEGs detected at this point (Table 4).The plant-pathogen interaction pathway and flavonoid biosynthesis were significantly enriched at all four time points.A total of 207 genes were mapped in the plant-pathogen interaction pathway.A total of 61 genes were mapped in the flavonoid biosynthesis pathway.In addition, the results indicated a large number of DEGs were enriched in the plant hormone signal transduction pathways.NPR1, TGA and PR-1 in the salicylic acid pathway were found to be upregulated in Pp-R and downregulated in Pp−S at 6 hpi.JAZ was found to be upregulated in Pp−R and downregulated in Pp-S at three time points.JAR1 and BAK1 were found to be downregulated in Pp−R and upregulated in Pp−S at three time points (Supplementary Figure S1).

KEGG Pathway Analysis of the DEGs
The Kyoto Encyclopedia for Genes and Genomes (KEGG) enrichment analysis was conducted to identify the main metabolic pathways associated with plant defense in the SX6907 infected by P. pachyrhizi.The DEGs were mapped to canonical pathways in the KEGG database.Based on the number of enriched DEGs, we identified the top 20 pathways (Figure 7).A total of 379 DEGs were mapped with the metabolic pathways related to plant defense at 6 hpi, which is a key time point for the most DEGs detected at this point (Table 4).The plant-pathogen interaction pathway and flavonoid biosynthesis were significantly enriched at all four time points.A total of 207 genes were mapped in the plant-pathogen interaction pathway.A total of 61 genes were mapped in the flavonoid biosynthesis pathway.In addition, the results indicated a large number of DEGs were enriched in the plant hormone signal transduction pathways.NPR1, TGA and PR-1 in the salicylic acid pathway were found to be upregulated in Pp-R and downregulated in Pp−S at 6 hpi.JAZ was found to be upregulated in Pp−R and downregulated in Pp-S at three time points.JAR1 and BAK1 were found to be downregulated in Pp−R and upregulated in Pp−S at three time points (Supplementary Figure S1).

Co-Expression Network Analysis of Soybean Resistance to P. pachyrhizi
In order to enhance the comprehension of the correlation between gene expression and the susceptibility or resistance of leaves to P. pachyrhizi and to ascertain the genes that are exclusively linked to resistance to P. pachyrhizi, the co−expression networks of 24 RNA-Seq datasets comprising leaf samples at four different time points with three replicates were analyzed using weighted gene co−expression network analysis (WGCNA).This analysis resulted in the identification of 24,386 DEGs.A WGCNA model was established using the chosen power values, resulting in the division of 24,386 genes into 33 modules.The module sizes ranged from 30 to 4436, with the exception of the "Grey" module, which comprised genes that could not be assigned to any module and lacked reference significance (Figure 9A).The "blue" and "brown" modules were highly associated with an individual infected stage between SX6907 and tianlong1, and the two modules were selected as the particular interest module for further analysis (Figure 9B, Supplementary Figure S2).A total of 451 genes were clustered in the "blue" module, which were differentially expressed between SX6907 and tianlong1 at 6 hpi and 24 hpi.These were categorized into three functional groups (BP, CC, and MF).The "blue" module in the BP group was mainly enriched in RNA processing (GO:0006396), plant−type hypersensitive response (GO:0009626), host-programmed cell death induced by symbiont (GO:0034050), cell death (GO:0008219), negative regulation of cell growth (GO:0030308), and programmed cell death (GO:0012501).The genes in the MF group were mainly enriched in the structural constituent of the cytoskeleton (GO:0005200) and RNA binding (GO:0003723).The genes in the CC group were significantly enriched in an extrinsic component of the plasma membrane (GO:0019897) and filiform apparatus (GO:0043680) (Supplementary Figure S3A).According to theKEGG pathway analysis, our results demonstrated that these genes were mainly involved in the plant-pathogen interaction (ko04626) and sphingolipid metabolism (ko00600) (Supplementary Figure S3B).A total of 432 genes were clustered in the "brown" module, which was differentially expressed between SX6907 and tianlong1 after infection.The "brown" module in the BP group was mainly enriched in the negative regulation of cell growth (GO:0030308), negative regulation of growth (GO:0045926), regulation of cell growth (GO:0001558), and cell activation involved in immune response (GO:0002263).The genes in the MF group were mainly enriched in 1,3−beta−D−glucan synthase activity (GO:0003843), sterol 24−C−methyltransferase activity (GO:0003838) and ubiquitin modification-dependent histone binding (GO:0061649).The genes in the CC group were significantly enriched in the filiform apparatus (GO:0043680) (Supplementary Figure S3C).The KEGG pathway analysis of our results demonstrated that these genes were mainly involved in the plant-pathogen interaction (ko04626) and sphingolipid metabolism (ko00600) (Supplementary Figure S3D).

Ca 2+ Signaling Pathway
The influx of Ca 2+ from the extracellular sources serves as a primary means of signal transduction, with calcium acting as a crucial second messenger ion in eukaryotic organisms [30].In plants, Ca 2+ signaling plays a central role in both the pattern-triggered immune (PTI) and effector-triggered immune (ETI) responses, eliciting characteristic cytoplasmic Ca 2+ elevations in response to the potential pathogens common to both forms of immunity [31].Our investigation revealed a differential expression of certain genes involved in the regulation of the Ca 2+ signaling pathway in SX6907 as compared to tianlong1 under P. pachyrhizi infection (Supplementary Figure S4 and Table S4), including CNGCs (Glyma.09G168600,etc.), CDPKs (Glyma.02G291300,etc.), and CaM/CML (Glyma.04G194800,etc.).PRR signaling is mediated by the CNGC families [32].The manifestation of a "defence no death" (DND) phenotype has been observed as a result of loss-of-function mutations in CNGC2 or CNGC4 [33].This present study has revealed the upregulation of certain CNGC genes at 6 h post-inoculation in the resistant material.Previous research has demonstrated that the disruption of CaM/CML gene expression or loss of CaM/CML function in mutated plants can significantly impair immune responses [34][35][36].
Our findings indicate that a greater number of CaM/CML genes were upregulated in the resistant material in response to P. pachyrhizi infection (Supplementary Table S4).Recent research has presented persuasive proof of the participation of CDPKs in the majority of immune signaling occurrences [37].Our investigation revealed the downregulation of a specific CDPK gene (Glyma.02G291300) in the resistant material following infection with P. pachyrhizi.This leads us to hypothesize that the calcium-signaling pathway and its interaction network may be crucial in soybean rust resistance and warrant further examination.

MAPK Signaling Pathway
Plants have not developed adaptive immunity mechanisms similar to those found in animals but have instead evolved a complex system with multiple layers to defend against invading pathogens.The initial response involves the recognition of pathogens by cell-surface pattern-recognition receptors (PRRs), which is known as PAMP−triggered immunity (PTI) [38].Upon activation of FLS2 and EFR, the MAPK signaling pathway is triggered, leading to the activation of the defense genes responsible for producing antimicrobial compounds.MAPK cascades play a crucial role in regulating various biological processes, including immunity, in plants [39].The present study utilized RNA−seq analysis to investigate the alterations in gene expression within the MAPK signaling pathway, specifically FLS2, MEKK1, BAK1, EFR, and WRKY, in response to P. pachyrhizi infection in soybean.The recognition of flagellin via FLS2 is widely recognized as a crucial component of plant innate immunity [40].FLS2 (Glyma.03G166300)was found to be upregulated in resistant material, while Glyma.08G083300(LRR receptor-like serine/threonine-protein kinase FLS2) was downregulated following infection with P. pachyrhizi.Following the perception of flagellin by FLS2, the initial signaling event is likely due to the recruitment of BAK1 into the flagellin receptor complex [41].Our investigation revealed that the expression of a BAK1 gene (Glyma.17G178700)was specific to susceptible materials, suggesting its potential involvement in the recognition and regulation of pathogen invasion.In Arabidopsis thaliana, the MEKK1-MKK1/MKK2-MPK4 mitogen-activated protein (MAP) kinase cascade has been shown to suppress cell death and immune responses [42].Notably, we observed upregulation of one MEKK1 gene in the resistant materials following 6 h of infection with P. pachyrhizi.The defense mechanisms of plants against pathogen attacks necessitate significant transcriptional reprogramming, which is under the regulation of WRKY transcription factors.Previous research has established the involvement of OsWRKY80, a positive regulator, in rice's defense response against the sheath blight pathogen R. solani [43].Our findings indicate that six WRKY transcription factors were upregulated in resistant material at 6 h post-infection with P. pachyrhizi, followed by a downregulation at 10 days post−infection.Additionally, seven WRKY transcription factors were exclusively upregulated in resistant material at 6 h post−infection.At 6 h post−inoculation, only two WRKY transcription factors exhibited upregulation exclusively in susceptible material.These transcription factors are believed to have significant involvement in the process of PTI and effector−triggered immunity during the interaction between soybean and P. pachyrhizi.These discoveries offer novel insights into the potential molecular mechanisms underlying soybean resistance to Asian soybean rust.

NLR Genes
ETI is the secondary response to pathogen invasion.Pathogens have the capacity to suppress PTI by injecting effector proteins into plant cells through secretion systems [44].To counteract the effect of effector proteins, plants have developed specific resistance mechanisms mediated by "R" genes [45].These genes encode five distinct types of R proteins that recognize effectors from various pathogens, with the largest type encoding the NLR protein.
NLRs play a crucial role in the recognition of pathogen effectors and the initiation of immune responses [46].This present study revealed that 45 NLR genes exhibited differential expression in the plant-pathogen interaction pathway between infected and non−infected resistant and susceptible materials (Supplementary Table S5).Specifically, 28 NLR genes were upregulated in the resistant materials, while 17 NLR genes were upregulated in the susceptible materials.Additionally, 12 NLRs were exclusively expressed in the resistant materials.Upon analysis of the identified resistance locus for Phakopsora pachyrhizi, it was observed that the Rpp1 locus region contained a cluster of NLR genes, whose protein sequences exhibited homology with BED finger NLR resistance protein in black cottonwood and putative disease resistance protein RPM1 in castor bean [47].Additionally, the Rpp3 region in Hyuuga was found to comprise 7 NLR genes [48].Furthermore, the NLR resistance genes in the Rpp4 region displayed a high degree of homology with the downy mildew resistance genes found in lettuce [10].Therefore, a lot of evidence showed that SBR resistance genes in soybean might be found in the NLR gene family; these differentially expressed NLR genes in SX6907 may be candidates for resistance to Asian soybean rust.

Flavonoid Biosynthesis and Disease Resistance
The flavonoid, a secondary metabolite, plays a crucial role in various plant functions, particularly in its ability to combat pathogens, which has garnered increasing attention [49].In tomato roots, flavonoids contribute to the systemic defense response against the gray mold induced by oligogalacturonide, while the isoflavones are essential in R gene−mediated resistance to soybean blight [50].The downregulation of the isoflavone synthetase gene via RNA interference in soybean roots results in a significant reduction of isoflavone accumulation by 95%, rendering the plant more susceptible to pathogens [51].In this investigation, the differential expression of 38 genes involved in the flavonoid biosynthesis pathway was examined in response to P. pachyrhizi infection in both resistant and susceptible materials (Supplementary Table S6).Notably, the upregulation of one specific dihydroflavonol 4−reductase gene (Glyma.17G173200)was observed in susceptible materials.Dihydroflavonol 4−reductase (DFR) is a crucial enzyme in the flavonoid pathway, responsible for the synthesis of anthocyanins, catechins, and procyanidins [52].Previous research has demonstrated that the overexpression of the DFR gene can enhance anthracnose resistance in tobacco and promote the accumulation of flavonoids in resistant leaves [53].This finding diverges somewhat from our own results, prompting further investigation into the relationship between DFR and plant disease resistance.Isoflavone reductase (IFR) is an enzyme that plays a role in the biosynthesis of isoflavonoids in plants.Previous research has demonstrated that the constitutive expression of GmIFR enhances resistance against P. sojae [54], while the induction of a rice isoflavone reductase-like gene, OsIRL, is triggered by rice blast fungal elicitor [55].Our study revealed that four isoflavone reductase genes (Glyma.11G070500,Glyma.11G070600,Glyma.06G317000, and Glyma.11G070200)exhibited an upregulated expression in the resistant materials.These genes have been implicated in the resistance to soybean rust.Caffeoyl−CoA 3-O-methyltransferase (CCoAOMT) is a key enzyme involved in the synthesis of lignin in plants [56].Previous research has demonstrated that correcting multiallelic mutations in a caffeoyl-CoA O-methyltransferase gene led to increased accumulation of resistance metabolites in Russet Burbank potato [57].Additionally, a ZmCCoAOMT2 gene has been found to confer quantitative resistance to southern leaf blight and gray leaf spot [58].In our study, we observed the upregulation of three CCoAOMT2 genes in the resistant materials.The functional significance of the flavonoid biosynthesis and disease-resistant genes identified here needs further investigation.

Plant Genotypes, Inoculations, and Experimental Conditions
SX6907 is a highly resistant germplasm and has an immune response against P. pachyrhizi isolate SS4 [7].The high-yield variety Tianlong 1, developed by the Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences (OCRI), is susceptible to the P. pachyrhizi isolate SS4.Soybean seeds were grown in the greenhouse at 24-26 • C and with a photoperiod of 18/6 h.After 14 days, the fully expanded primary leaves of the seedlings were collected and used for inoculation of the P. pachyrhizi isolate SS4.The inoculation method was carried out as described by Chen et al. [21]: The urediniospores were routinely multiplied on Tianlong 1 leaves and collected in tubes.Tween 20 (0.01% v/v) was used to adjust the urediniospore suspension to 10 5 urediniospores per ml, followed by placing the detached leaves on a wet filter paper-padded plate.The upper surface of the blade was covered with a layer of filter paper.Five to six leaves were placed on each plate.Each leaf was inoculated with four drops of urediniospore suspension at 5 µL per drop.On the first night after inoculation, the leaves were stored at 24 • C and 70% relative humidity under a 12/12 h photoperiod in a growth chamber.Each day, 1-2 mL of water was added to keep the filter paper moist (Figure 10).Then, the soybean leaves were harvested at four growth stages after infection: 0 h, 6 h, 24 h, and 10 d, respectively.Three biological replicates were collected for each sample.
v/v) was used to adjust the urediniospore suspension to 10 urediniospores per ml, followed by placing the detached leaves on a wet filter paper-padded plate.The upper surface of the blade was covered with a layer of filter paper.Five to six leaves were placed on each plate.Each leaf was inoculated with four drops of urediniospore suspension at 5 μL per drop.On the first night after inoculation, the leaves were stored at 24 °C and 70% relative humidity under a 12/12 h photoperiod in a growth chamber.Each day, 1-2 mL of water was added to keep the filter paper moist (Figure 10).Then, the soybean leaves were harvested at four growth stages after infection: 0 h, 6 h, 24 h, and 10 d, respectively.Three biological replicates were collected for each sample.

Observation by Laser Confocal Microscopy
A total of three inoculated soybean leaves were harvested at 6 h, 24 h, and 10 d after inoculation for histopathological analysis.Soybean leaf segments of 2-3 cm 2 were cut from the center of each inoculated leaf.A solution of ethanol/trichloromethane (3:1, v/v) containing 0.15% (w/v) trichloroacetic acid was used to fix and decolorize leaf sections for three to five days.A saturated chloral hydrate solution was used to clear the specimens until the leaf tissues were translucent.For microscopy, cleared leaf segments were stained with 0.1% wheat germ agglutinin (WGA) with Alexa Fluor488 (Thermo Fisher Scientific, Waltham, MA, USA) for 5-10 min, washed with distilled water, and viewed under differential interference contrast optics and UV light.Microscopy examinations were completed with a Zeiss confocal microscope.

Observation by Laser Confocal Microscopy
A total of three inoculated soybean leaves were harvested at 6 h, 24 h, and 10 d after inoculation for histopathological analysis.Soybean leaf segments of 2-3 cm 2 were cut from the center of each inoculated leaf.A solution of ethanol/trichloromethane (3:1, v/v) containing 0.15% (w/v) trichloroacetic acid was used to fix and decolorize leaf sections for three to five days.A saturated chloral hydrate solution was used to clear the specimens until the leaf tissues were translucent.For microscopy, cleared leaf segments were stained with 0.1% wheat germ agglutinin (WGA) with Alexa Fluor488 (Thermo Fisher Scientific, Waltham, MA, USA) for 5-10 min, washed with distilled water, and viewed under differential interference contrast optics and UV light.Microscopy examinations were completed with a Zeiss confocal microscope.

Library Preparation and Illumina Hiseq Sequencing
With 1 g of total RNA, the RNA-seq transcriptome libraries were prepared using Illumina's TruSeqTM RNA sample preparation kit (San Diego, CA, USA).After completing a polyA selection by oligo(dT) beads, the messenger RNA was isolated and fragmented with a fragmentation buffer.Illumina's protocol was followed for the cDNA synthesis, end repair, A-base addition, and ligation of the Illumina-indexed adapters.The cDNA libraries were then size-selected for 200-300 bp target fragments on the 2% Low Range Ultra Agarose, followed by 15 PCR cycles using the Phusion DNA polymerase (NEB).Paired-end libraries were quantified by TBS380 and then sequenced by the Illumina NovaSeq 6000 (150 bp × 2, Shanghai BIOZERON Co., Ltd., Shanghai, China).

Differential Expression Analysis and Functional Enrichment
We calculated the expression levels of each gene using the fragments per kilobase of exon per million mapped reads (FRKM) method to identify the DEGs between the two samples.The R statistical package edgeR (Empirical Analysis of Digital Gene Expression in R, http://www.bioconductor.org/packages/release/bioc/html/edgeR.html/,accessed on 20 March 2020) was used for the differential expression analysis.The DEGs between two samples were selected according to the following criteria: logarithmic fold changes greater than 2 and a false discovery rate (FDR) of less than 0.05.Goatools (https://github.com/tanghaibao/Goatools, accessed on 20 March 2020) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do,accessed on 20 March 2020) were used to analyze the GO functional enrichment and KEGG pathway analysis of the differentially expressed genes.When the Bonferroni-corrected p-value for the DEGs was less than 0.05, they were significantly enriched in GO terms and metabolic pathways.

WGCNA Construct
For the WGCNA, Pearson's correlation coefficient absolute value is used to measure the gene co-expression in the network.This absolute value is then raised to a power to generate the adjacency matrix.Based on the adjacency matrix, the topological overlap distance is then clustered based on the average linkage hierarchical clustering.Our modules were defined using the cutreeDynamic function with a minimum module size of 30 genes.The mergeCloseModules function was used to merge the modules with eigengene distances of less than 0.25.Biologically significant edges were determined for each module by empirically setting the weight value to 0.10.The biological significance of each module was explored through correlation analyses between the module eigengenes and measured agronomic traits.As a rule of thumb, a gene whose edges comprise over 15% of a module is considered the hub gene.

Conclusions
This study presents the first large-scale comparative transcriptomic profiling of resistant and susceptible soybean genotypes in response to the invasion of P. pachyrhizi.The results revealed a complex and massive gene network response, providing insight into the mechanisms directing resistance to P. pachyrhizi in soybeans.The analysis suggests that the NB-ARC domain-containing disease-resistant genes and cysteine-rich receptor-like protein kinase genes play an essential role in understanding pathogen invasion through a downstream resistance mechanism.P. pachyrhizi infection was associated with diverse plant defense response TF families, such as ERFs, WRKY, MYB, and bHLH.Furthermore, the genes involved in the Ca 2+ signaling pathway, MAPK signaling pathway, and flavonoid biosynthesis were upregulated in the resistant materials, thereby inducing systemic resistance.Moreover, the co-expression network analysis (WGCNA) identified that two modules were mainly enriched in RNA processing, plant-type hypersensitive responses, negative regulation of cell growth, and programmed cell death.In summary, this study provided an important resource to mine disease-resistant genes and will presumably provide a new theoretical basis for the optimization and innovation of genetic breeding for soybean crops.

Figure 2 .
Figure 2. The relative gene expression of 10 randomly selected genes examined by quantitative realtime PCR and RNA-seq.

Figure 2 .
Figure 2. The relative gene expression of 10 randomly selected genes examined by quantitative real-time PCR and RNA-seq.

Figure 3 .
Figure 3. Whole genomic gene expression for SX6907 and tianlong1 inoculation with P. pachyrhizi at 0 h, 6 h, 24 h, and 10 d.The y-axis represents the gene number.The blue box means fragments per kilo bases per million reads (FPKM) > 30, the orange box means FPKM 5-30, and the gray box means FPKM 1-5.

Figure 4 .
Figure 4. Number of differentially expressed genes in soybeans infected by Asian soybean rust.(A) The number of differently expressed genes in SX6907 compared with tianlong1 at four time points (0 h, 6 h, 24 h, and 10 d).The y-axis represents the DEGs number.Number of DEGs was counted with the criteria p < 0.05 and log2 (fold change) > 1. (B) Venn diagram comparison of differentially expressed genes (assigned by p < 0.05 and log2 (fold change) > 1) at four infection stages (0 h, 6 h, 24 h, and 10 dpi) in SX6907 compared with tianlong1.Red highlighted numbers represent the amount of upregulated DEGs, and blue highlighted numbers represent the amount of downregulated DEGs.

Figure 3 .
Figure 3. Whole genomic gene expression for SX6907 and tianlong1 inoculation with P. pachyrhizi at 0 h, 6 h, 24 h, and 10 d.The y-axis represents the gene number.The blue box means fragments per kilo bases per million reads (FPKM) > 30, the orange box means FPKM 5-30, and the gray box means FPKM 1-5.

Figure 3 .
Figure 3. Whole genomic gene expression for SX6907 and tianlong1 inoculation with P. pachyrhizi at 0 h, 6 h, 24 h, and 10 d.The y-axis represents the gene number.The blue box means fragments per kilo bases per million reads (FPKM) > 30, the orange box means FPKM 5-30, and the gray box means FPKM 1-5.

Figure 4 .
Figure 4. Number of differentially expressed genes in soybeans infected by Asian soybean rust.(A) The number of differently expressed genes in SX6907 compared with tianlong1 at four time points (0 h, 6 h, 24 h, and 10 d).The y-axis represents the DEGs number.Number of DEGs was counted with the criteria p < 0.05 and log2 (fold change) > 1. (B) Venn diagram comparison of differentially expressed genes (assigned by p < 0.05 and log2 (fold change) > 1) at four infection stages (0 h, 6 h, 24 h, and 10 dpi) in SX6907 compared with tianlong1.Red highlighted numbers represent the amount of upregulated DEGs, and blue highlighted numbers represent the amount of downregulated DEGs.

Figure 4 .
Figure 4. Number of differentially expressed genes in soybeans infected by Asian soybean rust.(A) The number of differently expressed genes in SX6907 compared with tianlong1 at four time points (0 h, 6 h, 24 h, and 10 d).The y-axis represents the DEGs number.Number of DEGs was counted with the criteria p < 0.05 and log2 (fold change) > 1. (B) Venn diagram comparison of differentially expressed genes (assigned by p < 0.05 and log2 (fold change) > 1) at four infection stages (0 h, 6 h, 24 h, and 10 dpi) in SX6907 compared with tianlong1.Red highlighted numbers represent the amount of upregulated DEGs, and blue highlighted numbers represent the amount of downregulated DEGs.

Figure 5 .
Figure 5. Heat map showing hierarchical cluster analysis of 74 expressed genes across all samples.Gradient scale represents expression levels, with red showing the highest expression to greenwith the lowest expression.

Figure 5 .
Figure 5. Heat map showing hierarchical cluster analysis of 74 expressed genes across all samples.Gradient scale represents expression levels, with red showing the highest expression to greenwith the lowest expression.

Figure 6 .
Figure 6.Heat map showing expression profiles of core defense genes that could be putatively involved in resistance development against P. pachyrhizi in Pp−R.Gradient scale shows Z−scores of DEGs where red represents the most induced expression and green depicts the highest repression.Abbreviations: CRK−cysteine-rich receptor-like protein kinase; GH−indole-3-acetic acid−amido synthetase; OST−organic solute transporter; R gene−NB-ARC domain-containing disease-resistant gene.

Figure 6 .
Figure 6.Heat map showing expression profiles of core defense genes that could be putatively involved in resistance development against P. pachyrhizi in Pp−R.Gradient scale shows Z−scores of DEGs where red represents the most induced expression and green depicts the highest repression.Abbreviations: CRK−cysteine-rich receptor-like protein kinase; GH−indole-3-acetic acid−amido synthetase; OST−organic solute transporter; R gene−NB-ARC domain-containing disease-resistant gene.

Figure 7 .
Figure 7.The KEGG pathway enrichment analysis of DEGs at four time points to P. pachyrhizi infection.

Figure 7 .
Figure 7.The KEGG pathway enrichment analysis of DEGs at four time points to P. pachyrhizi infection.

Figure 8 .
Figure 8. Analysis of differentially expressed transcription factors in four time points.(A).The ratios of differentially expressed transcription factor families.(B).Expression levels are shown for the WRKY, MYB, C2H2, and AP2−EREBP transcription factors in resistant material and susceptible material.FPKM values are represented by color gradient.Too high = red brick.

Figure 8 .
Figure 8. Analysis of differentially expressed transcription factors in four time points.(A).The ratios of differentially expressed transcription factor families.(B).Expression levels are shown for the WRKY, MYB, C2H2, and AP2−EREBP transcription factors in resistant material and susceptible material.FPKM values are represented by color gradient.Too high = red brick.

Figure 9 .
Figure 9. Weighted gene co-expression network analysis (WGCNA) of differentially expressed genes.(A) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure (1−TOM).The upper part of the figure is the gene cluster Tree constructed by dissTOM matrix constructed by weighted correlation coefficients; the lower part of the figure is divided into the distribution of genes in each module, the same color represents the same module; the Dynamic TreeCut color is the module identified by dynamicTreeCut method.(B) Hierarchical clustering analysis of all differentially expressed genes (DEGs) in particular interest modules.Blue color represents downregulation of genes.Red indicates upregulation of genes.

Figure 10 .
Figure 10.A simple diagram of the inoculation process.

Figure 10 .
Figure 10.A simple diagram of the inoculation process.

Table 1 .
Summary statistics of RNA-seq data and mapping results.

Table 3 .
Gene ontology enrichment of DEGs (biological process) in SX6907 compared with tianlong1 under P. pachyrhizi infection at 0 hpi, 6 hpi, 24 hpi, and 10 dpi.Number of DEGs contained in respective categories is given as a p-value.Only categories top 10 and p < 0.05 are displayed.

Table 4 .
DEGs enriched in the pathways related to plant defense at 6 hpi, 24 hpi, and 48 dpi.

Table 4 .
DEGs enriched in the pathways related to plant defense at 6 hpi, 24 hpi, and 48 dpi.