Genome-Wide Analyses Revealed Remarkable Heterogeneity in Pathogenicity Determinants, Antimicrobial Compounds, and CRISPR-Cas Systems of Complex Phytopathogenic Genus Pectobacterium

Abstract

The Pectobacterium genus comprises pectolytic enterobacteria defined as the causal agents of soft rot, blackleg, and aerial stem rot diseases of potato and economically important crops. In this study, we undertook extensive genome-wide comparative analyses of twelve species that conform the Pectobacterium genus. Bioinformatics approaches outlined a low nucleotide identity of P. parmentieri and P. wasabiae with other species, while P. carotovorum subsp. odoriferum was shown to harbor numerous pseudogenes, which suggests low coding capacity and genomic degradation. The genome atlases allowed for distinguishing distinct DNA structures and highlighted suspicious high transcription zones. The analyses unveiled a noteworthy heterogeneity in the pathogenicity determinants. Specifically, phytotoxins, polysaccharides, iron uptake systems, and the type secretion systems III-V were observed in just some species. Likewise, a comparison of gene clusters encoding antimicrobial compounds put in evidence for high conservation of carotovoricin, whereas a few species possessed the phenazine, carbapenem, and carocins. Moreover, three clustered regularly interspaced short palindromic repeats-Cas (CRISPR-Cas) systems: I-E, I-F, and III-A were identified. Surrounding some CRISPR-Cas regions, different toxin and antitoxin systems were found, which suggests bacterial suicide in the case of an immune system failure. Multiple whole-genome alignments shed light on to the presence of a novel cellobiose phosphotransferase system (PTS) exclusive to P. parmenteri, and an unreported T5SS conserved in almost all species. Several regions that were associated with virulence, microbe antagonism, and adaptive immune systems were predicted within genomic islands, which underscored the essential role that horizontal gene transfer has imparted in the dynamic evolution and speciation of Pectobacterium species. Overall, the results decipher the different strategies that each species has developed to infect their hosts, outcompete for food resources, and defend against bacteriophages. Our investigation provides novel genetic insights that will assist in understanding the pathogenic lifestyle of Pectobacterium, a genus that jeopardizes the agriculture sustainability of important crops worldwide.

Keywords: CRISPR-Cas; Pectobacterium; antimicrobial compounds; comparative genomics; dynamic evolution; horizontal gene transfer; pathogenicity determinants.

Figure 1

Pairwise heatmap of Pectobacterium species based on average nucleotide identity (ANI) and digital DNA-DNA hybridization (dDDH). Both ANI and dDDH are represented as percentage values among the twelve Pectobacterium species and two Dickeya species (D. zeae EC1 and D. solani IPO 2222). The upper diagonal displays ANI data, whereas the lower diagonal depicts the in silico dDDH data. The cut-off values for species delineation are 95 and 70% for ANI and dDDH, respectively.

Figure 2

Concatenated phylogenetic tree illustrating the relationships among all Pectobacterium species. A total of 59 strains from all known Pectobacterium species were selected to generate Neighbor-Joining tree based on 12 housekeeping genes: dnaA, dnaN, dnaX, fusA, gapA, gyrA, gyrB, recA, recN, rpoA, rpoD, and rpoS. The alignment of individual genes was performed using MUSCLE, and a concatenated tree was subsequently generated after aligning all 12 genes (total length ~ 18,269 bp). The Maximum Composite Likelihood method was applied to determine the evolutionary distances, with each node being supported by a bootstrap of 1,000 replicates to assess reliability. Bootstrap values are displayed next to the branches. The branches were colored to highlight the distinct strains clustered in each species. The dendrogram was drawn to scale, and all positions containing gaps and missing data were eliminated. DNA sequences of D. zeae EC1 and D. solani IPO 2222 were included as out-groups to root the tree. The evolutionary analysis was developed in MEGA X. Colored square and diamond shapes indicate the complete and draft genomes included in the study, respectively. Two main clades were identified and marked with dark red and green rectangles.

Figure 3

