Transcriptome-Based Identification of Differently Expressed Genes from Xanthomonas oryzae pv. oryzae Strains Exhibiting Different Virulence in Rice Varieties

Xanthomonas oryzae pv. oryzae (Xoo) causes bacterial blight (BB) in rice (Oryza sativa L.). In this study, we investigated the genome-wide transcription patterns of two Xoo strains (KACC10331 and HB1009), which showed different virulence patterns against eight rice cultivars, including IRBB21 (carrying Xa21). In total, 743 genes showed a significant change (p-value < 0.001 in t-tests) in their mRNA expression levels in the HB1009 (K3a race) strain compared with the Xoo KACC10331 strain (K1 race). Among them, four remarkably enriched GO terms, DNA binding, transposition, cellular nitrogen compound metabolic process, and cellular macromolecule metabolic process, were identified in the upregulated genes. In addition, the expression of 44 genes was considerably higher (log2 fold changes > 2) in the HB1009 (K3a race) strain than in the Xoo KACC10331 (K1 race) strain. Furthermore, 13 and 12 genes involved in hypersensitive response and pathogenicity (hrp) and two-component regulatory systems (TCSs), respectively, were upregulated in the HB1009 (K3a race) strain compared with the Xoo KACC10331 (K1 race) strain, which we determined using either quantitative real-time PCR analysis or next-generation RNA sequencing. These results will be helpful to improve our understanding of Xoo and to gain a better insight into the Xoo–rice interactions.


Introduction
Bacterial blight (BB) has become one of the most devastating diseases of rice and an important problem affecting rice production worldwide. Some bactericides have been developed to control the disease, but none is highly effective. Thus, genetic basis breeding for BB resistance has been found to be the most effective method for the control of BB [1]. At least 38 BB resistance (R) genes in rice have been identified [2,3], of which six (Xa2, Xa4, Xa7, Xa30, Xa33, and Xa38) have been physically mapped and eight (Xa1, xa5, xa13, Xa21, Xa26/Xa3, Xa10, Xa23 and Xa27) have been cloned [2][3][4][5][6][7][8][9][10][11]. Among them, the Xa21 gene was identified in the wild species Oryza longistaminata and shown to be highly effective against BB races of South and Southeast Asia [12]. Pattern recognition receptors (PRRs) are proteins of critical importance in detecting invading pathogens in plants and animals.

Figure 2.
Total number of R-seq reads from each Xanthomonas oryzae pv. oryzae strain library and mapped reads to the Xanthomonas oryzae pv. oryzae KACC10331 genome.

