Comparative Transcriptome Analysis of Resistant and Susceptible Tomato Lines in Response to Infection by Xanthomonas perforans Race T3

Bacterial spot, incited by several Xanthomonas sp., is a serious disease in tomato (Solanum lycopersicum L.). Although genetics of resistance has been widely investigated, the interactions between the pathogen and tomato plants remain unclear. In this study, tanscriptomes of X. perforans race T3 infected tomato lines were compared to those of controls. An average of 7 million reads were generated with approximately 21,526 genes mapped in each sample post-inoculation at 6 h (6 HPI) and 6 days (6 DPI) using RNA-sequencing technology. Overall, the numbers of differentially expressed genes (DEGs) were higher in the resistant tomato line PI 114490 than in the susceptible line OH 88119, and the numbers of DEGs were higher at 6 DPI than at 6 HPI. Fewer genes (78 in PI 114490 and 15 in OH 88119) were up-regulated and most DEGs were down-regulated, suggesting that the inducible defense response might not be fully activated at 6 HPI. Accumulation expression levels of 326 co-up regulated genes in both tomato lines at 6 DPI might be involved in basal defense, while the specific and strongly induced genes at 6 DPI might be correlated with the resistance in PI 114490. Most DEGs were involved in plant hormone signal transduction, plant-pathogen interaction and phenylalanine metabolism, and the genes significantly up-regulated in PI 114490 at 6 DPI were associated with defense response pathways. DEGs containing NBS-LRR domain or defense-related WRKY transcription factors were also identified. The results will provide a valuable resource for understanding the interactions between X. perforans and tomato plants.


