ABSTRACT
Soft-rot Enterobacteriaceae (SRE), typified by Pectobacterium and Dickeya genera, are phytopathogenic bacteria inflicting soft-rot disease in crops worldwide. By combining genomic information from 100 SRE with whole-transcriptome data sets, we identified novel genomic and transcriptional associations among key pathogenicity themes in this group. Comparative genomics revealed solid linkage between the type I secretion system (T1SS) and the carotovoricin bacteriophage (Ctv) conserved in 96.7% of Pectobacterium genomes. Moreover, their coactivation during infection indicates a novel functional association involving T1SS and Ctv. Another bacteriophage-borne genomic region, mostly confined to less than 10% of Pectobacterium strains, was found, presumably comprising a novel lineage-specific prophage in the genus. We also detected the transcriptional coregulation of a previously predicted toxin/immunity pair (WHH and SMI1_KNR4 families), along with the type VI secretion system (T6SS), which includes hcp and/or vgrG genes, suggesting a role in disease development as T6SS-dependent effectors. Further, we showed that another predicted T6SS-dependent endonuclease (AHH family) exhibited toxicity in ectopic expression assays, indicating antibacterial activity. Additionally, we report the striking conservation of the group 4 capsule (GFC) cluster in 100 SRE strains which consistently features adjacently conserved serotype-specific gene arrays comprising a previously unknown organization in GFC clusters. Also, extensive sequence variations found in gfcA orthologs suggest a serotype-specific role in the GfcABCD machinery.
IMPORTANCE Despite the considerable loss inflicted on important crops yearly by Pectobacterium and Dickeya diseases, investigations on key virulence and interbacterial competition assets relying on extensive comparative genomics are still surprisingly lacking for these genera. Such approaches become more powerful over time, underpinned by the growing amount of genomic information in public databases. In particular, our findings point to new functional associations among well-known genomic themes enabling alternative means of neutralizing SRE diseases through disruption of pivotal virulence programs. By elucidating novel transcriptional and genomic associations, this study adds valuable information on virulence candidates that could be decisive in molecular applications in the near future. The utilization of 100 genomes of Pectobacterium and Dickeya strains in this study is unprecedented for comparative analyses in these taxa, and it provides novel insights on the biology of economically important plant pathogens.
INTRODUCTION
Conflicts among microbes and against their hosts unfold under intense selective pressure, resulting in an abundant inventory of colonization traits in unicellular organisms from all domains of life (1–6). In prokaryotic genomes specifically, the allocation of gene inventories is recognized by frequent conservation within functionally clustered units. Such a pattern is determined by the strong selection imposed on the organization of prokaryotic genomes (7). The typical simultaneous occurrence of essential processes, such as transcription, translation, and protein localization in the prokaryotic cytoplasm, drives selective pressure on genome organization (8). In this scenario, the occurrence of clustered units constitutes a cohesive solution to facilitate cotranscription of functionally associated genes (9–12). In fact, this feature has been broadly exploited through guilt-by-association methods in several studies, enabling the elucidation of novel gene functions (9, 13, 14).
Bacteria recently reclassified in the Pectobacteriaceae family (15), responsible for causing soft rot-blackleg/wilt disease, are commonly designated, based on the previous taxonomic classification, as soft-rot Enterobacteriaceae (SRE). This group is mainly represented by the Pectobacterium and Dickeya genera (16). These microbes have attracted attention for having great impact on global food security, as they infect a considerable range of plant hosts (16, 17). Hence, roughly 100 SRE genome sequencing projects are publicly available, either completed or in progress (18–22). Underpinned by such a wealth of public data, recent transcriptome-based reports have depicted important facets of SRE biology. For example, the role of small RNAs (sRNAs) in the adaptive response of Pectobacterium atrosepticum exposed to nutrient-deficient environments has been investigated, revealing 68 regulated sRNAs in these conditions and leading to the discovery of nine novel sRNAs (23). In P. atrosepticum, the impact on the transcription of 26% of the genome upon deletion of expI revealed the critical role of quorum-sensing regulation for disease development (24). Similarly, transcriptome analyses demonstrated the impact of 32 isolated stress conditions, mimicking those found during infection, on Dickeya dadantii strain 3937 regulatory patterns, providing a detailed landscape on environmental triggers for gene expression (25). Furthermore, an investigation of the role played by the D. dadantii strain 3937 PecS global regulator during early colonization of leaf tissues uncovered more than 600 genes in its regulon (26). Thus, as a general rule, transcriptome-based approaches are powerful means to explore key strategies utilized by plant pathogens and plants’ critical physiological programs, which can optimize plant genetic enhancement, resulting in positive long-term impacts on food security (27–29).
The deployment of several exoenzyme families specializing in breaking plant cell wall is one of the most conspicuous strategies presented by SRE; hence, it is one of the most exhaustively investigated (30). The activity of plant cell wall-degrading enzymes (PCWDE) causes the release of by-products which can be taken up as nutrients by the bacterial cell (16). Some of these PCWDE include cellulases, pectate lyases (PLs), and pectin lyases (PNLs), among others (16, 31). Another important asset for disease development is the ability to biosynthesize high-molecular-weight polysaccharides, which may either be secreted extracellularly (EPS, or exopolysaccharide) or remain attached to the bacterial cell surface (e.g., lipopolysaccharide [LPS] and capsular polysaccharides [CPS]) (32).
Horizontal gene transfer is a powerful driver in the acquisition of novel functions by bacterial pathogens (33). In this context, genomically integrated bacteriophages (prophages) and toxin/antitoxin systems frequently exported by bacterial secretion systems (e.g., type I, III, and VI secretion systems) typify pathogenically important horizontally acquired islands (HAIs) (34, 35). The expression of HAIs is also a ubiquitous strategy to induce pathogenesis and to succeed in interbacterial competition (36). Specifically, the bacteriophage-like type VI secretion system (T6SS) has been implicated as a crucial ecological asset of several Gram-negative bacteria either as a virulence or bacterial competition agent (37).
In this article, regions of interest were gleaned by analyzing an original transcriptomic data set obtained from Pectobacterium carotovorum subsp. brasiliense strain PBR 1692 (here referred to as Pcb1692) during disease development. The aim was to survey the transcriptional activation of critical genomic regions required for virulence or interbacterial competition in Pcb1692 and assess their conservation in 100 SRE genomes. We report here the in planta coactivation of the carotovoricin homolog (Ctv) prophage in Pcb1692 along with a T1SS module immediately upstream. Comparative analysis added supporting evidence by unveiling an exquisite genomic conservation of the T1SS+Ctv block in Pectobacterium genomes. These results presumably point to a systemic association between these two themes in the Pectobacterium genus in which T1SS may export Ctv elements through the periplasm. Extensive gene neighborhood and protein domain architecture analyses, combined with large-scale sequencing, also shed light on the strong topological and transcriptional association of a previously predicted toxin/immunity pair (WHH- and SMI1_KNR4-containing gene products) with the T6SS machinery. Further evidence uncovered the infection-induced upregulation in Pcb1692 of a highly conserved group 4 capsule (GFC or G4C) biosynthesis cluster in SRE. The analyses also demonstrated that the GFC cluster may be the only capsule production region conserved in Pcb1692. Sequence analyses of gene products encoded by the gfcA locus in a large number of organisms unveiled high sequence variation, suggesting they have a role in GFC machinery as serotype-specific membrane proteins.
RESULTS AND DISCUSSION
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 (log2 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.
Mapping summary of RNA-Seq reads on Pcb1692 reference genome
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 log2 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-log2 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.
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 log2 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.
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.
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.
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).
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.
Conclusions.In this article, we presented an integrative approach that takes advantage of extensive public information to produce a reliable background for genome-wide studies in SRE organisms, including 100 strains from Pectobacterium and Dickeya genera. By combining this platform with an original transcriptome data set, we shed light on strong genomic associations among virulence and interbacterial competition determinants in the SRE group, supported by coordinated transcriptional upregulation of these elements in Pcb1692 throughout disease development in potato tubers. Collectively, these findings provide strong evidence for gene associations and/or infection-induced transcriptional demand within important pathogenicity themes, including PCWDE, T6SS, Ctv, and the capsule biosynthesis machinery, which may pave the way for biotechnological applications in the near future.
MATERIALS AND METHODS
Culture media, growth conditions, and total RNA extraction.Wild-type Pcb1692 was grown on a nutrient plate at 37°C for 16 h. Overnight cultures were prepared by inoculating a single colony into 10 ml Luria-Bertani (LB) broth at 37°C for 16 h with constant shaking at 200 rpm. To obtain potato inoculated with Pcb1692, healthy potato tubers (Solanum tuberosum cv. Mondial) were inoculated with the Pcb1692 (optical density at 600 nm [OD600] of 1) wild-type strain as previously described (117). The experiments were performed using three biological replicates, with three potato tubers per replicate. Macerated potato tissue was scooped out at 24 and 72 hpi and homogenized in double-distilled water. Bacterial cells were recovered by grinding the scooped, macerated potato tissues in 20 ml of double-distilled water using an autoclaved pestle and mortar. Starch material was removed by centrifuging ground tissue at 10,000 rpm for 1 min. Potato-inoculated and in vitro-cultured bacterial cells were stabilized using RNAprotect reagent (Qiagen, USA) according to the manufacturer’s instructions. Total RNA from in vitro-grown and potato-inoculated bacteria was extracted as previously described using the RNeasy minikit (Qiagen, Hilden, Germany).
Total RNA quality and cDNA library construction.The concentration and purity of each extracted total RNA sample was evaluated using spectrophotometric analysis (NanoDrop ND-1000; NanoDrop Technologies, Wilmington, DE) at a ratio of 230 to 260 nm. Using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.), total RNA sample concentration, RNA integrity number (RIN), and 28S/18S ratio were determined. Aliquots (200 ng) of total RNA from Pcb1692 extracted from in vitro-grown cells (16 h) and infected potato tubers at 24 and 72 hpi was used to prepare cDNA libraries. The Illumina sequencing service was provided by BGI Co., Ltd. (China). TruSeq RNA sample prep kit v2 (Illumina, USA) was used to construct the cDNA libraries by following the manufacturer’s protocol. In summary, the mRNA was cleaved into small fragments, followed by the synthesis of first-strand cDNA with random hexamer-primed reverse transcription. RNase H and DNA polymerase I were used to synthesize second-strand cDNA. The double-stranded cDNA was subjected to terminal modification by addition of adenosine and ligated with adapters. Adaptor-ligated fragments with suitable sizes were selected and enriched by PCR using the PureLink PCR purification kit (Invitrogen, USA) to create the cDNA libraries for sequencing. Paired-end sequencing (PE91) was performed on an Illumina HiSeq 2000 sequencing platform.
Read mapping, gene expression, and statistical analysis.Quality assessment of the raw data were performed by fastqc software (https://www.bioinformatics.babraham.ac.uk/projects/fastqc), and then low-quality regions were trimmed by Trimmomatic v 0.36 (100). The trimmed RNA-Seq reads were mapped to the Pcb1692 reference genome (GCF_000173135.1) using hisat2 v 2.1.0 (101). Raw read counts were performed in the R environment (https://www.r-project.org/) by the featureCounts package (102), and subsequent statistical analysis of differential expression was performed by the edgeR package (103). The threshold parameters used to assign differential expression were an FDR of <0.01 (104), absolute log2 fold change of >1 (upregulated), or log2 fold change of <−1 (downregulated). Gene transcriptional profiles were graphically rendered by Gitools (105).
Orthology analysis, domain architecture detection, and gene neighborhood screenings.Orthology relationships between protein sequences were assessed using the OrthoMCL (inflation, 1.5) pipeline (106), which takes tabular results from a BLASTP search under specific parameters (options: -seg yes -outfmt 6 -num_threads 3 -num_alignments 100000 -evalue 1e-05) as input. The sequences are clustered in orthologous groups labeled with OG numeric tags (e.g., OG_1 and OG_2). In parallel, all sequences were characterized using HMMER3 (107) supported by the Pfam database in order to generate in-house predictions of conserved domain architectures through hidden Markov model (HMM) profiles (108). Orthology and domain architecture information are then combined with genomic coordinates from each genome by custom Perl scripts in order to generate annotated gene neighborhood screenings.
Ectopic expression of predicted T6SS-dependent toxins.Removal of the rat FABP1 gene from pTrc99A (AddGene) removed some cut sites from the multiple-cloning site (MCS). KnpI and BamHI cut sites were reintroduced into the plasmid by incorporation of the cut sites into the 5′ region of the PCR primers (H5907_F/R) (Table 2). PCR amplification of PCBA_RS05790 (encoding a hypothetical protein upstream of WHH-containing nuclease) incorporated XbaI, KpnI, and BamHI upstream of the open reading frame (ORF) and BamHI and SalI downstream of the ORF. The PCR product was subcloned into pJET1.2/blunt and excised with XbaI/SalI restriction digestion. pTrc99A (without rat FABP1) was digested with XbaI/SalI and ligated with the excised PCR product. The gene was removed using BamHI digest and the linear plasmid religated. In this manner, BamHI and KnpI were reintroduced into the MCS to facilitate cloning of effector genes. The resulting plasmid is referred to as pTrc100, as its MCS differs slightly from that of pTrc99A, but the rest of the plasmid remains unchanged.
Primers used in this study
Effector genes (PCBA_RS05785, WHH-containing nuclease, SacI/PstI; PCBA_RS05775, D123-containing protein, KpnI/PstI digest; PCBA_RS18045, phospholipase, SacI/PstI digest; and PCBA_RS22965, AHH-containing nuclease, SacI/PstI digest) were ligated into arabinose-inducible expression plasmid pCH450. Immunity genes (PCBA_RS08690, AHHi, KpnI/PstI digest) were ligated into isopropyl-β-d-thiogalactopyranoside (IPTG)-inducible pTrc100 and transformed into E. coli DH5α. Overnight cultures were grown in 0.1% glucose to repress expression from the arabinose-inducible promoter. Cultures were washed in 10 mM MgSO4 adjusted to an OD600 of 0.05 in fresh LB broth supplemented with 5 μg/ml tetracycline, 50 μg/ml ampicillin, and 1 mM IPTG as required. After 30 min, expression from pCH450 was induced with 0.2% l-arabinose. Where necessary, expression from pTrc100 with IPTG was induced immediately. Optical density (at 600 nm) was measured hourly.
Collinearity analysis, prophage origin predictions, and capsule region inspection.In order to predict colinear conservation in SRE genomes, we first obtained genomic annotation and protein sequence data from all Pectobacterium and Dickeya strains available in the RefSeq database (https://www.ncbi.nlm.nih.gov/refseq) (see Table S7 in the supplemental material) until December 2017, totaling 100 strains, including Pcb1692. We then performed the synteny analysis using MCScanX (109), allowing a minimum match size of two genes to call syntenic blocks (option -s 2). MCScanX takes the tabular results of a nonstringent BLASTP search (86) as the input (option -evalue 1). The total number of positive matches of gene products (from one organism) against all other SRE protein data sets resulting from both BLASTP and MCScanX were parsed by in-house Perl scripts (https://www.perl.org/). Graphical ideograms were scripted in Circos (110), compiling syntenic regions in the genomes as ribbon links between chromosomes, as well as BLASTP and MCScanX hit counts from pairwise comparisons previously described. Prophage origin predictions were conducted by using PhiSpy software (111) in parallel with BLASTP search support (options -qcov_hsp_perc 40 -evalue 1e-05) against Phast (112) and NCBI bacteriophage sequence databases. In parallel, we performed genome-wide synteny predictions on Pcb1692 and P. carotovorum subsp. carotovorum strain BCS2 using recently predicted prophage regions published previously (35). The capsule regions were initially surveyed using HMM-profile search (113) based on the curated database of capsule-related domains made public by Rendueles et al. (83). The Pcb1692 genome next was programmatically scanned for segments harboring at least three consecutive matches carrying the above-mentioned capsule-related domains with two nonmatches gap, allowing for greater sensitivity. These segments were then manually inspected to identify detectable functionality in the adjacent genes that could support the existence of a capsule region.
Sequence search, alignment, and membrane topology predictions of GfcA-related entries.PCBA_RS09180-related sequences were obtained using the HHblits online tool with default parameters (87). HHblits positive hits were filtered by at least 40 matched amino acid columns in the HMM-HMM alignment, representing ∼40% of the query (PCBA_RS09180) sequence (Table S6.5). All above-threshold entries were retrieved. Out of 50 above-threshold entries, 27 were supported by publicly available genome-wide data, making them suitable for gene neighborhood screening. Genomic data for these 27 structures (complete genome/scaffold/contig) were then obtained from the RefSeq database for gene neighborhood screening. The CDART database was inspected to obtain extensive representativeness of YjbE, YjbF, Caps_synth_GfcC, and YjbH (Pfam entries PF11106, PF11102, PF06251, and PF06082). (88). Sequence alignments were performed using Clustal Omega (114). Prediction of transmembrane segments was conducted by using TMpred (89) and TMHMM (90) online servers in parallel. Signal peptides were predicted by SignalP 4.1 server (66) and Signal-BLAST (67). Detection of putative outer membrane proteins was conducted by using HHomp (115).
STRING network integration.In order to integrate Pcb1692 with the STRING database (77), orthology annotations with Escherichia coli strain K-12 and Pseudomonas aeruginosa strain PAO1 were used. The Pcb1692 entries conserving orthology with the model organisms were supplied to the STRING correlational database. This provided association between annotated orthologs with the respective Gene Ontology (GO) terms (116) for each sequence.
Accession number(s).Data determined in the course of this work have been deposited in NCBI’s Gene Expression Omnibus (GEO) and are accessible through GEO accession number GSE102557.
ACKNOWLEDGMENTS
This research study was funded by the National Research Foundation (NRF), South Africa, through Competitive Funding for Rated Researchers (CFRR 98993), NRF Bioinformatics and Functional Genomics (BFG 93685), and the NRF Research Technology and Transfer Fund (RTF 98654). D.Y.S. received an NRF BFG postdoctoral fellowship. D.B.-R. received a University of Pretoria postdoctoral fellowship. N.M. was funded by NRF Freestanding, PSA, and University of Pretoria Postgraduate bursaries. C.K.T. was funded by a University of Pretoria bursary. Any opinion, findings, conclusions, or recommendations expressed in this material are those of the author(s), and the NRF does not accept any liability in this regard.
Author responsibilities were as follows: conception or design of the study, D.B.-R., C.K.T., and L.N.M.; acquisition, analysis, or interpretation of the data, D.B.-R., C.K.T., N.M., D.Y.S., S.K., and L.N.M.; writing, D.B.-R., N.M., and L.N.M.
FOOTNOTES
- Received 24 August 2018.
- Accepted 29 October 2018.
- Accepted manuscript posted online 9 November 2018.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.02050-18.
- Copyright © 2019 American Society for Microbiology.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.↵
- 112.↵
- 113.↵
- 114.↵
- 115.↵
- 116.↵
- 117.↵