G. max cv. Kwangan was used for this study. Soybean seeds were sterilized with chlorine gas (100 ml of 5.25% NaOCl supplemented with 3 ml of 12 N HCl) in a fume-hood overnight and sowed on soils (Nongwoo Bio, Suwon, Korea). Soybean plants grew under long-day conditions (25°C, 16 h light/8 h dark, 120 μmol/m²/s) in a walk-in growth chamber.

The second trifoliate leaves of 3-week-old soybean plants were inoculated with PssB728a. We sprayed PssB728a suspension (OD₆₀₀ = 0.01) in 10 mM MgSO₄ containing 0.01% Silwet L-77 to the abaxial and adaxial side of leaves and kept the infected plants in the dark conditions with high humidity during the indicated periods. At least eight leaf discs were taken from infected leaves, and the number of bacteria in the leaf disc was counted using the typical serial dilution method. To analyze the epiphytic growth of PssB728a in plants, leaf discs of infected leaves with PssB728a were washed thoroughly with 10 mM MgSO_4, in which bacterial growth was counted. The washed leaves were then soaked in 70% ethanol during several seconds for surface sterilization and were ground in 10 mM MgSO₄ to measure the endophytic growth of PssB728a (Lee et al., 2012).

Trifoliate leaf tissues of soybean plants (~1 g) were harvested on the 5 days after spray-inoculation with 10 mM MgSO₄ (2 of the control group) and PssB728a (OD₆₀₀ = 0.01) (2 of treatment group) for the extraction of total RNA. RNA purification using RNeasy spin column (Qiagen, Hilden, Germany) was then conducted following the manufacturer’s instruction. After confirming the RNA quality (Bioanalyzer, Agilent, Palo Alto, CA, USA), the total RNA was subject to construct mRNA sequencing libraries using a TruSeq Stranded Total RNA LT Sample Prep Kit (Plant) according to the manufacturer’s instruction (Illumina, San Diego, CA, USA) (Martin and Wang, 2011). Each library from two biological replicates per treatment (total 4 libraries) was sequenced for 101 bp paired-end reads with a HiSeq 2500 platform (Illumina) serviced by Macrogen in Seoul, Korea.

The RNA-seq data were processed by TopHat2 and Bowtie2 (Langmead and Salzberg, 2012; Trapnell et al., 2009). After quality control of raw sequence reads, provided by Phred quality score and FastQC, adaptor sequences were trimmed using the Trimmomatic program to generate trimmed data (Bolger et al., 2014). The trimmed reads were aligned and assembled with the G. max v2.1 gene model (https://www.ncbi.nlm.nih.gov/assembly/GCA_000004515.4/) using the HISAT2 and StringTie (Kim et al., 2019; Pertea et al., 2015, 2016). The raw read counts, FPKM (fragments per kilobase of transcript per million mapping reads), TPM (transcripts per million), and simple description of the mapped transcripts to known soybean genes were presented in Supplementary Table 1.

Followed by the normalization by trimmed mean of M-values in edgeR, differentially expressed genes (DEGs) in infected plants were determined by two statistical outputs (|fold change (FC)| ≥ 2 and exactTest raw. P < 0.05) compared with the mock-inoculated plants (Chen et al., 2016). FC (infected samples/non-infected samples) and related statistical data of each DEG were listed in Supplementary Table 2.

The DEGs were used for gene ontology (GO) enrichment analysis using Protein Analysis Through Evolutionary Relationships (PANTHER) (Mi et al., 2019; Thomas et al., 2003) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Kanehisa et al., 2007). The GO terms and KEGG pathways with P < 0.01 and Q < 0.01 were considered to be significantly enriched in DEGs. The list of enriched GO terms and KEGG pathways related to each DEG were shown in Supplementary Tables 3 and 4. We can share detailed analysis procedures under a request.

PssB728a exploits two different habitats (leaf surface and apoplast) as a phytopathogenic bacterium. The epiphytic PssB728a population acts as an inoculum to infect host plants successfully (Hirano and Upper, 2000). The epiphytic bacterial growth was counted every 2 days after a spray-inoculation of PssB728a. The number of bacteria on the leaf surface increased 2 days post inoculation (dpi), and the bacterial growth was saturated on the 4th day and 6th day (4 and 6 dpi) with populations of 10⁶ colony forming unit (cfu) per leaf disc (Fig. 1A, left). Consistently, the endophytic bacterial population reached > 10⁸ cfu per leaf disc at 4 dpi (Fig. 1A, right), and the lesion increased in size as time progressed (Fig. 1B). The bacterial growth assay reveals that PssB728a causes a compatible interaction with soybean cv. Kwangan.

