In previous study, we identified two Pectobacterium brasiliense isolates, SM and DQ, with distinct pathogenicity levels. To comprehensively characterize the phenotypic traits of these strains, their virulence was evaluated using three distinct infection modes. When the stem base is injected, SM inoculated plants exhibited rapid withering and desiccating, indicating higher virulence compared to DQ in inducing aerial stem rot (Figure 1A). Using soil drenching approach, we found that SM had higher pathogenicity compared to DQ in inducing potato blackleg (Figure 1B). Plants inoculated with SM strain showed stem base darkening and wilting of the leaves within 2 weeks. The maceration areas observed on potato slices inoculated with SM are consistently larger compared to those inoculated with DQ (Figure 1C).

To further confirm that this virulence difference is not specific to one potato cultivar., we conducted pathogenicity assays using stem injection on multiple cultivated varieties and lab-preserved foreign resources (Figures 1D,E). Regarding the main domestic varieties, QS9 (Qingshu No. 9), YS108 (Yunshu No. 108), and L7 (Longshu No. 7) are relatively less susceptible than other main varieties. The pathogenicity index of SM on QS9 only reached 50% on the thirteenth day after infection. The pathogenicity index of SM on LS6 (Lishu No. 6), DF3 (Dafeng No. 3) and SY2 (Shiyan No. 2) reached 50% on the third day of infection (Supplementary Figure S1). Among them, JZ1 (Jingzhang No. 1) has the fastest onset, and the pathogenicity index of SM is greater than 60% on the third day of infection. However, the pathogenicity of DQ is very low among all varieties, with the pathogenicity index not exceeding 15%. For foreign varieties, the A292 variety had strong disease resistance, and the pathogenicity index of SM and DQ did not exceed 20% thirteenth days post inoculation (Figure 1E). In addition, the variety RA517-69 is the most susceptible to the disease, and the pathogenicity index of SM reached 50% 1 day post inoculation (Supplementary Figure S2). To rule out that this pathogenicity difference was due to the different growth rates of the bacteria, we measured the growth curves of the two bacterial strains (Supplementary Figure S3). The results showed that the growth rates of the two bacterial strains were similar, and both entered a growth plateau after 16 h. In summary, SM and DQ exhibit significant pathogenicity variation assayed on multiple potato varieties.

To gain a deep insight into the genome features of these two strains, we sequenced the entire genome of SM and DQ (Table 1). The genomes of SM and DQ consist of the single circular chromosome that are 4,953,627 and 4,764,461 bp in size with no plasmid. The average GC content in the whole genome of SM is 52.03%, which is similar to that of DQ (52.13%), SX309 (52.18%), and 1692 (52%). A total of 4,390 and 4,331 Open Reading Frames (ORFs) were predicted in the SM and DQ genomes, respectively (Table 1). In the genome of SM, in addition to 4,285 protein-coding genes (CDSs), the chromosome also contains 127 RNA genes, including 77 tRNA genes. In the genome of DQ, 4,227 protein-coding genes (CDSs) and 123 RNA genes containing 76 tRNA genes were predicted (Supplementary Table S1). Through analysis with the Pfam and SignalP databases, we classified 3,018 (67.77%), 2,984 (67.51%), 305 (6.85%) and 282 (6,38%) genes with different Pfam domains and signal peptides, respectively (Supplementary Table S1). Among them, 4,254 (95.53%) and 4,142 (93.71%) of the predicted genes of SM and DQ were assigned to KEGG categories, respectively (Supplementary Table S1).

To determine the phylogenetic status of SM and DQ, the phylogenetic tree of Pectobacterium spp. and Dickeya spp. was constructed based on 3,860 bp concatenated sequences of gapA, gyrA, icdA, mdh, proA and rpoS using MEGA11 software with Maximum Likelihood method and Tamura-Nei model (Supplementary Figure S4). The bootstrap values of 1,000 replicates display the significance of each branch. The MLSA-based phylogenetic tree indicated that strains SM and DQ belong to Pectobacterium brasiliense and have close distance with SX309 and the model strain 1692 as the near neighbors.