BLAST matrix between and within total proteomes of Pectobacterium genus. Pairwise protein comparison was performed using BLAST and protein coding sequences were compared among the genomes. A BLAST hit was considered significant when 50% of the alignment revealed identical matches and covered at least 50% Alignment. If two protein sequences were similar based on cut-off value, they were grouped within the same protein family. The color scale saturation changes from dark green to light green indicated the degree of homology between the proteomes, whereas the color scale from dark red to light red showed the homologous hits within the proteome itself (internal paralogs at the bottom row of the matrix). Twelve genomes of different Pectobacterium species were compared in the analysis.

Figure 4

Heat maps featuring the relationships across the twelve Pectobacterium species based on (A) codon usage and (B) amino acid usage. The identified gene data was plotted and used to construct heat maps in R with heatmap.2 function. The species were reorganized according to the profiles of codon and amino acid usage. Dendrograms at the top and the left side of both heat maps were drawn to illustrate the species clustering pattern. The row and column colors next to both dendrograms highlighted the different clusters observed after analysis. Histograms indicated the scale colors based upon the percentages of codon usage and amino acid frequencies.

Figure 5

Pan and core genome analysis among the Pectobacterium species. (A) Pan-Core genome plot. The pan and core genome were calculated using BLAST with a cutoff of 50% identity and 50% coverage. The total number of genes, conformed for pan and core genomes, are represented with green and orange squares, respectively. The bars in the plot represent the number of new genes (dark grey bars) or new gene families (light grey bars) found after the addition of each new genome. Therefore, the size of the pan-genome increased with each new genome. While the number of core genome decrease with the addition of each new genome. (B) Pan-genome tree. The dendrogram illustrates the grouping among the Pectobacterium species based upon shared gene families (core genome) defined in the pan and core genome analysis. Two main clades were defined, and the main branches of each clade are highlighted in green and violet colors. Shaded color rectangles were added manually to emphasize how the species were clustered according to the core genome.

Figure 6

Circular visualization of predicted Genomic Islands (GIs) among the 12 Pectobacterium species. Starting from the top left row to the bottom right—illustrated the circular GIs plots for the species: (A) P. atrosepticum SCRI1043, (B) P. carotovorum subsp. carotovorum PCC21, (C) P. carotovorum subsp. brasiliense BC1 (D) P. carotovorum subsp. odoriferum BC S7, (E) P. carotovorum subsp. actinidiae KKH3, (F) P. aroidearum PC1, (G) P. parmentieri SCC3193, (H) P. wasabiae CFBP 3304, (I) P. betavasculorum NCPPB 2795, (J) P. polaris NIBIO1392, (K) P. peruviense IFB5232 and (L) C. Pectobacterium maceratum PB69. The interactive visualization of distinct islands across the genomes are shown with colored blocks according to the description provided in the predictor tool: IslandPick, based on genome comparison (green); IslandPath-DIMOB, based on associated GIs features such as tRNAs, transposon elements, integrases and sequence bias (blue); SIGI-HMM, based on the codon usage bias with a Hiddden Markov model criterion (orange); Islander, based on the mapping of GIs in regards with the recurrent use of tRNA and tmRNA genes (turquoise); and, the integrated results of four tools (dark red). Antimicrobial resistance genes as well as pathogen-associated genes are also displayed as pink and yellow circular glyphs, respectively. Species with incomplete genomes, including P. c. subsp. actinidiae, P. betavasculorum, P. peruviense, and C. P. maceratum, the contigs were ordered against the closest relative reference genome (light green outer circle). The references were selected based upon the ANI and dDDH analysis; consequently, P. atrosepticum was used as the reference for P. betavasculorum and P. peruviense, P. c. subsp. brasiliense as the reference for P. c. subsp. actinidiae, and P. polaris as the reference for C. P. maceratum. Gray areas represented unaligned contigs.

Figure 7

Comparison between the horizontal acquired island 2 (HAI2) of Pectobacterium atrosepticum SCRI043 with the genomic island 6 (GI-6) of the newly reclassified species P. peruviense. (A) Genomic Islands (GIs) plots illustrating the location of HAI2 in P. atrosepticum and P. peruviense. The plots displayed the predicted GIs based on five approaches (from the inner to the outermost): Islander (turquoise), IslandPick (green), SIGI-HMM (orange), IslandPath-DIMOB (blue), and the integrated methods of the four tools (red). Pink and yellow dots indicated antimicrobial resistance and pathogen-associated genes. The olive green color in each island was manually added to highlight the location of the HAI2 among the three genomes. Moreover, the callout lines in each genome point out the approximate length size of HAI2. (B) Linear genomic representation of HAI2 in the three genomes. The arrows depict the genome organization within each island. In addition, the orientation of the arrows indicated forward or reversed the position of each gene. Arrows Color highlighted different genetic functions within the island as well as gene clusters. Dark yellow shaded regions represent areas with high identity between the analyzed genomes with cut-off values of 100% and 71%.

