Transcriptome and Comparative Genomics Analyses Reveal New Functional Insights on Key Determinants of Pathogenesis and Interbacterial Competition in Pectobacterium and Dickeya spp.

Transcriptome sequencing of Pcb1692 during in planta infection.Pcb1692 has been reported as one of the most aggressive Pectobacterium species known to date (17, 38, 39). Aiming to examine the transcriptional landscape of Pcb1692 during infection in potato tubers, a whole-transcriptome data set was generated, including samples harvested 24 and 72 h postinfection (hpi) along with an in vitro control (see Materials and Methods). The data set comprised 14 to 20 million transcriptome sequencing (RNA-Seq) paired-end reads for each stage, with 97% of the mapped reads in each sample uniquely mapped to the reference genome, implying good overall quality (Table 1). Subsequent analyses identified 1,743 protein-coding genes (43.5% of total annotated in Pcb1692) under infection-induced regulation (log₂ fold change of >1 or <−1; false discovery rate [FDR], <0.01) in the wild-type strain in at least one time range (see Table S1.1 in the supplemental material). Importantly, a recent study depicted D. dadantii’s transcriptome during infection on Arabidopsis thaliana, compared at 6 and 24 hpi, as being able to detect 13.5% of its protein-coding genes (575 out of 4,244) under regulation (up or down) (26). We take advantage of this recently published gene expression experiment, featuring a closely related wild-type SRE strain infecting a non-crop host, to comparatively examine up/downregulation in our original gene expression data set obtained from Pcb1692.

SRE’s most well-recognized strategy to accomplish plant tissue invasion involves disruption of plant cell wall by deploying a variety of pectinase-, cellulase-, and protease-related protein families (16). PLs and PNL are pectinases that participate in this process by breaking α-1-4-linkages of homogalacturonan backbone in the pectin molecule via a transelimination mechanism (40). The data set presented here emphasizes the role of PL in Pcb1692 infection by the significant activation of seven out of eight genes, occurring simultaneously during the first 24 hpi. This encompasses orthologs of well-documented genes (pelABCINZ, PCBA_RS04055, PCBA_RS04060, PCBA_RS04065, PCBA_RS10065, PCBA_RS03200, and PCBA_RS04070) and one unannotated gene (PCBA_RS18630) (Table S1.1). These findings are similar to those detected in D. dadantii, in which all nine PLs (pelABCDEILNZ) were activated while infecting A. thaliana (26). The pectin lyase gene pnl (PCBA_RS19200), which encodes a major PCWDE (24), is massively activated (>5 log₂ fold change), ranking among the top 1% most upregulated genes 24 hpi in Pcb1692 (Table S1.2). This finding suggests a specific high transcriptional demand of pnl in Pcb1692 during early disease development. Conversely, in the D. dadantii transcriptome, the pnl ortholog (Dda3937_03551) surprisingly shows no detectable regulation (Table S1.3) (26). Together, these results strengthen the notion of diversity in PCWDE pools transcriptionally activated by pectinolytic pathogens to surpass various obstacles imposed by hosts in distinct plant tissues (41). Whereas for some PCWDE families (e.g., PLs) the upregulation is basal, regardless of host characteristics, other PCWDE families (e.g., PNL) are transcriptionally activated under specific host-imposed cues.

Transcriptional profile and conservation of type VI secretion system and associated effectors in SRE.The T6SS is a bacteriophage-like structure that depends on the expression of 13 core genes (tss cluster) to assemble a complex that spans the cell envelope. In addition to the tss cluster, a hemolysin-coregulated protein (Hcp) forms the T6SS tube structure, which is propelled across the membranes with a piercing structure on its tip, consisting primarily of two proteins (VgrG/PAAR), into the target cell (42, 43). Furthermore, it has also been demonstrated that this piercing structure is able to accommodate independent proteins which function as effectors (34). Our data set underscores the importance of T6SS by showing the overwhelming activation of all 13 genes in the tss cluster (>3-log₂ fold change in transcriptional variation during the first 24 hpi). Moreover, a total of four Hcp secretion islands (HSIs), one adjacent to the tss cluster (HSI-1) and another three in different genomic contexts (HSI-2, -3, and -4), were also upregulated for at least one time point. These HSI exhibit an overall activation of two PAAR-coding and four vgrG genes (Table S2.1). In comparison, wild-type D. dadantii displays similar upregulation of most tss and four hcp genes out of six conserved in the genome upon infection on A. thaliana (Table S1.3). However, none of the VgrG- or PAAR-encoding genes (five and three, respectively) are transcriptionally regulated in D. dadantii in the same experiment (26). Notably, the role of VgrG-PAAR as a delivery mechanism for Rhs effectors into target cells is established in D. dadantii facing in vitro competition (44). Thus, although structural components of T6SS sheath and tube are similarly regulated when infecting distant hosts in Pcb1692 and D. dadantii, genes encoding the piercing tip components may be activated under specific cues (26) (Table S1.3).

