Transcriptome Analysis Identified Gene Regulation Networks in Soybean Leaves Perturbed by the Coronatine Toxin

The non-host specific Pseudomonas syringae phytotoxin Coronatine (COR) causes chlorosis and promotes toxicity by inducing physiological changes in plants. We performed transcriptome analysis to better understand plants' transcriptional and metabolic response to COR. Toward this end, mock-treated and COR-treated soybean plants were analyzed by RNA-Seq. A total of 4,545 genes were differentially expressed between the two treatments, of which 2,170 were up-regulated whereas 2,375 were down-regulated in COR treated samples. Gene annotation and pathway analysis conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases revealed that the differential genes were involved in photosynthesis, jasmonic acid (JA) synthesis, signal transduction, and phenylpropane metabolism. This study will provide new insights into COR mediated responses and extend our understanding of COR function in plants.


INTRODUCTION
Soybean [Glycine max (L.) Merr.] is one of the most important economic crops in the world but is vulnerable to various pests and diseases which drastically affect yield and quality. Soybean speck caused by Pseudomonas syringae pv. glycinea (Pgy) is a common disease, threatening soybean yield worldwide (Keen and Buzzell, 1991;Wrather et al., 1997). Pathogenic bacteria inject type three effectors (TTEs) through the type three secretion system (TTSS) in order to inhibit the plant immune response and cause virulence in compatible hosts. Besides effector proteins, pathogens also secrete phytotoxins to increase pathogenicity on plants. Phytotoxins promote pathogen growth by modifying host cell processes and by altering host metabolism (Galán, 2009). Coronatine (COR) is a non-host specific phytotoxin secreted by Pseudomonas syringae pathovars atropurpurea, glycinea, maculicola, morsprunorum, and tomato (Bender et al., 1999). Coronatine can be detected in healthy tissues adjacent to the bacterial lesions, indicating that COR could move systemically (Zhao et al., 2001). The main symptom caused by COR is chlorosis, however, COR does not destroy organelles or membrane systems (Palmer and Bender, 1995). In addition to causing chlorosis in plants, the compound can reopen stomata, which can promote bacterial invasion, proliferation, and development of plant disease symptoms (Geng et al., 2012).
Coronatine (C 18 H 25 NO 4 ) has a relative molecular weight of 319 g/mol and is a structural and functional analog of JA-Ile. It is formed by linking the α-amino acid Coronamic Acid (CMA) with the polyketone Coronafacic Acid (CFA) through an amide linkage (Brooks et al., 2004). Coronafacic Acid structure partially mimics the plant hormone jasmonic acid (JAs), and similar to Jasmonate, its activity is mediated by binding to the Coranatine Insensitive1 (COI1) receptor (Feys et al., 1994). The affinity of COR for COI1 is more than 1,000 times that of JA-Ile (Katsir et al., 2008). Coronatine induces JA biosynthesis and inhibits SA signaling, thereby dampening SA-dependent disease resistance (Zheng et al., 2012). The treatment of COR in plants also induces ethylene production and the release rate of ethylene is proportional to the concentration of COR in tobacco leaves (Ferguson and Mitchell, 1985;Kenyon and Turner, 1992). Although CMA is structurally similar to the ethylene precursor 1-aminocyclopropene 1-carboxylic acid (ACC), it does not induce the production of ethylene responsive genes (Uppalapati et al., 2005). Furthermore, treatment with CFA, CMA, MeJA, or ACC does not induce chlorosis, indicating that COR has additional functions than those attributed to either CFA or CMA (Uppalapati et al., 2005).
High-throughput RNA-Seq is a powerful tool to study gene expression changes in plants . In order to better understand the function of COR in soybean, we treated soybean plant leaves with COR which was followed by transcriptome analysis. The biological annotation of differentially expressed genes (DEGs) under treatment of COR provided new insights into COR mediated plant responses on soybeans, thereby extending our understanding of COR function in plants.