Figure 8

Genomes comparison across the Pectobacterium genus. The BLAST ring image generated using BRIGS displayed the genome comparison of 12 genomes within genus Pectobacterium. From the innermost origin of each lane, the image features were characterized: the axis with the size of genome (kbp), the GC content (black line), the GC skew (purple color indicates G’s), and afterward each color ring depicting each Pectobacterium species genome, namely, P. atrosepticum SCRI1043 (NC_004547), P. carotovorum subsp. carotovorum PCC21 (NC_018525), P. c. subs. brasiliense BC1 (NZ_CP009769), P. c. subsp. odoriferum BC S7 (CP009678), P. c. subsp. actinidiae KKH3 (NZ_JRMH00000000), P. aroidearum PC1 (NC_012917), P. parmentieri SCC3193 (NC_017845), P. wasabiae CFBP 3304 (NZ_CP015750), P. betavasculorum NCPPB2795 (NZ_JQHM00000000), P. polaris NIBIO1392 (NZ_CP017482), P. peruviense IFB5232 (NZ_LXFV00000000), and Candidatus Pectobacterium maceratum (NZ_PDVY00000000). The outermost ring identified the different loci of most relevant pathogenicity determinants and antimicrobial biosynthetic clusters. The complete genome of P. atrospeticum served as a template for the generation of the BRIG image based upon the nucleotide BLAST analysis. Blue arrows highlighted genomic regions that exhibited divergences when compared to the reference genome. Key acronyms (red legends) used to name the clusters or single virulence genes in the outermost ring are detailed: T1SS-MRP (Type I secretion system—multi-repeat protein), T1SS-prtDEF (Type I secretion system—protein transporters PrtDEF), T2SS (Type II secretion system), T3SS (Type III secretion system), T4SS (Type IV secretion system), T5SSa (Type V secretion system-hecAB), T5SSb (Unreported type V secretion system—described in this study), LPS (lipopolysaccharide), cps (capsular polysaccharide), ECA (Enterobacteria Common Antigen), EPS (Exopolysaccharide O-antigen), Iron-hmsHFRS (hemin storage), Iron- entAFECD (enterobactin), Iron-fusBACD (ferredoxin uptake), Iron-hasRADEF (heme acquisition), Iron-cbrABCD (achromobactin), 3H2B (3-hydroxy-2-butanone pathway), Flp/Tad (Fimbrial low-molecular-weight protein/Tight adherence protein), IV pilus (Type IV pilus), pilW (Pilus assembly protein PilW), pilABC (pilus production proteins PilABC), cfa (coronofacic acid), NRPS-syr (Non-Ribosomal Peptide-Synthetase like syringomycin), ars (arsenic resistance), Ctv (carotovoricin), ehp (phenazine antibiotic biosynthesis), colicin (colicin-like bacteriocin), Svx (virulence protein homologue to Xanthomonas campestris), SaxA1 (aliphatic isothiocyanate resistance protein SaxA), RplY (ribosomal protein RplY), Cit1 (citrate transporter), and Nip (Necrosis-inducing protein), and appA (phytase).

Figure 9

Comparison of the genetic organization of type III secretion system (T3SS) among the Pectobacterium species. The arrow position represented forward/reverse gene orientation. Arrow color signified specific gene composition within the T3SS. Gene names were provided at the top and bottom of the linear graph; the locus-tag was used for hypothetical or genes with no names. A pairwise alignment between the linear sequences was rendered based upon BLAST algorithm with cut-off values from 64% to 100%. Regions with higher nucleotide identity were displayed with a shaded grey. A light-yellow blur square highlighted the variable region among the genomes that interrupted the type III secretion gene cluster. Expanded legend entry acronyms are provided: hrc (hypersensitive response and pathogenicity conserved genes), hrp (hypersensitive response and pathogenicity or hairpin proteins), VgrG (valine-glycine repeat protein G), PAAR (proline-alanine-alanine-arginine repeat protein), dspE (disease-specific effector protein E), dspF (disease-specific chaperone protein F), hecB (hemolysin activation protein HecB), and hecA (hemolysin/hemagglutinin-like protein HecA).