Interestingly, the second major HSI in the Pcb1692 genome (HSI-2) spans ∼6.6 kb, harboring eight genes under positive regulation for at least one time point during infection (Fig. 1A). Of these, six (excluding only hcp and vgrG) are uncharacterized/unannotated (Table S2.1). Aiming to garner deeper functional insight on this region, we combined the analysis of domain architectures conserved in protein sequences with contextual genomic information from 100 SRE organisms. First, conserved domain analysis confirmed VgrG (PCBA_RS05800; Pfam entry PF05954) and PAAR gene products (PCBA_RS05770; Pfam entry PF05488) in this region, implying that the cluster is able to source both the Hcp tube and the piercing tip (VgrG/PAAR) assembly (Fig. 1B). Further, two uncharacterized genes (PCBA_RS05780 and PCBA_RS05785) flanked by the genes encoding VgrG and PAAR caught our attention for containing domain architectures SMI1_KNR4 and WHH (Pfam entry PF14414 and PF09346), respectively, in their products. Interestingly, the association between the superfamilies SUKH and HNH/ENDOVII, which encompass SMI1_KNR4 and WHH families, respectively, has been recognized as a recurrent theme in bacterial genomes (14). This duplet (WHH-SMI1_KNR4) was described as an immunity/toxin pair associated with contact-dependent growth inhibition (CDI) systems (e.g., T5SS and T6SS), although evidence actually linking these genes to the secretion systems remains scarce (14, 45–47). The WHH family was originally identified as a restriction endonuclease derived from the HNH domain (Pfam entry PF01844; CL0263) (48). Our analyses revealed that members of the WHH family are present in only 18 out of 100 SRE strains. In 3 of these 18 strains (Pectobacterium parmentieri strains CFIA1002 and RNS08421a and P. carotovorum subsp. brasiliense CFIA1033), the WHH-containing gene was found to be duplicated (Table S2.2). Out of the 18 genomes conserving WHH-containing genes, tight association to SMI1_KNR4 family members immediately downstream in 13 genomes was found, which corroborates the original report for this duplet (Fig. 1B and Table S2.2) (14). Analysis of the 13 genomes conserving this association (WHH-SMI1_KNR4) in SRE determined that one was not suitable for contextual genomic inspection in this specific region due to incompleteness of genome assembly (namely, P. betavasculorum strain NCPPB 2793). By assessing the remaining suitable structures conserving the WHH-SMI1_KNR4 duplet, we found 83% (10 out of 12) linkage with upstream HCP, and PAAR-encoding genes were mostly downstream (Fig. 1B). A duplication in the immunity gene encoding SMI1_KNR4 was also found in one strain, which corroborates the original report of this family (14) (Fig. 1B). These results imply a strong association of the WHH-SMI1_KNR4 duplet with HSIs in SRE genomes, such as P. carotovorum subsp. brasiliense (five strains, including Pcb1692), D. dadantii (two strains), Dickeya zeae (one strain), Dickeya chrysanthemi (one strain), and P. betavasculorum (one strain) (Fig. 1B and Table S2.2). These observations, taken together with the coordinated upregulation at the transcriptional level of WHH-SMI1_KNR4-encoding genes, along with all other HSI-2 neighboring elements in Pcb1692, suggest the T6SS-dependent secretion of these genes (Fig. 1). In this context, we provide a report of coordinated transcriptional regulation of a SUKH-1/HNH-ENDOVII system along with the surrounding HSI, in which the WHH-containing gene product may work as a T6SS effector recruited during infection.

• Open in new tab
• Download powerpoint

FIG 1

Transcriptional variation in planta of T6SS-related genes in Pcb1692 and linkage between the endonuclease WHH and HSI in SRE genomes. (A) Each row on the heatmap represents a gene labeled after its (i) respective encoded protein name, (ii) functional domain, or (iii) locus tag given annotation availability. The two columns of the heatmaps represent pairwise comparisons between samples; thus, each cell shows the variation in transcription found in each comparison, illustrated by a color scale ranging from −5 (dark gray) to 5 (red), representing the log₂ fold change in expression. For each gene label, statistical significance of regulation is represented with the following colors: red, upregulation; blue, downregulation; gray, no support for regulation. Column labels: Cwt, control wild type; Wt-24hpi/Wt-72hpi, wild type at 24/72 hpi in potato tubers. Four segments in the heatmap separate T6SS-related clusters, type VI secretion (Tss) and Hcp secretion island (HSI), located in different genomic regions, namely, Tss+HSI-1, HSI-2, HSI-3, and HSI-4. (B) Gene neighborhood of WHH-SMI1_KNR4 duplets in the 10 SRE strains conserved within HSIs. Among the 100 SRE strains analyzed, those exhibiting linkage between WHH-SMI1_KNR4- and Hcp/PAAR-encoding genes are represented. Each gene arrangement is associated with its respective strain(s), represented on the left side by 4-letter abbreviations: Dzea, Dickeya zeae; Pcb, Pectobacterium carotovorum subsp. brasiliense; Ddad, Dickeya dadantii; Pbet, Pectobacterium betavasculorum; Dchr, Dickeya chrysanthemi. For each species, the respective number of strains conserving that particular gene arrangement is represented in gray ellipses on the left. The WHH-SMI1_KNR4 duplets (Pfam PF14414 and PF09346, respectively) are highlighted in light blue. Hcp (T6SS_HCP; Pfam PF05638), VgrG (Phage_GPD; Pfam PF05954), and PAAR (PAAR_motif; Pfam PF05488) genes are highlighted in dark blue. Syntenically conserved blocks across the different gene arrangements are linked by light-brown areas. Duplicated genes in a single arrangement are linked by green shapes. Successfully clustered products in orthologous groups for which no conserved domains were detected are highlighted in dark blue.