Plant Material and Treatment
Soybean [Glycine max (L.) Merrill cv. Williams 82] was grown under a 12 h light /12 h dark period with light and dark temperatures of 25 and 20 • C, respectively. Four-week-old soybean plants were used in the study. Coronatine purchased from Sigma-Aldrich (St. Louis, MO) was dissolved in methanol to make a stock solution (600 µM, stored at −20 • C) which was subsequently diluted with water to a final concentration of 3 µM. After the addition of 0.004% Silwet L-77, the working solution was sprayed on soybean leaves and leaf samples were collected 10 h post spray. For control treatment, only water containing 0.004% Silwet L-77 solution was sprayed onto the plant leaves. Harvested samples were immediately frozen in liquid nitrogen and stored at −80 • C.

Library Construction and RNA-Seq
Each treatment (COR or mock) was carried out in triplicate. Total RNA was extracted from leaves of COR-treated and control plants using Trizol reagent (Invitrogen, Carlsbad, CA, USA; Rio et al., 2010), and the concentration of RNA was quantified using a NanoDrop-2000 nucleic acid spectrophotometer (Thermo Fisher Scientific, Wilmington, DE). Nucleic acid integrity was assessed by resolving RNA on a 1% (w/v) agarose gel. Quality control and library construction were entrusted to Huada Gene Technology Company (Wuhan, Hubei, China). In brief, from a pool of total RNA, the mRNA fraction was enriched with the help of Oligo-dT magnetic beads (Vazyme Biotech Company, Jiangsu, China). This was followed by RNA fragmentation and first strand cDNA synthesis using the random hexamer primers (Illumina, San Diego, CA, USA). Following second strand cDNA synthesis, the fragments were ligated with adapters and PCR amplified to construct cDNA libraries. BGISEQ-500 platform from Huada Gene Technology Co., Ltd. (Wuhan, Hubei, China) was used to generate raw library reads. The base quaility of the raw data was the ASCII (American standard code for information interchange) value of each character in the fourth line of the FASTQ format minus 64. Reads where more than half of the component bases had a quality score below 15 was considered a low-quality read and removed by Soapnuk. Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT) was used to align the filtered data to the reference genome (GCF_000004515.4_Glycine_max_v2.0; Kim et al., 2015). The RNA-seq data has been uploaded to NCBI gene expression omnibus server (accession number: PRJNA690873, Biosample ID: SAMN17266875, COR SRA ID: SRR13389957, and Mock SRA ID: SRR13389956).

