Konjac corms with symptoms of soft rot were collected from Qujing City, Yunnan Province and Ankang City, Shaanxi Province, China in 2019 and 2020 (Supplementary Table 1). Infected tissues were cut into small pieces and were sterilized in 75% (v/v) ethanol for 30 s, followed by three times rinses with sterile distilled water. The tissues were then mashed and diluted in sterile distilled water. A volume of 200 μl bacterial suspension from each dilution was spread on nutrient agar (NA) medium and incubated at 28°C for 48 h. A single colony was then picked and subcultured in nutrient broth (NB) medium. Pure cultures were obtained through successive streaking on NA.

For konjac slice assay, bacterial strains were grown overnight, washed and resuspended into sterile water with an OD₆₀₀ = 0.2. The tubers were sliced 0.5 cm thick and placed in a plastic food container containing wet paper tissues. Tuber slices were then inoculated with 20 μl of bacterial solution and incubated at 28°C for 48 h. Sterile water was used as a negative control. For in vivo assay, bacterial suspension was prepared with a concentration as mentioned above. Stems of 6-months-old konjac seedlings were then inoculated using a pin-prick inoculation method under greenhouse conditions. Similarly, sterile water was used as a negative control. The assay was repeated twice independently.

Pectinolytic activity assay for bacterial isolates was performed on crystal violet pectate (CVP) medium (Hélias et al., 2012). Briefly, pure bacterial suspensions were diluted into 10²–10³ CFU/ml. A volume of 100 μl dilution was then plated on a CVP medium and incubated at 28°C for 48 h.

Eleven pathogenic strains (per isolate per sample) were selected and then sequenced. Genomic DNA extraction, sequence library construction and sequencing were conducted at the Beijing Novogene Bioinformatics Technology Co., Ltd. Briefly, total DNA from each isolate at the exponential stage was extracted and assessed using agarose gel electrophoresis and Qubit 2.0 Fluorometer (Thermo Scientific). After library construction, sequencing was performed using a Nanopore PromethION platform and an Illumina NovaSeq PE150 platform. For most bacterial isolates, a hybrid assembly pipeline was conducted by Unicycler v0.4.8 using both Illumina reads and long reads (Wick et al., 2017). The genome of strain QJ003 was assembled by Raven v1.5.1 (Vaser and Šikić, 2021) and further polished by Pilon v1.24 (Walker et al., 2014). The quality of all genome assemblies was assessed by BUSCO v5.2.2 based on the dataset bacteria_odb10 (Manni et al., 2021). Prokka v1.14.5 was used for genome annotation (Seemann, 2014).

Pairwise average nucleotide identity (ANI) values were calculated by a Python module pyani v0.2.11 using ANIm method.¹ isDDH was estimated via the Genome-to-Genome Distance Calculator 3.0 web server using the recommended formula 2² (Meier-Kolthoff et al., 2021).

In addition to our 11 new assemblies, 53 complete genome sequences of Pectobacterium were retrieved from GenBank in May 2021.³ The accession numbers and other basic information for these downloaded genomes are provided (Supplementary Table 2). After confirming that all 11 P. aroidearum assembly sequences began with the gene dnaA, bacterial synteny was evaluated and visualized by multiple whole-genome alignment using progressiveMauve (Darling et al., 2010). Dot plot between two genomes was generated by D-Genies web server⁴ (Cabanettes and Klopp, 2018). A comparative genomic ring plot was generated using BLAST Ring Image Generator (BRIG; Alikhan et al., 2011).

Parsimony trees were inferred with kSNP v3.1.2 based on pan-genome SNPs (Gardner et al., 2015). The optimum k-mer size used for the Pectobacterium species and P. aroidearum isolates was 21 and 19, respectively, which was determined by subcommand Kchooser.

A maximum-likelihood (ML) tree for the P. aroidearum strains was inferred by FastTree 2.1.11 using generalized time-reversible (GTR) models based on core-gene alignment (Price et al., 2010). The Newick tree files were visualized using the online program iTOL v5 (Letunic and Bork, 2021) and MEGA11 (Tamura et al., 2021).

The Roary pipeline was used to infer the pangenome and a gene presence/absence matrix of Pectobacterium spp. and P. aroidearum, respectively (Page et al., 2015). The parameters for minimum blastp percentage identity (−i) and core gene (−cd) were adjusted to 90% and 100%, respectively.