To search for genetic features that might explain the virulence difference between these two strains, we performed comparative genomic analysis using the whole genome sequencing data. The genome alignment circular map was made with BLAST+ using BRIG (BLAST Ring Image Generator; Alikhan et al., 2011; Figure 2A). The map unveiled an exceptionally elevated level of homogeneity, exceeding 90% similarity between the two genomes. We calculated the ANI and DDH values of SM, DQ and 1692 (Supplementary Table S2). The findings indicate that the Average Nucleotide Identity (ANI) and DNA–DNA Hybridization (DDH) values between strains SM and DQ, when compared to each other, were approximately 97.27% and 86%, respectively. The DDH values calculated between strains SM and 1692, as well as between DQ and 1692, were found to be 79.2% and 80.5%, respectively. These values exceed the generally accepted species threshold of 70%, indicating a significant genetic relatedness between the strains. Additionally, the ANI values were 96.42% for SM and 96.09% for DQ compared to 1692, both clearly exceeding the threshold of 95% used for species demarcation. These results suggested that strains SM and DQ clustered closely and were at the same position in classification.

To assess the evolutionary distance between strains SM and DQ at the subspecies level, the whole genome sequences of SM and DQ were aligned with their closest sequenced relatives, SX309 and the type strain 1692 using Mauve (Figure 2B). The alignment revealed a high degree of similarity between the two genomes, with a pronounced collinearity. Notably, both SM and DQ genomes exhibited greater similarity to SX309 than to the type strain 1692. This observation supports the relationships elucidated in the phylogenetic analysis. In comparison to SX309 and 1692, there is no obvious gene inversion or deletion of large regions in SM and DQ genome, but a large area of local collinear block (LCB) inversion occurred. This pairwise comparison supports the analysis results that strains SM and DQ belong to the same subspecies.

In order to identify the specific genes with respect to pathogenicity in SM and DQ, we performed a pan-genomic comparative analysis of the two genomes (Figure 2C). The statistical analysis unveiled a comprehensive set of 4,357 pan-gene clusters, comprising 3,497 core orthologous gene clusters (80.26%) and 860 dispensable gene clusters (19.74%). Within the 860 dispensable gene clusters, 479 (10.99% of all pan gene groups) specific gene clusters unique to the SM genome were absent in the DQ genome. Conversely, the DQ genome harbors 381 unique gene clusters, and the functions of some of these unique genes remain currently unknown. Specific genes associated with KEGG functional categories were stated. Only SM has one specific gene in the category of cell motility, signal transduction, energy metabolism and nucleotide metabolism. The specific genes in SM and DQ strains are largely distributed in membrane transport and carbohydrate metabolism (Supplementary Table S3). In summary, the core-pan gene analysis indicated that the SM strain possesses a greater number of unique gene groups.

The production of plant cell wall-degrading enzymes and secretion of these enzymes are required for Pectobacterium brasiliense infection. To identify whether SM and DQ strains have different plant cell wall-degrading enzymes and secretion systems which are related to pathogenicity, we compared the plant cell wall degrading enzymes and secretion systems genes in SM and DQ. A total of 42 and 41 related genes encoding pectinases, cellulases, and proteases were identified in the SM and DQ genomes, respectively (Supplementary Table S4). Concisely, the SM genome encompasses 18 genes that encode pectin-degrading enzymes. Notably, these pectin degradation function genes are conserved in both the SM and DQ genomes. Likewise, 9 protease-encoding genes were identified in both the SM and DQ genomes. Furthermore, the SM genome features 15 genes associated with cellulose degradation, encompassing cellulases, endoglucanase, beta-glucosidase, and cellulose synthase operon proteins. Among these genes, a majority are conserved in both genomes. However, one gene encoding endoglucanase was found to be absent in DQ. This gene belongs to glycosyl hydrolase 12 family, and it functions in polysaccharide catabolic process which may contribute to pathogenicity. In general, there is little significant difference in the number of cell wall-degrading enzymes between the two bacteria except one endoglucanase S gene which is encoded specifically in SM.