Gene Annotation and Classification
Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) was used to evaluate gene expression levels, and DEGseq R package (Anders and Huber, 2010) was used to determine the DEGs between the two samples. The screening threshold for differential genes was | log 2 Fold Change | ≥ 1, FDR < 0.01, FPKM ≥ 1.
Based on Gene Ontology (GO) (http://www.geneontology. org/) analysis, DEGs were divided into three major functional categories including molecular function, cellular component, and biological process. These genes were further classified

Quantitative Real Time PCR (qRT-PCR) Validation
In order to assess the reliability of the data obtained from RNA-Seq, we measured the relative fold change of five genes by qRT-PCR analysis (Yu et al., 2014). RNA was extracted from frozen leaf tissue using Trizol reagent. In order to determine RNA integrity, the RNA samples were resolved on a 1% (w/v) agarose gel. For first strand cDNA synthesis, the PrimeScript TM RT (TaKaRa, Mountain View, CA), and the gDNA Eraser (TaKaRa, Mountain View, CA) kits were used. For qPCR analysis, the TB Green TM Premix Ex Taq TM (TaKaRa, Mountain View, CA) kit was used under the following set of conditions. A pre-denaturation step of 95 • C for 30 s, 40 cycles with each cycle employing a denaturation temperature of 95 • C for 5 s, and an annealing/extension step with a temperature of 60 • C for 30 s followed by melt curve analysis. For each treatment, we analyzed a total of three biological replicates, each of which included three technical repeats. The 2 − CT method was used for the relative quantification of the five genes. The primers used for the qRT-PCR analysis are listed in Table 1. The GLYMA_20G141600 gene that encodes for polyubiquitin served as an internal control.

Overview of the Soybean Transcriptome
Transcriptome sequencing was performed on samples of COR and mock treated soybean leaves. A high proportion of clean reads (>96%) were obtained after filtration of the raw reads and more than 90% of the total reads had a quality value of 30, suggesting the sequencing reads were of high quality. Hierarchical Indexing for Spliced Alignment of Transcripts FIGURE 1 | Volcano plot for DEGs. The X-axis represents the log 2 -transformed expression changes between COR-treated and mock samples, and the Y-axis represents the -log10-transformed p-value denoting the significance of differential expression. The red dots, blue dots, and gray dots denote up-regulated, down-regulated, and non-differential genes, respectively.
(HISAT) was used to align the filtered data to the reference genome, and the mapping result statistics are shown in Table 2 ( Kim et al., 2015). Gene expression based on fragments per kilobase per million mapped reads (FPKM), generated sequence information for 41,290 and 40,505 reference genes from the COR and mock samples, respectively. Most of the expressed genes were distributed between 10 and 100 FPKM ( Table 3). Transcriptome sequencing analysis was carried out on soybean leaves from the COR and the mock groups. The overall transcriptome data analysis showed that the transcriptome sequencing data was of high quality and could meet the basic requirements necessary for downstream analysis.

Analysis of Differentially Expressed Genes
Our results showed that the COR-vs.-mock comparison yielded 4,531 DEGs. In order to summarize the results obtained from our analysis, and to make visualization of the regulatory patterns more intuitive, we have drawn a volcanic map (Figure 1). The volcano diagram clearly and intuitively shows the expression profile of differential genes in soybeans under different treatments. Although COR-treatment resulted in a greater number of down-regulated genes (2,370), most of the genes with large fold differences were up-regulated (2,161; Figure 1).

Gene Ontology Enrichment Analysis
The GO database was used to functionally annotate DEGs into the following three major categories: Biological process, Cellular component, and Molecular function. There were 4,545 significant DEGs between COR and mock treated samples. Twelve GO subclasses belonged to biological processes sub-category, 12 subclasses were found within the cellular components category while 9 subclasses were enriched in the molecular function category (Figure 2). Within the biological process category, cellular processes, and metabolic processes accounted for a relatively high proportion of hits. In the cell components category, cell, organelle, membrane, and membrane part were the four subclasses with higher percentage representation whereas within the molecular function category, binding, and catalytic activity were the two majorly represented subclasses.

KEGG Annotation of DEGs
The KEGG database integrates gene catalogs obtained from a completely sequenced genome with higher-level system functions at the cell, species, and ecosystem level. The integrative function of the database transforms the expressed gene data into gene networks. By mapping DEGs to the reference genes in the KEGG database, we identified 2,100 DEGs involved in 131 KEGG pathways. Enrichment analysis showed that KEGG enrichment with a Q-value <0.05 was reliable, and on the basis of that threshold, a total of 30 KEGGs were found to be significantly enriched (Figure 3, Supplementary Table 3). The KEGG pathway with the largest number of enriched genes was related to "plant hormone signal transduction" (ko04075), with 181 genes annotated in this pathway. Although there were only 26 genes in the "Photosynthesis-antenna proteins" (ko00196) pathway, the enrichment factor reached 0.764, which was the highest of all annotated pathways. Analysis of the enriched KEGG pathway showed that the vast majority of KEGG genes belonged to the "metabolic branch" and were related to photosynthesis, hormone transduction, and secondary metabolism (Figure 3,  Supplementary Table 3).

Chlorophyll and Photosynthesis Related Genes Were Down-Regulated by COR Treatment
Coronatine treatment resulted in the down-regulation of a large number of genes involved in chloroplast Frontiers in Sustainable Food Systems | www.frontiersin.org  Table 3), indicating that COR had a significant influence on plant health. Chlorophyll a/b binding protein was significantly down-regulated, whereas genes related to magnesium removal, such as magnesium dechelatase, chlorophyllase, uroporphyrin-III C-methyltransferase, pheophorbidase, pheophorbide an oxygenase, and glucuronosyltransferase were significantly upregulated (Supplementary Table 3). Furthermore, COR induced the expression of chlorophyllase, an enzyme that degrades chlorophyll (Takamiya et al., 2000;Xie et al., 2008). In the three photosynthetic pathways: Photosynthesis-antenna proteins (ko00196), Carbon fixation in photosynthetic organisms (ko00710), Photosynthesis (ko00195), there were 129 genes with significant gene expression changes (Supplementary Table 1). Out of these, only 15 genes were up-regulated while 114 genes were down-regulated, indicating that photosynthesis was significantly downregulated upon COR treatment.