Analysis of antibacterial activity of type VI secretion system-related toxins.Aiming to investigate the potential of HSI-borne genes from Pcb1692 to have antibacterial activity, four relevant genes were cloned into an expression vector and ectopically expressed in Escherichia coli, as described by Koskiniemi et al. in 2013 (44). To this end, the previously described WHH-containing gene product (PCBA_RS05785) was selected alongside three other putative toxins (PCBA_RS05775, PCBA_RS18045, and PCBA_RS22965) for this analysis. PCBA_RS05775 belongs to the cell division cycle protein 123 family (D123; Pfam entry PF07065), which seems not to be accompanied by an adjacent immunity-encoding gene (Table S2.3). The PCBA_RS18045 gene is located in HSI-1, within the Pcb1692 tss gene cluster, and possesses a PAAR domain followed by an alpha/beta hydrolase fold with a GxSxG catalytic lipase motif (Pfam entry PF12697), characteristic of phospholipases (49). This phospholipase family member is followed by a putative immunity gene downstream that carries an Ank domain (Pfam entry PF00023). The PCBA_RS22965 gene has a conserved HNH/EndoVII fold derivative known as AHH (Pfam entry PF14412). Similar to WHH family members, the AHH proteins were predicted to be functionally associated with CDI systems in the original report of the family (14). By scanning 100 SRE genomes, we found that 23 of these conserve AHH-encoding genes, out of which one (Dickeya zeae strain CSLRW192) harbors two copies of the gene (Table S2.3). In most SRE genomes harboring members of the AHH family, no detectable domain in the downstream immunity protein can be observed. In Pcb1692, similar to six other SRE genomes conserving AHH-encoding genes, a contig break occurs up to five genes prior to this AHH member, hindering direct gene neighborhood-based conclusions. Nonetheless, in other strains from P. carotovorum subsp. brasiliense, P. carotovorum subsp. carotovorum, and P. atrosepticum species conserving 70% to 95% of overall genomic synteny with Pcb1692 (Fig. S1), the AHH genes are solidly neighbored by hcp, vgrG, and a DcrB domain (Pfam entry PF08786)-encoding gene, suggesting a role as a T6SS-associated effector (Table S2.3). Such an arrangement precisely matches the previously assessed HSI-4 region in Pcb1692 (Fig. 1A). In addition to the original report of the AHH family, these results suggest that the AHH member plays a role in Pcb1692 infection as a T6SS-associated nuclease.

The results show that no cell death occurred due to expression of three of the four effectors between 60 and 240 min (Fig. 2A). However, expression of the AHH effector caused a reduction in E. coli growth, indicating that it has a toxic effect on E. coli. This growth inhibition could be counteracted by coexpression of the toxin and immunity gene (AHH+i) in E. coli, suggesting that the previously observed reduced growth was due to the expression of the AHH nuclease (Fig. 2B). Thus far, no experimental data exist for the AHH and WHH nucleases, as well as the D123 protein, as T6SS effectors. This report elucidates the antibacterial activity of an AHH-containing protein.

• Open in new tab
• Download powerpoint

FIG 2

Growth of E. coli DH5α expressing effectors from pCH450 induced with 0.2% l-arabinose after 30 min of growth or effector and immunity protein from pCH450 and pTrc100, respectively. (A) Expression of phospholipase, WHH nuclease, and D123 effectors in E. coli shows growth comparable to that of the empty vector control (p450). Individual expression of these effectors does not contribute to bacterial growth inhibition under the conditions analyzed. (B) Expression of AHH nuclease decreases the growth rate of E. coli. Expression of both AHH and its immunity protein (AHH+i) negates the toxic effect of the effector alone; however, after 240 min a plateau in growth is reached due to protein overexpression.

Several T6SS effectors targeting eukaryotes, bacteria, or both have been identified thus far (25, 50, 51). Here, we have demonstrated that one of the four effectors identified inhibits growth of E. coli cells. It is possible that, for the other effectors tested, the in vitro conditions used were not conducive to growth inhibition. We may also need to consider that our hypothesis that the T6SS of Pcb1692 is mainly antibacterial is not entirely correct. Bernal et al. (52) have noted that, thus far, no plant-targeting effectors have been identified as T6SS dependent. As the D123 protein has no reported function and lacks a downstream immunity protein, it could be a candidate for a plant-targeted effector. Since some effectors have interkingdom activity, it may be necessary to reevaluate whether the effectors we have identified as putative antibacterial effectors have key functions within the eukaryotic host. Moreover, it must also be taken into consideration that T6SS could have functions other than contact-dependent antibacterial competition, which may favor conservation of effectors involved in other biological programs that will enforce successful host colonization.