The SM and DQ genomes harbored multiple secretion systems, including the type II secretion system (T2SS), type III secretion system (T3SS), and type VI secretion system (T6SS), which were closely related to bacterial pathogenicity (Figure 3). Both SM and DQ lacked Type IV secretion system and Type V secretion system. SM had one gene lapB specifically for Type I secretion system. The T2SS genes were conserved between SM and DQ, except that pulF was absent in DQ. Many plant pathogenic bacteria employed T3SS to inject a diverse array of toxic proteins into plant cells to facilitate their infection. An hrp/hrc gene cluster containing 12 genes was found in the SM and DQ genomes, and these genes were highly conserved in the two genomes. T6SS was widespread in many Gram-negative bacteria and deliver toxic effector proteins to neighboring bacteria or host cells. In this study, the prevalence of T6SS was observed in both the SM and DQ genomes, with the SM genome containing 31 T6SS-associated genes and the DQ genome containing 29 such genes. There were 4 and 5 vgrG genes and 13 and 10 hcp genes in the genomes of SM and DQ, respectively. These genes encoded extracellular structural components of the secretion machinery and specific effectors (Supplementary Table S5). The copy number of vgrG and hcp genes appeared to be different between two subspecies. In summary, SM possesses more secretion system gene cluster genes compared to DQ. Notably, SM harbors two genes specifically related to type I secretion system and type II secretion system.

To determine the differences in gene expression between SM and DQ strains, transcriptome sequencing was performed for two strains under two distinct conditions: cell grown in nutrient broth (in vitro) and cells recovered from infected potato stems (in vivo). RNA-Seq analysis was validated by qRT-PCR with three biological replicates (Supplementary Figure S5). Our analysis revealed distinct gene expression patterns in both strains under different conditions. In SM, a comparison between in vivo and in vitro conditions unveiled 951 upregulated and 758 downregulated genes (Supplementary Table S6). Similarly, DQ exhibited 940 upregulated and 837 downregulated genes in the potato host environment compared to in vitro culture (Supplementary Table S6). Moreover, a comparison of conserved genes between SM and DQ showed differential expression patterns (Figure 4A). Under in vitro conditions, 565 genes displayed differential expression, with 341 genes upregulated and 224 genes downregulated. In the potato host environment, 260 genes exhibited significant differential expression, comprising 112 upregulated and 148 downregulated genes (Supplementary Table S6). Among these, 167 differentially expressed genes were found to overlap between the two conditions, with 73 upregulated and 94 downregulated genes. This suggests that the expression levels of these genes were constitutively different and independent of growth conditions. We observed a large number of strain-specific different expressed genes under in vitro conditions, while a small portion of strain-specific genes changed their expression in host environment (Figure 4A). This observation implied that biochemical and physiological activities of these two strains were distinct when grown in NB medium, although they were phylogenetically similar. Conversely, when infecting plants, these two strains tend to employ similar virulent mechanisms or strategies for host colonization.

Comparing to DQ strain, SM exhibited significant upregulation of T2SS related genes in plant host (Figure 4B). Similar results were also observed for T3SS genes (Figure 4B). Genes related to bacterial motility were consistently upregulated in both conditions, but more significant in vitro, in which genes, such as fliC, flg1, fliD and flgL, exhibited the higher expression levels in NB medium compared to plant host (Figure 4B). Additionally, KEGG pathway analysis revealed that flagellar associated genes were enriched (Supplementary Figure S6). Taken together, our data suggest that SM strain has higher expression of genes related to T2SS, T3SS and motility, which are important virulence traits for bacterial infection.

Given the fact that the flagellar assembly pathway in SM was enriched in transcriptome analysis, we conducted an assay to compare swimming and swarming motility of two strains on NB plates. The results indicated that SM exhibited greater swimming and swarming motility compared to DQ, which is consistent with our transcriptome data (Figure 5A).

Similarly, T2SS upregulation may also increase the secretion of pectinases, cellulases, and proteases, which are critical enzymes for Pectobacterium pathogenicity. To test this hypothesis, we evaluated these enzyme activities on the assay plates. The results indicated that the high virulent strain SM exhibited higher activity of pectinase (Figure 5B), protease (Figure 5C), and cellulase (Figure 5D). In summary, SM exhibits greater swimming and swarming motility, higher activities of PCWDEs than DQ, which explains the higher virulence of SM strain.