For enrichment analysis, a COG functional category was assigned to the protein-coding genes of P. aroidearum QJ036 using eggNOG-mapper v2.1.4 (Cantalapiedra et al., 2021). Functional enrichment analysis was performed using Fisher’s exact test within the R environment (v4.0.0). p values for multiple comparisons were adjusted with the Benjamini and Hochberg method (Benjamini and Hochberg, 1995).

Protein secretion systems (types I, II, III, IV, V, and VI) were detected by uploading the protein sequences to TXSScan/MacSyFinder (Abby et al., 2014, 2016). SecRet6⁵ was used to predict the T6SS-associated proteins (Li et al., 2015). Domain identification and protein family classification was conducted by InterProScan (Blum et al., 2021). Multiple sequence alignment was then performed by T-Coffee (Di Tommaso et al., 2011) and visualized by ENDscript (Robert and Gouet, 2014).

PCWDE proteins of the genera Pectobacterium and Dickeya were downloaded from the UniProt database and used for building a reference database. BLAST+ (v2.6.0) was applied to find putative PCWDEs in 14 P. aroidearum genomes using an E-value threshold (E-value < 1e − 5). Only hits with at least 35% identity and 50% query coverage were kept. If a protein annotation from BLAST+ results was not consistent with that derived from eggNOG-mapper (v2.1.4), this protein annotation was examined manually.

Genomic islands (GIs) were predicted using the IslandViewer 4 webserver.⁶ IslandViewer 4 integrated three GI prediction methods (IslandPath-DIMOB, SIGI-HMM and IslandPick) as well as annotations of virulence factors, pathogen-associated genes, and antimicrobial resistance (AMR) genes.

To further determine the causal agent of konjac soft rot in China, A. konjac tubers with foul-smelling rot symptoms were collected from Qujing, Yunnan Province and Ankang, Shaanxi Province, China. After bacterial isolation and purification, a total of 11 isolates were shown to cause typical soft rot symptoms on konjac tubers in vitro (Figure 1A; Supplementary Figure 1). As expected, the pectolytic ability of these candidate pathogens was further confirmed based on the formation of deep cavities on crystal violet pectate (CVP) media (Figure 1B; Supplementary Figure 2). 16S rRNA gene sequencing and BLAST analysis showed that all 11 16S rRNA genes share at least 99% identity with that of P. aroidearum (NR_159926), which is consistent with previous findings (Xie et al., 2020). Since all candidate pathogens belonged to the same species based on 16S rRNA gene sequences and morphological similarity, only one strain OJ036 was selected for pathogenicity test in vivo. Stem rot and wilting symptoms were visible only within 2 days for seedlings inoculated with strain QJ036 while no symptoms were observed for the control group (Figure 1C). This represents strong evidence that the causal agent of konjac soft rot is P. aroidearum.

In addition, we did pathogenicity tests in vitro to explore the host range of P. aroidearum. Interestingly, P. aroidearum was found to cause typical rot symptoms on sweet potato (Ipomoea batatas (L.) Lam.), jícama (Pachyrhizus erosus (L.) Urb.), yacón (Smallanthus sonchifolius (Poepp.) H.Rob.) and taro (Colocasia esculenta (L.) Schott), which has not been reported (Figure 1D). Although in vivo pathogenicity tests are needed, these results suggest a broader host range of P. aroidearum.

To further confirm the taxonomic status and investigate the genomic diversity of konjac soft rot pathogens, whole-genome sequencing was performed for the 11 strains using a Nanopore PromethION platform and an Illumina NovaSeq platform with more than 100x coverage depth (Supplementary Table 3). After hybrid genome assembly using filtered data, one circular chromosome without gaps was obtained for each isolate. High scores (>99%) were achieved from a BUSCO assessment of genomic completeness, indicating the high quality of our assemblies (Supplementary Figure 3). The parsimony tree based on genome-wide single nucleotide polymorphisms (SNPs) clustered these Pectobacterium spp. into seven well-resolved clades, of which clade VII can be further clustered into three subclades (Figure 2). Notably, all konjac soft rot pathogens form a monophyletic group (clade IV) with P. aroidearum L6, P. carotovorum subsp. carotovorum PC1 and P. carotovorum subsp. carotovorum PCCS1, indicating that these strains belong to the same species (Figure 2). This result is further supported by average nucleotide identity (ANI) and digital DNA–DNA hybridization (dDDH) analyses, two widely used methods for the taxonomy of prokaryotes (Chun et al., 2018). Sequence comparisons between the trains in clade IV show that ANI is greater than 95% and dDDH greater than 70%, which are above the cut-off values for species delineation (Supplementary Table 4). Taken together, these results clearly indicate that all these strains in clade IV should be classified as P. aroidearum and that three strains (P. carotovorum subsp. carotovorum PC1, P. carotovorum subsp. carotovorum PCCS1 and P. carotovorum subsp. carotovorum PCC21) are incorrectly named.