Genes Related to ROS Generation Were Upregulated in Response to COR Treatment
When plants are under biological stress, rapid generation of reactive oxygen species (ROS) can be observed (Wojakowska et al., 2013;Gao et al., 2017). While ROS cannot prevent the development of macroscopic symptoms, its production limits the invasiveness of pathogens (Zou et al., 2005). Reactive oxygen species strengthens cell walls through oxidative crosslinking of glycoproteins, and acts as a second messenger in some cell signaling pathways (Lamb and Dixon, 1997). However, ROS also exhibits potential toxicity as its increased accumulation can accelerate membrane lipid peroxidation (Mittler, 2002;Mittler et al., 2004). In order to increase their resistance to oxidative damage, plants require high levels of antioxidant enzymes (Sudhakar et al., 2001). The antioxidant enzymes or ROS scavengers can reduce ROS levels (Keppler and Baker, 1989). Various ROS scavengers, including ascorbic acid peroxidase (PER), glutathione, superoxide dismutase (SOD), and catalase (CAT) maintain ROS homeostasis in different compartments of plant cells (Mittler et al., 2004). Most of the enzymes involved in ROS detoxification, such as glutathione S-transferase (GST), CAT, PER, and glutathione reductase (GR), were up-regulated in COR treated samples, while SOD-related enzymes were downregulated (Table 4).

COR Regulation of JA Biosynthesis and Related Hormone Pathways
Plant hormones not only play an essential role in growth and development, but they also contribute to how plants respond  to both biotic and abiotic stressors. For instance, Jasmonate Acid (JA) modulates developmental as well as innate immune responses. Jasmonate Acid belongs to a small class of lipidderived molecules and acts as a salicylic acid (SA) antagonist and an ethylene (ET) synergist (Gfeller et al., 2010). The first step in JA synthesis is accomplished by lipoxygenase (LOX)-catalyzed oxidation of linolenic acid (Vick and Zimmerman, 1984). In the JA synthesis pathway, 12 genes were up-regulated and only 2 genes were down-regulated. Among the genes encoding LOX, two genes were down-regulated, while five genes were up-regulated. In particular, the genes encoding allene oxide synthase (AOS) and allene oxide cyclase (AOC) were all upregulated, suggesting that COR induced the activation of the JA biosynthetic pathway (Supplementary Table 2). Coronatine treatment also affected the expression of ethylene (ET) synthesis genes (Supplementary Table 2). In this study, a total of 14 DEGs involved in ethylene synthesis were identified, of which 13 genes were up-regulated and 1 gene was down-regulated (Supplementary Table 2), indicating that COR could promote ethylene biosynthesis.

