Pathogenicity assays
Resistance and susceptible phenotypes of common bean genotypes BAT93 and JaloEEP558 were demonstrated with bacterial growth and symptom development after inoculation with X. phaseoli pv. phaseoli strain CFBP6546R. For both genotypes, bacterial population sizes decreased during the first day, then increased rapidly up to day 5 and stabilized over time, which is a typical dynamic for common bean colonization by Xanthomonas [51]. However, bacterial population sizes were significantly higher on JaloEEP558 than on BAT93 (p-value < 0.05) at 8 and 15 days post inoculation (Fig. 1A). Moreover, BAT93 presented almost no symptoms, while 39% to 50% of the leaf area were symptomatic on JaloEEP558 (Fig. 1B). Symptoms 15 days post inoculation (DPI) were characterized by the apparition of necrotic areas that were not detected at 8 DPI. These results show that JaloEEP558 is susceptible, while BAT93 is resistant to strain CFBP6546R, as previously described for X. phaseoli pv. phaseoli isolate W18 [50].
Whole genome sequencing and annotation
The genomes of BAT93 and JaloEEP558 were sequenced and annotated to serve as basis for the mapping of corresponding RNA-Seq reads. For each genotype, Illumina sequencing produced around 3.2 × 108 paired reads totalizing 6.4 × 104 Mbp. The resulting assembly consisted of 36,622 scaffolds totalizing 453.0 Mbp for BAT93, and 31,483 scaffolds totalizing 449.3 Mbp for JaloEEP558, with a coverage of approximately 100X for both genotypes (Table 1). These assemblies represented around 84% of the P. vulgaris v2.1 reference genome (537.2 Mbp) from the Andean genotype G19833 [18] or 77% of the estimated genome size (~587 Mbp), and 82% of the BAT93 reference genome (549,6 Mbp) published by Vlasova et al in 2016 [52]. Conversely, structural annotation predicted 33,275 and 32,914 protein coding genes for BAT93 and JaloEEP558, respectively, which are numbers higher than those found in the reference genomes of G19833 (27,433) or BAT93 (30,491). Reciprocal BLASTp (e-value ≤ 1 × 10-6) using all predicted genes showed that our BAT93 and JaloEEP558 genome sequences shared more homologs (27,208) together than with both G19833 and BAT93 reference genomes (Table S1). Thus, to avoid possible biases due to different sequencing and annotation methods, we used our own version of the BAT93 genome for RNA-Seq mapping and comparative transcriptomic analysis between BAT93 and JaloEEP558. To ensure the quality of our data, we used a subset of homologs between BAT93 and JaloEEP558 that had the same best hit after BLASTp on the G19833 reference genome. After removal of paralogs, we obtained 20,787 unique genes that we used for all the analyses presented hereafter (Table S2).
Analysis of differentially-expressed genes
To explore the responses of BAT93 and JaloEEP558 to strain CFBP6546R, transcriptomes of inoculated leaves were produced. In other Xanthomonas-plant pathosystems such as tomato [53], rice [54] or sweet orange [55, 56], bacterial effectors impacted plant transcriptomes around 24 to 48 hours after leaf infiltration. Here, we performed RNA-Seq analyses on inoculated leaves 48 hours after infiltration. Illumina sequencing led to a total of 139.4 and 138.2 million raw reads for BAT93 and JaloEEP558 respectively, with an average of 23.2 million raw reads per sample for BAT93 and 23.0 million raw reads per sample for JaloEEP558 (Table 1). After stringent quality check, data cleaning and mapping, we obtained an average of 20.2 million mapped reads on genes per sample (87% of raw reads) for BAT93 and 20.7 million mapped reads on genes per sample (90% of raw reads) for JaloEEP558. Water-inoculated and bacteria-inoculated plants formed distinct groups after principal component analysis, confirming the similarity of biological replicates within each condition and the similarity of both genotypes after water treatment (Fig. S1). Similar trends were observed using a Pearson correlation matrix, indicating that the bacteria had a significant effect on the transcriptomes of both genotypes. After bacterial inoculation, a total of 5,581 out of the 20,787 homologs were differentially expressed in at least one genotype compared to water inoculation (Table S3), using adjusted p-value < 0.05 and |log2FC| > 1.5 (see Materials and Methods). Comparison of RT-qPCR and RNA-Seq values on 10 genes presenting different patterns of expression in BAT93 and JaloEEP558 revealed a high correlation for these genes (Pearson r = 0.95), further confirming the reliability of RNA-Seq data (Fig. S2).
Global impact of X. phaseoli pv. phaseoli on the transcriptome of common bean
The susceptible genotype initiated a more intense and diverse biological response than the resistant genotype. Indeed, differential expression analysis identified 2,576 and 4,503 differentially-expressed genes (DEGs) in BAT93 and JaloEEP558 respectively, which represents 12% and 22% of the total homologs, respectively (Fig. 2). Three groups of DEGs could be defined (Fig. 3). First, the core transcriptome consisting of 1,482 genes simultaneously induced (758) or repressed (724) in both genotypes. Then, two specific transcriptomes, consisting of 1,094 genes specifically induced (291) or repressed (803) in BAT93, and of 3,021 genes specifically induced (1,367) or repressed (1,654) in JaloEEP558. Enrichment tests identified 83, 32 and 126 gene ontology (GO) terms enriched in the core, the BAT93-specific and the JaloEEP558-specific transcriptomes, respectively (Table S4). To remove redundancy and poorly-informative GO terms, we focused our analysis on biological process GO terms summarized using REVIGO [57]. This analysis highlighted 20, 11 and 29 GO terms enriched in the core, the BAT93-specific and the JaloEEP558-specific transcriptomes, respectively (Fig. 4). We hypothesized that the core transcriptome was representative of the general response of common bean to X. phaseoli pv. phaseoli, while the differences observed between the two specific transcriptomes likely reflected the differences in resistance or susceptibility observed between both genotypes faced to X. phaseoli pv. phaseoli.
Resistance was linked to downregulation of photosynthesis while susceptibility was linked to downregulation of defenses. In BAT93, the most significantly enriched GO terms were related to photosynthesis (Table S4), with a large majority of repressed genes (Fig. 4), suggesting that the resistant genotype strongly suppresses production of primary energy. On the other hand, the most significantly enriched GO terms in JaloEEP558 were related to cell death (Table S4), with a majority of repressed genes annotated as R genes from the NLR family. Together with the fact that genes belonging to the GO term “defense responses” were predominantly repressed in this genotype (Fig. 4), this suggests that global defense responses are suppressed in the susceptible genotype. Interestingly, “signaling” was enriched in the specific transcriptomes of both BAT93 and JaloEEP558 as well as in the core transcriptome. DEGs belonging to “signaling” were mostly upregulated in the core transcriptome, while specific transcriptomes comprised both upregulated and downregulated genes. This suggests that common signaling pathways were induced in both genotypes, while other were specifically up- or downregulated in one genotype or another.
GO analysis done using up- and downregulated genes separately highlighted different GO terms than when using the whole dataset of DEGs (Table S5). It pointed out an upregulation of RNA metabolism and gene expression, as well as in nitrogen compound metabolic process in JaloEEP558, while the later was downregulated in BAT93.
Detailed differences between resistant and susceptible genotypes
MapMan visualization gave a detailed overview of the differences between the resistant and susceptible genotypes (Fig. 5). In accordance with the GO analysis, Wilcoxon Rank Sum Test obtained with MapMan showed that genes related to photosynthesis were enriched in the BAT93-specific transcriptome, while defense genes related to biotic stress and signaling (RLK and PR genes) were enriched in the JaloEEP558-specific transcriptome (Table S6). Additionally, genes related to cell wall modification, ethylene signaling pathway and fatty acid metabolism were specifically enriched in JaloEEP558. To analyze the differences between both genotypes in more detail, we generated lists of genes from different classes (see Materials and Methods) and performed a KEGG analysis of hormonal signaling pathways (Table S7, Fig. S3).
NLR and RLK genes
We observed a specific repression of 30 NLR genes in the susceptible genotype while only one was specifically repressed in BAT93 (Table S7). Additionally, two NLR genes were specifically induced in JaloEEP558, while one (Phvul.006G056500) was both induced in BAT93 and repressed in JaloEEP558, thus appearing as a good candidate for being involved in CBB resistance. The RLK family presented a more complex pattern than NLRs, with a large number of RLK genes specifically induced or repressed in each genotype. More RLKs were repressed than induced in both genotypes, with 11 and 40 specifically induced and 34 and 125 specifically repressed in BAT93 and JaloEEP558, respectively (Table S7). As a result, more RLKs were differentially-expressed in the susceptible genotype than in the resistant genotype.
Kinases
A majority of kinases other than RLKs was induced in the core transcriptome (15/17, Table S7). In particular, four genes encoding calcium-dependent protein kinases (CDPKs) and calcineurin B‐like protein-interacting protein kinases (CIPKs) were induced in both genotypes, indicating that genes linked to calcium signaling were induced in common bean during the interaction with X. phaseoli pv. phaseoli whatever the outcome (resistance or disease). In the JaloEEP558-specific transcriptome, a large modulation of the expression of kinases occurred, which followed the same trend of induction (20/29) than what was observed in the core transcriptome. On the other hand, BAT93-specific transcriptome was less impacted and a majority of kinases (10/12) was repressed.
Transcription factors
A larger induction of transcription factors (TFs) was specifically observed in JaloEEP558 (137/268) compared to BAT93 (28/91, Table S7). TFs from the APETALA2/ethylene response factor (AP2/ERF) superfamily presented the most differences between both genotypes (Fig. S3). Indeed, much more AP2/ERF genes were specifically induced in JaloEEP558 (32) than in BAT93 (4). In the meantime, less AP2/ERF genes were repressed in JaloEEP558 (6) than in BAT93 (8). This indicates a specific induction of AP2/ERF genes in the susceptible genotype 48 h post-inoculation. Similar patterns were observed for the MYB family, although with smaller differences between both genotypes compared to AP2/ERFs (Fig. S3). Other major TF families such as WRKY and GRAS were more impacted (i.e. both more induced and more repressed) in JaloEEP558 than in BAT93.
Lipid metabolism
Lipid metabolism was more impacted in the susceptible genotype (72 specific DEGs) than in the resistant genotype (23 specific DEGs, Table S7). Significantly, a large majority of genes involved in fatty acid synthesis (20/21) was specifically repressed in JaloEEP558 (Fig. 5A, Table S7). Some genes encoding proteins involved in lipid degradation, such as Alpha/beta hydrolases were also repressed in JaloEEP558, while others were upregulated, indicating that lipid degradation was impacted in different ways in the susceptible genotype [58]. On the other hand, fatty acid desaturases were mostly repressed in BAT93. Notably, the homolog of fatty acid desaturase 8 (FAD8, Phvul.006G068600) was both induced in JaloEEP558 and repressed in BAT93. These enzymes could play a role in defense activation [59, 60], hormone synthesis [61], increased membrane fluidity and influence membrane properties of chloroplasts [62]. Here, repression of these genes in the resistant genotype could be linked to the induction of plant defenses.
Photosynthesis and sugar metabolism
A global downregulation of photosynthesis and sugar metabolism occured in the resistant compared to the susceptible genotype (Fig. 5A). In particular, genes linked to the biosynthesis of photosystems I and II, the NAD(P)H deshydrogenase and the ATP synthase where largely repressed in BAT93 (Fig. 5B). By contrast, in JaloEEP558 three genes linked to the photosynthetic electron transport chain were induced. This result indicated a specific repression of the photosynthetic electron transport chain in the resistant genotype. More generally, the whole reactions taking place in the chloroplast appeared down-regulated in the resistant genotype, including the Calvin cycle, tetrapyrrole synthesis, photorespiration, the ascrobate and glutathione redox pathways (Fig. 5A) Sucrose and starch biosynthesis and degradation, and glygolysis were also specifically repressed in BAT93 (Fig. 5A).
Cell wall
Susceptibility was linked to the induction of cell wall modification and repression of cellulose biosynthesis while resistance was linked to the repression of genes involved in cell wall modification. Indeed, most genes involved in cell wall modification and degradation were induced in JaloEEP558 and repressed in BAT93. This was the case of xyloglucan endotransglycosylases, expansins, pectin methylesterase inhibitors, and glycosyl hydrolases (Fig. 5A and Table S7). On the other hand, genes involved in cellulose biosynthesis such as cellulose synthases were repressed in JaloEEP558.
Hormone signal transduction
Resistance and susceptibility were marked by opposite hormone signal transduction pathways. Indeed, resistance was linked to upregulation of SA signaling, while susceptibility was linked to upregulation of the ethylene pathway and downregulation of the SA and cytokinin pathways (Fig. S4). Several genes involved in resistance signaling were specifically repressed in JaloEEP558 such as the homolog of PATHOGENESIS-RELATED 1 (PR1, Phvul.006G197100) that was also specifically induced in BAT93 [63]. PR1 being a marker of SA accumulation, these results suggest that upregulation of the SA pathway occurred in the resistant genotype while it was suppressed in a susceptible context. Consistently, cytokinin signaling was suppressed in JaloEEP558. This is highlighted by the strong induction of a homolog of genes encoding type A Arabidopsis response regulators (A-ARR), which are negative regulators of cytokinin signaling and SA-dependent basal immunity [64]. As described in the MapMan enrichment analysis, many genes involved in ethylene signaling were specifically induced in the susceptible genotype, including AP2/ERF genes (Table S6, Fig. S3). This trend is confirmed by the induction of homologs from ETHYLENE RESPONSE (ETR) and ERF1/2 genes (Fig. S4). Genes involved in other hormonal signaling pathways (auxin, brassinosteroids, gibberellins or abscisic acid) also showed contrasting expression profiles between resistant and sensitive genotypes, although the differences observed were not strong enough to be conclusive.
Candidate genes for resistance to CBB
To search for genes putatively responsible for the resistance observed in BAT93, we have cross-checked different pieces of information that may suggest the involvement of certain genes in resistance to CBB. The different criteria we searched for were (i) that there was a large difference of log2FC values between the susceptible and resistant genotypes, (ii) that the gene was induced in one genotype and repressed in the other, and (iii) that the gene colocalized with a locus of resistance to CBB.
Genes with high log2FC difference
Interestingly, 201 DEGs had a large difference of expression (|Dlog2FC| > 5) between BAT93 and JaloEEP558 (Table S8). Among those genes, 148 (74%) were more expressed in JaloEEP558 than in BAT93. This proportion was higher than what was observed for the whole dataset (52%), and this difference was even higher when focusing on the top 50 genes with the largest |Dlog2FC|, among which 46 had a higher log2FC in JaloEEP558 than in BAT93 (Table S3). This indicates that the greatest differences observed between the two genotypes were due to genes strongly induced in JaloEEP558 and/or strongly repressed in BAT93. Of these 148 genes, the Heat Shock Protein (HSP) family was the most represented with 13 members, four of which belonged to the top 10 genes more expressed in JaloEEP558 than in BAT93. Significantly, three genes from the Lateral Organ Boundaries (LOB) family belonged to the list of 148 genes. Two of these LOB genes (Phvul.007G195100 and Phvul.008G257400) were ranked second and 15th of the most induced genes in JaloEEP558, suggesting that the induction of LOB genes is important for susceptibility. On the other hand, 53 genes had a higher log2FC in BAT93 than in JaloEEP558 (Table S8). A large proportion of these genes belonged to families related to resistance such as RLKs (7), NLRs (2) and PRs (2), indicating that a classical resistance response occurred in BAT93.
Genes with opposite reactions to X. phaseoli pv. phaseoli
Only 16 DEGs were simultaneously induced in one genotype and repressed in the other one (Table S9). Among those, 11 were both repressed in BAT93 and induced in JaloEEP558, suggesting that they contribute to either suppressing defenses and/or promoting disease (Fig. 3). Of those, six had also a |Dlog2FC| > 5, which encoded different proteins including a plasmodesmata-located protein (Phvul.001G229800), an HSP (Phvul.001G039700), a Lipid-Transfer Protein (Phvul.008G137100), a kinase (Phvul.008G081300), a pectin methylesterase inhibitor (Phvul.002G318500) and a sulfite exporter from the TauE/SafE family (Phvul.001G061000) (Table 2, Fig. 6). On the other hand, five genes were both induced in BAT93 and repressed in JaloEEP558, suggesting that they positively regulate resistance to CBB. Four of these genes had also a |Dlog2FC| > 5, including a placenta-specific 8 (PLAC8) family gene (Phvul.003G265800), an S-ribonuclease binding protein (SBP) family gene (Phvul.009G119200), a PR gene (Phvul.006G197100), and an ovate family gene (Phvul.009G057100).
Genes within QRLs
We positioned the four CBB QRLs previously identified in BAT93 on the bean reference genome and extracted the DEGs from specific transcriptomes located within these putative QRLs (Table S10). A total of 600 DEGs were retrieved in these four CBB QRLs, among which 24 (4.0%) had a |Dlog2FC| > 5, which is similar to the ratio observed in the whole transcriptomic data where 201 out of 5,581 DEGs (3.6%) had a |Dlog2FC| > 5 (Fig. 6). In addition, 12% and 21% of the total genes were differently expressed in the QRLs of BAT93 and JaloEEP558, respectively, which is not different from what was observed for the whole transcriptome (Fig. 2). Thus, no specific pattern of expression could be observed in the QRLs compared to the rest of the common bean genome. Among the genes in QRLs with |Dlog2FC| > 5, the three most repressed genes in JaloEEP558 were one NLR gene (Phvul.005G014200) and two RLK genes (Phvul.007G051300 and Phvul.007G030300). Thus, these three genes represent good candidates for CBB resistance. The previously described HSP (Phvul.002G231500, Phvul.002G212000 and Phvul.009G080200) and LOB (Phvul.005G049700) genes specifically induced in JaloEEP558 or repressed in BAT93 were also retrieved in these QRLs.
Only one gene presented the three characteristics of being induced in BAT93 while repressed in JaloEEP558, having a |Dlog2FC| > 5 and locating in a QRL to CBB. This gene (Phvul.009G057100) is located within the QRL associated to marker D0157 on chromosome 9. It belongs to the ovate family, which is involved in plant growth regulation and can suppress elongation [65–67].