Overall genomic features of P. aroidearum strains, including genome size, GC (guanine-cytosine) content, number of protein-coding sequences (CDS), and number of RNA genes, were quantified (Table 1). The length of P. aroidearum genomes ranges from 4,865,541 bp (QJ315) to 5,057,072 bp (QJ003). GC content varies between 51.6% and 51.9%. The highest number of putative CDS was observed in AK042 (4,469), whereas the lowest number of CDS was found in QJ311 (4,277). In addition, the number of tRNA, rRNA and CRISPRs (Clustered Regularly Interspaced Short Palindromic Repeats) is the same among our 11 genome assemblies. We found no evidence of plasmids in P. aroidearum.

For synteny analysis, the program Mauve was used to perform whole-genome multiple alignment to detect genomic rearrangement. Interestingly, a couple of large-scale chromosomal inversions were observed among these 11 P. aroidearum isolates (Figure 3A). These strains could be clustered into three groups based on genomic rearrangement, which is further supported by the phylogenetic analysis using genome-wide SNPs (Figure 3A) and BLAST comparisons of the whole genome (Supplementary Figure 4). Except for strain QJ311, the grouping of the other strains is consistent with their geographical distribution (Supplementary Table 1). To further confirm the above results, a total of three strains (QJ002, QJ036 and AK042) selected from each group were used for pairwise genome alignment. As expected, dot plots also revealed the genomic inversions between QJ036 and either QJ002 or AK042 (Figure 3B).

The major mechanisms underlying chromosomal rearrangements are recombinational exchanges between homologous sequences such as ribosomal operons and mobile genetic elements (MEGs) including transposons, insertion sequence (IS) elements, and prophages (Raeside et al., 2014). The origins of the chromosomal inversion were investigated by examining the sequences bordering each inversion. The inversion that occurred between QJ036 and QJ002 resulted from homologous recombination of ribosome operons which are located exactly at the border regions of this inversion (Figure 3A). However, the reason for the inversion between QJ036 and AK042 is unknown.

The pangenomes of plant pathogens are often associated with species- or strain-specific virulence, host specificity or adaptive potential, or evolutionary history (Amir et al., 2020). To investigate the genomic plasticity of genus Pectobacterium, the 64 complete genomes from 9 Pectobacterium spp. including our 11 assemblies were used to identify the core- and pan-genomes. The number of core genes was 2,228, which only accounts for 11.31% of the pangenome, which contained a total of 19,698 genes (Supplementary Figure 5A). As more genomes were added, the pangenome trend showed a gradual expansion, implying an open pangenome of Pectobacterium strains (Supplementary Figure 5B).

Similarly, a pangenome of 6,630 genes was identified from the 14 P. aroidearum genomes, with a core genome of 3,575 genes (53.92% of the pangenome; Figure 4A). The number of accessory genes of each strain varied from 701 (PC1) to 898 (QJ003; Figure 4B). The pangenome fitted cumulative curve showed that with the addition of more genomes, the number of core genes remained relatively stable while the number of total genes continued to increase, indicating that the pangenome is still open (Figure 4C). In addition, a maximum likelihood (ML) tree based on core gene alignment is congruent with the matrix representing the presence and absence of core and accessory genes (Figure 4D).

Functional enrichment analysis was performed to uncover the biological roles of core and accessory genes using the QJ036 genome as an example. We classified all genes in the QJ036 genome into three groups including genus-core genes (2228), species-core genes (1347) and accessory genes (729). Enrichment analysis based on the COG (Clusters of Orthologous Groups of proteins) database showed that genus-core genes are overrepresented in categories associated with essential life activity such as energy production and conversion, the transport and metabolism of amino acid, nucleotide and lipid, and translation, ribosomal structure and biogenesis (Table 2). The species-core genes are enriched for virulence-associated categories, including the transport and metabolism of carbohydrate and inorganic ion, cell motility, and intracellular trafficking and secretion (Table 2). Interestingly, the category of replication, recombination and repair is the only enriched COG for the accessory genes (Table 2). Further investigation listed all genes belonging to this category and found that a majority of these genes were annotated as recombinase/integrase and mobile genetic elements such as transposons, insertion sequence (IS) elements (Supplementary Table 5).