To investigate the role of the cellulase in bacterial pathogenicity, we analyzed transcription levels of cellulase genes including glycoside hydrolases, β-glucosidases and cellulose synthases in both strains under in vivo and in vitro conditions. The heatmap showed that the expressions of most glycoside hydrolases were induced in plant in both strains, except GH5-1 and GH8. It is noteworthy that the endoglucanase S (GH12) gene, which was absent in DQ strain, exhibited significant high expression level in SM strain in vitro, with a fold change level exceeding 6 (Figure 6A). This result indicated that endoglucanase S (GH12) may play a crucial role in the process of bacterial infection and be a factor that determined the higher pathogenicity in SM.

To test this hypothesis, we cloned GH12 with its own promoter into a multi-copy plasmid to overexpress GH12 in DQ. The potato slices maceration assay demonstrated that DQ overexpressing GH12 (GH12OE) had higher pathogenicity compared to the wild type strain (Figure 6B). However, GH12OE failed to fully complement the virulence to SM level, implying that other SM specific factors also contributed to its virulence. To further elucidate the functional properties of GH12, we constructed knock-out and complementary expression strains for the GH12 gene and proceeded to evaluate their pathogenicity. The potato slices inoculation assay demonstrated that, compared to potato slices inoculated with the wild-type SM strain and those with the GH12-complemented expression strain (CΔGH12), that inoculated with the GH12 knock-out (ΔGH12) SM strain exhibited smaller lesion areas following maceration (Figure 6C). This observation suggests that the absence of the GH12 gene results in a reduced pathogenicity of SM. Taken together, we show that endoglucanase GH12, which is absent in DQ, significantly contributes to bacterial pathogenicity, likely via accelerated plant cell wall degradation.

Bacterial strains, primers, plasmids and gene sequences are listed in Supplementary Tables S7, S8. Pectobacterium brasiliense SM and DQ were isolated from typical blackleg symptomatic potato stems in Sanming, Fujian Province, and Erdos, Inner Mongolia, respectively, in 2019. These strains are routinely cultured in Difco Nutrient Broth (NB) medium or on Difco Nutrient Agar (NA) at a temperature of 30°C. Escherichia coli DH5α carrying pH2-GH12 plasmid was cultured at 37°C in Luria Broth medium. Kanamycin (50 μg mL⁻¹) was supplemented when needed. Strains were stored in the medium supplemented with 50% glycerol at −80°C.

Three inoculation approaches are conducted in our study. For aerial stem rot, 100 μL of bacterial suspension (10⁸ CFU/mL) was injected into the stem base of potato plants. Inoculated plants were placed in a greenhouse at 21°C with 80% humidity. The plant phenotypes were observed and photographed continuously for 1 week after infection for the calculation of pathogenicity index. The pathogenicity assay of potato blackleg was conducted through a bacterial soil drench experiment. Each pot was filled with 500 mL of bacterial suspension (10⁸ CFU/mL) and the incidence rate was monitored over a two-week period. For potato chip experiment, potato tubers were washed in clean water, rinsed in alcohol, and then sliced and dried on a clean bench. Three potato chips in each replicate were inoculated with 100 μL of bacterial suspension (OD600 = 0.1 or 0.3) by pipetting and incubated at 30°C for 24 h. Maceration area of the lesions were measured using ImageJ. All experiments were conducted in three replicates.

Multilocus sequence analysis (MLSA) was used to conduct phylogenetic relationships among 6 housekeeping genes including gapA, gyrA, icdA, mdh, proA and rpoS from 16 strains of Pectobacterium spp. and 8 strains of Dickeya spp. Amino acid sequences were aligned using MUSCLE. Subsequently, the six housekeeping gene sequences were concatenated in the same order. The evolutionary history was inferred by using the Maximum Likelihood method from MEGA11 software (Tamura et al., 2021), and 1,000 bootstrap replicates were included in a heuristic search by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated by Tamura-Nei model (Ortí et al., 1997).

