MPMI PhytoFrontiers Phytobiomes all journals
RESEARCHOpen Access icon OPENOpen Access license

Transcriptional Profiling of Three Pseudomonas syringae pv. actinidiae Biovars Reveals Different Responses to Apoplast-Like Conditions Related to Strain Virulence on the Host

    Affiliations
    Authors and Affiliations
    • Elodie Vandelle1
    • Teresa Colombo2
    • Alice Regaiolo1
    • Vanessa Maurizio1
    • Tommaso Libardi1
    • Maria-Rita Puttilli1
    • Davide Danzi1
    • Annalisa Polverari1
    1. 1Department of Biotechnology, University of Verona, Verona, 37134, Italy
    2. 2National Research Council of Italy (CNR), Institute of Molecular Biology and Pathology (IBPM) c/o Department of Biochemical Sciences “A. Rossi Fanelli”, “Sapienza” University of Rome, Rome, 00185, Italy

    Published Online:https://doi.org/10.1094/MPMI-09-20-0248-R

    Abstract

    Pseudomonas syringae pv. actinidiae is a phytopathogen that causes devastating bacterial canker in kiwifruit. Among five biovars defined by genetic, biochemical, and virulence traits, P. syringae pv. actinidiae biovar 3 (Psa3) is the most aggressive and is responsible for the most recent reported outbreaks; however, the molecular basis of its heightened virulence is unclear. Therefore, we designed the first P. syringae multistrain whole-genome microarray, encompassing biovars Psa1, Psa2, and Psa3 and the well-established model P. syringae pv. tomato, and analyzed early bacterial responses to an apoplast-like minimal medium. Transcriptomic profiling revealed i) the strong activation in Psa3 of all hypersensitive reaction and pathogenicity (hrp) and hrp conserved (hrc) cluster genes, encoding components of the type III secretion system required for bacterial pathogenicity and involved in responses to environmental signals; ii) potential repression of the hrp/hrc cluster in Psa2; and iii) activation of flagellum-dependent cell motility and chemotaxis genes in Psa1. The detailed investigation of three gene families encoding upstream regulatory proteins (histidine kinases, their cognate response regulators, and proteins with diguanylate cyclase or phosphodiesterase domains) indicated that cyclic di-GMP may be a key regulator of virulence in P. syringae pv. actinidiae biovars. The gene expression data were supported by the quantification of biofilm formation. Our findings suggest that diverse early responses to the host apoplast, even among bacteria belonging to the same pathovar, can lead to different virulence strategies and may explain the differing outcomes of infections. Based on our detailed structural analysis of hrp operons, we also propose a revision of hrp cluster organization and operon regulation in P. syringae.

    Copyright © 2021 The Author(s). This is an open access article distributed under the CC BY-NC-ND 4.0 International license.

    Pseudomonas syringae pv. actinidiae is a Gram-negative pathogenic bacterium that causes bacterial canker in kiwifruit. In the last 10 years, severe outbreaks of the disease worldwide have caused huge economic losses, particularly in Italy, New Zealand, and China, which are among the largest producers. Following an epiphytic phase of multiplication, P. syringae pv. actinidiae can enter kiwifruit plants via natural openings such as stomata and hydathodes, by mechanical wounding, or by pollen dissemination (Vanneste et al. 2011). Symptoms include blossom necrosis and brown to black leaf spots often surrounded by a chlorotic halo. From primary infection sites, bacteria can move systemically through the twigs and trunk, generating extensive cankers that produce white to orange exudates (Scortichini et al. 2012).

    Five different biovars of P. syringae pv. actinidiae have been defined based on genetic, biochemical, and biological traits, including virulence and toxin production. Biovar 1 (here described as Psa1), which can synthesize phaseolotoxin, includes the original P. syringae pv. actinidiae strains found in Japan in 1984. Biovar 2 (here described as Psa2) has only been found in Korea (1992 to 1997), and produces the toxin coronatine. Biovar 3 (here described as Psa3) was responsible for the severe outbreaks of bacterial canker in Italy, New Zealand, Chile, and China between 2008 and 2010, and is the most virulent despite the apparent absence of toxins (Vanneste 2017). Biovar 4, identified in Australia and New Zealand, was initially classified as a low-virulence P. syringae pv. actinidiae strain but was subsequently reclassified as a different pathovar: P. syringae pv. actinidifoliorum (Cunty et al. 2015; Patel et al. 2014). Finally, biovars 5 and 6 were found in small areas in Japan in 2012 and 2015, respectively, and are poorly characterized (Fujikawa and Sawada 2016, 2019). Given the major differences in biovar behavior and the severe economic losses associated with recent Psa3 outbreaks, it is increasingly important to gain insight into the evolution and pathogenicity of P. syringae pv. actinidiae, and this has elevated the species to the status of a model organism for the analysis of bacterial infections in woody plants. In particular, we need to understand the virulence of the most aggressive strains and the triggers that allow endophytic P. syringae pv. actinidiae populations to become virulent, which could facilitate the development of new control strategies.

    Effector proteins are essential bacterial pathogenicity factors exported into plant cells by the specialized type III secretion system (TTSS). They subvert host cellular processes to support bacterial proliferation, thus promoting disease. TTSS components are encoded by hypersensitive reaction and pathogenicity (hrp) and hrp conserved (hrc) genes (Alfano et al. 2000; Arnold et al. 2003; Shindo et al. 2016), which are present in almost all phytopathogenic bacteria, including P. syringae pv. actinidiae (Marcelletti et al. 2011; McCann et al. 2013). The expression of hrp and hrc genes is one of the most important events during infection, and is induced following contact with plant tissues or in minimal hrp-inducing medium (HIM) (Aldon et al. 2000; Galán et al. 2014; Rahme et al. 1992; Rico and Preston 2008). The latter is widely used to study the behavior of bacteria in contact with host cells because it reproduces apoplast-like conditions, including minimal levels of nutrients and glucose or fructose as available carbon sources, and a low pH (Rico and Preston 2008).

    As a first step in the search for Psa3 virulence determinants, genomic analysis revealed genes encoding four TTSS-related effector proteins (HopH1, HopZ5, HopAM1-2, and HopAA1-2) that are missing in other biovars and, therefore, may account for the aggressiveness of Psa3 (Marcelletti et al. 2011; McCann et al. 2013). However, this is not sufficient evidence to confirm the molecular basis of differential virulence because it is not clear whether or not these genes are expressed. Therefore, transcriptomic analysis is required to determine gene expression profiles under different conditions, revealing how the bacteria adapt to interactions with host cells or new growth conditions, as well as shedding light on the regulatory networks that govern such adaptive responses. HIM is often used in this context to trigger the expression of TTSS genes (Rico and Preston 2008).

    To identify P. syringae pv. actinidiae genes and pathways associated with the activation of virulence and strain-dependent aggressiveness, we designed a custom multistrain microarray covering the pan-genome of Psa1, Psa2, and Psa3 as well as the model organism P. syringae pv. tomato DC3000. We carried out a comprehensive transcriptomic analysis of these biovars, including two different Psa3 strains representing the two recent outbreaks that occurred simultaneously in New Zealand (ICMP18884/V-13) and Italy (CRAFRU8.43), to identify differentially expressed genes following the inoculation of minimal medium (HIM) compared with King’s B (KB) rich medium, the latter representing standard in vitro growth conditions. Although the minimal medium is not an exact replica of the apoplast, there are indeed important reasons why it should be regarded as an acceptable substitute depending on experimental purpose. In addition to its use for many years to study the induction of virulence in phytopathogenic bacteria (Huynh et al. 1989; Kim et al. 2009; Rahme et al. 1992)), results generated using minimal medium have not been contradicted by more recent transcriptomic studies in planta, suggesting that this approach generates reliable data, although in planta studies do reveal a broader set of mechanisms (McAtee et al. 2018; Nobori et al. 2018). Moreover, whereas in planta studies provide authentic readouts, the data generally represent the bulk analysis of RNA extracted from pathogen populations, which may display a significant degree of heterogeneity in gene expression and regulation within the host environment (Preston 2017). In the context of our study, which aimed to compare different strains of the same pathovar, as much variation caused by the environment as possible needed to be eliminated. Thus, the minimal medium was ideal, guaranteeing homogeneous bacterial populations. Attempting the same analysis in planta would have generated too much noise and the differences between strains would be far more difficult to discern.

    Hypotheses based on the transcriptomic data were explored using a bioinformatics pipeline to compare i) the structural features of specific P. syringae pv. actinidiae genomic regions and ii) the whole-genome identification and transcriptomic characterization of three bacterial gene families acting as virulence activators or TTSS regulators: histidine kinases (HKs), their response regulators, and proteins involved in cyclic (c)-di-GMP metabolism and signaling. The results highlighted structural and functional differences that might account for the differential aggressiveness and diverse virulence mechanisms in the three P. syringae pv. actinidiae biovars, and also questioned widely accepted paradigms of bacterial pathogenicity by indicating the ability of even biovars within the same P. syringae pathovar to deploy different responses and pathogenicity mechanisms following contact with the host. Thus, our study provides new insight into the molecular basis of virulence in bacterial pathogens that infect tree crops.

    RESULTS

    Multistrain microarray design.

    To find transcriptomic differences among different P. syringae pv. actinidiae biovars, we designed a multistrain microarray containing the whole set of annotated sequences from four P. syringae pv. actinidiae strains belonging to the three main biovars best characterized at the time: NCPPB3739 (Psa1), ICMP19073 (Psa2), ICMP18884 (Psa3), and CRA-FRU8.43 (Psa3). The P. syringae pv. tomato DC3000 genome was included as a control to distinguish general and P. syringae pv. actinidiae-specific transcriptomic features. Finally, we also included annotated transcripts from three integrative conjugative element (ICE) sequences. Annotated transcripts on the multistrain microarray (representing all P. syringae strains and ICE sequences) were collected as cDNA sequences from different public sources (Supplementary Table S1) in July 2015. Additional genome annotations representing the CRA-FRU 8.43 strain were obtained from Prof. Giuseppe Firrao (University of Udine).

    We carried out in silico comparisons to identify sublists of shared and strain-specific sequences from the cDNA libraries, allowing us to prepare a nonredundant cDNA collection comprising the unique multistrain set of annotated protein coding sequences. These BLAST comparisons identified a subset of genes shared by all strains (core genome) as well as subsets of partially shared and strain-specific genes (dispensable genome) The multistrain protein sequence dataset (20,561 sequences) was uploaded to the Agilent eArray web-based application for the design of a five-genome P. syringae high-density microarray for single-strain and cross-strain comparative transcriptomics (Supplementary Fig. S1). Extended datasets for the sequences and bioinformatics pipeline described herein are available on request.

    Assignment of proteins to orthologous clusters for cross-strain comparison.

    To define groups of orthologous sequences, we used reciprocal best hits (Chervitz et al. 1998; Mushegian et al. 1998; Remm et al. 2001; Rubin et al. 2000; Wheelan et al. 1999) to detect orthologous genes in our Pseudomonas dataset based on significant reciprocal similarity between amino acidic sequences. We used BLASTP (Altschul et al. 1990) to perform all pairwise sequence alignments between any two proteins encoded in the five Pseudomonas genomes and three ICE sequences, and recorded the best hit for each protein.

    Close cross-strain investigation of the strain-specific and multistrain sets of protein coding genes highlighted residual redundancy that failed to capture obvious relationships due to minimal differences in the annotated protein sequences. To remove this confounding redundancy and improve the performance of the microarray in cross-strain gene expression comparisons, we further grouped the proteins into superclusters based on slightly relaxed values in the output of BLASTP all-against-all runs (E < 0.001; at least 94% identity over at least 90% of the aligned sequences). The above procedure assigned the full and redundant collection of protein sequences collected for different P. syringae strains from different annotation sources (Supplementary Fig. S1) to a nonredundant set of 10,839 protein groups.

    The identified protein groups were then compared across P. syringae strains to identify sets of protein-coding genes shared by all strains (core genome) as well as partially shared and strain-specific genes (dispensable genome) (Tettelin et al. 2005). This analysis defined a core P. syringae pv. actinidiae or P. syringae pv. tomato genome of 3,551 genes, which expanded to 4,206 genes when we only considered the sequences common to P. syringae pv. actinidiae biovars (Fig. 1A). This broadly agreed with a previous report showing that P. syringae pv. actinidiae includes 3,916 single-copy orthologs shared by four biovars: Psa1, Psa2, Psa3, and Psa5 (Zhao et al. 2019). Finally, the comparison of strains CRAFRU8.43 and ICMP18884 revealed a core Psa3 genome of 6,103 unique supercluster IDs and strain-specific accessory genomes of 307 and 73 unique IDs, respectively (Fig. 1B). This broadly agrees with the 363 clade-specific genes previously identified in Psa3 strains (Zhao et al. 2019).

    Fig. 1.

    Fig. 1. Core and dispensable proteomes of the Pseudomonas syringae pv. actinidiae (Psa) strains. Venn diagrams show the common and unique proteins in A, the five P. syringae pv. actinidiae strains considered in this study and B, in the strains CRAFRU8.43 and ICMP18884 representing Psa3. Venn diagrams were generated using Draw Venn Diagram based on superclusters to avoid protein redundancy among strains.

    Download as PowerPoint

    General transcriptional profiles of P. syringae pv. actinidiae biovars and P. syringae pv. tomato under apoplast-like conditions.

    Next, the four P. syringae pv. actinidiae strain were cultured in rich medium (KB) or minimal medium (HIM) to identify genes modulated specifically under apoplast-like conditions. Samples were harvested after 4 or 8 h of incubation for microarray analysis. The hierarchical clustering of samples based on the similarity of gene expression profiles resulted in dendrograms with high reproducibility among the biological replicates (Supplementary Fig. S2A). The first level of clustering showed a clear separation between the pathovars P. syringae pv. actinidiae and P. syringae pv. tomato, whereas the second level separated the Psa3 strains (CRAFRU8.43 and ICMP18884) from Psa1 (J35) and Psa2 (KN.2). Finally, the third level separated Psa1 and Psa2. The growth conditions (rich or minimal medium) and sampling time points (4 or 8 h) played only a minor role in the clustering of P. syringae pv. actinidiae strains, but the separation of rich and minimal medium was clearly observed for P. syringae pv. tomato. Principal component analysis confirmed the clustering revealed by the dendrogram (Supplementary Fig. S2B and C). Overall, these results indicated that microarray data extrapolated from the experiment can identify genes with statistically significant differential expression among different pathovars or biovars and different media.

    Preliminary analysis revealed a smaller number of differentially expressed genes (DEGs) in the Psa3 strain CRAFRU8.43 (914 DEGs) than in the others (approximately 1,500 DEGs) when comparing growth on minimal and rich medium at both time points, with a false discovery rate (FDR) < 0.05 and an absolute log2 fold change (log2FC) > |1| (Fig. 2A; Supplementary Datasets S1 to S10). The removal of genes modulated in the same manner at both time points confirmed this trend, revealing 648 DEGs in CRAFRU8.43 compared with approximately 1,000 to 1,400 DEGs for the others (Supplementary Fig. S3). Overall, of those genes that were differentially expressed, approximately 50% of genes were upregulated and approximately 50% were downregulated in all strains. When considering the time points individually, the two Psa3 strains (CRAFRU8.43 and ICMP18884) behaved in a similar manner, with most DEGs detected at 8 h (Supplementary Datasets S5 to S8), whereas Psa1 (Supplementary Datasets S1 and S2), Psa2 (Supplementary Datasets S3 and S4), and P. syringae pv. tomato (Supplementary Datasets S9 and S10) featured a higher number of DEGs at 4 h. Although this indicates that Psa3 responds slowly to apoplast-like conditions and shows a later response peak than the other strains, more detailed analysis (described below) indicated that a rapid response to nutrient depletion in the Psa3 strains may have masked the modulation of gene expression at 4 h. Finally, a comparison of DEGs between 4 and 8 h for each strain individually revealed that approximately 30% of DEGs were commonly modulated at both time points in Psa3 (CRAFRU8.43 and ICMP18884) as well as P. syringae pv. tomato, whereas this proportion dropped to 10 to 15% for Psa1 and Psa2, suggesting that these latter strains respond to minimal medium by regulating different metabolic processes progressively over time (Fig. 2B).

    Fig. 2.

    Fig. 2. Transcriptome profiles of different Pseudomonas syringae pv. actinidiae (Psa) strains grown under apoplast-like conditions. Pto = P. syringae pv. tomato. A, Number of differentially expressed genes (DEGs) in the different P. syringae pv. actinidiae strains incubated for 4 or 8 h in hrp-inducing minimal medium compared with King’s B rich medium. B, Number of DEGs commonly regulated at both time points in the different P. syringae pv. actinidiae strains expressed as a percentage of the total number of DEGs in each strain.

    Download as PowerPoint

    Comparison of DEGs among bacterial strains.

    Next, we compared transcriptional data at both time points (without redundancy) to avoid any “time effect” on gene expression. As stated above, strain CRAFRU8.43 appeared less responsive to minimal medium than the other strains. The comparison of CRAFRU8.43 and ICMP18884 showed that almost 80% of DEGs in CRAFRU8.43 (495 of 648) were also modulated in ICMP18884, suggesting a similar response to growth conditions but involving fewer genes (Supplementary Fig. S4). Therefore, to focus on the differences between P. syringae pv. actinidiae biovars, we grouped the two Psa3 strains for further analysis.

    Genes induced in minimal medium were compared between strains to identify common and unique upregulated genes. The Venn diagram in Figure 3A (and Supplementary Dataset S11) shows that 150 genes were commonly induced in all strains, including genes encoding products related to starvation (e.g., a glycerol uptake facilitator, the phosphate starvation-inducible protein PslF, and members of the sodium-solute symporter protein family). Furthermore, 56 of these genes corresponded to hypothetical proteins, the identification of which may provide further insights into the response of P. syringae to apoplast-like conditions. Finally, we observed the common strong upregulation (log2FC > 5) of a gene encoding a C1 family papain-like cysteine protease. Only 51 genes were induced in all P. syringae pv. actinidiae strains but not in P. syringae pv. tomato, suggesting that there was no clear pathovar-specific reaction to minimal medium (Fig. 3A; Supplementary Dataset S11). P. syringae pv. actinidiae-specific induced genes included the copper resistance and binding protein-encoding gene copZ and several genes related to the type VI secretion system, found to be significantly induced in Psa3 in planta (McAtee et al. 2018). There were many similarities between Psa3 and P. syringae pv. tomato, with 126 commonly induced genes not regulated in Psa1 and Psa2, including 41 encoding TTSS components and others encoding flagellar components or enzymes involved in pyoverdine synthesis, which are important aspects of virulence (Supplementary Dataset S11). Additional P. syringae pv. tomato-specific genes were involved in coronatine synthesis and Psa3-specific genes encoded several peptidases. This suggests that P. syringae pv. tomato and Psa3 share some common mechanisms for the induction of virulence but also trigger specific virulence pathways. Whereas most strains featured 250 to 300 uniquely upregulated genes (representing approximately 30% of the induced genes in each strain), only 89 such genes were found in Psa1, suggesting a less-specific response (Supplementary Dataset S11). These genes encoded several proteins related to c-di-GMP signaling or chemotaxis. In Psa2, genes related to classical virulence mechanisms were not upregulated but the unique genes included many related to malonate or sarcosine metabolism as well as sugar transport (Supplementary Dataset S11). Finally, 92 Psa2 genes corresponded to hypothetical proteins, the function of which will provide insight into the behavior of this strain in the apoplast. Interestingly, although the Psa2 genome contains genes for the complete coronatine biosynthesis pathway, none of these genes was induced under our experimental conditions.

    Fig. 3.

    Fig. 3. Comparison of genes modulated in hrp-inducing medium in different Pseudomonas syringae strains. Psa = P. syringae pv. actinidiae, Pto = P. syringae pv. tomato. Venn diagrams show common and unique differentially expressed genes A, upregulated or B, downregulated in the different strains. Venn diagrams were generated using Draw Venn Diagram.

    Download as PowerPoint

    Genes suppressed in minimal medium were compared between strains to identify common and unique downregulated genes (Supplementary Dataset S12). The Venn diagram in Figure 3B revealed a similar distribution to the upregulated genes discussed above. We identified 120 genes repressed in all strains, including many encoding ribosomal proteins or with other roles in translation. Only 57 genes were repressed specifically in P. syringae pv. actinidiae, with no enrichment for particular functions. Among the 71 genes commonly downregulated in Psa3 and P. syringae pv. tomato, we observed the gene encoding Lon protease, a negative regulator of the TTSS, consistent with the upregulation of TTSS-related genes in these strains. In most strains, we detected approximately 200 uniquely downregulated genes. Again, the exception was Psa1, where the number dropped to 64, supporting the lower specificity of its response.

    Biovar-specific responses to apoplast-like conditions.

    To identify the pathways potentially triggered in each strain by the perception of minimal medium, we used gene ontology (GO) analysis to reveal the functionally enriched categories among the genes upregulated in minimal medium (Fig. 4). Only three categories were commonly enriched in all strains (transcription, styrene catabolism, and vesicle organization), whereas the most enriched categories in Psa3 and P. syringae pv. tomato (defined as –log10[1/adjusted P value] ≥ 15) were related to pathogenesis, the TTSS, and flagellum-dependent cell motility. Another three categories were solely enriched in Psa3 (iron transmembrane transport, glycerol-3-phosphate transmembrane transport, and fumarate metabolism), whereas sporulation and organic hydroxyl compound transport functions were solely enriched in P. syringae pv. tomato. In contrast, the GO terms that were most enriched in Psa1 included categories related to flagellum-dependent cell motility, chemotaxis, and signal transduction, whereas less-enriched categories related to polyhydroxybutyrate biosynthesis and disaccharide metabolism were solely enriched in Psa1. GO enrichment analysis in Psa2 revealed an overall lower degree of enrichment compared with other strains (–log10[1/adjusted P value] = 3 to 4) and highlighted functional categories related to signal transduction and the metabolism of carbohydrates, alditol phosphate, branched amino acids (valine and leucine), thioester, cyanate, xylulose, acetyl-CoA, and olefin.

    Fig. 4.

    Fig. 4. Gene ontology (GO) categories enriched among upregulated genes in different Pseudomonas syringae strains under apoplast-like conditions. Psa = P. syringae pv. actinidiae, Pto = P. syringae pv. tomato. Data show clusters of GO categories overrepresented in differentially expressed genes that are upregulated in different P. syringae pv. actinidiae strains in hrp-inducing minimal medium compared with rich King’s B medium, relative to the whole genome. Cluster enrichment was determined using Blast2GO Pro.

    Download as PowerPoint

    When we applied the same strategy to identify enriched GO categories among the genes downregulated in minimal medium, the functions of the most enriched categories in each strain (Supplementary Fig. S5) and the functions that were enriched solely in a particular strain (Supplementary Fig. S6) were more diverse, pointing to a profound reprogramming of global metabolism. The downregulated genes in all strains were enriched for functions related to translation, ATP synthesis, mitochondrial transport, and glycine betaine biosynthesis, suggesting the general reallocation of energy reserves (Supplementary Fig. S5A). Some categories were uniquely enriched among the downregulated genes in P. syringae pv. actinidiae, including those related to xanthine metabolism and ribosome biogenesis (Supplementary Fig. S5B). Others were shared between Psa3 and P. syringae pv. tomato, including functions related to ion homeostasis, citrate transport, lysine, proline and γ-aminobutyric acid catabolism, folic acid biosynthesis, and heat responses (Supplementary Fig. S5C). Several enriched GO categories were shared by at least two strains (Supplementary Fig. S5D). GO categories that were enriched among the downregulated genes uniquely in a given strain included the removal of superoxide radicals in Psa1; branched-amino acid transport, carbohydrate utilization and sucrose metabolism, lipid A and lipopolysaccharide core region biosynthesis, and aromatic amino acid biosynthesis in Psa2; and bacteriocin immunity, guanosine-containing compound biosynthesis, and superoxide metabolism in Psa3 (Supplementary Fig. S6).

    Interestingly, sample clustering based on enriched GO categories among the upregulated genes indicated that Psa3 and P. syringae pv. tomato shared similar traits, with a clear induction of pathogenesis-related processes, in contrast to Psa1 and Psa2 (Supplementary Fig. S7A). However, sample clustering based on enriched GO categories among downregulated genes showed greater similarity between Psa1 and Psa3 on the one hand, and Psa2 and P. syringae pv. tomato on the other (Supplementary Fig. S7B). This confirms that transcriptome profiling is needed to fully understand the similarities and differences among bacterial strains, and indicates that the behavior of Psa1 is intermediate between the mild Psa2 and virulent Psa3 strains.

    Another interesting observation concerned the hrp/hrc cluster, which contains 26 open reading frames organized in seven operons (Alfano et al. 2000). The analysis of gene expression within the hrp/hrc cluster (Fig. 5A) showed that all genes, except those encoding the upstream regulators HrpR and HrpS, were strongly upregulated in Psa3 and also in P. syringae pv. tomato (Supplementary Fig. S8). In contrast, only a few genes (hrpA1, hrpZ, hrpF, and hrpO) were induced in Psa1, and the expression level was lower than in Psa3 (log2FC approximately 4 in Psa3 and approximately 2 in Psa1). Only hrpT was induced in Psa2, in common with the other strains.

    Fig. 5.

    Fig. 5. Modulation of hrp/hrc genes and type III secretion system (TTSS) regulators in Pseudomonas syringae pv. actinidiae (Psa) strains incubated in minimal medium. A, Graphical representation of hrp/hrc gene expression in the different P. syringae pv. actinidiae strains. Each bar represents the log2 fold change of expression for the corresponding gene. Color blocks indicate the organization of hrp/hrc genes in operons. B, Hierarchical clustering of the absolute expression level of the hrp/hrc genes and the main TTSS regulators in the different P. syringae pv. actinidiae strains in minimal (hrp-inducing) and rich (King’s B) medium at 4 and 8 h. Clusters were generated using MeV with normalized fluorescence values. Blue triangles indicate the different clusters. C, Graphical representation of the trend in gene expression for clusters 1 and 2 in the different P. syringae pv. actinidiae strains.

    Download as PowerPoint

    To confirm the differential response of the TTSS in the P. syringae pv. actinidiae biovars, we transformed them with a reporter system in which the expression of green fluorescent protein (GFP) is driven by the hrpA1 promoter (phrpA1-GFP) (Vandelle et al. 2017). The modified strains were used to monitor hrpA1 promoter activation as a marker of TTSS activity in minimal medium (Supplementary Fig. S9). In agreement with the microarray data, hrpA1 promoter activity was strong in Psa3, already significant after incubation for 1 h in minimal medium, whereas the induction was slower and approximately 50% weaker in Psa1, and no signal was detected in Psa2 over a period of 6 h, confirming that the TTSS is not induced in this strain under these conditions. Moreover, the addition of kiwifruit leaf extract to the minimal medium did not trigger Psa2 TTSS expression, whereas it significantly enhanced hrpA1 promoter activity in both Psa1 and Psa3 (Supplementary Fig. S9).

    Biovar-dependent responses of the TTSS.

    Given the importance of the TTSS in bacterial pathogenicity, and based on the differential expression of hrp and hrc genes in the P. syringae pv. actinidiae biovars described above, we investigated the correlation between hrp and hrc gene expression and virulence in more detail. Hierarchical clustering of hrp and hrc gene expression profiles across all samples using microarray normalized fluorescence values, including genes encoding the major regulators AefR, RhpR/RhpS, RhpC/RhpP, RpoN, AlgU, CsvR/CsvS, and Lon (Xie et al. 2019), revealed two main clusters. The first included all modulated hrp and hrc genes except hrpT, hrpV, and hrcN whereas the second included the major TTSS regulators (Fig. 5B). In cluster 1, Psa3 clearly responded more strongly to minimal medium compared with rich medium, and a similar but weaker response was also observed for Psa1, whereas Psa2 showed no significant modulation in response to minimal medium. In contrast, all P. syringae pv. actinidiae biovars showed identical responses to minimal medium in the second cluster at both time points (Fig. 5C). Surprisingly, all of the positive regulators of TTSS were significantly induced in Psa2 after only 4 h of incubation in minimal medium, whereas only some of them were induced in Psa1 and Psa3 (Supplementary Fig. S10). In contrast, the negative regulator Lon protease was repressed solely in Psa3. The same hierarchical clustering approach showed that the expression level of hrp/hrc genes and TTSS regulators was higher in Psa3 grown in rich medium at 4 h compared with other rich medium samples but declined at 8 h. This indicates that at least some hrp and hrc genes and regulators are already induced by the overnight incubation of Psa3 in rich medium prior to the switch to minimal medium, potentially explaining the apparent weaker response of Psa3 at 4 h.

    The conserved structure of the hrp/hrc cluster in different P. syringae pv. actinidiae biovars.

    The hrp/hrc cluster expression profiles revealed that the genes significantly upregulated in Psa1 almost invariably corresponded to the first genes of each operon, based on the previously described operon structure (Alfano et al. 2000), with hrpO as the only exception. This was mirrored by the expression profiles in Psa3, where the abovementioned genes were expressed more strongly than the others. These data support the general observation that gene expression in operons declines with distance from the transcriptional start, due to the progressively increasing likelihood of ribosomal detachment (Lim et al. 2011). Although the general proximal–distal expression gradient supports the reliability of our microarray data, hrpO displayed a higher log2FC value than the other genes in the hrpJ operon (hrcN, hrcQ, hrcV, and hrpJ) and the hrp box associated with the hrpP gene is >300 bp upstream, not as close to the transcriptional start site as other hrp boxes in the cluster (Supplementary Fig. S11). Therefore, the role of the hrp box in the regulation of hrpP expression is unclear. We propose a hypothesis in which hrpO is not the last gene of the hrpJ operon, as currently assumed, but the first gene of the hrpU operon, which lies immediately downstream. We first tested our hypothesis in silico by operon prediction using the P. syringae pv. actinidiae ICMP18884 genome as a template, which suggested the presence of a single operon spanning all of the genes from hrpJ to hrcU (Supplementary Fig. S12). Given these inconclusive data, we tested for the presence of polycistronic transcripts in the Psa3 strain CRAFRU8.43. First, we detected the presence of individual transcripts for hrcN, hrpO, and hrpP, and their expression only in response to minimal medium confirmed our microarray data and ensured the efficient retrotranscription (Supplementary Fig. S13). We then designed primers that spanned pairs of genes (hrcN-hrpO and hrpO-hrpP) and, using the same cDNA as for individual gene-transcript-level analysis, we detected amplicons only when using the hrpO and hrpP primers, thus demonstrating that hrpO and hrpP produce a polycistronic transcript and are part of a single operon (Fig. 6A). This confirms that hrpO is the first gene of the hrpU operon and lies just upstream of hrpP, and further implies that the hrp box previously assumed to regulate hrpP may actually regulate hrpO, even though it is found within the hrpO coding sequence (Fig. 6B). The lack of expression in Psa1 and Psa2 made it impossible to test for polycistronic transcripts in these biovars but the same operon structure was confirmed by the presence of polycistronic transcripts in P. syringae pv. tomato (Supplementary Fig. S14). These data suggest that the organization of operons hrpJ and hrcU may be conserved among P. syringae strains and that the presence of an upstream hrp box does not necessarily specify the position of an operon.

    Fig. 6.

    Fig. 6. Proposed revision of hrpJ and hrcU operon organization in Pseudomonas syringae pv. actinidiae biovar 3 (Psa3). A, Reverse-transcription PCR analysis of polycistronic transcripts containing the hrpN, hrpO, and hrpP sequences in Psa3 (CRAFRU8.43). Primers were designed to amplify polycistronic transcripts encoded by the hrpN-hrpO and hrpO-hrpP genes. Total RNA extracted from CRAFRU8.43 cells cultured for 4 or 8 h in minimal medium (hrp-inducing medium [HIM]) was reverse transcribed and amplified by PCR. Genomic DNA (g) was used as positive control. B, Coding sequence of hrpO including the hrp box (underlined). The nucleotides highlighted in red correspond to the hrp box consensus motif (GGACC-N15-CCAC-N2-A).

    Download as PowerPoint

    The complete sequences of the hrp/hrc clusters from each P. syringae pv. actinidiae biovar were aligned to identify genomic differences that might explain the strong induction of hrp/hrc gene expression in Psa3, the weak induction in Psa1, and the absence of a response in Psa2. The hrp/hrc clusters of Psa1 and Psa3 were closely related, with few nucleotide substitutions, whereas weaker similarity was observed between Psa1 or Psa3 and Psa2, in particular in the region covering the hrcU operon (Supplementary Fig. S15A). To ensure that sequence variations in the hrcU operon did not reflect poor genome annotation, we further aligned the hrp/hrc clusters of two additional Psa2 genomes (ICMP19071 and ICMP19072, which showed high identity with ICMP19073), confirming that the sequence of the Psa2 hrcU operon diverges from the other two biovars (Supplementary Fig. S15B). However, the nucleotide variations between Psa2 and Psa3 translated to 98 to 100% identity at the protein level, suggesting that most of the nucleotide changes resulted in synonymous codons. Overall, these data indicate that the modulation of hrp/hrc cluster genes in different P. syringae pv. actinidiae biovars is probably related to different signaling pathways and regulators rather than differences in the hrp/hrc genome sequence.

    Common and biovar-specific expression profiles of key signaling proteins.

    In bacteria, signal transduction is largely mediated by members of the HK family, which phosphorylate downstream response regulator (RR) proteins as part of a two-component system (Gao and Stock 2009). In some cases, the HK and RR domains are present on the same hybrid histidine kinase (HHK) protein (Buschiazzo and Trajtenberg 2019). Therefore, we focused on the role of these key signaling proteins in the response of P. syringae pv. actinidiae biovars to the apoplast-like conditions imposed by minimal medium.

    Putative HKs were identified using MIST3, retrieving 65, 61, and 66 proteins containing at least one putative HK domain in the genomes of Psa1, Psa2, and Psa3, respectively. Additional candidates were identified by using BLASTP to screen all HKs from each biovar against the genomes of the other two biovars, and by screening all three biovars with the HKs identified in P. syringae pv. tomato (Lavín et al. 2007). Finally, a whole-genome BLASTN search using probes corresponding to HK sequences confirmed, at the nucleotide level, the presence or absence of further putative HKs. Overall, the analysis identified 72 common candidates (21 HHKs and 51 HKs) in all P. syringae pv. actinidiae strains, as well as 3 (2 HHKs and o1 HK) unique to Psa2, although the 2 HHKs were also represented in the Psa3 genome as pseudogenes (Supplementary Fig. S16). We found that 29 of the genes were differentially expressed (24 induced and 5 repressed), including 16 that were modulated in Psa1 (15 induced and 1 repressed), 22 that were modulated in Psa2 (17 induced and 5 repressed), and only 11 that were modulated in Psa3 (8 induced and 3 repressed). Of the 75 identified HKs and HHKs, 37 contained one or more transmembrane domains (maximum = 13) and 8 included a signal peptide. Various other domains were also present in these proteins (e.g., GAF, PP2C, 2CSK_N, dCache1, CHASE3, KdpD, Usp, KinB, PHY, CheW, Hypoth_Ymh, and RsbRD), although the most widely represented domains were PAS/PAC (27 proteins) and HAMP (19 proteins).

    Signal perception requires the presence of sensors at the onset of environmental changes; thus, it is tempting to speculate that basal expression levels of HK or HHK proteins (i.e., in rich medium) may also account for the specific response of different biovars to minimal medium. The hierarchical clustering of HK or HHK expression profiles using normalized fluorescence intensities, representing cells grown in rich medium, yielded one cluster containing Psa1 and Psa2 at both time points, and a second cluster containing Psa3 at both time points. These results suggest that HK or HHK expression profiles in Psa1 and Psa2 are similar in rich medium (at least for the first 8 h). In contrast, Psa1 and Psa3 clustered together and were separated from Psa2 in minimal medium, further supporting the unique behavior of Psa2 under apoplast-like conditions.

    We identified eight major gene clusters with specific expression profiles in the different biovars (Fig. 7A). Cluster 1 included HK or HHK genes that were constitutively expressed at high levels in Psa3 and at low levels in Psa2, and induced by minimal medium in Psa1 (Fig. 7B). Cluster 7 included HK or HHK genes that showed constitutively higher expression levels in Psa3 in rich or minimal medium compared with the other two biovars (Fig. 7C). The corresponding proteins, whose expression was not modulated by switching from rich to minimal medium, are interesting candidates for external signal perception in the most virulent biovar, potentially accounting for the rapid induction of virulence given their constitutively high expression level. Cluster 7 also included one gene that was expressed at particularly high levels in Psa3 under all conditions; however, it was annotated as a pseudogene due to the presence of a small deletion producing a slightly truncated product compared with its counterpart in Psa2 (WP_020312785.1). The presence of transcripts encoded by this pseudogene was confirmed by reverse-transcription (RT)-PCR (Supplementary Fig. S17).

    Fig. 7.

    Fig. 7. Hierarchical clustering of the expression profiles of histidine kinases in Pseudomonas syringae pv. actinidiae. A, Gene and sample trees obtained from the hierarchical clustering of the absolute expression level of all members of the histidine kinase (HK) family identified in this study, including hybrid histidine kinases. Normalized expression values based on microarray data were used for hierarchical clustering based on Pearson’s distance metric. The color scale represents higher (yellow) or lower (blue) expression levels with respect to the median transcript abundance of each gene across all samples. Psa3 is represented by the strain ICMP18884. Gene and sample clustering was performed with MeV using a distance threshold adjustment of 1.5. Blue triangles indicate different clusters of HK genes showing unique expression trends among P. syringae pv. actinidiae (Psa) biovars depending on the time and growth conditions. B, Detailed description of cluster 1 and the corresponding graphical representations of expression trends among P. syringae pv. actinidiae biovars. Cluster 1 = genes expressed at higher levels in Psa3 than the other strains in both media, and induced in Psa1 in minimal medium. C, Detailed description of cluster 7 and the corresponding graphical representations of expression trends among P. syringae pv. actinidiae biovars. Cluster 7 = genes expressed at higher levels in Psa3 than other strains in both media and not modulated in the other strains. Samples are not clustered in B and C and, therefore, the order differs from the main cluster shown in A.

    Download as PowerPoint

    The procedure described above for HK and HHK genes was repeated for the identification of RR candidates. This led to the identification of 74 common RR genes, together with 1 uniquely found in Psa2 and 1 found in Psa2 and Psa3 but not Psa1 (Supplementary Fig. S18). Any sequences containing HK domains were excluded because they were already counted among the HHK genes discussed above. Only one of the RR candidates contained a transmembrane domain but various other domains were identified (e.g., PAS/PAC, HPT, CheW, AAA ATPase, HTH_8 (Fis-type), GGDEF, EAL, HD-GYP, CheB_methylest, CHASE3, LyTR, PP2C_SIG, TPR, and ANTAR), although the most widely represented domains were Trans_reg_C (20 proteins) and HTH_LuxR (12 proteins). This was consistent with the function of RRs, which often act as direct transcriptional regulators. As already observed for the HK and HHK genes, the RR genes were mainly upregulated in response to minimal medium, with only three genes repressed (one in both Psa1 and Psa3, one in Psa2 alone, and one in Psa3 alone). Although Psa3 is the most virulent strain, only 12 RR genes were upregulated in response to minimal medium (mainly at 8 h), compared with 20 in Psa1 and 21 in Psa2 (mainly at 4 h). As discussed above for HK and HHK genes, the RR genes involved in rapid external signal transduction may be already expressed in Psa3 before the switch to minimal medium. The analysis of normalized fluorescence values revealed RR genes that were not modulated in Psa3 but were expressed at higher levels compared with the other biovars (Fig. 8). The RR genes formed 10 distinct clusters, with cluster 1 comprising those expressed at higher levels in Psa3 than the other biovars, even in rich medium (Fig. 8A). Interestingly, some of these genes were slightly upregulated in Psa1 grown in minimal medium but expression levels remained low in Psa2 under all conditions (Fig. 8B). Within this cluster, the gene encoding protein WP_017682831.1 was constitutively expressed at a high level in Psa3. This gene is not present in the Psa1 genome, and is found close to one of the HHK pseudogenes in Psa3, in a chromosomal region that has been newly acquired by the most recent biovars (Psa2 and Psa3).

    Fig. 8.

    Fig. 8. Hierarchical clustering of the expression profiles of response regulators in Pseudomonas syringae pv. actinidiae. A, Gene and sample trees obtained from the hierarchical clustering of the absolute expression level of all members of the response regulator family identified in this study. Normalized expression values based on microarray data were used for hierarchical clustering based on Pearson’s distance metric. The color scale represents higher (yellow) or lower (blue) expression levels with respect to the median transcript abundance of each gene across all samples. P. syringae pv. actinidiae biovar 3 (Psa3) is represented by the strain ICMP18884. Gene and sample clustering was performed with MeV using a distance threshold adjustment of 1.5. Blue triangles indicate different clusters of response regulator members showing unique expression trends among P. syringae pv. actinidiae biovars depending on the time and growth conditions. B, Detailed description of cluster 1 and the corresponding graphical representations of expression trends among P. syringae pv. actinidiae biovars. Cluster 1 = genes expressed at higher levels in Psa3 than other strains in both media and not modulated in the other strains. Samples are not clustered in B and, therefore, the order differs from the main cluster shown in A.

    Download as PowerPoint

    HKs are also regulated at the posttranslational level by the second messenger c-di-GMP, which plays a key role in bacterial signal transduction (Cheng et al. 2019). This molecule is produced by diguanylate cyclases, which possess a canonical diguanylate cyclase domain (GGDEF), and is degraded by phosphodiesterases (PDEs), which may possess EAL or HD-GYP domains (Richter et al. 2019). Such proteins comprise the second-most common group of two-component system output domains, establishing a link between the two major systems in bacterial signaling. HHKs are activated by c-di-GMP, which binds to a specific pseudo-Rec domain (Dubey et al. 2020). Therefore, we carried out a whole-genome screen for GGDEF, EAL, and HD-GYP domain-containing genes and analyzed their expression as above. We identified 46 genes potentially involved in c-di-GMP metabolism, 17 encoding GGDEF proteins, 4 encoding EAL proteins, 17 encoding proteins with both domains, and 7 encoding proteins with an HD-GYP domain (Supplementary Fig. S19). Of these genes, 40 were common to all strains, and the remaining 6 were missing in at least one biovar. In total, 19 of the encoded proteins contained at least one transmembrane domain and 10 also contained a PAS/PAC domain. Other domains were present in the GGDEF/EAL proteins, including pseudo-Rec (see above), dCache1, CHASE, MASE2, GAF, MHYT, and HAMP. In contrast, the HD-HYP proteins contained NTP_transf_2, GlnD_UR_UTase and ACT, PolyA_Pol and PolyA_pol_RNAbd, and tRNA_edit domains.

    Five of the genes encoding GGDEF or EAL proteins were commonly modulated in the different P. syringae pv. actinidiae biovars in response to minimal medium: three upregulated (WP_003378933.1, WP_017683475.1, and WP_017683804.1) and two downregulated (WP_017683933.1 and WP_017684025.1). However, most of the genes related to c-di-GMP signaling were modulated in a biovar-specific manner. Three of four GGDEF genes modulated in Psa3 were downregulated, whereas two GGDEF genes were upregulated in Psa1 and Psa2 (although different genes in each biovar). Only two EAL genes were modulated in response to minimal conditions and were uniquely upregulated in Psa1. These data point to a trend of declining c-di-GMP levels in Psa1 and Psa3 (suppression of synthesis or promotion of degradation) whereas the upregulation of genes involved in c-di-GMP synthesis tended to maintain steady levels in Psa2. Hierarchical clustering of normalized fluorescence values revealed six gene clusters with transcript profiles differing among the three biovars (Fig. 9A). Cluster 1 contained seven genes with constitutively higher expression levels in Psa3 compared with the other biovars. Cluster 3 contained four genes with constitutively lower expression levels in Psa3 compared with the other biovars, and these genes were induced in response to minimal medium in Psa1 but not modulated by the medium in Psa2 (Fig. 9B). The two genes from cluster 1 showing the highest expression levels in Psa3 encoded one EAL protein and one HD-GYP protein, both involved in the degradation of c-di-GMP, whereas the gene in cluster 3 most strongly repressed in Psa3 encoded a GGDEF protein involved in c-di-GMP synthesis (Fig. 9B). The behavior of the P. syringae pv. actinidiae biovars was reflected at the level of sample clustering, with Psa1 and Psa2 forming one cluster in rich medium (indicating similar expression of genes involved in c-di-GMP metabolism) but resolving to different clusters following the switch to minimal medium, whereas Psa3 formed a single cluster encompassing both media (indicating the limited effect of growth conditions on the regulation of these genes). Taken together, these data indicate that Psa3 has a lower constitutive level of c-di-GMP than the other biovars. This was confirmed using the pCdrA-gfpC reporter system previously described for c-di-GMP level analysis (Cerna-Vargas et al. 2019). Indeed, the fluorescence signal, proportional to c-di-GMP level in bacterial cells, was significantly lower in Psa3 compared with the other biovars (Fig. 10A and B), demonstrating that basal c-di-GMP levels are actually lower in Psa3. On the other hand, the high levels present in Psa1 are likely depleted following exposure to apoplast-like conditions, thanks to the induction of fosfodiesterase activities, whereas those in Psa2 are unaffected.

    Fig. 9.

    Fig. 9. Hierarchical clustering of the expression profiles of cyclic di-GMP-related enzymes in Pseudomonas syringae pv. actinidiae. A, Gene and sample trees obtained from the hierarchical clustering of the absolute expression level of all genes containing GGDEF, EAL, and HD-GYP domains identified in this study. Normalized expression values based on microarray data were used for hierarchical cluster analysis based on Pearson’s distance metric. The color scale represents higher (yellow) or lower (blue) expression levels with respect to the median transcript abundance of each gene across all samples. P. syringae pv. actinidiae biovar 3 (Psa3) is represented by the strain ICMP18884. Gene and sample clustering was performed with MeV using a distance threshold adjustment of 1.5. Blue triangles indicate different clusters of genes containing GGDEF, EAL, or HD-GYP domains showing unique expression trends among P. syringae pv. actinidiae biovars depending on the time and growth conditions. B, Detailed description of cluster 1 and the corresponding graphical representations of expression trends among P. syringae pv. actinidiae biovars. Cluster 1 = genes expressed at higher levels in Psa3 compared with other strains in both media and not modulated in other strains. C, Detailed description of cluster 3 and the corresponding graphical representations of expression trends among P. syringae pv. actinidiae biovars. Cluster 3 = genes expressed at lower levels in Psa3 than other strains in both media and upregulated in Psa1. Samples are not clustered in B and C and, therefore, the order differs from the main cluster shown in A.

    Download as PowerPoint
    Fig. 10.

    Fig. 10. Levels of cyclic (c)-di-GMP, biofilm formation, and swarming motility in the different Pseudomonas syringae pv. actinidiae (Psa) biovars. A and B, Fluorescence emission by P. syringae pv. actinidiae strains belonging to different biovars harboring the c-di-GMP biosensor plasmid pCdrA::gfp previously described by Cerna-Vargas et al. (2019). A, Transformed strains were grown on solid King’s B (KB) medium and pictures were taken with a fluorescence microscope after 48 h. Experiments, including three technical replicates, were repeated twice with the same results. B, Transformed strains were washed in liquid hrp-inducing medium (HIM) for 15 min and fluorescence intensity was measured in microtiter plates with a fluorimeter. Shown are means of replicates and standard errors from three technical replicates. The experiment was repeated twice with similar results; ** indicates P value < 0.01 with Student’s t test. C, Biofilm formation was quantified in Psa1, Psa2, and Psa3 cell suspensions cultured in 96-well plates for 24 h at 28°C in KB medium or HIM and stained with crystal violet. Absorbance was measured at an optical density at 590 nm. Wells containing no cells were used as negative controls. Different letters indicate a statistically significant difference (P < 0.05). D, Swarming motility of different P. syringae pv. actinidiae biovars was evaluated in solid KB medium containing 0.4% agar. Pictures were taken after 48 h of incubation. The experiment was repeated twice with similar results.

    Download as PowerPoint

    Inverse relationship between TTSS induction and biofilm production.

    Biofilm formation is positively regulated by c-di-GMP (Wang et al. 2019). Our expression data and c-di-GMP level assessment indicated that Psa2 may produce more c-di-GMP than the other biovars, explaining the absence of TTSS induction in this biovar. Therefore, we tested for biofilm formation in the three strains, anticipating that Psa2 would readily form biofilms but the other biovars would not. As expected, we found that Psa1 and Psa3 were poor biofilm producers in both media, reaching maximum crystal violet absorbance values of 0.8 and 0.4 after 24 h in rich and minimal medium, respectively. Psa2 formed biofilms more readily, achieving absorbance values of 1.6 in rich medium (twofold higher) and 1.2 in minimal medium (threefold higher) over the same period (Fig. 10C). Differences were already visible after only 4 h, with absorbance values of 0.6 (Psa2) and 0.2 (Psa1 and Psa3) in both media. By contrast, in line with the inverse regulation of biofilm formation and swarming motility in Pseudomonas spp. (Engl et al. 2014; Fishman et al. 2018), Psa2 showed a very low swarming capacity on KB agar compared with Psa1 and Psa3 (Fig. 10D).

    DISCUSSION

    P. syringae causes frequent epidemic diseases in herbaceous and woody crops and is among the best-studied bacterial phytopathogens (Lamichhane et al. 2014, 2015). P. syringae pv. actinidiae has recently emerged as a model organism for the analysis of bacterial pathogens affecting woody species due to the severe impact of recent disease outbreaks (Morris et al. 2019). However, it is unclear why certain P. syringae pv. actinidiae strains are so virulent, particularly those representing biovar Psa3 (Vanneste 2017). Therefore, we investigated the aggressiveness of the main three P. syringae pv. actinidiae biovars using a multistrain microarray, which remains a useful approach despite the advance of RNA sequencing technology because it provides a direct visual signature for the comparison of multiple samples. Intrinsic redundancy due to cross-hybridization was avoided by whole-genome comparisons, allowing the microarray to be used simultaneously for transcriptomic analysis and the tentative functional annotation of genes not yet represented in the available genome sequences of poorly characterized strains. Our single multistrain microarray also provided a common reference name for all P. syringae pv. actinidiae strain genes and orthologs from the control pathovar P. syringae pv. tomato.

    P. syringae pv. actinidiae biovars and P. syringae pv. tomato share a common response to nutrient deficiency.

    The microarray data revealed the common upregulation of many genes involved in transcriptional regulation as a core mechanism shared among the P. syringae pv. actinidiae and P. syringae pv. tomato strains. The modulated genes included those encoding σ factors (which control the transcription of broad sets of target genes) and proteins related to phosphate starvation, including a PhoH-like protein, the transcriptional regulator PhoP, the phosphate starvation-inducible protein PsiF, and the phosphonate metabolism protein PhnM. Phosphate limitation has a pleiotropic effect on bacterial physiology, triggering the degradation of polyphosphate, the accumulation of phosphate-free membrane lipids, and the accumulation of polyhydroxyalkanoate (PHA) in regulatory granules (Hokamura et al. 2015; Romano et al. 2015). The P. syringae pv. actinidiae biovars and P. syringae pv. tomato induced two main genes related to PHA synthesis, encoding the phasins PhaI and PhaF. However, the functional category related to PHA appeared enriched only in Psa1, which induced three other PHA-related genes despite the generally low number of genes uniquely upregulated in this biovar. We also observed the strong common induction of a MaoC-like hydratase (MaoC), an (R)-hydratase that links β-oxidation to the PHA biosynthesis pathways (Fiedler et al. 2002; Fukui and Doi 1997; Fukui et al. 1998; Wang et al. 2013). PHAs provide intracellular carbon and energy reserves under nutrient-limited conditions with excess carbon (Hokamura et al. 2015); thus, this pathway may provide a common adaptive response allowing P. syringae to overcome nutrient limitation within the apoplast.

    The stringent response triggers the downregulation of genes encoding translational components and the simultaneous upregulation of genes involved in amino acid biosynthesis and transport (Potrykus and Cashel 2008). This reflects the combined actions of the DnaK suppressor protein DksA and the nucleotide alarmones (p)ppGpp on RNA polymerase (Kim et al. 2018). Accordingly, we observed the induction of two genes encoding DnaK suppressor proteins but we did not observe the modulation of genes related to stress-induced alarmone production (RelA and SpoT). All strains also induced an operon encoding the serine protein kinase PrkA, a protein containing a von Willebrand Factor type A domain, and a SpoVR-like protein (Vercruysse et al. 2011) that has been implicated in host–pathogen interactions (Seshadri et al. 2015).

    One of the most strongly induced genes in all strains (including P. syringae pv. tomato), with a log2FC > 5, encoded a C1 family papain-like cysteine protease, which potentially regulates the activity of virulence factors. Secreted proteases can help pathogens avoid recognition by the plant immune systems, and many TTSS‐dependent effector proteins in P. syringae pv. tomato (e.g., AvrRpt2, AvrPphB, HopPtoN, and HopX) are clan CA proteases, which share the same papain‐like three-dimensional structure, order of catalytic residues, and other conserved features (Nimchuk et al. 2007; Rawlings et al. 2006; Shindo and Van der Hoorn 2008). However, the candidate we identified is unlikely to be an effector because it was also strongly induced in Psa2, in which the TTSS is inactive. Given the importance of proteases for protein quality control systems (Figaj et al. 2019), this protease may be a stress-response protein that allows P. syringae pv. actinidiae and P. syringae pv. tomato to adapt in the hostile host environment. Alternatively, the protease may be active in the apoplast (for instance, playing a role in host immunity evasion), as reported for the P. syringae pv. tomato alkaline protease AprA (Pel et al. 2014).

    P. syringae pv. actinidiae biovars also show unique responses to apoplast-like conditions.

    The microarray data also revealed many biovar-specific responses to nutrient restriction. In Psa1 and P. syringae pv. tomato, we observed a significant upregulation of chemotaxis-related genes. This contrasts with a previous study showing the repression of such genes in P. syringae pv. syringae B728a in the apoplast environment, probably reflecting differences in the media used in these experiments. Two chemotaxis pathways (che1 and che2) are required for the complete fitness of P. syringae pv. tomato but primarily function during the epiphytic phase, because the corresponding mutants cannot enter the apoplast (Clarke et al. 2016). However, this does not rule out the possibility that chemotaxis is also required for virulence, allowing bacteria to move within the apoplast to colonize the host more efficiently. The main mechanism associated with chemotactic control is flagellum-related motility, which was induced in P. syringae pv. tomato, Psa1, and Psa3 but not Psa2. Among the motile strains, Psa1 showed the most significant enrichment in flagellum-related genes, and the HKs WP_017683751.1 and WP_017682456.1 (present in all strains) were specifically induced only in this biovar. Both HKs are located in operons also containing genes related to flagellum assembly and chemotaxis. The specific induction of these HKs may account for Psa1 motility in response to minimal medium. In Psa2, the inability to induce genes encoding flagellum components may reflect the upregulation of σ factor AlgU, a repressor of flagellar gene expression in planta (Markel et al. 2016). This may avoid host immune responses because flagellum components act as pathogen-associated molecular patterns.

    In Psa3 (and P. syringae pv. tomato), we also observed the induction of genes related to iron transport and chelation; in particular, several genes involved in pyoverdine synthesis. Iron is required for many key metabolic functions in bacteria, including the tricarboxylic acid cycle, electron transport chain, and DNA synthesis (Earhart 1996), and it also induces several virulence genes in P. syringae pv. tomato (Kim et al. 2009). However, iron levels in the apoplast do not limit P. syringae pv. tomato growth (O’Leary et al. 2016) and PvdS-regulated iron-scavenging systems are not required for P. syringae pv. tomato pathogenesis (Jones and Wildermuth 2011), suggesting that iron scavenging does not account for the aggressiveness of Psa3. Indeed, as proposed by McAtee et al. (2018), iron transport-related genes, together with phosphate acquisition-related genes described above, may account for Psa3 adaptation to nutrient depletion. The induction of anmK, encoding an anhydro-1,6-muramic acid kinase involved in cell wall metabolism, may contribute to the virulence of Psa3 given its similar role in P. syringae pv. tomato (O’Malley et al. 2019).

    In Psa2, we observed the induction of genes related to branched-chain amino acid metabolism, including those required for the synthesis of core metabolic precursors such as pyruvate, acetyl-CoA, and oxaloacetate. Branched-chain amino acids are important nutrients but also act as signals that induce virulence gene expression which, in some Gram-negative bacteria, is regulated via the leucine-responsive regulatory protein Lrp (Baek et al. 2009; Kaiser and Heinrichs 2018). The Psa2 lrp gene was not modulated by minimal medium, suggesting that branched-chain amino acids regulate virulence via other pathways, as shown for Xanthomonas oryzae pv. oryzae (T. Li et al. 2019). Leucine and valine are also precursors for the synthesis of diffusible signal factors (DSFs), which act as quorum-sensing molecules. Several genes related to leucine and malonate metabolism were upregulated in Psa2, including crotonase/enoyl-CoA hydratase, malonate decarboxylase subunit γ, the malonate transporter subunit MadL, and a malonyl CoA-acyl carrier protein (ACP) transacylase. This suggests the induction of malonyl-ACP synthesis, which may feed into the fatty acid synthase (FAS) elongation cycle to generate DSF precursors (Zhou et al. 2017). Accordingly, we also observed the upregulation of the FAS gene fabG encoding β-ketoacyl-ACP reductase specifically in Psa2. Given that quorum-sensing signals have not yet been identified in P. syringae pv. actinidiae, our data provide tantalizing evidence that DSFs may fulfill this role in Psa2 (Cellini et al. 2020; Javvadi et al. 2018; Patel et al. 2014).

    The Psa2 genome is unique among the P. syringae pv. actinidiae biovars because it carries a complete metabolic pathway for the synthesis of coronatine (McCann et al. 2013) and these toxins have been synthesized in optimized medium at 18°C (Han et al. 2003). Interestingly, we did not observe the modulation of coronatine-related gene expression under our apoplast-like conditions in Psa2, although these genes were modulated in P. syringae pv. tomato. However, the corR gene, encoding a positive regulator of coronatine biosynthesis, was strongly induced by minimal medium in Psa2 but not in P. syringae pv. tomato. In a previous study, the P. syringae pv. tomato corR gene was slightly induced in both Hoitink-Sinden medium supplemented with sucrose and hrp-derepressing medium (Sreedharan et al. 2006), suggesting that induction was below the limit of detection in our experiments. Conversely, the expression of corS, encoding the cognate HK, was significantly lower in Psa2 than P. syringae pv. tomato. Both elements of the two-component CorRS system are required for coronatine biosynthesis (Rangaswamy and Bender 2000); therefore, the low expression of corS in Psa2 may be insufficient to switch on coronatine metabolism. The differential expression of corR and corS may reflect their relative positions in the corRS operon, as described for certain structured operons in other bacteria (K. Li et al. 2019; Sharma et al. 1984). Furthermore, the operon may be regulated by external factors, as reported for the temperature-dependent regulation of coronatine synthesis in the soybean pathogen P. syringae pv. glycinea PG4180, whereas coronatine production in P. syringae pv. tomato is temperature independent (Palmer and Bender 1993; Smirnova et al. 2008). Given the 99% sequence identity between the CorS proteins in Psa2 and P. syringae pv. glycinea, it is tempting to hypothesize that coronatine synthesis in Psa2 is also temperature dependent, and that our experimental conditions (28°C) do not favor this metabolic pathway.

    Genomic analysis reveals new features of the TTSS-related pathogenicity island.

    The TTSS plays a key role in Pseudomonas–host interactions, and an activating signal is required to induce the secretion of virulence effectors (Galán et al. 2014) that suppress the host immune system and allow the bacteria to establish a niche in the apoplast (Xin et al. 2018). Although the signals are poorly understood, there is compelling evidence that TTSS activation occurs upon contact with target cells (Aldon et al. 2000; Rahme et al. 1992) and that the hrp/hrc gene cluster is controlled by multiple physiological and environmental factors that are replicated by minimal medium (Rico and Preston 2008). Therefore, we anticipated the significant induction of hrp and hrc genes in P. syringae pv. tomato, as reported previously (Nobori et al. 2018, 2019; Rico and Preston 2008). We also observed the same phenomenon in Psa3, whereas the hrp and hrc genes in Psa1 were induced weakly and those in Psa2 were not induced at all. These data suggest that different biovars of the same pathovar respond in different ways to the apoplast-like environment and, therefore, may rely on different pathogenicity strategies.

    A comparison of the complete 23,721-bp hrp/hrc cluster revealed high cross-strain conservation, except for a slight divergence within the hrcU operon in Psa2 compared with Psa1 and Psa3; however, this did not appear to affect the promoter regions or gene products. The hrp/hrc cluster showed 99.94% identity between Psa1 and Psa3 and 99.46% identity between Psa2 and Psa1/Psa3, not including the gap in Psa2 corresponding to the hypervariable intergenic region. Focusing on the hrcU operon, the identity between Psa2 and the other biovars decreased to 97.6%. Overall, the hrp/hrc cluster showed a canonical structure delimitated by hrpR and hrpK in all biovars, and a similar structure in P. syringae pv. tomato (Alfano et al. 2000). However, the detection of a polycistronic transcript featuring hrpO and hrpP suggests that the current organizational model of the hrpU and hrpJ operons should be revised, and the presence of the hrp box within the coding sequence of hrpO, the first gene in the operon, suggests a novel mode of operon regulation that may be more flexible (Ishihama 2012). Our data suggest that the hrp box within the hrpO coding sequence may not necessarily regulate hrpP but could facilitate the regulation of hrpO and downstream genes by the σ factor HrpL.

    Another key feature of TTSS-containing pathogenicity islands is the presence of an exchangeable effector locus (EEL) downstream of the hrp cluster (Alfano et al. 2000). Our Psa2 and Psa3 strains featured the canonical EEL structure (including tRNALeu, queA, and tgt sequences) although two of the three Psa2 genomes (ICMP 19071 and ICMP 19072) contained a phage insertion within the EEL, supporting the observation that tRNA genes are often found at phage integration sites (Langille et al. 2010). Multiple rearrangements have been reported in Psa1 and Psa3 biovars (Poulter et al. 2018) and, accordingly, we found that the tRNALeu, queA, and tgt sequences were dispersed to other chromosomal locations in Psa1 (ICMP9617 and ICMP9853). Such rearrangements may be selectively neutral but the genome structure may affect transcriptional regulation (Poulter et al. 2018). Given that the EEL plays a role in bacterial fitness (Alfano et al. 2000), it is tempting to speculate that the lower virulence of Psa1 compared with Psa3 may, in part, reflect the loss of the canonical EEL structure.

    The modulation of positive TTSS regulators does not correlate with the induction of hrp and hrc genes in different P. syringae pv. actinidiae biovars.

    The three P. syringae pv. actinidiae biovars contained similar pathogenicity islands yet the hrp/hrc gene cluster in Psa2 was not induced. This primarily depends on HrpL, a master regulator of bacterial TTSS genes (Waite et al. 2017). The corRS operon contains an HrpL-dependent promoter and, accordingly, coronatine production is abolished in P. syringae pv. tomato hrpL mutants (Fouts et al. 2002; Sreedharan et al. 2006). Therefore, the absence of hrpL expression in Psa2 may contribute to the lack of coronatine-related gene expression described above. However, CorR is a positive regulator of the hrp/hrc regulon in P. syringae (Sreedharan et al. 2006) and a P. syringae pv. tomato corR mutant showed a reduction and delay in the expression of hrpL and an impaired hypersensitive response on Nicotiana benthamiana. One hypothesis to explain the absence of CorRS signaling in Psa2 is the suboptimal temperature, which might prevent the induction of hrpL; however, the hrpA1 promoter was not activated in Psa2 even at 24 or 18°C, thus refuting this hypothesis (Supplementary Fig. S20).

    Intriguingly, many upstream regulators of the TTSS were strongly induced in Psa2 after 4 h in minimal medium (log2FC > 1.5, adjusted P value < 0.05) but modulation in Psa1 was negligible and modulation in Psa3 was only observed after 8 h. Our analysis included AefR (Deng et al. 2009); AlgU, which regulates hrpL and hrpRS directly (Markel et al. 2016); and CsvR, an HK that cooperates with CsvS to induce hrpRS (Fishman et al. 2018). Moreover, the hrpA1 gene is required for full expression of both secreted proteins (HrpW) and components of the Hrp secretion machinery (HrcC and HrcJ) in minimal medium and in planta (Wei et al. 2000). Not only was the hrpA1 gene not induced in Psa2 but also the transcript was up to eightfold less abundant compared with the other biovars in minimal medium and up to threefold less abundant in rich medium, and was expressed at significantly lower levels than all other hrp genes located in the same operon, except hrpZ1. These data suggest that hrp/hrc cluster genes may be repressed in Psa2, due to the low expression of hrpA1. Surprisingly, hrpT was upregulated in Psa2 even though it should be dependent on HrpL binding to the hrp box controlling the hrpC operon, suggesting an alternative and independent form of transcriptional regulation. However, HrpT inhibits the expression of hrpL indirectly (Ortiz-Martín et al. 2010; Preston et al. 1998); thus, the induction of hrpT may also explain the repression of the TTSS in Psa2. Finally, it is worth mentioning that the absence of TTSS induction in Psa2 is not due to the lack of plant-derived signals, known to regulate positively TTSS expression in several phytopatogenic bacteria belonging to the P. syringae complex (Haapalainen et al. 2009; Rahme et al. 1992; Tang et al. 2006). Indeed, the addition of kiwifruit leaf extract to the minimal medium did not trigger hrpA1 promoter activity in Psa2 carrying the reporter system, whereas it strongly upregulates fluorescence increase in the other two biovars, in particular the most aggressive Psa3 (Supplementary Fig. S9). This not only confirms that plant leaf extract mimics reliably in planta conditions to study plant-regulated virulence mechanisms but also demonstrates that the lack of TTSS induction or, alternatively, TTSS repression in Psa2 occurs independently of plant signals perceived by the bacteria.

    The induction of TTSS regulators in Psa3 did not fit with the strong induction of hrp/hrc cluster genes already observed at 4 h in this biovar. However, the strong expression of these genes in rich medium indicated that they were already upregulated during overnight growth of the inoculum. Indeed, after >20 h, the rich medium is depleted of nutrients and behaves similarly to minimal medium, which may trigger TTSS-related responses to nutrient deficiency in Psa3. This is supported by the declining TTSS-related transcript levels in rich medium after 8 h. On the other hand, the massive induction of these genes in minimal medium agrees with the perception of nutrient-independent signals in the apoplast such as acidic pH triggering a more profound TTSS response. The ability of Psa3 to integrate multiple external signals may contribute to its virulence, enabling cells to rapidly deploy their resources for an attack on the host plant. The upregulation of anmK specifically in Psa3 was also notable because a nonsense mutation in anmK was recently shown to impair the TTSS in a P. syringae pv. tomato gacA mutant background . Moreover, the alignment of hrpL promoters from the different biovars revealed a single mutation in Psa3 that created a putative supplementary binding site for the transcription factor LysR (Supplementary Fig. S21). It is unclear whether this influences hrpL expression and TTSS activation in Psa3 but the presence of such a mutation possibly influencing the expression of a master regulator provides new insight into the TTSS-related virulence of Psa3.

    The role of two-component systems and c-di-GMP in biovar-dependent responses.

    Pathogenic bacteria use sophisticated two-component signaling systems in order to adapt to environmental conditions within the host and to trigger virulence effectors (Gotoh et al. 2010; Wang and Qian 2019). We investigated the potential role of two-component systems in the virulence of Psa3 but found no evidence of HK, HHK, or RR genes that were uniquely present or specifically upregulated in this biovar. However, the comparison of transcript levels revealed some HK, HHK, and RR genes that were constitutively expressed at higher levels in Psa3 (i.e., in both rich and minimal medium) but not modulated in Psa1 or Psa2, or that were induced in Psa1 but not Psa2 when switched to minimal medium. In particular, two HHK genes that are missing from Psa1 and annotated as pseudogenes in Psa3 were nevertheless strongly expressed in Psa3. One of them (WP_020312785.1) featured a 5-bp deletion in Psa3 leading to a frameshift and premature stop codon compared with its ortholog in Psa2 (Supplementary Fig. S22A and B) but the promoters were conserved in both biovars. The promoter sequence conservation and high level of expression indicates the products are likely to be functional. Apparent pseudogenes have previously been shown to produce noncoding RNAs with a role in gene expression (Pink et al. 2011) or even active products corresponding to truncated proteins (Kandouz et al. 2004; Zhang et al. 2006). The corresponding gene in Psa2 encodes an HHK similar to the Escherichia coli EvgS protein, which perceives mild acidic conditions (pH 5.5) via a cytoplasmic linker (Eguchi and Utsumi 2014). The translation of the Psa3 pseudogene in silico indicated two putative truncated proteins, the first corresponding to the periplasmic region (including the PBDb and PAS domains) and the second also containing the cytoplasmic region, including three key catalytic domains (HK, HPT, and REC), similar to the pH-sensing cytoplasmic portion of EvgS (Supplementary Fig. S22C and D). In E. coli, EvgS activates EvgA, an RR containing a helix-turn-helix LuxR DNA-binding motif (Ma et al. 2004). Accordingly, we found that a similar RR gene (WP_017682831.1) was also strongly expressed in Psa3 and was located close to the HHK-like pseudogene WP_020312785.1 discussed above, within a genomic segment absent from Psa1. Therefore, our data suggest that a truncated HHK derived from this apparent pseudogene may interact with the RR to transduce a pH-dependent signal for TTSS activation in the apoplast. This agrees with evidence that pseudogenes can be called upon as a reservoir of genetic information to deal with environmental stress or conditions that promote mutation (Harrison and Gerstein 2002; Harrison et al. 2002).

    The transduction of signals is facilitated by second messengers such as c-di-GMP and (p)ppGpp (Chatnaparat et al. 2015; Pérez-Mendoza et al. 2014). The synthesis of c-di-GMP requires a diguanylate cyclase containing a GGDEF domain whereas excess c-di-GMP is degraded by PDEs containing an EAL or HD-GYP domain. P. syringae pv. tomato produces one diguanylate cyclase and one PDE that modulate c-di-GMP levels to control virulence (Aragón et al. 2015; Engl et al. 2014). Genetic manipulation of the c-di-GMP content has shown that the accumulation of this signal promotes biofilm production while inhibiting motility and the TTSS (Wang et al. 2019). We found 46 proteins with GGDEF or EAL/HD-GYP domains potentially involved in c-di-GMP metabolism and signaling. Interestingly, GGDEF-containing proteins were mostly repressed in Psa3 but not the other biovars, whereas some putative PDEs were constitutively expressed at higher levels in Psa3 than the other biovars. Two PDE genes (WP_017682825 and WP_017682824) were located just downstream of the HHK pseudogene and cognate RR gene discussed above, suggesting that this region of the genome, missing in Psa1, is a key determinant of pathogenicity. In proteins with both GGDEF and EAL domains, one of the domains is usually inactive and fulfils a regulatory function, or a third regulatory domain is present to disjoin the activity of the GGDEF and EAL domains (Chou and Galperin 2015; Galperin 2010). We identified two hybrid GGDEF/EAL proteins that were strongly expressed in Psa3. One of them (WP_017684829) also included a GAF domain, suggesting that other small molecules such as c-di-GMP or cAMP can bind to the GAF domain and modulate c-di-GMP synthesis or degradation (Tamayo et al. 2007). The other (WP_017683620) also contained PAS/PAC and REC domains, suggesting that it may function as an RR, thus linking the activation of an RR directly to the modulation of c-di-GMP by promoting or repressing diguanylate cyclase or PDE activity (Tamayo et al. 2007). These data suggest that Psa3 tends to reduce c-di-GMP levels by suppressing its synthesis or accelerating its degradation. Given the increased virulence of Pseudomonas mutants overexpressing the PDE BifA (Aragón et al. 2015), the aggressiveness of Psa3 may also involve the hrp/hrc-dependent depletion of c-di-GMP leading to the induction of TTSS-dependent signaling (Wang et al. 2019; Xie et al. 2019).

    In Psa2, the higher c-di-GMP content may account for the inactive TTSS and overall lower virulence. The level of c-di-GMP is regulated by the chemoreceptor PscA, and P. syringae pv. tomato pscA mutants accumulate this second messenger and favor biofilm formation over swarming motility (Cerna-Vargas et al. 2019). Although pscA was upregulated in Psa2 in minimal medium at 4 h, it was still expressed at higher levels in Psa3 in both media at both time points (log2FC > 2; adjusted P value < 0.05). This low level of pscA expression may contribute to the higher c-di-GMP levels in Psa2 which, in turn, would suppress motility-related genes such as flhA, fliE, and fliN (Wang et al. 2019). Indeed, we found that these genes were expressed at lower levels in Psa2 than Psa3 (log2FC > 1.5, adjusted P value < 0.05) and, accordingly, Psa2 showed a lower swarming motility capacity compared with the other P. syringae pv. actinidiae biovars. The higher c-di-GMP level in Psa2 would also explain the greater propensity of this biovar for biofilm formation in both media compared with the other biovars.

    Taken together, our data suggest that c-di-GMP is key marker of aggressiveness in different P. syringae pv. actinidiae biovars, and acts as a switch between sessile behavior (biofilm formation) and virulence (TTSS activation). Psa3 appears to constitutively display, and likely maintain, lower c-di-GMP levels by the constitutive expression of PDE genes and the suppression of diguanylate cyclases, whereas the expression of diguanylate cyclases in Psa2 leads to the repression of TTSS responses by c-di-GMP. In line with the intermediate phenotype we observed, Psa1, which displays a higher c-di-GMP level compared with Psa3, suppresses c-di-GMP-degrading enzymes containing an EAL domain when transferred to minimal medium, thus requiring more time than Psa3 to deplete any reserves of this second messenger, leading to the weaker induction of TTSS responses compared with Psa3.

    Conclusion.

    Our multistrain transcriptome profiling strategy allowed us to identify subtle differences between closely related bacterial strains under different conditions and to analyze genes such as corR and corS in Psa2 that have yet to be annotated in particular strains. More importantly, we were able to confirm the expression of pseudogenes in Psa3 that may have been excluded from analysis based solely on single-strain gene catalogs. We focused on TTSS effectors which play a key role in the fate of bacterial interactions in susceptible and resistant plants (Dillon et al. 2019; Lee et al. 2019). Comparative genomics has previously identified strain-specific effectors that were assumed to confer virulence such as those present in Psa3 but not Psa1 (Marcelletti et al. 2011; McCann et al. 2013). However, this study highlights the importance of transcriptional analysis to confirm whether such effectors, and the corresponding TTSS components, are expressed and correlated with virulence traits. Differences in the timing of TTSS induction and secondary virulence mechanisms may provide an advantage to certain strains, increasing their aggressiveness. Plants can block the translocation of TTSS effectors although the mechanism is unclear; thus, the ability of bacteria to respond quickly to the apoplast environment would confer an advantage, allowing them to overcome these defenses (Crabill et al. 2010).

    The identification of factors that determine virulence such as specific two-component systems and corresponding signaling pathways and regulatory networks is necessary for the development of strategies to reduce the virulence of plant pathogens. Our work shows how P. syringae strains belonging to the same pathovar, infecting the same host with different degrees of aggressiveness, may display profoundly different responses to the apoplast, thus accounting for biovar diversity in terms of virulence. P. syringae pv. actinidiae is not only a devastating kiwifruit pathogen but also a very important model to study the regulation of bacterial TTSS activity. In particular, we present evidence that c-di-GMP may be a key factor controlling the aggressiveness of P. syringae pv. actinidiae biovars, and further investigations should focus on the direct targets of c-di-GMP that control the biofilm-TTSS switch.

    MATERIALS AND METHODS

    Microarray design and fabrication.

    Whole-genome sequences representing four P. syringae pv. actinidiae strains and P. syringae pv. tomato, as well as three ICE sequences, were included in the design of an Agilent custom high-density microarray chip (Agilent Design ID: 078853; GEO accession: GPL27505). At the time, complete genome sequences were available for Psa3 ICMP18884 and P. syringae pv. tomato DC3000, and partial sequences were available for the others. The sequences and gene annotation data were collected in July 2015 from the NCBI GenBank and RefSeq databases and the University of Udine (CRAFRU8.43 annotations from Prof. Giuseppe Firrao). If multiple genome annotations were available for a given strain, a unique reference transcriptome was created by merging the data to avoid artificial redundancy at the strain level, giving the following priority for each transcript: RefSeq > GenBank > University of Udine (Pruitt et al. 2005). Nonredundant sequences annotated for each strain were then grouped by collapsing identical sequences into one representative to create a multistrain, nonredundant collection of protein coding sequences (CDSs) covering all annotated transcripts for all strains (n = 20,561 total CDSs). The Agilent eArray web tool was used to design probes based on genomic sequences representing each gene family using the following parameters: method = Tm Matching Methodology, probe length = 60 bp, number of probes per target = 3, transcriptome details = “use target file as transcriptome”, and probe design = “design without 3′ bias”. Gene sequences were annotated using the Blast2GO suite v4.1.9 (Conesa et al. 2005) after mapping against the NCBI nonredundant database with BLASTX v2.6.0 (Altschul et al. 1990). The eArray design application yielded 18,598 high-performance probes (lowest possible target ambiguity) interrogating 20,554 protein sequences, with 14,457 unambiguous probes matching a unique target. The remaining set of probes targeting a 60-mer string shared by two or more transcripts (here dubbed “ambiguous probes”, n = 4,141) reflected residual transcript redundancy retained in our comprehensive microarray design originating from both biological (multiple strains sharing related sequences) and artificial (multiple collections of genome annotation data) sources. On the other hand, we favored full representation of annotated sequences to minimize hybridization biases due to strain-specific variation in transcript sequences and to conservatively treat inherent uncertainty of annotated sets of transcripts.

    The design procedure for the multistrain microarray is summarized in Supplementary Figure S1.

    Bacterial strains and growth media.

    Single colonies of the Psa3 strains CRA-FRU 8.43 and ICMP18884/V-13, the Psa1 strain J35, the Psa2 strain KN.2, and the P. syringae pv. tomato strain DC3000 grown on rich solid medium (KB agar) were inoculated into KB medium and incubated overnight at 28°C, with shaking at 200 rpm. When cells reached the late log phase, they were collected by centrifugation (5,000 × g, 10 min, room temperature), washed three times in liquid KB medium (King et al. 1954) or HIM (Rico and Preston 2008), then resuspended at a final density of 2 × 108 CFU/ml in HIM or 2 × 107 CFU/ml in KB and incubated for 4 or 8 h at 28°C, with shaking at 200 rpm. The cells were harvested by centrifugation as above and the pellets (approximately 2.4 × 109 cells) were stored at −20°C. The procedure was carried out twice to obtain three independent biological replicates of each sample for microarray analysis.

    RNA extraction and microarray chip hybridization.

    RNA was extracted from each sample using the Spectrum Plant Total RNA Kit (Sigma-Aldrich) and quantified using a Nanodrop spectrophotometer (Thermo Fischer Scientific). RNA quality was evaluated using an Agilent RNA 6000 Nano Kit Bioanalyser. RNA was processed and labeled for microarray analysis using the One-Color Microarray-Based Gene Expression Analysis Low-Input Quick Amp WT Labeling kit (Agilent Technologies), according to the manufacturer’s instructions.

    Microarray data analysis.

    The fluorescence intensity for each probe was measured using an Agilent G4900DA SureScan Microarray Scanner System with Agilent Scan Control software, and data were extrapolated using Agilent Feature Extraction software. Fluorescence intensities were calculated by robust multiarray averaging, including adjustment for the background intensity, log2 transformation, and quantile normalization. Normalized data were used to identify DEGs with threshold values of P < 0.05 and log2FC values > |1|. DEGs were compared across different strains or conditions using the online software Calculate and Draw Custom Venn Diagrams. Annotation of gene sequence GO terms and GO functional enrichment analysis was carried out using Blast2GO v4.1.9 (FDR < 0.01).

    Gene sequence clustering for interstrain comparisons.

    A nonredundant dataset of amino acid sequences (n = 17,047) representing all proteins annotated from RefSeq or GenBank collected in this study (Supplementary Fig. S1) was created by taking each identical sequence only once. To allow interstrain comparisons of protein sets, we grouped the different proteins annotated for each strain of interest from different sources (Supplementary Fig. S1) based on sequence similarity and then analyzed the presence or absence of the given protein cluster across Pseudomonas strains.

    Grouping of proteins was performed in two steps based on results of all-against-all alignments of amino acid sequences performed with the BLASTP software, version 2.7.1 (Altschul et al. 1990).

    In the first grouping step, any pair of protein sequences scoring as reciprocal best hits in the BLASTP output was grouped in the same cluster. To assess the quality of this procedure, we compared clusters of proteins called in our analysis against a set of nonredundant RefSeq protein records annotated for P. syringae (marked by protein identifier WP_XXXXXXX.1, where XXXXXXX is a unique numeric code) and found that all such instances annotated in multiple strains were correctly assigned to the same cluster of proteins.

    Next, we further grouped clusters of protein sequences marked by high sequence similarity among cluster members (i.e., E < 0.001; at least 94% identity over at least 90% of the aligned sequences in the BLASTP output) in superclusters.

    The above procedure assigned the full and redundant collection of protein sequences (n = 56,042 FASTA sequences) collected for different P. syringae strains from different annotation sources (Supplementary Fig. S1) to a total of 10,839 protein groups. These protein groups were mapped to the 18,598 probes of the P. syringae multistrain microarray to allow comparative transcriptomic investigation.

    Identification of HK, HHK, and RR sequences and genes related to c-di-GMP metabolism.

    HK, HHK, and RR genes were retrieved using MiST v3.0 (Gumerov et al. 2020). Genes containing GGDEF or EAL domains were retrieved by using BLASTP to search annotated proteins in the P. syringae pv. actinidiae genomes with representative domain sequences from the RCSB Protein Data Bank: protein O35014 for the EAL domain and protein B8GZM2 for the GGDEF domain. Genes containing an HD-GYP domain were similarly identified by screening with domain ID COG2206 from the Conserved Domains and Protein Classification databank. Proteins that were not common to all three P. syringae pv. actinidiae biovars were used as BLASTP queries to search the proteins annotated in other P. syringae pv. actinidiae genomes to confirm biovar specificity, including biovars not included on the microarray (ICMP9853 for Psa1, ICMP19071 and ICMP19072 for Psa2, and CRAFRU8.43 for Psa3). A similar BLASTP search against the P. syringae pv. actinidiae genomes was carried out for all HK, HHK, and RR genes previously identified in P. syringae pv. tomato (Lavín et al. 2007). Finally, for proteins still absent from one or more P. syringae pv. actinidiae strains, the probes matching the corresponding genes were used as BLASTN queries against the complete P. syringae pv. actinidiae genomes to account for potential errors in genome annotation and the presence of pseudogenes, thus verifying the presence or absence of the corresponding genes at the nucleotide level.

    Phylogenetic analysis.

    For each protein family, an unrooted phylogenetic tree was built using the phylogeny.fr “a la carte” pipeline as previously described (Ariani et al. 2017), with slight modifications. Phylogenetic trees were built using protein sequences and “alignment curation” was avoided in the pipeline due to the high variability among protein sequences. Results were downloaded in Newick format for tree representation using the Interactive Tree of Life (iTOL) v3.2.4. SMART (Letunic et al. 2002; Schultz et al. 1998) was used to identify domain architecture and protein motifs with default parameters. The dataset was uploaded for the visualization of protein classification (HK, HHK, GGDEF, EAL, and HD-GYP) and the presence of other domains using iTOL.

    Hierarchical clustering.

    Normalized fluorescence intensities from microarray experiments were imported as data matrices into MeV (Howe et al. 2011). The data were adjusted as median center genes or rows and clustered using the hierarchical clustering module. Gene and sample trees were clustered with optimized gene and sample leaf orders using Pearson correlation and average linkage clustering. The trees were subsequently cut into clusters using a distance threshold (0.5 to 1) empirically adjusted to highlight the most relevant features of the trees.

    Comparison of hrp cluster sequences.

    The complete sequences of the hrp cluster (hrpRhrpK) were retrieved for P. syringae pv. actinidiae strains ICMP9617, ICMP19073 and ICMP18884 in the Pseudomonas Genome Database (Winsor et al. 2016) and multiple sequences were aligned using Geneious v11.0.3.

    Colony-based c-di-GMP reporter assays.

    Fluorescence intensity analyses using the c-di-GMP biosensor pCdrA::gfpC were carried out according to methods described previously by Cerna-Vargas et al. (2019), with slight modifications. Briefly, the reporter plasmid pCdrA::gfpC was transformed into Psa1, Psa2, and Psa3 strains by electroporation. Strains growing on selective medium containing antibiotics were further confirmed as positive using colony PCR. For the analyses on solid medium, overnight bacterial cultures of the positive strains grown in KB medium were adjusted to an optical density at 600 nm (OD600) of 0.8 and 20-µl drops were spotted on KB agar plates. After 48 h of incubation at 28°C, fluorescence intensity was visualized with a Leica DFC7000 T microscope (Leica Microsystems, Wetzlar, Germany). Fluorescence was visualized employing a GFP filter set (emission/excitation filter: 535/485 nm). Pictures were taken using Leica Application Suite software V4.13. Alternatively, for fluorescence measurement in liquid medium, overnight bacterial cultures of the positive strains grown in KB medium were washed twice with HIM and resuspended in HIM at an OD600 of 1, and 200-µl aliquots were inoculated in microtiter plates. Fluorescence intensity was measured using the fluorometer plate reader Tecan Infinite M Plex (Tecan Trading AG, Switzerland).

    Biofilm assay.

    Biofilm formation was quantified in Psa1 (J35), Psa2 (KN.2), and Psa3 (CRAFRU8.43) as previously described (O’Toole 2011), with slight modifications. Bacterial cells were grown overnight in liquid KB medium before washing and resuspending the cells in KB medium or HIM as described above. The OD was adjusted to OD600 = 0.2 (KB medium) or OD600 = 1 (HIM) and 100-μl aliquots of each culture were incubated at 28°C in a transparent 96-multiwell plate for 24 h, with shaking at 180 rpm. Cell-free media were used as negative controls. At the end of the incubation period, the microtiter plates were washed with tap water and 150 µl of 0.1% crystal violet was added to each well and incubated at room temperature for 30 min. The plates were rinsed three to four times by submerging them in water, then blotted carefully on a stack of paper towels to remove excess cells and dye before drying for 30 min. To solubilize the crystal violet, we added 125 µl of 30% acetic acid to each well and incubated them at room temperature for 30 min. Absorbance was measured in a plate reader at 590 nm using 30% acetic acid in water as the blank.

    Swarming motility assay.

    P. syringae pv. actinidiae strains belonging to different biovars were grown at 28°C for 24 h on KB agar. Cells were resuspended in KB medium to an OD600 of 0.8. Bacterial suspension (5 µl) was spotted onto soft nutrient broth yeast extract agar (0.4% [wt/vol] agar). Plates were incubated at 28°C under dark conditions. Swarm colonies were photographed after 48 h of incubation.

    RT-PCR.

    Total RNA was treated with TURBO DNase (Thermo Fisher Scientific) and first-strand cDNA was synthesized with SuperScript III Reverse Transcription (Invitrogen) according to the manufacturer’s instructions. The cDNA was amplified by PCR using GoTaq G2 DNA Polymerase (Promega Corp.) with the following parameters: initial denaturation at 95°C for 2 min; followed by 30 cycles of denaturation at 95°C for 30 s, annealing at primer-specific temperature (Supplementary Table S2) for 30 s, and extension at 72°C for primer-specific duration (Supplementary Table S2); followed by a final fill-in reaction at 72°C for 7 min. Genomic DNA was used as a positive control. The PCR products were separated by 2% agarose gel electrophoresis in 1× Tris-acetate-EDTA buffer. The primer sets are listed in Supplementary Table S2.

    ACKNOWLEDGMENTS

    We thank M. Scortichini for providing the P. syringae pv. actinidiae strains, G. Firrao for providing genomic information about strain CRAFRU8.43, N. Vitulo for providing bioinformatics support, and Regione Veneto for funding.

    AUTHOR-RECOMMENDED INTERNET RESOURCES

    Blast2GO Pro: https://www.blast2go.com

    Conserved Domains and Protein Classification databank: https://structure.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml

    Draw Venn Diagram: http://bioinformatics.psb.ugent.be/webtools/Venn

    eArray web tool: https://earray.chem.agilent.com/earray

    Geneious v11.0.3: https://www.geneious.com

    Interactive Tree of Life v3.2.4: https://itol.embl.de

    NCBI nonredundant database: https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins

    phylogeny.fr database: http://www.phylogeny.fr

    RCSB Protein Data Bank: http://www.rcsb.org

    Ref-Seq databases: https://www.ncbi.nlm.nih.gov

    The author(s) declare no conflict of interest.

    LITERATURE CITED

    The author(s) declare no conflict of interest.

    Funding: Regione del Veneto “Progetto di innovazione per la difesa della pianta del kiwi e per la valorizzazione dei suoi frutti” and LR number 32 del 9/08/1999–art. 4 “Ricerca di interesse regionale e sperimentazione” (CUP) number H16D14000090002.

    Current address for A. Regaiolo: Johannes Gutenberg-Universität Mainz, Institut für Mikrobiologie und Weinforschung, Mainz, Germany.