RNA-seq Reads Mapping and GO Term Enrichment of Differentially Expressed Genes
To compare the transcriptome profiles of the Xoo KACC10331 (K1 race) and HB1009 (K3a race) strains, RNA sequencing libraries were constructed. The two libraries generated about 16.1 and 19.6 million reads, which were mapped to the X. oryzae pv. oryzae KACC10331 genome sequence in the NCBI GenBank database (NC_006834.1) with 96.9% and 96.4% matched reads to NCBI annotated gene regions ( Figure 2). There were 4065 protein-coding genes (CDS) among 4281 annotated genes in the genome of KACC10331 by genome analysis [24]. A total of 4062 and 3995 protein-coding genes were detected (ě5 mapped reads) in the Xoo KACC10331 (K1 race) and HB1009 (K3a race) strains, respectively ( Figure 2 were detected (≥5 mapped reads) in the Xoo KACC10331 (K1 race) and HB1009 (K3a race) strains, respectively ( Figure 2).

Figure 1.
Virulence assay (A) and lesion length (B) of rice varieties inoculated with two Xanthomonas oryzae pv. oryzae strains (K1 and K3a race). Asterisks represent statistically significant differences relative to a susceptible rice cultivar (IR24) (paired, two-tailed Student's t test, ** p-value < 0.01).  To identify DEGs, we used the R-package DEGseq [25]. Genes with significant differences in expression level between the two strains were evaluated using the criterion of a p-value less than 0.001 in a t-test, resulting in 743 DEGs when the HB1009 (K3a race) strain was compared with the Xoo KACC10331 strain (K1 race) (Table S1). To further analyze these DEGs, the transcripts were categorized according to their annotated function with respect to biological processes, molecular functions, and cellular components, on the basis of the blast and GO term annotation using Blast2GO software [26]. The biological processes mediated by these DEGs were primarily associated with cellular processes, metabolic processes, and single organism processes, among others ( Figure 3, Table S2). The most major molecular functions were binding and catalytic activity. Cell was the most majority in the cellular component category. Most of these categories contained a larger numbers of genes with typically upregulated expression than those of genes with typically downregulated expression ( Figure 3, Table S2). In the GO term enrichment analysis using Blast2GO (Fisher's exact test), four remarkably enriched GO terms, DNA binding, transposition, cellular nitrogen compound metabolic process, and cellular macromolecule metabolic process, were identified in the upregulated genes. However, there were no GO terms enriched in the genes with downregulated expression (Table S3). Thus, the RNA-seq and GO term enrichment analysis results suggested that some of these genes were predominantly expressed in HB1009 (K3a race), since the expression of 44 genes (a total of 147 enriched genes) was considerably higher (log2 fold changes > 2) in the HB1009 (K3a race) strain than in the Xoo KACC10331 (K1 race) strain (Table S4). To identify DEGs, we used the R-package DEGseq [25]. Genes with significant differences in expression level between the two strains were evaluated using the criterion of a p-value less than 0.001 in a t-test, resulting in 743 DEGs when the HB1009 (K3a race) strain was compared with the Xoo KACC10331 strain (K1 race) (Table S1). To further analyze these DEGs, the transcripts were categorized according to their annotated function with respect to biological processes, molecular functions, and cellular components, on the basis of the blast and GO term annotation using Blast2GO software [26]. The biological processes mediated by these DEGs were primarily associated with cellular processes, metabolic processes, and single organism processes, among others ( Figure 3, Table S2). The most major molecular functions were binding and catalytic activity. Cell was the most majority in the cellular component category. Most of these categories contained a larger numbers of genes with typically upregulated expression than those of genes with typically downregulated expression ( Figure 3, Table S2). In the GO term enrichment analysis using Blast2GO (Fisher's exact test), four remarkably enriched GO terms, DNA binding, transposition, cellular nitrogen compound metabolic process, and cellular macromolecule metabolic process, were identified in the upregulated genes. However, there were no GO terms enriched in the genes with downregulated expression (Table S3). Thus, the RNA-seq and GO term enrichment analysis results suggested that some of these genes were predominantly expressed in HB1009 (K3a race), since the expression of 44 genes (a total of 147 enriched genes) was considerably higher (log2 fold changes > 2) in the HB1009 (K3a race) strain than in the Xoo KACC10331 (K1 race) strain (Table S4).

Differential Expression of Two-Component Systems between the Xoo K1 Race and K3a Race Strains
To cope with various environmental conditions, including limited-nutrient niches and various toxic substances produced by the host, bacteria have evolved a mechanism involving two-component regulatory systems (TCSs). TCSs are composed of histidine kinases (HKs) and response regulators (RRs). In response to environmental stimuli, the HK phosphorylates the cognate RR, which then regulates gene expression [32]. In this study, among the 48 TCSs, 15 TCSs showed significantly different expression levels in the HB1009 (K3a race) strain compared with Xoo KACC10331 (K1 race). Twelve TCSs (XOO0336, XOO0519, XOO1105, XOO1207, XOO2227, XOO2228, XOO3659, XOO3709, XOO3710, XOO3871, XOO3875, and XOO4008) and three TCSs (XOO0423, XOO2322, and XOO2798) were upregulated and downregulated, respectively, in the HB1009 (K3a race) strain compared with the Xoo KACC10331 (K1 race) strain ( Figure 5, Table 2, Table S1). In addition, three (XOO1105, XOO2227, and XOO4008) out of 12 significantly upregulated TCSs were over-represented in the GO term enrichment analysis (Table S4).    Quorum sensing (QS) is the process in which bacteria monitor their own population density by sensing the levels of extracellular signal molecules that are released by the microorganisms, and the QS systems might enable pathogens to overcome the host defense mechanisms [33]. RpfC and RpfG serve as a two-component system for the perception and transduction of the extracellular diffusible signal factor (DSF) family signal-mediated QS [34][35][36]. The expression pattern of the RpfC/RpfG two-component regulatory system (XOO2870 and XOO2871) associated with QS was different between the HB1009 and KACC10331 strains (Table S1). Furthermore, Ryan et al. [37] reported that several domain proteins, including GGDEF, EAL, and HD-GYP, are involved in a network of signal transduction systems for responding to different environmental factors to modulate the level of the second messenger cyclic di-GMP in X. campestris pv. campestris. In this study, a TCS regulatory protein with a GGDEF (XOO2787) or with both GGDEF and EAL domains (XOO0520) had upregulated (0.59 and 1.26 log2 FC, respectively) expression in the HB1009 strain. In contrast, a TCS regulatory protein with an HD-GYP domain (XOO2798) had significantly downregulated expression (´2.2 log2 FC) in the HB1009 strain compared with the KACC10331 strain. In addition, XOO2860 (encoding cyclic di-GMP phosphodiesterase A) had upregulated expression in the HB1009 strain compared with the KACC10331 strain (Table S2).
Rice lines carrying Xa21 (encodes a leucine-rich repeat receptor-like kinase) are able to induce an effective defense response to multiple strains of the bacterial Xoo pathogens [38]. In addition, several genes that are required for activation of XA21-mediated immunity (rax) have been identified in Xoo [13][14][15]. Furthermore, Lee et al. [39] and Zhang et al. [31] reported that Xoo requires a regulatory TCS called RaxRH to regulate expression of 10 rax genes, including raxA, raxB, raxC, raxST, raxP, raxQ, raxR, raxR2, raxH, and raxH2. Our results revealed that all of these rax genes were not expressed at significantly different levels between the HB1009 and KACC10331 strains according to the DEGseq analysis (p-value > 0.001, FALSE), except for raxR2 (XOO0423) ( Table 2, Table S2). However, some of these genes, including raxP (XOO3397), raxR (XOO3535), raxB (XOO3543), raxA (XOO3544), and raxST (XOO3546), showed different expression patterns between the two strains in the qRT-PCR analysis ( Figure 5). As described above, the expression of raxR2 (XOO0423) associated with TCSs was significantly downregulated in the HB1009 (K3a race) strain ( Table 2, Table S1), which showed reduced virulence against IRBB21 (carrying Xa-21) (Figure 1). In contrast, raxR (XOO3535) was upregulated in the HB1009 (K3a race) strain, although the raxR gene did not show significantly different expression levels in the DEGseq analysis (p-value > 0.001, FALSE) ( Table 2, Table S2). The findings of the current study are consistent with those of Lee et al. [39] who reported that a response regulator encoded by phoP (raxR2; XOO0423) is upregulated in a raxR gene knockout mutant strain. The PhoPQ-regulated protein (XOO1731) not only is required for AvrXa21 activity but also regulates various cellular activities as a regulator of virulence in Salmonella and other species [40][41][42]. Upregulation of the phoPQ-regulated protein was detected in the HB1009 strain, although this gene did not show significantly different expression levels in the DEGseq analysis (p-value > 0.001, FALSE).

Discussion
In the past several years, a large number of Xoo pathogenesis-related genes have been isolated and characterized, as many researchers have tried to elucidate the mechanisms of the Xoo-rice interaction. However, many aspects of the mechanisms of the Xoo-rice interaction are still not clearly understood. For example, although several genes involved in the activation of XA21-mediated immunity have been identified, a key gene required to elicit rice-resistant protein Xa21 expression has not been identified [43]. As described earlier, control of BB involves the introduction of resistance genes to confer major gene-specific resistance against some pathogen races. However, this type of resistance frequently results in rapid changes in the pathogen diversity, and new races of the pathogen are able to overcome the deployed resistance [44,45].
In this study, we presented the comparison of genome-wide transcriptional patterns between two Xoo strains and used transcriptome profiling to identify genes that might be involved in the different pathotypes and that may be specific to the Xoo race. GO term enrichment analysis identified four remarkably enriched GO terms in the upregulated genes, and the expression of 44 genes (a total of 147 enriched genes) was considerably higher (log2 fold changes > 2) in the HB1009 (K3a race) strain than in the Xoo KACC10331 (K1 race) strain (Table S3). Among the 44 genes, 13 and seven genes were identified as transposase and hypothetical genes, respectively. In addition, six genes, identified as transcription regulators, showed a higher level of expression in the HB1009 (K3a race) strain than in the Xoo KACC10331 (K1 race) strain. Although the additional research is needed to determine the biological function of these genes, these results suggest that these enriched genes might be involved in host-specific immunity responses or facilitate virulence processes in the pathogen. Xoo virulence and the regulatory genes required for pathogenicity are commonly found in pathogenicity islands (PAIs) that encode for a type III secretion system (TTSS) assembled from hrp gene products [46][47][48]. In this study, several genes associated with the expression of hrp genes were significantly upregulated in HB1009. This result suggests that the expression levels of these hrp genes may differ for these two strains in rice varieties, since hrp genes were shown to be essential for bacterial pathogenicity in susceptible hosts and hypersensitive reaction (HR) induction in host and non-host plants in the plant pathogen [47,49]. TAL effectors also play key roles in host immune responses or virulence processes in the pathogen, and the genomic sequences revealed that Philippine (PXO99A), Japanese (MAFF311018), and Korean (KACC10331) strains contain 19, 17, and 15 TAL effector genes, respectively. In addition, Southern hybridization analysis revealed that the HB1009 strain contains 13 TAL effector genes in its chromosome (data not shown). Due to the repetitive structure of the TAL effector coding sequences, NGS-based transcriptome analysis is not sufficient to evaluate their expression. Although additional approaches such as site-directed mutagenesis of these TAL effectors and evaluation of host immunity responses are needed to facilitate the understanding of the function of these TAL effectors, transcriptome analysis revealed that the expression levels of these genes were slightly upregulated in the HB1009 strain; however, these genes did not show significantly different expression levels in the DEGseq analysis (p-value > 0.001, FALSE).
Zhang et al. [31] reported that the bacterial motilities were significantly different among the Xoo strains (C5, China race 5; P2 and P6, Philippines race 2 and 6, respectively), which showed different virulence patterns against various rice cultivars. This study also demonstrated that the genes encoding Hrp proteins and T3 effectors were significantly downregulated in C5, whereas no significantly different virulence pattern was observed among those strains in a virulence assay with various rice cultivars. Thus, they hypothesized that strong motility might compensate for the weaker expression of Hrp proteins in C5, which allows C5 to exhibit similar virulence levels with the other two Xoo strains.
We also conducted a bacterial motility assay with the KACC10331 and HB1009 strains on semi-solid swarm plates. However, there was no significant difference in bacterial motility between the two strains (data not shown). These results suggest that the different virulence patterns of the two Xoo strains might not be due to differences in bacterial motility but, rather, might be due to differently expressed genes, including hrp and TCSs, as well as an enriched gene set according to the enrichment analysis between the two Xoo strains.
Although further studies are required to elucidate the relationship between the differently expressed genes and the Xoo-rice interaction, our findings should facilitate the identification of genes that are involved in this interaction, especially elicitors for the expression of the rice-resistance genes such as avrXa21, and will provide valuable insights to aid in developing future strategies to control BB caused by newly evolved strains such as HB1009 (K3a).

Virulence Assay
Seeds of the rice varieties were sown in a seedling nursery, and 30-day-old seedlings were transplanted. Leaves of 50-day-old (at the tillering stage) rice varieties were artificially infected with the Xoo strains (10 9 cells¨mL´1) by clipping the leaf tips with sterile scissors. The infected rice plants were then grown and maintained under greenhouse conditions (25-30˝C, 60% relative humidity). Distilled water was used as the control treatment. Twenty-one days after inoculation, the means values and standard deviations (SDs) of lesion lengths were calculated for each triplicate set of experiments.

Total RNA Isolation, Illumina Sequencing, and Data Analysis
Total RNA was extracted from a stationary phase culture (OD 600nm 0.8) using an RNeasy mini kit according to the manufacturer's instructions (Qiagen, Hilden, Germany). Total RNA was further treated using RNase-free DNase set (Qiagen, Hilden, Germany) following manufacturer's instructions to discard DNA contamination. Next-generation sequencing (NGS) library preparation and sequencing were performed with 1 µg of each total RNA using the Illumina HiSeq2000 platform in accordance with the manufacturer's protocol (Illumina, Inc., San Diego, CA, USA). The raw data were deposited to the National Center for Biotechnology Information (NCBI) Sequence Reads Archive (SRA) with accession numbers SRP066131 (Xoo KACC10331) and SRP066133 (Xoo HB1009). Sequencing reads were processed by SolexaQA software [50] to control the quality of raw data. The trimmed reads were mapped to the Xoo KACC 10331 genome (NCBI accession no. NC_006834) using Bowtie aligner (http://bowtie-bio.sourceforge.net/index.shtml). The expression level was quantified by the reads per kilobase per million mapped reads (RPKM). DEGseq [25] was applied to identify differentially expressed genes (DEGs) with significance defined as a p-value less than 0.001. Gene Ontology (GO) term annotation and enrichment analysis of the upregulated and downregulated genes were carried out using Blast2GO software [26].

Quantitative Real-Time RT-PCR Assay
To validate the RNA-seq data, a subset of differentially expressed genes (DEGs) was verified by quantitative real-time RT-PCR (qRT-PCR). An independent set of cell cultures and total RNA samples from the two Xoo strains were prepared following the same protocol as for the Illumina analysis. The sequence of each gene was obtained from the Xoo 10331 database (http://www.ncbi.nlm.nih.gov) and used for designing primers by IDT (Integrated DNA Technologies, Coralville, IA, USA, https://eu.idtdna.com/Primerquest/Home/Index). RNA samples from three independent replicates were treated with DNase before cDNA synthesis. Quantitative real-time PCR analysis was performed using a RotorGene 6000 system (Qiagen, Hilden, Germany) in 25-µL reactions containing 12.5 µL SensiFAST SYBR No-ROX kit (Bioline, Sydney, Australia), 5 pmol of each primer (Table S5), and 25 ng of cDNA template. Reaction conditions were as follows: 3 min of holding at 95˝C followed by 40 cycles of 95˝C for 5 s, 60˝C for 10 s, and 72˝C for 15 s. The gene expression levels (arbitrary units) were normalized on the basis of transcript amounts of 16S RNA as an internal reference. The relative expression level of each gene is defined as ∆Ct = Ct target´C t 16S to represent the difference between the transcript abundance of genes examined and the transcript abundance of 16S RNA.