Bacterial secretion systems and plant cell wall degrading enzymes (PCWDEs) play a key role in the interaction between soft rot bacteria and host plants. We mined the 14 P. aroidearum genomes to identify and compare these determinants of pathogenicity. All six secretion systems (types I–VI) were detected using the program TXSScan. Type I secretion system (T1SS), type II secretion system (T2SS), type III secretion system (T3SS) and type VI secretion system (T6SS) were conserved in all tested P. aroidearum strains (Figure 5A). However, type IV secretion system (T4SS), a versatile secretion system that is involved in protein translation, bacterial conjugation and DNA uptake/release, is not detected in several P. aroidearum genomes (Figure 5A). Notably, P. aroidearum L6 has two copies of each T4SS subtypes (typeG and typeT) that are classified according to different mating pair formation complexes (MPF; Guglielmini et al., 2013). For T6SS, although each P. aroidearum strain has a conserved T6SS cluster, the number of predicted effectors and immunity proteins varies considerably (Supplementary Table 6).

Although the presence of type V secretion system (T5SS) in Pectobacterium spp. is still inconclusive (Li et al., 2018; Van Gijsegem et al., 2021), two subtypes of T5SS were identified by TXSScan in our assemblies. To further confirm the prediction, sequence analysis was performed using putative T5SS proteins from the QJ036 strain. In two-partner secretion (TPS or type 5b), the passenger domain (TpsA protein) and the β-domain (TpsB protein) are encoded by two separate genes which are frequently, but not always, encoded in an operon (Wells and Henderson, 2013). Domain analysis based on the InterPro database showed that each putative T5bSS protein (QJ036_01923, QJ036_02058 and QJ036_00084) identified by TXSScan contains a signal sequence, two polypeptide-transport-associated (POTRA) domains and a shlB domain (β-domain; Figure 5B). In addition, these proteins displayed high similarity in primary and secondary structure with the known TpsB protein CdiB of Escherichia coli (Figure 5C). These results strongly indicate that these three proteins are members of TpsB family. To find the TpsA proteins in P. aroidearum, TpsA protein fhaB (Uniprot: P12255) of Bordetella pertussis was used as a query to search against protein database of QJ036 using BLASTP. Domain analysis of the three BLASTP hits (QJ036_01921, QJ036_02059, and QJ036_00083) shows that each contains a signal sequence, a haemagg_act domain and a filamentous hemagglutinin repeat region (Figure 5B). As expected for a type 5b, these three TpsB proteins are very close to their corresponding TpsA proteins according to the locus tags. These results therefore provide evidence for the existence of three copies of T5bSS in P. aroidearum QJ036. However, the presence of T5aSS is not well supported by our sequence analysis results.

PCWDEs seem to be very conserved in the P. aroidearum strains examined in this study. In total, 21 PCWDEs genes were found in each strain including two cellulase genes (celV, celS), one oligogalacturonate lyase gene (ogl), eight pectate lyase genes (pel1, pel2, pel3, pelX, pelY, pelW, pelE, pelL), two pectin acetylesterase genes (paeX, paeY), one pectin lyase gene (pnl), two pectinesterase genes (pemA, pemB), three polygalacturonase genes (pehX, pehK, pehN, pehA), two rhamnogalacturonan lyase genes (rhiE, rhiN) and one protease gene (prtC; Supplementary Table 7).

Genomic islands (GIs) are clusters of consecutive genes likely acquired via horizontal gene transfer (HGT), which may facilitate microbial adaptation by disproportionately encoding factors involved in virulence or antimicrobial resistance (Bertelli et al., 2019). Here, GIs of P. aroidearum isolates were identified using the online webserver IslandViewer 4 (Bertelli et al., 2017). The location of GIs in each genome was visualized, which revealed that P. aroidearum isolates with similar GIs distribution patterns tended to be from the same sampling site (Supplementary Figure 6). In addition, the number of GI genes in each genome ranged from 454 (PC1) to 866 (AK042), accounting for about 10%–20% of the total genes of a genome (Figure 6A). However, none of these genes were annotated as virulence/resistance factors, and the majority of these genes were annotated as hypothetical proteins. To further explore the gene content of GIs, we collected all GI-derived genes with a gene ID and calculated the occurrence number of each gene (Supplementary Table 8). The proteins encoded by the top 11 genes include prophage integrase (IntA, IntS), tyrosine recombinase XerC, transporter (YflS, CitN), major exported protein HcpA, DNA-binding transcriptional repressor YiaJ, Acetyl-CoA:oxalate CoA-transferase YfdE, ribose import permease RbsC, Tyrocidine synthase 3 TycC and DNA topoisomerase 3 TopB (Figure 6B).