COR Induced the Expression of Genes Involved in the Production of Stress Associated Proteins and Plant Secondary Metabolites
The accumulation of secondary metabolites constitutes an important component of plant stress response (Dixon et al., 2002). Some secondary metabolites play important defense roles during pathogen invasion (Erb and Kliebenstein, 2020). Vegetable Storage Protein (VSP) can be rapidly synthesized or degraded and plays an important role in the process of plant stress adaptation (Staswick et al., 1991). In our study VSP was up-regulated upon COR treatment (Table 4). Glucosinolates (GS) are also defense-related secondary metabolites, and related genes were also up-regulated in the presence of COR (Table 4). Coronatine also induced the expression of key enzymes in the phenylpropane pathway such as phenylalanineammonialyase (PAL), chalcone synthase (CHS), and phenylalanineammonialyase (PAL) ( Table 4,  Supplementary Table 3). These results suggest that COR may affect plant defense responses by regulating plant secondary metabolism, which is known to play important roles during plant pathogenesis.

qRT-PCR Validation of RNA-Seq Results
A total of 5 DEGs (GLYMA_18G057900, GLYMA_05G039900, GLYMA_17G023100, GLYMA_04G166900, and GLYMA_18G072200) identified through RNA-Seq analysis were selected for qRT-PCR verification. Our comparison results (Figure 4) showed that the relative expression levels of these genes in RNA-seq and qRT-PCR followed the same trend, suggesting that the results of RNA-seq were reliable.