Transcriptome sequencing libraries were made with total RNA from soybean leaves collected on day 5 after inoculation with 10 mM MgSO₄ (non-infected sample, NI) and PssB728a (infected sample, IF). Among 60,506 known soybean transcripts, 31,608 genes, except 28,898 not detected in at least one sample, were used in the transcriptome analysis. The correlation matrix of Pearson’s coefficient among four samples (NI_1, NI_2, IF_1, and IF_2) indicated that two independent biological replicates were highly reproducible (Fig. 2A). In addition, multidimensional scaling analysis showed that NIs and IFs were separately clustered, and components I and II reflected a 97.4% and 2.1% difference among samples, respectively (data not shown). These statistical analyses show that the transcriptome-seq data can represent the transcriptional reprogramming occurring in infected soybean leaves with PssB728a at 5 dpi.

In a volcano plot, the up-regulated and down-regulated genes with |FC| ≥ 2 and P < 0.05 were displayed on the right (brown dots) and the left (blue dots) sides, respectively (Fig. 2B). About 8,900 up-regulated genes (14.7% of known soybean transcripts) and 7,780 down-regulated genes (12.9%) were assigned as DEGs in infected leaves with PssB728a (Fig. 2C). Hierarchical clustering of expression levels of DEGs (up to 3,000 of 16,680 genes) indicated that two biological replicates for non-infected plants and infected plants had also been done reproducibly (Fig. 2D). Therefore, we suggest that DEGs’ information obtained from the transcriptome analysis helps us to understand the transcriptional changes at 5 dpi in susceptible soybean plants to PssB728a infection. Whole lists of DEGs were presented in Supplementary Table 2.

The cellular and physiological processes of plants are altered upon pathogen infection. Molecular networks composed of complex interactions between/among genes, proteins, and metabolic compounds are involved especially in plant immune signaling, transcription of immune-related genes, and biosynthesis of phytohormones and secondary metabolites, as the outcome of a plant-microbe interaction (Delplace et al., 2022). A compatible interaction between host and pathogen leads to successful infection and disease (Glazebrook, 2005). Different from an incompatible interaction accompanying a hypersensitive response due to intensified immune responses, the compatible interaction entails interruption of plant immunity by effectors and toxins derived from pathogens, resulting in suppression of immune responses during pathogenesis (Chisholm et al., 2006; Dangl and Jones, 2001; Dodds and Rathjen, 2010; Jones and Dangl, 2006). To analyze which biological processes and molecular functions of DEGs were implicated in the compatible interaction between cv. Kwangan and PssB728a, the GO enrichment analysis was accomplished, which explains the associations between gene products and GO terms using the PANTHER classification system.

First, all up-regulated and down-regulated genes were subjected to GO analysis (GO biological process) to point out biological processes in which DEGs are involved. The top 10 categories exhibiting the lowest P-values were listed with the fold enrichments determined by the ratio of the number of observed genes versus that of expected genes (Fig. 3A and B). The up-regulated DEGs were enriched in ‘small molecule metabolic process (GO:0044281)’, ‘child terms of organonitrogen compound metabolic process (GO:1901564) and cellular amide metabolic process (GO:0043603)’, and ‘generation of precursor metabolites and energy (GO:0006091)’ (Fig. 3A). On the other hand, certain kinds of genes involved in ‘regulation of metabolic process (GO:0019222)’, ‘macromolecule metabolic process (GO:0043170)’, and ‘protein modification process (GO:0006464)’ were significantly down-regulated by PssB728a at 5 dpi (Fig. 3B). In addition, large sets of genes associated with RNA biological process appeared to be influenced by pathogen infection, because the expression of those genes involved in RNA biosynthesis, RNA metabolic process, mRNA metabolic process, and regulation of RNA biosynthetic process decreased upon bacterial infection (Supplementary Table 3). In the same line with it, the number of observed genes related to transcription was significantly less in the infected plants compared with that of expected genes of the soybean genome (Fig. 3C). Supporting this, the expression of 99 genes which are implicated in the metabolism of purines (heterocyclic aromatic organic compound consisting of two rings fused, adenine and guanine here), essential building block components for DNA and RNA, was altered by PssB728a infection (Supplementary Table 4). Second, we also analyzed the molecular function of DEGs’ products (GO molecular function and GO-slim molecular function). Both up- and down-regulated DEGs were significantly enriched in 62 GO terms, including transferase activity, catalytic activity, kinase activity, and small molecule binding (Supplementary Table 3). These results suggest that primary and secondary metabolisms are mainly affected by PssB728a infection, whereas transcriptional activity decreased at 5 dpi.