INTRODUCTION
Bacterial spot caused by at least four distinct species of Xanthomonas (X. euvesicatoria, X. vesicatoria, X. perforans, and X. gardneri), severely affects marketability of both fresh-market and processing tomato Stall et al., 2009). Based on their virulence on a group of tomato genotypes, the causative agents of bacterial spot can be classified into five physiological races T1-T5 (Jones et al., 2005). It is difficult to manage the disease once the outbreak occurs (Yang et al., 2007;Stall et al., 2009). Identification and characterization of the host resistance as well as regulatory genes can provide insights into understanding the mechanisms of host resistance and the use of host resistance to develop new cultivars with durable resistance.
The resistance to X. perforans race T3 in tomato can be classified into two types, hypersensitive response (HR) and field resistance. HR is usually conditioned by single dominant genes, while field resistance is generally controlled by multiple genes (Yang, 2013). To date, four major dominant genes, Rx4 in Solanum pimpinellifolium accession PI 128216 (Robbins et al., 2009;Pei et al., 2012), Xv3 in an unimproved breeding line Hawaii 7981 (Wang et al., 2011), Rx LA1589 in S. pimpinellifolium accession LA 1589 (Sun et al., 2011a), and RXopJ4 in S. pennellii accession LA 716 (Sharlach et al., 2013), conferring HR to X. perforans race T3 have been identified and mapped. A recent study suggests that Rx4, Xv3, and Rx LA1589 are possibly the same gene (Zhao et al., 2015). Although these genes also confer field resistance to race T3, they need modifiers and may be affected by genetic background (Scott et al., 2001;Robbins et al., 2009;Wang et al., 2011;Sharlach et al., 2013). Therefore, the tomato lines with HR usually show partial resistance to race T3 in the field. On the contrary, the S. lycopersicum var. cerasiforme accession PI 114490 shows no HR but high resistance to race T3 conditioned by at least five quantitative trait loci (QTL) in the field (Sun et al., 2011b(Sun et al., , 2014Scott et al., 2015). In addition to identification and genetic mapping of genes or QTLs, efforts have also been making on discovery of genes participating in a complex molecular network of regulation during the time-course of HR to race T3 in the unimproved tomato line Hawaii 7981 using microarray analysis approach (Gibly et al., 2004;Balaji et al., 2007), and identification of genes differentially expressed in the resistant line PI 114490 and a susceptible line OH 88119 during the timecourse of the race T3 infection using cDNA-AFLP techniques (Du et al., 2014). Comparison of genes differentially expressed in PI 114490 and OH 88119 provides us some information to understand the mechanism of resistance to bacterial spot race T3 in tomato during the process of symptom development at the 3, 4, and 5 days post spray-inoculation stages (Du et al., 2014). However, both microarray and cDNA-AFLP can only identify part of the genes involved in resistance or defense processes due to their low throughput limitation. To fully unravel the mechanisms of field resistance to race T3 of bacterial spot, it is necessary to identify more genes differentially expressed during different infection times in tomato.
With the availability of a high-quality tomato genome sequence and second-generation sequencing, RNA-seq technology has rapidly become a popular tool for genomewide expression profiling, providing the potential to better understanding the comprehensive host-pathogen interactions. Many studies have, therefore, tried to elucidate a nearly complete picture of inducible defense response pathways using RNA-seq analysis (Kim et al., 2011;Li et al., 2012;Gao et al., 2013;Tan et al., 2015). Through comparing the different RNA-seq data about hosts and pathogens interactions, similar or specific sets of genes activated in different plant pathosystems have been discovered. Meanwhile, using the expression pattern regulated by pathogen infection, it is also possible to identify the known defense gene family member involved in the tomato immune response. For example, some plant R genes, such as pepper Bs3 and rice Xa27 belonging to NBS-LRR gene family, are transcriptionally activated by corresponding transcription activator-like effector (TALE) protein (Gu et al., 2005;Röemer et al., 2007). Taking advantage of the TALE-induced expression pattern, a Xanthomonas TAL-effector activated resistance gene Bs4C has been identified using RNA-seq (Strauss et al., 2012). Therefore, RNA-seq technology can be used as a tool to isolate pathogen induced R gene, which could avoid laborious positional cloning and investigate the expression pattern of different defense signaling genes family members as well as identify novel defense related genes.
In the current study, it is first time using RNA-seq to investigate transcript dynamics of the field resistant line PI 114490 and the susceptible line OH 88119 in response to X. perforans race T3 at 6 h and 6 days post-infection. Largescale expression profiling identified both common and specific differentially expressed genes in two tomato lines at different time points after inoculation. In addition, 8 induced genes were validated by quantitative real-time PCR (qRT-PCR) experiments. Our results may provide a subset of potential candidate defenserelated genes in tomato-X. perforans interaction, which will help in better understanding the molecular basis of field resistance to bacterial spot in tomato.

Plant Materials and X. perforans Race T3 Inoculation
Two tomato lines, a cherry tomato line PI 114490 with high level of field resistance and an elite processing breeding line OH 88119 with susceptibility to X. perforans race T3 (Sun et al., 2011b), were used in this study. For inoculation with X. perforans race T3, tomato seedlings were grown in a growth room at 28 • C/25 • C (light/dark), with a 14-h photoperiod. The race T3 strain Xv829 was cultured on yeast, dextrose, calcium carbonate (YDC) agar medium (Lelliot and Stead, 1987) at 28 • C for 48-72 h prior to inoculation. Bacterial cultures were diluted to a concentration of approximately 3 × 10 8 colony forming units per ml using sterile solution containing 10 mM MgSO 4 ·7H 2 O and 0.025% (v/v) Silwet L77. Six-week-old plants were spray-inoculated with the bacterial suspension, and plants sprayed with the sterile solution containing 10 mM MgSO 4 ·7H 2 O and 0.025% (v/v) Silwet L77 were used as controls (mock-treatment). The plants were misted with water twice a day (9:00 a.m. and 5:00 p.m.) from 1 day before inoculation to 6 days after inoculation to increase humidity and prolong leaf wetness for disease development.

RNA Isolation, RNA-seq Library Preparation, and Sequencing
Total RNA was isolated from leaf samples collected at started point (mock-treatment), 6 h-post-inoculation (HPI) and 6 day-post-inoculation (DPI) using the TRIzol reagent (Invitrogen, USA) following to the manufacturer's protocol. The concentration of total RNA was determined by NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, USA), and the RNA integrity value (RIN) was checked using RNA 6000 Pico LabChip of Agilent 2100 Bioanalyzer (Agilent, USA).
Three libraries (mock-treatment, race T3 treatment at 6 HPI and 6 DPI) for each tomato line were constructed for RNA-seq analysis. Two hundred micro-grams total RNA was incubated with 10 units DNase I (Ambion, USA) at 37 • C for 1 h, and then nuclease-free water was added to bring the sample volume to 250 µl. The MicroPoly(A)Purist ™ Kit (Ambion) was used to enrich mRNA according to the manufacture's protocol. The enriched mRNA was dissolved in 100 µl RNA storage solution and the final concentration was determined by NanoDrop 2000 Spectrophotometer. Double-stranded cDNA was synthesized from the mRNA according to the method described previously (Ng et al., 2005) with some modifications. First-strand cDNA was synthesized from 10 µg mRNA with Gsul-oligo dT primers using 1000 units of superscript II reverse transcriptase (Invitrogen). After incubation at 42 • C for 1 h, the 5 ′ -CAP structure of mRNA was oxidized by NalO 4 (Sigma, USA) and ligated to biotin hydrazide, which was used to select complete mRNA/cDNA by binding Dynal M280 beads (Invitrogen). The second-strand cDNA synthesis was subsequently performed using Ex Taq polymerase (TaKaRa, Japan). The polyA and 5 ′ adaptor were removed by GsuI digestion. The double-strand cDNA was ultrasonically fractioned with a cDNA size fractionation column at a range of 300-500 bp, and then purified using Agencourt Ampure beads (Agencourt, USA

Analysis of Illumina Reads and Identification of Differentially Expressed Genes
Raw reads were firstly cleaned by removing adaptor sequences and unreliable low quality sequences (reads with unknown sequences "N"). The filtered clean reads were then mapped to S. lycopersicum reference genome sequence (SL2.40 version) and related annotation was obtained from the International Tomato Annotation Group (ITAG) Release 2.3 Predicted CDS (https://www.sgn.cornell.edu/) using Bowtie2 with default settings (Langmead et al., 2009). Clean reads mapped to reference genome sequence with multiple genes were filtered. Reads per kilo bases per million uniquely mapped reads (RPKM) values per sample for all genes were calculated based on both the total number of reads that mapped to the tomato reference CDS database and corresponding genes length (Mortazavi et al., 2008). Differentially expressed genes (DEG) were identified by DEGseq package using the MA-plot-based method with Random Sampling model (MARS) (Wang et al., 2010). The q-value was adjusted using the method described in Benjamini and Hochberg (1995). "False Discovery Rate (FDR) ≤ 0.001 and the absolute value of Log2fold-change ≥ 2" were used as the threshold to determine the significance of gene expression difference.

Gene Ontology and KEGG Pathway Analysis
Gene ontology analysis for DEGs was performed using GoPipe (Chen et al., 2005) through BLASTP against Swiss-Prot and TrEMBL database. The metabolic pathway was constructed based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database using BBH (bi-directional best hit) method (Kanehisa et al., 2010). The KO number was obtained for each protein and used for constructing metabolic pathways.

RT-PCR and qRT-PCR Validation
To verify the results of RNA-seq analysis, 8 differentially expressed genes (Receptor like kinase, Polyphenol oxidase, Harpin-induced protein-like, LRR receptor-like serine/threonine-protein kinase, 1-aminocyclopropane-1carboxylate oxidase, Calmodulin 2, Universal stress protein family protein, and WRKY transcription factor 16) were selected for semi-quantitative RT-PCR and qRT-PCR analysis. First strand cDNA was synthesized by M-MLV reverse transcriptase (TaKaRa) with the same batch of RNA samples used in RNA-seq analysis. The specificity of primers for each gene (Table S1) was confirmed by analyzing PCR products on agarose gel and melting curve during real-time PCR. RT-PCR and qRT-PCR was conducted according to the previous protocols (Pei et al., 2012;Du et al., 2014).

Analysis of RNA-seq Libraries and Reads Alignment
Approximately 43 million Illumina raw reads were generated from six libraries with an average of 7 million 100 bp raw reads in each individual library. The raw reads were deposited in the National Center for Biotechnology Information Sequence Read Archive under the accession number SRP065597. In the six libraries, 97.9-99.3% of raw reads were clean. Of the clean reads, 46.5-77.3% could be mapped uniquely to one location within the tomato reference genome, 8.4-40.3% could be aligned to highly conserved domain sequences shared by different annotated genes, and 9.6-21.3% failed to map to the tomato genome sequence ( Table 1). The uniquely mapped clean reads matched 21,168-21,993 with an average of 21,526 genes (approximately 62%) to the tomato annotated genes for the six libraries ( Table 1). Approximately 20,000 genes were mapped at both mock-and T3treatment samples (Image S1). The number of specific matched genes was higher in the PT6d library (890) than in other five libraries (Image S1).
The genes RPKM values were calculated only for those uniquely mapped reads. The number of genes with RPKM values of 0 varied from 12,734 to 13,559, with an average of 13,200 in the six libraries (Table S2). Of which 10,117 were common in all libraries, indicating a very low level or no expression of those genes in the six samples. The numbers of genes with RPKM value of 0 were slightly higher in the three libraries of the susceptible line OH 88119 than the corresponding samples in the resistant line PI 114490. Based on the expression analysis using RT-PCR FIGURE 1 | Numbers of differentially expressed genes (DEGs) identified in tomato lines PI 114490 and OH 88119 at 6 h and 6 days post-inoculation of Xanthomonas perforans race T3. The DEGs were identified using FDR ≤ 0.001 and the absolute value of Log 2 fold-change ≥ 2 as the threshold for the significance of genes expression difference. PM and OM: PI 114490 and OH 88119 respectively mock-inoculated with the sterile solution containing 10 mM MgSO 4 ·7H 2 O and 0.025% (v/v) Silwet L77. PT and OT: PI 114490 and OH 88119 respectively inoculated with T3. and qRT-PCR, genes with RPKM value less than 2.28 needed more amount of cDNA and genes with RPLM value greater than 10.7 needed less amount of cDNA for obtaining clear bands (data not shown). Therefore, genes with RPKM values of 0-3 were considered as at low level of expression, while genes with RPKM values beyond 10 were considered to be expressed at high level. The number of highly expressed genes was less in PT6h library compared to OT6h, but more in PT6d than OT6d (Table S2).

Differentially Expressed Genes in
Response to X. perforans Race T3 More stringent criteria with smaller FDR (≤0.001) and higher fold-change value (absolute value of Log 2 fold-change ≥ 2) were used to analyze the global transcriptome dynamics of both resistant line PI 114490 and susceptible line OH 88119 in response to race T3 during 6 h and 6 days post-inoculation. Overall, the number of differentially expressed genes (DEGs) was higher in PI 114490 than in OH 88119 at both infection stages, and the number of DEGs was considerably higher at 6 DPI than at 6 HPI in both tomato lines (Figure 1). At the 6 HPI, there were 591 and 299 genes were differentially expressed in PI 114490 and OH 88119, respectively. Most DEGs showed down-regulations in both tomato lines, whereas only 78 and 15 genes were upregulated in PI 114490 and OH 88119, respectively (Figure 1). However, the numbers of DEGs dramatically increased to 1656 and 1154 at 6 DPI in PI 114490 and OH 88119, respectively, and the numbers of up-regulated genes were obviously greater than that of down-regulated genes (Figure 1). Consideration of numbers and fold changes of these DEGs, resistant tomato line PI 114490 had a higher frequency and stronger fold changes of DEGs at 6 DPI, compared to PI 114490 at 6 HPI and OH 88119 at both time points (Figure 1).

Genes with Differentially Expressed Patterns in Two Tomato Lines during Infection with Race T3
In order to gain comprehensive insights into understanding molecular mechanisms of resistance induced by X. perforans race T3, the overlap differentially expressed patterns of DEGs between PI 114490 and OH 88119 at two infection stages were analyzed using a Venn diagram (Figure 2). For tomato line PI 114490, 282 genes (8 up-regulated and 274 down-regulated) were co-modulated in both 6 HPI and 6 DPI compared to mocktreatment (Figure 2A), while 86 genes (6 up-regulated and 80 down-regulated) were co-modulated in both 6 HPI and 6 DPI compared to mock-treatment in line OH 88119 ( Figure 2B). Comparison between the two tomato lines discovered 51 and 344 DEGs co-regulated at 6 HPI and 6 DPI, respectively (Figures 2C,D). Among the 51 co-regulated DEGs at 6 HPI, 50 were co-down-regulated, and only 1 (Solyc03g114530.2.1) encoding Strictosidine synthase family protein was co-upregulated. The expression of Strictosidine synthase (Str), the key gene of terpenoid indole alkaloid (TIA) biosynthetic pathway, has been found to be induced by fungal elicitor (Pasquali et al., 1992) and salinity stress in Catharanthus roseus (Dutta et al., 2013). Go classification analysis of the 50 down-regulated genes, showed that the dominant subcategories in biological process were response to stimulus (Table S3). Among the 344 DEGs co-regulated at 6 DPI, 326 genes were commonly up-regulated, and 18 genes were commonly down-regulated ( Figure 2D). In the face of a large number of commonly up-regulated genes, particular emphasis was placed on the highly up-regulated genes resulting in 32 DEGs with the Log 2 fold-change ≥ 5 in both lines at 6 DPI ( Table 2), of which a large proportion (40.6%) were assigned to Go term "response to stimulus." For example, marker genes of stress signaling such as pathogen-related genes (PR) and osmotin-like protein were intensely induced, and PR proteins are inducible by infection with various types of pathogens in many plant families (van Loon et al., 2006). In addition to these stimulus genes, other genes also have been clarified related to defense response. Clustering analysis of DEGs among the four comparison groups resulted in 37 clusters (Figure 3). Despite there were substantial number of DEGs, none of these genes was commonly up-regulated at four comparison groups. Cluster 10 was the largest with 711 DEGs, followed by cluster 8 representing 223 common up-regulated genes in PI 114490 and OH 88119 at 6 DPI. Common up-regulated genes in cluster 8 indicated that both PI114490 and OH88119 used some same sets of genes to response to race T3 infection but the induced intensity of transcriptional levels were different between two genotypes. Cluster 10 and cluster 20 consisting of the DEGs were only upregulated (compared to mock-inoculated) in PI 114490 and OH 88119 at 6 DPI response to race T3. These specific DEGs in PI 114490 from cluster 10 ( Table 3) were promising candidates for revealing the resistance response mechanism. Cluster 32 comprised only one unknown function gene Solyc02g084850.2.1 TABLE 2 | RPKM for common highly up-regulated (Log 2 fold-change ≥ 5) genes in tomato lines P I114490 and OH 88119 at 6 days post-inoculation (DPI). with strongly down-regulated in both PI 114490 and OH 88119 at two infection stages ( Table 3). Cluster 27 contained 7 genes strongly down-regulated at 6 HPI but strongly upregulated at 6 DPI (

Gene Ontology Classification Analysis of DEGs
Go enrichment categorized for 52.5, 54.2, 50.5, and 51.2% of DEGs into functional groups from PT6h, OT6h, PT6d, and OT6d, respectively. More assigned function of DEGs was covered in biological process and molecular function categories than in cellular component categories (Image S2). The dominant subcategories of both PI 114490 and OH 88119 at 6 HPI and 6 DPI in each main category were cellular process, cell and catalytic activity, respectively (Image S2). Despite these enrichment subcategories were similar in both tomato lines, the individual genes contributing to the common enriched subcategories were substantial diversified between the two tomato lines.

Metabolic Pathway by KEGG Analysis of DEGs
To characterize the pathway enrichment of the identified DEGs, gene classification was performed on the basis of KEGG analysis (Kanehisa et al., 2010). Only significant pathway categories among four comparisons were selected and listed in Table 4. Genes involved in photosynthesis and oxidative phosphorylation pathways were mainly up-regulated in OH 88119 at 6 DPI, but down-regulated in PI114490 at 6 HPI and 6 DPI ( Table 4). It was obvious that genes associated with defense response pathways (labeled in bold in Table 4) were significantly upregulated in PI114490 at 6 DPI. The defense-related pathways with most representation by DEGs were plant hormone signal transduction, plant-pathogen interaction and phenylalanine metabolism. A total of 16 defense response genes in KEGG pathway "plant-pathogen interaction" (KO: 04626) exhibited significant differences between PI 114490 and OH 88119 in response to T3 infection ( Table 5). Of which, 5 and 3 were downregulated at 6 HPI in PI 114490 and OH 88119, respectively. However, the numbers of DEGs increased in both tomato lines at 6 DPI. A total of 15 (14 up-regulated and 1 down-regulated) in PI 114490 and 8 (6 up-regulated and 2 down-regulated) gens in OH 88119 was detected. It was noteworthy that 5 DEGs (Solyc05g050350.1.1, Solyc10g079420.1.1, Solyc11g071740.1.1, Solyc09g014990.2.1, and Solyc12g009220.1.1) had the similar expression pattern between two tomato lines at both 6 HPI and 6 DPI ( Table 5).

Validation of RNA-seq Data for Selected Genes by RT-PCR and qRT-PCR
To evaluate the validity of RNA-seq analysis, transcriptional levels of selected 8 DEGs representing a wide range of expression levels and patterns were determined by RT-PCR and qRT-PCR analysis. Based on the functional annotation in literatures, majority of these DEGs were associated with massive defense response processes including ethylene biosynthesis (Solyc07g049530.2.1), defenserelated enzymes (Solyc02g072470.2.1, Solyc07g008600.1.1, Solyc01g057000.2.1, Solyc02g036480.1.1, Solyc04g072070.2.1 and Solyc10g081170.1.1), and second metabolite biosynthesis (Solyc08g074630.1.1). In general, the expression data provided by qRT-PCR were almost consistent with profiles detected by RNA-seq at all time-points (Figure 4), confirming the trends of up-or down-regulation of all the analyzed genes. The differences in the magnitudes of changes observed between the qRT-PCR and RNA-seq results might be caused by different algorithms.
It was validated by qRT-PCR that four genes Solyc02g072470.2.1, Solyc02g036480.1.1, Solyc08g074630.1.1, and Solyc07g008600.1.1 were steadily up-regulated during two time-points infection in both PI 114490 and OH 88119, differing in their induced strength. Solyc07g049530.2.1 and Solyc04g072070.2.1 initially displayed down-regulated at 6 HPI but significantly reinforced at 6DPI in both tomato lines, and the up-regulated fold-change was higher in PI 114490. The expression pattern of Solyc01g057000.2.1 was opposite between two tomato lines, which was up-regulated in the resistant line PI 114490 while down-regulated in the susceptible line OH 88119 (Figure 4).

DISSCUSSION
The S. lycopersicum var. cerasiforme accession PI 114490 has been considered as a durable source for resistance to bacterial spot due to its high level of quantitative field resistance to four races T1-T4 (Yang, 2013). There is no or few lesions on PI 114490, but the bacterial population of race T3 in its leaves is not significantly different from that in the susceptible line OH 88119 (Sun et al., 2014). So considering on the aspect of phenotype, PI 114490 shows higher tolerant response to avoid the Our previous study using cDNA-AFLP techniques has already identified DEGs between the two tomato lines (Du et al., 2014). Therefore, we started with one biological replication for RNA-seq analysis. An average of 7 million sequence reads were obtained from RNA-seq and the clean reads could match approximately 62% of reference genes (Table 1), which was closed to previous tomato RNA-seq data with 20 million sequence reads (Tang et al., 2013), indicating that the sequencing depth was sufficient for the transcriptome coverage. The well correlation between qRT-PCR and RNA-seq data (Figure 4) also demonstrated the reliability of RNA-seq analysis. Of course, more replicates should provide more reliable information for understanding the genes involved in the process of resistance/defense to the pathogen of bacterial spot race T3 in tomato.
Inducible Defense Response to X. perforans Interfered at 6 HPI but Activated at 6 DPI in Both Resistant and Susceptible Tomato Lines Despite PI 114490 and OH 88119 have different responses to race T3 infection (Sun et al., 2014), lots of common DEGs were regulated in both two tomato lines. At 6 HPI, 50 common downregulated genes existed in both two tomato lines (Figure 2C) suggested that the inducible defense response pathways were not completely activated at 6 HPI. Meanwhile, down-regulation of the minority of pathogen-related genes in cluster 27 (Figure 3) indicated that at least some defense response had been interfered by pathogen infection as early as 6 HPI.   Among the 344 common DEGs at 6 DPI, a majority of 326 genes were commonly up-regulated ( Figure 2D) in both tomato lines. The dominant subcategories in biological process were cellular process, macromolecule metabolism and response to stimulus (data not shown). Six up-regulated genes in plantpathogen interaction pathway were shared by the two tomato lines at 6 DPI ( Table 5). Previous studies have proven that the function of these regulated genes plays important role in bacterial spot defense response (Cui et al., 2015). Therefore, defense response of tomato lines was activated 6 days after infection by race T3. Accumulation expression levels of those co-regulated genes were likely to increase basal defense response and improve the disease resistance to race T3 infection in both resistant and susceptible tomato. So the susceptible line OH 88119 also succeeded to properly activate the expression of many defense-related genes.

Resistant and Susceptible Tomato Lines Displayed Differentially Expressed Patterns in Photosynthesis Pathway
Plant-pathogen interaction response usually alters the expressional level of genes associated with photosynthesis (Major et al., 2010). Many photosynthesis-related genes including chlorophyII a/b binding proteins and photosynthetic reaction center proteins were strongly repressed in susceptible hybrid poplar leaves at 9 days inoculated with Melampsora medusa (Miranda et al., 2007). However, our RNA-seq data showed that race T3 pathogen caused opposite impact on photosynthesis between resistance and susceptible tomato lines. Genes involved in photosynthesis pathway were mainly up-regulated in the susceptible tomato line OH 88119 at 6 DPI, but down-regulated in the resistant line PI 114490 at 6 HPI and 6 DPI ( Table 4). To the best of our knowledge, the reason of this divergent response between resistant and susceptible plants in photosynthesis pathway has not previously been known in plant-pathogen interactions.

Analysis of Differentially Expressed Genes Involved in Plant Immune Response Pathways
Plants possess different sophisticated defense strategies which are usually accompanied by transcriptional changes to protect themselves against pathogen attacks. The plant immune response pathway can be mainly divided into two branches (Jones and Dangl, 2006). The first branch and initial defense response is usually activated by pathogen-associated molecular patterns (PAMP) whose presence is recognized by plant plasma membrane-localized pattern recognition receptors (PRRs), which results in PAMP-triggered immunity (PTI). The second branch is activated by nucleotide-binding leucine-rich repeat (NB-LRR) resistance gene, which directly or indirectly interacts with virulence factor, and results in activation of effector-triggered immunity (ETI). Since the plant-pathogen interaction pathways (KEGG: 04626) summarized the genes involved in defense network of PTI and ETI, so we focused on the transcript dynamics of this KEGG pathways and enumerated the DEGs imparting the plant immune response ( Table 5).
Several PAMPs have been identified from bacteria. Flagellin is one of important and well-characterized PAMPs that can be recognized by the LRR receptor-like serine/threonine-protein kinase (FLS2) (Gómez-Gómez and Boller, 2000). In our study, RLK protein Solyc02g070890.2.1 with the highest amino acid similarity (54%) to FLS2 protein from Arabidopsis thaliana than other tomato RLK genes was down-regulated at 6 HPI, and slightly up-regulated at 6 DPI in PI 114490, but this gene was not differentially expressed in OH 88119 upon race T3 infection. RLK protein Solyc02g072470.2.1 was the highest induced member of RLK family in our RNA-seq data. qRT-PCR validated that Solyc02g072470.2.1 were significantly induced in both tomato lines at 6 DPI, which displayed more than 35and 18-fold-induction in PI 114490 and OH 88119, respectively (Figure 4). The highest induced member of LRR receptor-like serine/threonine-protein kinase was Solyc07g008600.1.1 in RNAseq analysis, which was verified up-regulated in both two tomato lines (Figure 4).
Recognition of PAMPs initiates downstream signaling pathways involving WRKY transcription factors to confer defense response against bacterial, fungal pathogens and nematodes (Asai et al., 2002;Bhattarai et al., 2010). Based on KEGG analysis, WRKY transcription factor Solyc09g014990.2.1 homologous to WRKY33 was down-regulated at 6 HPI but highly up-regulated at 6 DPI in both genotypes ( Table 5). Since members of WRKY families are thought to play important role in defense response, we further investigated the expression patterns of all 85 WRKY families genes annotated in the International Tomato Annotation Group (ITAG) Release 2.3 Predicted CDS. A total of 29 WRKY genes were differentially expressed. Among them, none was significantly up-regulated at 6 HPI but 17 (58.6%) were up-regulated in both tomato lines at 6 DPI ( Figure 5) (Huang et al., 2012). In the current study, none of these WRKY genes were detected at 6 HPI, but SlWRKY39/53/80/81 were significantly induced in both tomato lines, and SlWRKY8/23 were specific up-regulated in PI 114490 at 6 DPI. SlWRKY72b are transcriptionally up-regulated during disease resistance mediated by the R gene Mi-1, and virus-induced gene silencing of this gene in tomato resulted in a clear reduction of Mi-1-mediated resistance as well as basal defense against root-knot nematodes and potato aphids (Bhattarai et al., 2010). WRKY gene Solyc06g070990.2.1, the highest identity with SlWRKY72b than other WRKY protein in tomato, were significantly induced at 6 DPI in PI 114490 ( Figure 5).
Intracellular Ca 2+ influxes have also long been recognized as essential and early events downstream of multiple PAMP perception, leading to local and systemic acquired resistance (Lecourieux et al., 2005;Boudsocq et al., 2010). In the process of Ca 2+ signaling, some members of Calcium-dependent protein kinase (CDPK) family was immediately induced by flg22 and Avr-Cf9 interaction (Romeis et al., 2000;Boudsocq et al., 2010). Unexpectedly, none of CDPK genes was induced as early as 6 HPI in both tomato lines. Calmodulin (CALM or CaM), a major Ca 2+ sensor, play positive regulatory resistance roles in different host-pathogen interaction (Park et al., 2004;Choi et al., 2009). Only one CALM gene Solyc10g081170.1.1 was specifically up-regulated at 6 DPI in PI 114490 (Table 5), which was also validated by qRT-PCR analysis (Figure 4).
In the process of plant and pathogen interaction, bacterial type III effectors which are directly delivered into the host cells via the type III secretion system can bypass or disrupt host first line of PAMP-triggered immunity. But resistant plants contain R proteins belonging to the nucleotide -binding site-leucinerich repeat (NBS-LRR) family to directly or indirectly detect bacterial effectors, which will activate downstream signaling and lead to pathogen resistance (DeYoung and Innes, 2006). Previous papers have evidence that some R genes can be activated by the specific effectors, which are displayed by the accumulation of R genes expression levels. Therefore, we paid attention to those putative R genes containing NBS-LRR domain that could possibly be induced by the X. perforans infection. A total of 298 genes were annotated as NBS-LRR proteins for tomato in SGN database, and the RPKM values of most genes were low in the six sequenced samples, indicating the expression level of most putative R genes was low in tomato plant. The RNA-seq data showed that none of putative R genes were up-regulated in both tomato lines at 6 HPI. But the number of up-regulated putative R genes in PI 114490 at 6 DPI was more than that of in OH 881119 (Figure 6). Six putative R genes were commonly up-regulated in both tomato lines. Among these genes, Solyc07g056190.2.1 displayed the highest fold-change (Figure 6). In addition, there were 11 and 1 putative R genes were specifically up-regulated in P I114490 and OH88119, respectively (Figure 6).

Up-regulated Genes
Previously, we used cDNA-AFLP technique to identify differential expression genes in tomato lines PI 114490 and OH 88119 in response to race T3 infection at 3, 4, and 5 days. A total of 60 genes were commonly up-regulated at different levels in two tomato lines (Du et al., 2014). Among these genes, 14 genes were also up-regulated in this study (Table S5). Based on the functional annotation in literatures, majority of these genes were associated with massive defense response process including jasmonic acid biosynthesis pathways (Solyc03g122190.2.1 and Solyc04g079730.

CONCLUSION
In this study, high-throughput sequencing data was used to identify X. perforans responsive modules, which were enriched for differentially expressed genes between resistant line PI 114490 and susceptible line OH 88119 at different time point after sprayinoculation of race T3. The data provided a whole view on the variation tendency of genes involved in tomato-X. perforans interaction. A large number of candidate genes involved in tomato-X. perforans interaction were identified. Comparison of RNA-seq data with our previous cDNA-AFLP data revealed common up-regulated genes, and some of them were involved in jasmonic acid and ethylene biosynthesis.

AUTHOR CONTRIBUTIONS
Conceived and designed the experiments: YW, WY. Performed the experiments: HD. Analyzed the data: HD, JY, WY. Contributed reagents/materials/analysis tools: YW, JY, WY. Wrote the paper: HD, WY.  Table S3 | Common down-regulated genes associated with "response to stimulus" according to Gene Ontology classification in tomato lines PI 114490 and OH 88119 at 6 h post-inoculation (HPI). Table S4 | Gene ontology (GO) analysis of up-regulated genes for biological process in tomato lines PI 114490 and OH 88119 at 6 d post-inoculation.
Table S5 | Sets of common differentially expressed genes in RNA-seq analysis with previous cDNA-AFLP results.
Image S1 | Venn diagram showing the overlaps and specific expressed genes between mock samples and samples spray-inoculated with Xanthomonas perforans race T3 in tomato lines PI 114490 and OH 88119. PM and OM: PI 114490 and OH 88119 respective mock-treatment with the sterile solution containing 10 mM MgSO 4 .7H 2 O and 0.025%(v/v) Silwet L77. PT and OT: inoculation with bacterial spot race T3 in PI 114490 and OH 88119, respectively.
Image S2 | Distribution of the differentially expressed genes (DEGs) within Go secondary categories of molecular function, cellular component and biological process.