DISCUSSION
Bacterial COR has been shown to induce the activation of LOXs, which catalyze the synthesis of JA and MeJA. Through a FIGURE 4 | qRT-PCR validation of DEGs identified by RNA-seq. Five differentially expressed genes (GLYMA_18G057900, GLYMA_05G039900, GLYMA_17G023100, GLYMA_04G166900, and GLYMA_18G072200) identified through RNA-Seq analysis were selected for qRT-PCR validation, and the gene expression changes between COR-treated and mock samples detected by RNA-seq and qRT-PCR were compared. The white columns represent the mean fold change of the five genes through qRT-PCR analysis (error bars represent standard error of mean) while the gray bars represent fold changes obtained through RNA-Seq analysis. Pearson correlation coefficient used to measure the strength of the correlation between both methods showed a significant positive correlation (r = 0.974, p = 0.0053).
positive feedback loop, JA can also activate LOXs and promote JA accumulation even further (Bell and Mullet, 1991). Both COR and JA induce ethylene biosynthesis (Czapski and Saniewski, 1992;Uppalapati et al., 2005) which plays an important role in the development of chlorotic symptoms associated with soybean speck disease (Lund et al., 1998). During plantpathogen interactions, SA signaling activates resistance against biotrophic and hemi-biotrophic pathogens, whereas JA and ethylene (ET) activate resistance against necrotrophic pathogens while suppressing SA dependent signaling (Glazebrook, 2005). As bacteria expressing COR induce the production of JA and ET, we would expect a dampening of SA-dependent responses and enhanced growth of biotrophic pathogens. Coronatine-induced ethylene has also been shown to accelerate plant senescence and inhibit photosynthesis in a dose-dependent manner (Ueda and Kato, 1980;Stall and Hall, 1985;Kenyon and Turner, 1990;Bleecker and Kende, 2000;He et al., 2002;Uppalapati et al., 2005). Consistent with previous studies, we showed that in COR-treated soybean leaves, photosynthetic related genes were down-regulated, suggesting that COR is involved in inducing chlorotic symptoms through the inhibition of the photosynthetic machinery.
Reactive oxygen species generation is one of the initial responses to pathogen invasion (Lamb and Dixon, 1997). High levels of ROS are produced in plant cells undergoing stress. Superoxide anion (O 2− ), hydrogen peroxide (H 2 O 2 ), and hydroxyl radical (OH − ) are the main forms of ROS produced during photosynthesis (Karpinski et al., 2003;Apel and Hirt, 2004). Pseudomonas syringae infection is reported to result in COR-induced decomposition of chlorophyll, which in turn promotes ROS production (Mur et al., 2010). As high level of ROS leads to cell toxicity (Mittler, 2002), its overproduction usually results in the generation of free radical scavenging enzymes. GST is one such enzyme, which is known to be involved in regulating redox potential of infected tissues (Zhang et al., 2015). Consistent with these findings, COR has been reported to cause GST mRNA accumulation 4-9 h after bacterial inoculation (Greulichi et al., 1995). However, not all antioxidant enzymes are up-regulated upon ROS generation. As some studies show that COR inhibits the expression of thylakoid-localized Cu/Zn SOD (Ishiga et al., 2009). These results may explain why we observed a down-regulation of SOD expression and an upregulation of GSTs upon COR treatment.
In Arabidopsis, COR also induces the expression of genes related to phenylpropane metabolism (Attaran et al., 2014). Flavonoids are secondary metabolites that are synthesized via the phenylpropane pathway and are involved in imparting resistance against necrotrophic pathogens. Lamb and Dixon (1997) proposed that the isoflavone phytoalexin pathway may also be the main participant involved in the modulation of oxidative stress responses (Lamb and Dixon, 1997;Winkel-Shirley, 2001). The up-regulation of PAL and CHS indicated that COR targets the phenylpropane metabolic pathway and promotes the synthesis of flavonoids which in turn induce the production of ROS-scavenging enzymes. Of note, ROS scavengers may be beneficial to the host as well as the pathogen. To increase the resistance to oxidative damage induced upon bacterial infection and to improve plant health, plants need high levels of these antioxidants. However, COR can inhibit the early onset of pathogen-associated molecular patterns (PAMP)triggered immunity (PTI) and induce disease-related necrotic cell death at later stages of disease progression. While COR has not been shown to function as a ROS scavenger, tomato plants infected with COR-producing bacterial strains have reduced ROS levels through the activity of ROS scavenging enzymes (Ishiga et al., 2009). It is therefore not surprising that most ROS scavenger-related genes were up-regulated in our study.
In addition to ROS scavengers, we focused on other genes that were specifically regulated by COR. These included the vegetative storage proteins (VSP) and glucosinolate (GS). Vegetative storage proteins is known to be upregulated during insect herbivory and is specifically induced by JA (Benedetti et al., 1995). Other plant hormones below a certain concentration threshold do not modulate VSP levels in plant cells (Anderson et al., 1989). In Arabidopsis, VSP expression was not detected in the Arabidopsis coi1 mutant seedlings, and VSP in wild-type Arabidopsis seedlings was induced upon MeJA and COR treatments (Benedetti et al., 1995). Thioglucoside and its hydrolyzed metabolites are involved in multiple functions including antioxidant activity, chemical protection, and resistance to biotic and abiotic stressors (Guo et al., 2013). The ability of coi1 mutants to induce indole-GS through MeJA was severely impaired, suggesting that COI1 plays a positive role in inducing GS production (Mewis et al., 2006;Brader et al., 2007). It is worth noting that GS, VSP, and ethylene expression can be induced by both COR and JA, indicating that these genes may be induced by COR and JA in a similar manner. Interestingly, the ability of COR to induce chlorosis is also achieved through COI1 (Mecey et al., 2011). The above results indicate that COI1 plays a major role in the regulation of COR-induced responses in plants.
The phytotoxin COR plays an important role in promoting the pathogenicity of bacteria, and the virulent activity of bacteria lacking COR is significantly reduced (Geng et al., 2014). Our results showed that COR induced the expression of genes related to JA and ethylene synthesis, suggesting that COR may activate these two hormone pathways to regulate plant defenses. Coronatine also induced expression of Phenylalanine ammonia lyase (PAL), which is the key enzyme involved in the production of phenylpropane metabolites in plants. This indicated that COR could dampen plant defenses through perturbing phenylpropane metabolism. Our comprehensive transcriptional analysis of the interaction between COR and soybean through RNA-Seq would help better understand the non-host resistance (NHR) response in plants and would also provide a theoretical basis for the improvement of soybean resistance against nonadapted pathogens.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
XZ, BH, and XG designed the research strategy and analyzed the sequenced data. SS, ZZ, and TL conducted part of experiments. XZ, AJA, and XG wrote the manuscript. HW, ZL analyzed part of data. All authors contributed to the article and approved the submitted version.