Genomic conservation and transcriptional regulation of prophages during in planta infection.The impact of HAIs originating from prophages in bacterial pathogenic behavior is pervasive for both plant and animal pathogens (53–55). Importantly, selective advantage conferred by prophage-borne genes to plant pathogens has been associated with the activity of bacterial secretion systems, as reported in Pseudomonas syringae and Ralstonia solanacearum (35, 56, 57). In this context, by using a combination of relevant genome-wide approaches (see Materials and Methods), we detected two prominent prophage regions in Pcb1692 exhibiting significant activation at the transcriptional level during infection (Fig. 3A). The first region encompasses 20 genes in the Pcb1692 genome strongly upregulated in the first 24 hpi (Fig. 3A). Moreover, 90% of these genes are syntenic with genomic segments in at least 54/60 other Pectobacterium strains (Fig. 3B and Table S3.1). Interestingly, a genomic segment in P. carotovorum subsp. carotovorum similar to the one found in Pcb1692 was originally described as being of prophage origin, regarded as carotovoricin (Ctv) consisting of 19 genes (58) (Table S3.2). The report of this bacteriocin has been recently corroborated by in silico analyses in the same species, even extending the putative range of the cluster to 22 genes in P. carotovorum subsp. carotovorum (termed PecaPC1-p2) (Table S3.3) (35). In this context, we also detected in P. atrosepticum a segment of 12.7 kb containing 11 gene products highly similar to Ctv sequences (ECA_RS18415-ECA_RS18485). Intriguingly, this segment is located within a previously characterized bacteriophage spanning ∼36 kb, designated ECA41 (59). In addition, ECA41 contribution to virulence toward potato hosts was reported, although the mechanism of putative effectors’ action remains unknown (59). Curiously, those three reports describing Ctv, PecaPC1-p2, and ECA41 occurred without knowledge convergence based on the detectable homology linking the regions in P. atrosepticum and P. carotovorum subsp. carotovorum, which our analyses now reconcile under the same Ctv-homologous root. This preliminary observation, combined with interpretations of the above-mentioned reports, may provide a robust background for future studies on the role of carotovoricin in pathogenesis.

• Open in new tab
• Download powerpoint

FIG 3