Numerous genes encoding enzymes such as hydrolases, lyases, transferase, and metabolite interconversion enzymes were expressed highly upon bacterial infection (Fig. 4A). In contrast, a large number of genes encoding kinases, especially non-receptor serine/threonine protein kinases (about more than 150 genes) and protein modifying enzymes such as ubiquitin-protein ligases, expressed lower in infected soybean plants than mock-inoculated plants (Fig. 4B). In fact, KEGG pathway analysis revealed that the high population of DEGs was mostly associated with metabolic pathways (carbonate, secondary metabolite, starch and sucrose, amino sugar metabolism) and protein processing (Supplemental Table 4).

Transcriptional profiles of several crop plants after pathogen infections clarified particular metabolic and signaling pathways tightly related to compatible and incompatible interaction (Cregeen et al., 2015; Gyetvai et al., 2012; Meng et al., 2021; Sonah et al., 2016; van Esse et al., 2009; Wang et al., 2010).

Using 3,775 genes mapped in KEGG pathways (http://www.kegg.jp/kegg/pathway.html) among 16,688 DEGs, gene-set enrichments were analyzed to examine which pathways were activated or suppressed at 5 dpi in infected leaves (Kanehisa et al., 2022). Fig. 5 presents the top 15 pathways with P-values < 7.47E-45 and Q-values < 6.63E-44 as follows; metabolic pathways, biosynthesis of secondary metabolite, carbon metabolism, plant hormone signal transduction, MAPK signaling pathway, plant-pathogen interaction, photosynthesis, protein processing in the endoplasmic reticulum (ER), purine metabolism, pyruvate metabolism, etc. Many genes were mapped simultaneously into different pathways. For example, among 137 genes in the category of glycolysis [7], whose expressions were affected upon bacterial infection, 49 genes (about 36%) participate in the pyruvate metabolism [15]; 50% of genes (49 genes) in pyruvate metabolism work at glycolysis-related pathway as well.

Cellular detoxification of reactive oxygen species (ROS) formed from different metabolic pathways and upon various unfavorable stress conditions is an essential biological process in living organisms to maintain homeostasis (Huang et al., 2019). Upon bacterial infection, expression of the biosynthetic genes of antioxidant compounds appeared to be increased as expected. For instance, the expression of genes encoding enzymes involved in ascorbate biosynthesis, such as mannose-6-phosphate isomerase 1 (EC: 5.3.1.8, Entrez_Gene_ID: 100815644), phosphomannomutase (EC: 5.4.2.8, Entrez_Gene_ID: 732605), and mannose-1-phosphate guanylyltransferase 1-like (EC: 2.7.7.13, Entrez_Gene_ID: 100777120), increased during infection (Supplementary Table 4). In addition, numerous glutathione S-transferases, which function as one of the major detoxifying enzymes by conjugating tripeptide (γ-Glu-Cys-Gly) glutathione to various hydrophobic and electrophilic substrates (Liu et al., 2015a), were actively transcribed in the leaves infected with PssB728a (Supplementary Table 4).

Furthermore, a variety of genes related to an ER stress response, which occurs when the capacity of the ER to fold proteins becomes saturated and leads to activation of stress-adaptive response called unfolded protein response (UPR) (Bao and Howell, 2017; Ruberti et al., 2015), showed the altered transcription patterns under PssB728a infection. Three soybean orthologs of an ER stress sensor, the ER-resident transmembrane protein inositol-requiring enzyme 1 (serine/threonine-protein kinase/endoribonuclease IRE1, EC: 2.7.11.1 or EC: 3.1.26, Entrez_Gene_IDs: 100785904, 100788389, 102661527), reduced their mRNA expression upon bacterial infection (Supplemental Table 4). An ortholog of another ER-resident transmembrane protein ER stress sensor, PERK (eukaryotic translation initiation factor 2-alpha kinase 4/eIF-2-alpha kinase GCN2 isoform X1 (EC:2.7.11.1, Entrez_Gene_ID: 100797996), which has been characterized as an ER stress sensor in animal cells but not identified experimentally in plant cells (Park and Park, 2019), also showed the reduced expression in infected leaves. In contrast, several ER-stress responsive genes and genes whose products participate in the protein secretion process exhibited higher expression upon bacterial infection (Supplemental Table 4). An ER stress-inducing compound, tunicamycin enhanced the susceptibility to P. syringae pv. tomato DC3000 in Arabidopsis while gene expression of PR1 and FLS2-mediated PTI such as callose deposition increased (Chakraborty et al., 2017). ire1 and bzip60 mutants, which are deficient in IRE1, an ER stress sensor activating UPR and its splicing target bZIP-transcription factor, bZIP60, respectively, were more susceptible to P. syringae pv. maculicola (Moreno et al., 2012). Consistently, knock-down tobacco plants of IRE1 and bZIP60 via virus-induced gene silencing were more sensitive to the necrotic fungal pathogen, Alternaria alternata (Xu et al., 2019). The transcriptome profile of soybean plants infected with PssB728a suggests that ER-mediated biological process is significantly affected in a compatible interaction (Saijo et al., 2009).

Among the top 15 pathways (Fig. 5), we analyzed further details about plant hormone signal transduction [5], MAPK signaling pathway-plant [10], and plant-pathogen interaction [11]. Genes participating in auxin- and cytokinin-mediated signaling were up-regulated in soybean leaves under PssB728a infection, whereas ethylene (ET)-, jasmonic acid (JA)-, and salicylic acid (SA)-mediated signaling seemed to be down-regulated in PssB728a-infected leaves (Supplementary Fig. 1). In the hormone-biosynthetic aspects, transcript levels of genes participating in the last half of α-linolenic acid metabolic pathway for JA biosynthesis decreased by more than 2-folds in the infected plant (Supplementary Fig. 2). On the other side, soybean orthologs encoding three critical enzymes for the ET biosynthesis, S-ADENOSYLMETHIONINE SYNTHETASE (SAM), 1-AMINOCYCLOPROPANE-1-CARBOXYLATE SYNTHASE (ACS), and AMINOCYCLOPROPANECARBOXYLATE OXIDASE (ACO) were actively transcribed under PssB728a infection (Supplementary Table 4). Additionally, a putative ortholog of Arabidopsis ISOCHORISMATE SYNTHASE1 (ICS1), ICS was strongly transcribed in the infected leaves (EC:5.4.4.2, Entrez_Gene_ID: 100799025) albeit genes involved in shikimate pathway did not show any trends on transcriptional change under infection (Supplementary Table 4). Note that we do not have any experimental clues on the dynamics of these hormones in the infected plants. To measure hormone levels during infection will be necessary to understand the compatible interaction. Taken together, we assume that PssB728a successfully suppresses plant immune signaling activated by ET, JA, and SA and employs growth-related hormone signaling for their colonization.

MAPK signal cascades consisting of three sequential serine/threonine kinases (MAPKKK → MAPKK → MAPK) transmit environmental and developmental signals upon a range of stresses, including pathogen infection, heat, freezing, ROS, drought, and wounding (Meng and Zhang, 2013). Various hormones such as auxin, abscisic acid (ABA), JA, SA, and ET affect the MAPK signaling pathway (Jagodzik et al., 2018). Upon PssB728a infection, ABA receptors PYRABACTIN RESISTANCE 1 (PYR1)/PYRABACTIN RESISTANCE 1-LIKE (PYL) and its downstream signaling components type 2C protein phosphatases (PP2C), negative regulators of ABA signaling, displayed higher mRNA expression. In contrast, the expressions of MAPK signaling components were reduced, probably leading to a non-efficient and unsuccessful stress-adaptive response with disease increase (Supplementary Fig. 3). MAPK pathway relays signals from membrane-residing immune receptors to downstream components for defense activation in plants (Thulasi Devendrakumar et al., 2018). Most of the MAPKKKs, MAPKKs, and MAPKs detected in this study were down-regulated in the transcript levels on day 5 after PssB728a infection (Supplementary Fig. 3). Besides MAPK signaling cascades, genes encoding extracellular (co)receptors, CYCLIC NUCLEOTIDE GATED CHANNELs, and Ca²⁺-signaling proteins were weakly expressed in the infected leaves of soybean (Supplementary Fig. 4). However, a soybean ortholog of FLAGELLIN-SENSITIVE 2 (FLS2) showed increased mRNA expression in the soybean leaves infected with PssB728a. However, its downstream pathways, MEKK1-MKK4/5-MPK3/6 and MPK4, reduced their transcription. In addition, soybean homologous genes of well-identified defense regulators in Arabidopsis and tomato, such as EDS1, RIN4, SGT1, and Pti1/5/6, displayed reduced mRNA levels in the infected leaves (Supplementary Fig. 4). This KEGG pathway analysis suggests that the plant immune responses of soybean cv. Kwangan are compromised by PssB728a.

In this study, we developed a pathosystem between soybean cv. Kwangan and PssB278a, a causal pathogen of bacterial brown spot disease, and presented the transcriptional changes in soybean leaves on day 5 after PssB728a infection. Unexpectedly, mRNA expression of homologous genes to well-identified immune-related genes in model and other crop plants decreased or was not dramatically altered in soybean leaves on day 5 after infection. In contrast, primary and secondary metabolisms underwent a change in PssB728a-infected conditions. Thus, these analyses represent the transcriptional changes at the late phase of a compatible interaction and reflect the suppression of plant immunity in susceptible plants by virulent pathogen infection. Furthermore, this study identifies large number of candidates for future functional analyses and provides a foundation to understand molecular mechanisms implicated in the compatible interaction of a susceptible Korean soybean cultivar and bacterial brown spot disease.