The genomic DNA was extracted using the SDS method (Lim et al., 2016). The extracted DNA was assessed via agarose gel electrophoresis and quantified using the Qubit® 2.0 Fluorometer (Thermo Scientific). The complete genome of Pectobacterium brasiliense SM and Pectobacterium brasiliense DQ was sequenced at Beijing Novogene Bioinformatics Technology Co., Ltd. using PacBio Sequel platform and Illumina NovaSeq PE150. Single Molecule Real-Time (SMRT) sequencing library was constructed using SMRT bell TM Template kit, version 1.0, with an insert length of 10 kb. Long reads which longer than 6000 bp, were selected as seed sequences using the automatic error correction function in the SMRT portal, while the remaining shorter reads were aligned to these seed sequences using Blasr. The preliminary assembly results were further refined using the Variant Caller module in the SMRT Link software. The corrected assembly result serves as a reference sequence, which is then aligned with Illumina data using Burrows-Wheeler Aligner (BWA). The circularity of the chromosomal sequence was verified by assessing the overlap between the head and the tail. Then the initial site was further refined through BLAST with the DNAa database. Coding genes were retrieved using GeneMarkS program. tRNA genes were identified using tRNAscan (Lowe and Chan, 2016). Genes functions were predicted based on a BLASTP search against the Non-Redundant Protein Database. KEGG annotation results were obtained using the Diamond software. Furthermore, secretory proteins and secretion systems proteins were predicted using SignalP database (Bendtsen et al., 2004) and Swiss-Prot (The UniProt Consortium, 2023), respectively.

Average Nucleotide Identity (ANI) values were calculated for the genomes using online Average Nucleotide Identity Calculator (EzBioCloud; Yoon et al., 2017). In silico DNA–DNA hybridization (DDH) was calculated using the Genome-to-Genome Distance Calculator (GGDC; Auch et al., 2010). The synteny analysis results of the genomic sequence were obtained by Mauve (Darling et al., 2004) using SM genome as the reference genome. The genome alignments were performed with BLAST+ using BRIG (Alikhan et al., 2011). Orthology assessment was conducted through CD-HIT rapid clustering of similar proteins to obtain core genes and specific genes. Each annotated gene is classified into a pan-genome structure, including core (present in all genomes analyzed), or unique (present in only one genome) pan-genome fragment.

Total RNA samples were quantified using spectrophotometry. RNA quality assessment was performed using the Agilent 2,100 Bioanalyzer (Agilent Technologies, CA, United States) at Beijing Novogene Bioinformatics Technology Co., Ltd. To remove rRNA and purify mRNA from total RNA, probes were employed. Initially, cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase. DNA polymerase I and RNase H were used to synthesize the second strand of cDNA. PCR was then performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index primers. Subsequently, PCR products were purified using AMPure XP beads, and the library quality was assessed using the Agilent Bioanalyzer 2100 system. After cluster generation, libraries were prepared and sequenced on the Illumina Novaseq platform, generating 150 bp paired-end reads. The FPKM of each gene is calculated based on the length of the gene and the number of reads of the gene (Mortazavi et al., 2008). Differential expression analysis was performed using the DESeq R software package (1.18.0). The resulting p-values were adjusted using Benjamini and Hochberg’s method to control the false discovery rate (Benjamini and Hochberg, 1995). KOBAS software was utilized to identify the statistical enrichment of differentially expressed genes in the KEGG pathway.

Swimming and swarming motility were tested by inoculating bacterial suspension on nutrient agar plates solidified with 0.3 and 0.6% agar, respectively. The colony diameter was measured after 24 h of incubation. Detection of pectinase activity was assayed on M63 polygalactonate (PGA) plate. After incubation for 48 h, the culture medium was stained with 10% copper acetate aqueous solution. Degradation of PGA results in the appearance of a white halo around bacterial colonies, which is measured to assess total pectinase activity. The ability to produce cellulases was analyzed on M63 with 0.2% glycerol and 1% carboxymethylcellulose. After incubation for 48 h, the plates were flooded with 0.1% Congo red solution for 10 min and then washed with 1 M NaCl for 5 min. The diameter of the halo that appeared around the colonies indicated cellulase activity and was subsequently measured. Determination of protease activity on medium containing 1% skim milk. After 24 h of incubation, measure the diameter of the transparent halo around the colony, reflecting protease activity. All the plate assays mentioned above were inoculated with 2 μL of bacterial suspension containing 10⁸ CFU ml⁻¹. All experiments were performed three times with three replicates.