Figure 10

Schematic sequence organization of the type I-F, type I-E, and type III-A clustered regularly interspaced short palindromic repeats-Cas (CRISPR-Cas) across the genomes of the 12 Pectobacterium species. (A) CRISPR-Cas subtype I-F marked in orange. (B) CRISPR-Cas subtype I-E system marked in green. (C) CRISPR-Cas type III-A marked in yellow. Brown horizontal pentagons represented the CRISPR loci while the arrows depicted specific gene composition such as the cas or csm and other surrounding genes with either known or unknown function. The arrow and pentagon positions represented sequence orientation. Descriptions of gene names are provided inside the arrows or pentagons as well as at the top or bottom of each linear graph; genes without a name or hypothetical proteins labeled with locus-tag. The linear graphs were drawn manually and represent the scale of the respective gene lengths. Key acronyms used in the linear sequences: CRISPR (clustered regularly interspaced short palindromic repeats), cas, csy, and csm (genes encoding CRISPR-associated proteins), Lrp/AsnC (leucine responsive regulatory protein/asparagine synthase C transcription regulator family); YitT (membrane protein YitT), Amt (aspartate/tyrosine/aromatic aminotransferase), Pept (M48 family peptidase), Hcp (hemolysin-coregulated protein). VapC (toxin protein), VagC (antitoxin/virulence-associated protein), HigB (mRNA interferase toxin HigB), RelE/ParE (toxin-antitoxin system, mRNA interferase toxin RelE/toxin ParE), HicA (mRNA interferase toxin HicA), HicB (antitoxin HicB), pemK (mRNA interferase toxin pemK).

Figure 11

Circos plot description of main genetic features, relevant virulence, and antibiotic biosynthesis clusters among all Pectobacterium species. Seven circular layers have been rendered to show unique features based upon the comparative genomic analyses. From the outermost layer to the innermost layer: the area of each genome; sizes of each genome (axis in kb); gene function types across the species, name of each species indicated in 3rd layer—acronyms used in the plot starting with Pectobacterium atrospecticum (Patrosepticum), P. carotovorum subsp. carotovorum (Pcarotovorum), P. c. subsp. brasiliense (Pbrasiliense), P. c. subsp. odoriferum (Podoriferum), P. c. subsp. actinidiae (Pactinidiae), P. aroidearum (Paroidearum), P. parmentieri (Pparmentieri), P. wasabiae (Pwasabiae), P. betavasculorum (Pbetavasculorum), P. polaris (Ppolaris), P. peruviense (Pperuviense), and Candidatus Pectobacterium maceratum (Cmaceratum); virulence genes (pink lane) and key gene regulators (green lane), as indicated in the legend (5th layer); scattered dots of all genes, where the red and blue dots depict negative and positive gene orientation, respectively; and the inner layer showing the connection among homologues genes of all gene clusters covered in this study. The colors of the connector lanes of each cluster are labeled on the right side of the figure. Expanded legend entry acronyms are provided: T1SS (Type I secretion system), T2SS (Type II secretion system), T3SS (Type III secretion system), T4SS (Type IV secretion system), T5SS (Type V secretion system), T6SS (Type VI secretion system), Flp/Tad (Fimbrial low-molecular-weight protein/Tight adherence protein), cfa (coronofacic acid), NRPS (Non-Ribosomal Peptide-Synthetase) like syringomycin, 3H2B (3-hydroxy-2-butanone pathway), ECA (Enterobacteria Common Antigen), Cps-Poly (capsular polysaccharide), LOS/LPS (lipopolysaccharide/lipooligosaccharide), EPS—O-ag (Exopolysaccharide O-antigen), entAFECD (enterobactin), cbrABCD (achromobactin), fecIRABCD (ferric citrate uptake), hasRADEF (heme acquisition), hsmHFRS/pgaABCD (hemin storage), and fusABCD (ferredoxin uptake).