Transcriptional variation in Pcb1692 and genomic range of Ctv and PcbPr1 insertions in four Pectobacterium strains. (A) The heatmaps show transcriptional variation profiles of T1SS+Ctv block and serralysin/inhibitor genes (top) and from a 7-kb segment of PcbPr1 (bottom). The heatmap is represented as described for Fig. 1. For Ctv and PcbPr1, contig coordinates are represented in parentheses immediately below labels, which corresponds to panel B in this figure. Genes are preferentially represented by the domain architecture or by the respective orthologous group (OG_#). na, none of the previously described labels is applicable. (B) Drastically contrasting degree of conservation between Ctv and PcbPr1 is represented in four selected Pectobacterium strains. Genomic regions from P. carotovorum subsp. carotovorum strain PC1 (light purple) and BCS2 (light green), P. atrosepticum strain SCRI1043 (light orange), and Pcb1692 (light blue) are represented in circular ideograms. The next two consecutive inner radii show the degree of conservation of the regions, represented by bars. The height of the bars in each radius represents the frequencies of BLASTP (protein sequence similarity; purple/green color range) and MCScanX (syntenic conservation; red/blue color range) hits for each gene product in the respective strains (PC1, BCS2, SCRI1043, and Pcb1692) compared pairwise to those of the remaining 99 SRE strains. The color range key used in the bars represents low and high conservation in BLASTP and MCScanX comparisons. Genes with 30 or fewer positive hits (out of 99) are shown in purple and red, respectively, in each radius. Genes with more than 30 positive hits are shown in green and blue. The links binding genomic regions represent the insertion range of prophages for which conserved synteny could be detected (Ctv, orange; PcbPr1, pink). The green outer radius adjacent to the ideograms represents the range of BLASTP similarity to the original Ctv described by Yamada et al. (58) (Table S3.2). (C) Gene neighborhood screening of T1SS+Ctv block. Genes were separated into plant cell wall degradation (PCWD; light red), T1SS (blue), and Ctv (light green) groups or a group unrelated to these classes (gray). The genomic organization found in Pectobacterium genomes is to the right of the names for species conserving arrangements abbreviated as described in Table S8. To the left of each species name is the number of strains harboring the particular arrangement.

Furthermore, we observed that Ctv is neighbored by a complete set of T1SS genes in P. carotovorum subsp. brasiliense (PCBA_RS08610-PCBA_RS08620) immediately upstream (Fig. 3C and Table S3.4). In addition, conserved domain analyses revealed two genes in the Ctv cluster containing SLT (Pfam entry PF01464) and Gp37 (Pfam entry PF09646) architectures, both reportedly conserved in bacterial virulence factors (60, 61). Therefore, the presence of predicted effector functionalities in these genes, along with a known secretion system, prompted us to assess the degree of conservation of the T1SS+Ctv structure. The results support a remarkable conservation of this genomic architecture, displaying a prominent block of at least 13 genes in Pectobacterium genomes. In contrast, the T1SS+Ctv genomic association is consistently absent from all 39 Dickeya genomes analyzed (Table S3.4). The 13 genes in the T1SS+Ctv gene block include (i) two genes encoding a peptidase/inhibitor pair, (ii) three T1SS components, and (iii) eight Ctv genes under solid conservation in 96.7% of Pectobacterium genomes (59 out of 61). Importantly, this high conservation trend seems to occur regardless of the overall synteny between these genomes (Fig. S1). The T1SS is recognized for its relative simplicity, composed of three core proteins which form a tunnel-like structure in the inner membrane, enabling molecule transfer from cytosol to the extracellular space fueled by ATP hydrolysis (62). Although typically regarded as a signal-independent system, there is a known group of products exported in a T1SS-dependent manner that contains N-terminal signals. This group encompasses bacteriocins and microcins (63–65). Since Ctv constitutes a bacteriocin system, an additional layer of support to the T1SS-Ctv functional association could be added by searching for signal peptides within Ctv sequences. Given the generally poor conservation described in T1SS signal peptides, sensitive search setups were carried out (66, 67). The predictions from two different methods corroborate the possible T1SS/Ctv association by the combined detection of N-terminal signal peptides in five proteins encoded by genes in this region. In addition, some genes lacking known conserved domains, for which the protein sequences were clustered in orthologous groups (see Materials and Methods) OG_2978, OG_2823, and OG_2853, are collectively conserved in 96% to 100% of Pectobacterium genomes (Fig. 3C and Table S3.4). Curiously, these highly conserved genes in Pectobacterium encode small products (up to 109 amino acids) and could undertake a role in bacterial competition as small antimicrobial peptides. Here, we unraveled the strikingly dominant theme encompassing T1SS and Ctv clusters in Pectobacterium genomes tied by different sources of evidence. This linkage suggests either that Ctv is an addictive selfish element colonizing these genomes or that, similar to other reported T1SS/bacteriocin associations, this system was recruited in the Pectobacterium lineage to export Ctv-borne products either toward or through the cell membranes.

The second prophage, which apparently has not been described thus far, has low similarity to other SRE genomes, being mostly confined into P. carotovorum subsp. brasiliense strains (here referred to as PcbPr1) and one P. carotovorum subsp. carotovorum strain (Fig. 3B). BLASTP searches of PcbPr1 sequences resulted in 5 proteins out of 64 (7.8%) returning matches against at least 90 SRE species (out of 99) (Fig. 3B). In contrast, overall 82.3% of Pcb1692 protein-coding sequences (3,376 out of 4,099) successfully match at least 90 SRE (Fig. S2). The low conservation of PcbPr1 is also observable in terms of genomic organization, as 93.7% of the whole PcbPr1 structure (60 out of 64 genes) conserves synteny with less than 10% of all Pectobacterium genomes analyzed (Fig. 3B and Fig. S2). Despite the low degree of conservation, PcbPr1 harbors an internal segment spanning 7 kb, in which 88.2% of the genes (15 out of 17) are cohesively upregulated over the first 72 hpi (Fig. 3A and Table S4). Importantly, 12 out of 17 genes in this region encode proteins for which known conserved domains could not be detected. In addition, most of these gene products were clustered in small orthologous groups with sequences from other P. carotovorum subsp. brasiliense and P. carotovorum subsp. carotovorum isolates, implying lineage-specific origin (Fig. 3A). This observation supports the previous synteny analysis, reinforcing the possibility of PcbPr1 comprising a lineage-specific insertion in P. carotovorum subsp. brasiliense and P. carotovorum subsp. carotovorum. Notably, genomic conservation of functional lineage-specific prophages is observable in both animal and plant pathogens, as typified in studies using E. coli and R. solanacearum (68–70). Interestingly, it has been reported that two HAIs in R. solanacearum, namely, Rasop1 and Rasop2, showed no similarity to other phages previously sequenced in the species (35). In order to assess gene expression patterns of Rasop1 and Rasop2 during infection, we gleaned data from an existing whole-transcriptome experiment (71). Upregulation of 58.6% and 62.5% of Rasop1 (27 out of 46) and Rasop2 (20 out of 32) genes was detected upon infection of tomatoes (Table S5) (71). The evidence presented underscores lineage-specific insertion conservation as an important strategy for some organisms, probably providing a competitive advantage against close species not bearing the same traits. Therefore, PcbPr1 conservation and upregulation during infection might have implications in Pcb1692 success when facing in planta competition against other Pectobacterium spp. As one of the few relatively large lineage-specific genomic regions compared to those of other Pectobacterium spp., we speculate that PcbPr1 could be a key asset for the reported high level of virulence observed in P. carotovorum subsp. brasiliense.

Characterization of dedicated polysaccharide-biosynthetic clusters.The composition of serotype-specific polysaccharide gene clusters gives rise to biochemical diversity in bacterial cell surfaces (72). The fundamental role of cell surface in virulence has been well established in both animal and plant pathogens (32, 73–75). In SRE, however, aside from the establishment of an E. coli rffG homolog in P. atrosepticum as a probable player in enterobacterial common antigen biosynthesis (76), overall knowledge on this subject remains limited. In this context, our data set enabled detection of a 25-kb region in Pcb1692 functionally associated with polysaccharide biosynthesis, exhibiting significant coactivation of 17 out of 24 genes within the first 24 hpi in potato tubers (Fig. 4A). A homologous region in D. dadantii 3937 exhibits a different pattern while infecting A. thaliana leaves by regulating only 3/24 genes positively and 2/24 negatively (Table S1.3). Through orthology-based annotation of Pcb1692 genes compared with those of model organisms (see Materials and Methods) and subsequent integration with the STRING correlational database (77), most of the gene products were annotated into the LPS biosynthetic pathway (Fig. 4B). Out of the remaining eight unannotated entries, seven were successfully characterized by conserved domain inspection through HMM profiles (78, 79), revealing well-known domains in polysaccharide production, such as an integral membrane polysaccharide-specific transporter (PST; Pfam entry PF14667), an acetyltransferase (AT; Pfam entry PF13302), a cupin-like protein (Pfam entry PF05523), and four glycosyltransferases (GT; PFAM clans CL0110 and CL0111) (Fig. 4A and Table S6.1). The remaining unannotated gene (PCBA_RS09180) is related to gfcA-yjbE, described in E. coli, and will be discussed in detail in the next section. As a general rule, all eight of these entries unsurprisingly display high sequence variation, mostly comprising weakly conserved blocks among SRE (Fig. 4A). These lineage-specific genes presumably comprise a serotype-specific block in SRE polysaccharide biosynthesis.

• Open in new tab
• Download powerpoint

FIG 4

In planta transcriptional profile, genomic conservation, and annotation of dedicated polysaccharide clusters in Pcb1692. (A) Heatmap of GFC operon as described in the legend to Fig. 1. Asterisks next to genes indicate those for which no orthology relationship could be assessed; as such, they are not represented in the network. Next to each gene name, horizontal bars (percentage-stacked) display their respective conservation compared with SRE genomes of Dickeya (39 strains) and Pectobacterium (60 strains) genera. These bars include overall frequencies of BLASTP (sequence similarity) and MCScanX (syntenic conservation) searches of Pcb1692 gene products against 99 SRE complete data sets. The bars are divided in segments: syntenic with Dickeya spp. (blue), BLASTP positive against Dickeya spp. (red), syntenic with Pectobacterium spp. (gray), and BLASTP positive against Pectobacterium spp. (yellow). (B) Correlational network denoting the annotated gene association between successfully predicted Pcb1692 orthologs in E. coli and P. aeruginosa according to the STRING database (77) as well as showing respective GO term annotations. Thickness of network links are proportional to the combined scores obtained from the STRING algorithm. Blue lines represent new coregulation evidence found in our transcriptome data set between genes for which other categories of associations are known. Those entries not annotated in the three major GO terms shown in the network are highlighted in gray. (C) Transcriptional variation of waa (top) and wec (bottom) gene clusters represented as described for panel A.

By integrating our annotations with STRING networks, we next observed new coexpression (coordinately upregulated in our data set) associations between Pcb1692 orthologs during the course of infection (Fig. 4B and Table S6.2). The coexpression linking several genes encoding GTs and modification enzymes (e.g., rhlC, rfbA, glf, and wbpE) in our data set were undetected in STRING, indicating an uncommon genomic architecture in Pcb1692, detailed below (Fig. 4B). Some other associations are surprisingly novel, such as those among gfcBCD and wza-wzc genes, all reportedly participants of capsule biosynthesis (80, 81) (Fig. 4B). In this regard, functional studies have elucidated the role of Wza-Wzc as a transmembrane complex which spans the entire periplasm and which is required for capsule assembly (82). As for gfcABCD, although the exact functions remain unknown, these genes are regarded as the typifying theme of group 4 capsule production clusters (32, 81, 83). Therefore, the presence of wza-wzb-wzc and gfc orthologs indicates that Pcb1692 possesses, and significantly activates during infection, the genetic apparatus to produce GFC. The contrast between expression patterns found in experiments involving D. dadantii and Pcb1692 suggests that the transcriptional activation of the region is triggered by specific cues related to the host.

It is noteworthy that a recently proposed model for programmatic large-scale identification of GFC operons based on E. coli did not include internal serotype-specific genes (83). Indeed, a typical GFC operon in E. coli carries the gfc genes immediately upstream of wza-wzb-wzc paralogs (i.e., gfcE, etp, and etk) with no adjacent serotype-specific genes. However, our report supports the presence of lineage-specific genes flanked by wza-wzb-wzc and gfcBCD in Pcb1692, for which transcriptional upregulation is required upon infection, potentially representing a serotype-specific block (Fig. 4A). Moreover, the presence of poorly conserved gene organizations upstream of gfc genes seems to be a general feature in SRE genomes, as depicted by extensive gene neighborhood analysis (Fig. S3 and Table S6.3). These blocks in SRE may vary from 13 to 27 genes harboring diverse family compositions consistently flanked by wza-wzb-wzc upstream and gfcBCD downstream. Although blocks ranging from 16 to 20 serotype-specific genes seem to be preferred across SRE genomes, there is no direct association between any particular species and specific block sizes (Fig. S3 and Table S6.4). This unusual organization of the GFC lineage-specific region in SRE compared to the one in E. coli might explain why some coexpression associations (e.g., gnd and rhlC) in our data set were not present in the STRING database as mentioned above (Fig. 4B). The functional implications for SRE, and especially Pcb1692 harboring this lineage-specific gene array flanked by wza-wzb-wzc and gfcABCD, will be discussed in detail in another section.

Importantly, a recent study predicted that 40% of the bacterial lineages analyzed conserve more than one capsule system (83). This observation prompted us to systematically survey the Pcb1692 genome for other possible capsule region occurrences. Since LPS, EPS, and capsules share many gene families in their functional units, which greatly hinders fully automatic predictions, the regions in the Pcb1692 genome were manually inspected following an initial automatic screening (see Materials and Methods). The analysis allowed additional detection of two conspicuous polysaccharide clusters in Pcb1692, although neither of these regions conserves characteristic domains that typify capsule-related regions (Table S6.1). Thus, these findings indicate that GFC is the only known capsule group produced in Pcb1692. Nonetheless, these two additional regions carry eight and seven waa and wec orthologs, respectively, which are commonly associated with LPS biosynthesis and could be important assets during infection; therefore, they are worth investigating (84, 85). By analyzing their transcriptional patterns, an increased demand for GFC and wec gene transcription within the first 24 hpi was observed. In contrast, in the waa region, only 3 out of 16 genes were upregulated (Fig. 4C). Further, significant downregulation of two genes (waaF and hldD) and three borderline predictions of downregulation (waaG, waaQ, and waaC) (Table S6.1) indicated a general low demand for transcription of the waa region during infection. Curiously, between 24 and 72 hpi a slightly negative transcriptional modulation occurs in most of the genes in the three regions; however, it does not significantly impact the overall trend for upregulated genes (Fig. 4C). Unlike GFC or waa gene clusters, the wec region is remarkably conserved in nearly all SRE, encompassing 6 (out of 11) genes displaying infection-induced upregulation in Pcb1692 (Fig. 4C). These observations elucidate the unequal transcriptional demand for three distinct polysaccharide clusters during infection of potato by Pcb1692. While wec and especially gfc gene clusters seem to be consistently recruited at the transcriptional level, the transcription profile of the waa cluster is mostly flat, suggesting different functional demands during infection for these regions.

Analysis of gfcA-related sequences and genomic contexts.The inability to assess either orthology or domain conservation for PCBA_RS09180, in addition to the current lack of information regarding the function of gfcA, prompted an in-depth comparative investigation into this locus. Based on results presented above, 88% of the SRE genomes analyzed carry orphan (i.e., not clustered by OrthoMCL) genes upstream of gfcBCD (Table S6.3). In the remaining 22% of strains, we detected gene products from small clusters populated with sequences from three (OG_5444; Dickeya solani, D. dadantii, and D. chrysanthemi), two (OG_10199; D. solani and D. chrysanthemi), and one (OG_7594; D. dianthicola) species (Table S6.3). This clustering pattern in groups populated by a small number of species implies lineage-specific expansions. Furthermore, overall only 10% of SRE strains carry the YjbE domain (Pfam entry PF11106) originally described in E. coli. These include the two P. betavasculorum strains analyzed and six out of seven P. atrosepticum strains in the analysis. Hence, the presence of some PCBA_RS09180-related sequences in small clusters, represented by no more than three species (e.g., OG_5444 and OG_7594) due to the high level of sequence variation, constitutes a typical serotype-specific pattern.

In order to test these results more broadly, we next expanded the gene neighborhood analysis beyond the SRE group. Sequence-based searches were performed in order to retrieve a set of PCBA_RS09180-related proteins in bacteria. Unsurprisingly, BLASTP (86) search using PCBA_RS09180 sequence against NCBI (nonredundant protein database) was ineffective. On a preliminary level, this is consistent with the observations found in SRE that showed high sequence variation in PCBA_RS09180-related sequences. We then utilized a more sensitive approach (87) to obtain distantly related PCBA_RS09180 proteins. Within 50 positive matches, 27 were supported by publicly available genome-wide data and were suitable for gene neighborhood screening (Table S6.5). From these, three encode YjbE-containing products (Pfam entry PF11106), corroborating the relationship between PCBA_RS09180 and yjbE and/or gfcA. On the other hand, the remaining 24 gene products display no detectable domain (similar to PCBA_RS09180). The results revealed that even in evolutionarily distant organisms, 91% (22 out of 24) of the genes distantly related to that encoding PCBA_RS09180 are confined upstream of gfc and/or yjb homologous operons (Table S6.6). Thus, although this locus has been termed yjbE and/or gfcA and the respective domain described in E. coli is designated YjbE, this locus is in fact highly variable and the YjbE domain is weakly represented. An additional confirmation was obtained by consulting a large repository for domain architecture annotation (88). The presence of the YjbE domain in bacteria is at least 20-fold less frequent than that of YjbF, YjbH, and Caps_synth_GfcC (Table S6.7). Together, these results uncover the highly variable nature of the gfcA locus in a broad range of prokaryotes, which could be a serotype-specific player in the capsule biosynthesis machinery.

The gfcA-encoded products from 100 SRE strains were then aligned to assess possible conserved residues in their sequences. The results revealed a conspicuously conserved N-terminal segment harboring an MKKTLxxLxxxxAxxxxxxA motif (Fig. S4). Further, we conducted predictions for the secondary structure of PCBA_RS09180 sequence compared to those of GfcA and YjbE using two different methods (89, 90). Together, the analyses cohesively show two possible transmembrane sections in both PCBA_RS09180 and GfcA with apparent cytoplasmic orientation of the N-terminal region. Conversely, secondary-structure predictions for YjbE were inconclusive. Also, the prediction of an N-terminal signal peptide in PCBA_RS09180 and GfcA also strengthens the possibility of shared functionality. In addition, these findings indicate a closer functional relationship between PCBA_RS09180 and GfcA than PCBA_RS09180 and YjbE. Strikingly, further extensive signal peptide predictions made by two different methods in 100 gfcA sequences from SRE also returned 96% overlapping positive results in N-terminal regions. This result matches those for the conserved region found in the previous alignment, reinforcing their functional conservation (Fig. S4). Moreover, secondary-structure predictions in the same set of sequences supported that 80% of gfcA orthologs can conserve between one to three transmembrane regions, which is in accordance with the results from GfcA. Together, these results suggest that PCBA_ RS09180, and presumably SRE orthologs, could be directed to the inner membrane (similar to GfcA) to function as membrane proteins. This corroborates the current hypothesis of GfcABCD performing an auxiliary role in polysaccharide translocation across the membranes (91).

Analysis of GFC elements and wza-wzc system organization in Pcb1692.The currently accepted model for polysaccharide biosynthesis via the Wza-Wzc transmembrane conduit involves association of several functional modules, such as (i) a series of peripheral GTs and sugar modification enzymes, (ii) a Wzx-like (PST family) transmembrane flipping protein, (iii) an initiation sugar transferase, and (iv) a Wzy-like polymerase (32). The Wza-Wzc system is often associated with two highly similar systems that produce group 1 and 4 capsules (32). Curiously, the gene arrangement in SRE’s GFC region resembles the one from group 1 capsule described in other species carrying a serotype-specific block downstream of wza-wzb-wzc (Fig. S3) (92). Nonetheless, the genome-wide absence of the group 1 capsule-characteristic wzi gene suggests it could not be produced in the vast majority of SRE species, except for P. betavasculorum and P. carotovorum subsp. actinidiae (data not shown). In this context, the results for Pcb1692 were examined in light of comparative analyses with SRE and then superimposed onto the canonical Wza-Wzc model in order to propose an organization for the GFC machinery in P. carotovorum subsp. brasiliense. A total of eight GTs or modification enzymes, predicted to be soluble or weakly attached to the inner membrane, located in this region exhibit cohesive upregulation within the first 24 hpi in Pcb1692 (Fig. 5). These enzymes exhibit zero to two transmembrane (TM) sections, which is in accordance with the Wza-Wzc functional model (Fig. 5). The conserved PST within the GFC region presents 12 TM helices, which is described as the consensus number of TM sections for Wzx-like flippases (93). The functional domain detected in this sequence is also classified in the same Pfam clan (CL0222), which encompasses the domain found in the canonical Wzx (Pfam entry PF01943). Further, in Pcb1692, the Glycos_transf_4 (Pfam entry PF00953) domain found in the WbpL ortholog (PCBA_RS09080) is also conserved in WbaP/WecA sequences, which are the reported initiation transferase linked to Wzc in other bacteria (Table S6.1) (94, 95). This WbpL ortholog also features 11 TM helices solidly predicted, which matches the exact number found in WecA from E. coli (Fig. 5). As for the Wzy-like activity, however, the direct association to this Pcb1692 region is unclear. The WzyE activity in E. coli is described as being similar to that of members of the GT-B subfamily (96, 97). However, GT-B is the only unrepresented GT subfamily within the described GFC region in Pcb1692. Interestingly, a GT-C member found within the Pcb1692 GFC region conserves 10 TM helices, matching the number found in Pcb1692's direct ortholog of WzyE (PCBA_RS14975). The presence of periplasmic loops could also be predicted in this GT-C member, consistent with previously described Wzy proteins (Fig. 5) (98, 99).

• Open in new tab
• Download powerpoint

FIG 5

Schematic organization of GFC region organization in Pcb1692. At the top left, the color assignments of known or unknown protein functionalities are shown using the following abbreviations: MPA, membrane-periplasmic auxiliary protein; PTP, protein tyrosine phosphatase; OMA, outer membrane auxiliary protein; GT, glycosyltransferase. A second box, framed by continuous lines, shows a binary table with graphical representations for the presence (T, true) or absence (F, false) of five parameters analyzed in the protein sequences. Abbreviations: TM, transmembrane segments; OMP, outer membrane protein; 3DS, available three-dimensional structure; UR, significant infection-induced upregulation in Pcb1692; LC, low conservation among SRE (detailed below), as assessed by synteny and sequence similarity. The number of detected TMs is within the pink shape. Numbers with an asterisk indicate the prediction methods diverge, whereas the lack of an asterisk indicates these methods converge in the number of TMs predicted. LC is attributed to sequences exhibiting fewer than 40 and 70 positive hits, respectively, in protein similarity (BLASTP) and synteny-based (MCScanX) pairwise comparisons using 100 SRE strains (see Table S6.1 for details). The arrows indicate the transcriptional orientation (sense or antisense) of the genes within the region, and the curly braces indicate ranges of genes transcribed in a given orientation for which gene names are displayed (top to bottom) according to the genomic order. A dashed line above the region representation highlights the serotype-specific block found between wza-wzb-wzc and gfcABCD in Pcb1692.

By inspecting outer membrane (OM) propensity (see Materials and Methods) in sequences from the GFC region, only GfcD and Wza returned positive results. This opens a question on the localization of GfcB, which is generally regarded as an OM protein but surprisingly could not be predicted as one (Fig. 5) (PDB entry 2IN5). The only structurally resolved encoded product from the gfc operon currently published is GfcC, which is hypothesized to be a periplasmic protein (91). Together, the results presented here show that Pcb1692 conserves canonical elements of the Wza-Wzc system (i.e., wzx-wzy and wecA) separated from the GFC region (Fig. 4), which could support their encoded products as putative participants in the machinery. Additionally, evidence also indicates that PST, GT-4, and WbpL proteins encoded by genes within the serotype-specific block of the GFC region remarkably mimics structural features found in Wzx, Wzy, and WecA, respectively, suggesting they may undertake similar roles in the wza-wzc system (Fig. 5). Combined, these analyses expand the current view of GFC region gene organization by applying extensive contextual comparisons among SRE genomes, domain analysis, and membrane topology predictions. The presence of serotype-specific blocks flanked by wza-wzb-wzc and gfc genes in SRE sheds light on a novel genomic architecture that combines the known pattern from group 1 capsule regions with typical genes involved in group 4 capsule biosynthesis.