European Population of Pectobacterium punjabense: Genomic Diversity, Tuber Maceration Capacity and a Detection Tool for This Rarely Occurring Potato Pathogen

Enterobacteria belonging to the Pectobacterium and Dickeya genera are responsible for soft rot and blackleg diseases occurring in many crops around the world. Since 2016, the number of described species has more than doubled. However, some new species, such as Pectobacterium punjabense, are often poorly characterized, and little is known about their genomic and phenotypic variation. Here, we explored several European culture collections and identified seven strains of P. punjabense. All were collected from potato blackleg symptoms, sometimes from a long time ago, i.e., the IFB5596 strain isolated almost 25 years ago. We showed that this species remains rare, with less than 0.24% of P. punjabense strains identified among pectinolytic bacteria present in the surveyed collections. The analysis of the genomic diversity revealed the non-clonal character of P. punjabense species. Furthermore, the strains showed aggressiveness differences. Finally, a qPCR Taqman assay was developed for rapid and specific strain characterization and for use in diagnostic programs.


Introduction
Bacteria belonging to the Pectobacterium and Dickeya genera are causal agents of soft rots in a wide variety of host plants [1]. The agricultural and economic impact ranks these pathogens classified within the group of soft rot Pectobacteriaceae (SRP) [2] among the most studied phytopathogenic bacteria [3]. They produce extracellular enzymes degrading plant cell walls [4,5] causing the breakdown of plant tissues. Several of the SRP can generate blackleg symptoms, characterized by a blackening of the stem base [6]. The losses can be significant [7] and, to date, chemical, physical or biocontrol applications have shown limited efficacy. Therefore, prophylactic methods such seed sanitation and breeding potato cultivars to improve natural resistance are still highly recommended [8].
Since the wide-scale use of next generation sequencing (NGS) for strain characterization, it has become clear that the genetic diversity of SRP is very high [9]. Analysis of NGS data resulted in the elevation of different subspecies to the species level [10] and also in the description of new species. As a consequence, the number of species belonging to these two genera has increased from five Pectobacterium spp. and seven Dickeya spp. before 2016 to 18 Pectobacterium spp. and 12 Dickeya spp. until now [11,12]. The NGS analysis further resulted in the revision and refining of the taxonomic position of strains maintained in laboratory or reference collections, sometimes for a long time [10,13]. The isolates were collected from different backgrounds, including diseased plant samples and surface waters [14][15][16][17][18][19][20].
Historically, the isolates collected in temperate regions from blackleg-diseased plants were identified as Erwinia atroseptica, Erwinia carotovora subsp. wasabiae and Erwinia chrysanthemi [21][22][23]. On the basis of the current taxonomy, P. atrosepticum, P. brasiliense, P. parmentieri, P. versatile, P. polaris and P. carotovorum are the most prevalent Pectobacterium species responsible for soft rot and/or blackleg diseases [24][25][26][27][28][29][30]. It should be noted, however, that the ability of pectinolytic bacteria to grow as well as develop blackleg and soft rot symptoms depends on many environmental factors such as temperature and humidity, in a species-dependent way [31]. Concerning Dickeya spp., D. dianthicola and the clonal pathogen D. solani are still the only two species isolated from blackleg symptoms in European potato fields [32][33][34].
Our study focused on P. punjabense species. The type strain P. punjabense SS95 T was isolated in 2017 from blackleg symptoms sampled in a Pakistani potato field. The average nucleotide identity and in silico DNA-DNA hybridization values allowed one to distinguish this strain as belonging to a new species, taxonomically close to P. parmentieri and P. wasabiae. In addition, other genomic analyzes showed that several genes of P. punjabense SS95 T involved in siderophore transport and metabolite uptake are absent in P. parmentieri and P. wasabiae genomes [45]. Since this initial description, there are no records of the pathogen from other countries. We investigated the distribution, the genetic diversity and phenotypic variation of Pectobacterium punjabense in Europe by the analysis of strains deposited in culture collections. We also describe the development of a qPCR TaqMan assay for specific detection of P. punjabense.

Bacterial Strains
In order to assess the presence of P. punjabense in European potato soft rot Pectobacteriaceae collections, strains of the French RNS collection from inov3PT, the Dutch IPO collection from Wageningen University and Research and the Polish IFB collection from the Intercollegiate Faculty of Biotechnology were screened by sequencing of the housekeeping genes gapA and dnaX [46,47] or gyrA, rpoA, rpoS and recA for the Polish collection [14,48,49]. Four strains, namely, P9A19a, RNS08.28, IPO3715 and IFB5596 (Table 1), all originating from blackleg symptoms and isolated on CVP medium [50], showed gapA and dnaX sequences similar to those of the P. punjabense SS95 type strain. These strains, stored in glycerol at −80 • C, were plated on Tryptone Yeast extract (TY) rich medium and were used in all following steps of this work. Seven type strains of other Pectobacterium species, namely, P. polonicum DPMP315 T , P. wasabiae CFBP3304 T , P. parmentieri RNS08.42.1a T , P. betavasculorum NCPPB2795 T , P. zantedeschiae 9M T , P. peruviense IFB5232 T and P. atrosepticum CFBP1526 T , were also used as reference in genomic and phenotypic comparative assays. The draft genomes of the four P. punjabense candidates were sequenced using a MiSeq Paired-End Illumina technology. Then, paired-end reads of each strain were assembled using the de novo assembly tool of CLC genomics workbench 10.1.1 (https://digitalinsights. qiagen.com accessed on 12 February 2021). The draft genome sequences were deposited in GenBank to which accession numbers have been assigned. The genome sequences representative of the seven Pectobacterium species closest to P. punjabense, as well as those of P. punjabense SS95 T , were used as references in the comparative analyses (Table 1).

Multi-Locus Sequence Analysis (MLSA)
The sequences of the fourteen housekeeping genes (acnA, fusA, gapA, glyA, groEL, gyrB, mdh, mtlD, purA, recA, rplB, rpoD, rpoS, secY) of P. punjabense candidates and reference type strains, recovered from the corresponding genome sequences, were concatenated, aligned and used to build a phylogenetic tree using the maximum composite likelihood method with MEGA X [51]. Bootstrap values were calculated from 1000 replicate iterations. The average nucleotide identity (ANI) calculator of the EzBioCloud server [52] was used for genome comparison and to confirm the taxonomic assignation of each P. punjabense candidate genome. In addition, the Genome-to-Genome Distance Calculator of the Leibniz Institute DSMZ [53] was used to calculate the DNA-DNA hybridization (DDH) value in order to obtain a second robust parameter for each genome comparison. Then, read mapping of each P. punjabense candidate to the P. punjabense SS95 T reference genome was performed. Finally, the variant detection tool of CLC genomics workbench 10.1.1 (https://digitalinsights.qiagen.com accessed on 12 February 2021) was used to estimate the number of SNPs for each P. punjabense candidate, with the P. punjabense SS95 T genome as a reference.

BRIG Analysis
A sisual comparison of genome homology was performed using BRIG (BLAST Ring Image Generator) [54]. P. punjabense SS95 T was used as the reference genome and was compared to the other genomes of P. punjabense. The analysis was performed with the default settings.

Core Protein Phylogeny
A core phylogenetic analysis conducted from the protein sequences deduced from the concatenated core genes shared by all the genomes was compared with the clustering obtained in the multi-locus sequence analysis (MLSA). The phylogenomic analysis was performed with the PhyloPhlAn computational pipeline (https://huttenhower.sph. harvard.edu/phylophlan accessed on 18 January 2021), which uses the most conserved universal proteins from full proteomes to extract the phylogenetic signal [55]. The maximum likelihood tree was built based on 399 protein sequences. The core protein analysis was performed using the genomic sequences from the type and reference strains of the most closely related members of the genus Pectobacterium.

Carbon Source Utilization Profiles
P. punjabense strains (IFB5596 and 139A2, both co-isolated from the same sample in 1996, as well as RNS0828a, P9A19a and IPO3715) and strains of other closely related species such as P. atrosepticum (SCRI1043, ICMP1526 T ), P. betavasculorum (CFBP2122 T, CFBP5271), P. peruviense (IFB5232 T , IFB9229), P. parmentieri (SCC3193, IFB5322), P. wasabiae (CFBP3304 T , CFBP3308) and P. zantedeschiae (9M T and 2M) were analyzed for the usage of different carbon sources using GEN III plates, and the results were compared with BIOLOG utilities. The assays were performed according to the manufacturers' instructions.

Fatty Acid Methyl Esters (FAME) Composition
The whole-cell derived fatty acids were isolated directly from bacteria growing on TSA plates for 2 days at 28 • C. After saponification, methylation, extraction and washing according to the procedure described by Sasser et al. [56], the obtained fatty acid methyl esters (FAME) were analyzed using a gas chromatography-FID detector (Agilent 7820A) and an Ultra 2-HP capillary column (cross-linked 5% phenyl methyl silicone, 25 m, 0.22 mm id and 0.33 µm film thickness). A standard No. 1200-A solution (MIDI Inc., Newark, DE, USA) was used to calibrate the GC. Qualitative and quantitative analyses of the fatty acids were performed using the TSBA library (ver. 6.2B) and Sherlock Microbial Identification System software (MIDI Inc., Newark, DE, USA).

Tuber Maceration Assay
Assays were conducted in order to compare the tuber maceration capacity of P. punjabense to that of other Pectobacterium reference strains. In addition to P. punjabense P9A19a, RNS08.28, IPO3715 and SS95 T , three P. punjabense isolates RNS16-153-1A, RNS18-61 and RNS18-78 were also used in the tuber maceration assay. These three isolates were collected in 2016 or 2018 from blackleg symptoms during field surveys in France and were identified as P. punjabense based on dnaX and gapA housekeeping gene sequences ( Figure S1) in the course of this work. The tuber assay was adapted from previous works, as described in Blin et al. [34]. Inoculum suspensions of each of the tested strains were prepared from overnight cultures in TY broth and calibrated at 10 8 cfu/mL. Two pipette tips containing 10 µL inoculum suspension each were driven into the flesh of surfacedisinfected potato tubers (cv. Bintje). Eight tubers per strain were thus inoculated for a total of sixteen inoculation points per strain. A negative control with 10 µL of NaCl 0.8% was used.
Inoculated potato tubers were incubated at 25 • C in the dark in humid chambers with water-saturated air. After five days, tubers were cut across inoculation points and symptom severity was scored based on a six-class visual scale according to rot extension. To each class, a coefficient was assigned whose value increased with symptom severity (0, 0.2, 0.4, 0.6, 0.8 and 1). A pathological index (from 0 to 100) representative of the aggressiveness of each strain was calculated as follows: (1) Firstly, the data were analyzed using a Kruskal-Wallis non parametric rank test. Secondly, a pairwise post hoc Tukey test was performed to discriminate strains if significant differences were detected by the Kruskal-Wallis rank test.
Using these P. punjabense "specific regions", primers and probes were designed using Primer3 [57]. In silico specificity of candidates was tested by a local blast on all the Pectobacterium genomes of the study, as well as with the Standard Nucleotide BLAST on the NCBI Nucleotide collection database. A qPCR TaqMan assay specific for P. punjabense was developed on a LightCycler96 Instrument real-time PCR system (Roche Life Science). The amplicon has an expected size of 217 bp and is located in a gene coding for a hypothetical protein. Amplification reactions were performed with 10 µL of FastStart Essential DNA Probes Master 2x (reference 06402682001, Roche Life Science), 3 µL of 3.33 µM primers P. punjabense 2b-F TCCTTCAGCCAGAGAACCAG, 3 µL of 3.33 µM primers P. punjabense 2b-R AACAACAATACCGGCAAGTGG, 2 µL of 2 µM probe P. punjabense 2b TGCAGGC-CTTGTAACTCCGCT labelled with a 5 reporter dye (FAM-6-carboxyfluorescein (FAM)) and a 3 -end quencher dye (BHQ1), synthesized by Eurofins Genomics (Ebersberg, Germany). Finally, 2 µL of template DNA was added to a final volume of 20 µL. Cycling parameters were 95 • C for 10 min, 40 cycles of 95 • C for 10 s and 60 • C for 1 min. Taqman assays were performed on a large number of strains to test the specificity of the tool. Pure DNA samples were prepared using the MasterPure™ Complete DNA kit (Lucigen) and diluted at 1 ng/µL. A serial dilution of the P. punjabense SS95 T extracted DNA ranging from 2 ng to 2.5 fg was used to determine the sensitivity of the Taqman assay. Each quantity was tested in triplicate.

Identification and Diversity of P. punjabense Candidates
The MLSA phylogenic tree including strains tentatively identified as P. punjabense derived from European collections and type strains from phylogenetically related Pectobacterium species was built from the concatenated sequences of fourteen housekeeping genes extracted from genomes. All P. punjabense candidates grouped together with the type strain P. punjabense SS95. The resulting cluster presented a robust bootstrap value ( Figure 1). The core protein phylogeny confirmed the status of P9A19a, RNS08.28, IPO3715 and IFB5596 as P. punjabense strains (Figure 2).  The ANI and DDH values were calculated to consolidate the assignation of the European isolates P9A19a, RNS08.28, IPO3715 and IFB5596 to P. punjabense species. When the ANI and DDH values were calculated between P. punjabense genomes (Table 2), all values were higher than the accepted cut-off value for species delineation: 95 for ANI values and 70 for DDH [58,59]. The ANI and DDH values of P. punjabense isolates ranged from 98.6 to 100 and from 88.6 to 100, respectively. These infra-specific variations highlight the non-clonal character of P. punjabense strains. The ANI and DDH parameters also confirmed that P. polonicum DPMP315 is closely related to P. punjabense strains, with a mean ANI value calculated at 93.9 and a DDH of 54.9. Then, P. parmentieri RNS08.42.1a and P. wasabiae CFBP3304 are the two other strains closest to the P. punjabense strains with ANI values ranging between 91.0 and 91.5 and DDH values between 43.3 and 44.4. Table 2. Average nucleotide identity (ANI) and DNA-DNA hybridization (DDH) values for each genome comparison. Comparison between P. punjabense strains are marked in green.

Strain
Code (1) In addition, a variant detection analysis to estimate the number of SNPs between each P. punjabense strain and the type strain SS95 yielded approximately 35,000 SNPs for RNS08.28 and IPO3715. For IFB5596, approximately 33,000 SNPs were found. Surprisingly, no SNP was found for P9A19a, suggesting a clonal origin of these Pakistani and French isolates. Finally, pairwise ANI values, SNP analysis and also BRIG analysis used for genome comparison ( Figure S2) confirmed that the four European P. punjabense strains are genetically different.

Physiological and Structural Characteristics
The physiological and biochemical features of all the P. punjabense strains were very similar in the BIOLOG GEN III plates system ( Figure 3). None of the single carbon sources alone could discriminate P. punjabense from the other Pectobacterium species. Noticeably, all P. punjabense strains assimilated D-glucuronic acid and L-serine as the sole carbon source, while most other strains tested could not. Only P. punjabense, P. parmentieri and P. polonicum grew in the presence of the inorganic compound potassium tellurite. P. punjabense strains displayed a high sensitivity for different antibiotic compounds such as nalidixic acid, aztreonam, minocycline and formic acid compared with the other species. In particular, the two strains of P. betavasculorum were able to grow in the presence of each of these four antibiotic compounds.

Occurrence of P. Punjabense in Collections
The three bacterial collections used in our study gather SRP collected over many years from different sources, although strains from blackleg=diseased plants are overrepresented. For collection, a CVP medium has been used that allows the isolation of all SRP with largely the same efficiency [50]. We therefore could make an estimate on the isolation incidence of P. punjabense relative to other SRP strains (Table 3). For example, over the five sampling campaigns between 2015 and 2019, 1663 strains were deposited in the RNS collection, including only four P. punjabense strains, all collected from different fields. In the IPO collection, among 1012 SRP characterized between 1963 and 2020, only one P. punjabense strain was identified. From 2031 strains sampled in Poland between 1996 and 2014 (including 1500 strains isolated in 1996), four P. punjabense isolates were identified, but all were isolated from the same blackleg diseased plant, and they showed identical housekeeping genes sequences.

Aggressiveness
The P. punjabense strains varied in their tuber maceration capacity (Figure 4). For strains IPO3715, RNS08-28, SS95 T and RNS18-78, the pathological index was significantly lower than that of RNS 18-61, which had an index similar to that of the reference strain P. zantedeschiae 9M T , previously described as a strain with a high maceration capacity [18]. Strains P9A19a and RNS16-153-1a presented a pathological index of 71 and 79, respectively, indicating an intermediate aggressiveness level compared to the other strains of the panel. Finally, none of the P. punjabense strains tested showed a potential maceration capacity equivalent to the referent aggressive strain P. atrosepticum CFBP6276 under the conditions of the experiment.

Development of a qPCR TaqMan Assay Specific for P. Punjabense
The qPCR Taqman 2b assay gave a strong signal (Ct values ranging from 17.58 to 20.29) with DNA extracted from the eight P. punjabense strains (the five strains used for the genomic work and the three additional isolates included in the pathogenicity assay), but not with DNA from 59 different Pectobacterium and Dickeya species, all calibrated at 1 ng/µL (Table 4).
To determine the threshold level, a ten-fold serial dilution of P. punjabense SS95 T DNA ranging from 2 ng to 2.5 fg was tested. The detection threshold of the qPCR Taqman 2b was obtained at 20 fg of P. punjabense DNA, with a mean detection value of 34.52 Ct and a standard deviation of 0.98 ( Table 5). The analysis of raw curves generated by the LightCycler96 Instrument real-time PCR system allowed fixing the efficiency of the qPCR reaction at 1.99, which is very close to the expected theoretical value of 2.00.

Discussion
The recent development of next generation sequencing technologies considerably improved knowledge of bacterial genomics including those of the genus Pectobacterium. It became evident that this genus is highly diverse, and NGS analysis also led to improved identification of strains in culture collections and the description of several new Pectobacterium species. This is the case of P. punjabense species, whose type strain SS95 was isolated from blackleg diseased potatoes in Pakistan in 2017. Until now, SS95 was the only P. punjabense strain described [45]. By exploring strains of different European work collections on the basis of sequences of housekeeping genes homology, we uncovered seven strains assignable to P. punjabense, and the genomes of four of them were sequenced. MLSA ( Figure 1) and core phylogenic protein analysis (Figure 2), as well as ANI and DDH genomic parameters ( Table 2), confirmed that these four new strains indeed belong to P. punjabense. The P. punjabense strains could not be distinguished from other species tested on the basis of the fatty acid profile (Table S1).
The strains RNS08-28 and P9A19a were isolated in France in 2008 and 2015, respectively, while strain IPO3715 was isolated in the Netherlands in 2013. The oldest strain in our sample, IFB5596, isolated in 1996, comes from a Polish potato field. These data provide interesting information on the prevalence of this species. First, P. punjabense is present in Europe and is not limited to Pakistan, where it was initially described. Furthermore, the Polish strain IFB5596 indicates that P. punjabense has been present in Europe for at least 25 years but has always gone unnoticed due to a lack of specific detection tools. All strains identified as P. punjabense were isolated from blackleg diseased plants, but they remain very rare among isolated strains; in the three collections analyzed, the isolation frequency of P. punjabense strains varied from 0.05 to 0.24% of the total number of SRP strains isolated during the sampling period ( Table 3).
The genomic data of the P. punjabense strains also provided valuable information on the diversity within this species. In particular, the ANI comparisons between P. punjabense strains showed values ranging from 98.6 to 100.0 (Table 2), which is comparable to those observed in the closely related species P. parmentieri, with values ranging from 98.92 to 99.97 [60]. Surprisingly, the two P. punjabense strains SS95 T and P9A19a showed exactly the same genomic sequences (ANI of 100.0; DDH of 100.0; no SNP found) although they were isolated in two distant countries and in two different years. In addition, no seed potatoes of the cultivar from which strain P9A19a was isolated are exported to Pakistan. There is no information on the cultivar from which the Pakistani strain SS95 was isolated, but there is no trade of Pakistani seed potatoes into Europe. As a result, no available data link these two strains to a possible common event. If we exclude the comparison between SS95 T and P9A19a, the average ANI for the other P. punjabense comparisons is 98.8 ± 0.1 and the mean DDH value is 89.3 ± 0.8. The relatively high genomic diversity is consistent with the finding that the pathogen has been present for a long time in Europe [61]. This hypothesis is reinforced by the high number of SNPs observed when comparing the different strains of P. punjabense identified to the type strain (more than 30,000). Therefore, the situation of P. punjabense over time is clearly distinct from that of the emergence of the clonal pathogen D. solani in European potato fields [33], where less than 100 SNPs were observed between two D. solani genomes [62].
With our current knowledge, the question of the P. punjabense reservoir remains unresolved. Indeed, the low frequency of P. punjabense found among isolates collected from blackleg symptoms, combined with the genomic diversity observed, suggests that Ppunjabense, is present in other niches beside potato plants, but no reference is available to support this hypothesis. In recent surveys of water from French rivers, no P. punjabense have been isolated [63], while P. versatile [10] or P. aquaticum [17] were recovered. P. punjabense was also not found within the French Collection for Plant-associated Bacteria (CIRM-CFBP), which includes over 265 Pectobacterium strains isolated from 1944 to 2020 and which represent isolates from a high diversity of plants [28]. Furthermore, the low incidence of P. punjabense strains sampled from blackleg-diseased plants in comparison to other species of Pectobacterium, such as P. brasiliense, P. parmentieri or P. atrosepticum, indicates that the potato plant is probably not the preferred host for this species.
The P. punjabense isolates tested were able to macerate potato tubers, indicating that they can cause soft rot during storage of tubers. For most strains investigated in this study, the ability of P. punjabense to cause potato blackleg has not been determined yet, but in a field experiment with vacuum-infiltrated seed tubers, strain IPO3715 did not result in diseased plants, whereas in the same experiment P. brasiliense showed a high disease incidence of 80% [64]. Therefore, further investigations are needed to assess the host range of P. punjabense.
Testing phenotypic characteristics using BIOLOG plates revealed that the P. punjabense strains have a lower capacity to grow in the presence of most of the antibiotic compounds tested, compared to the other pectinolytic bacteria of the panel (Figure 3). These data suggest that in complex environments, such as field soil where microbial competition is intense and phytosanitary products are widely used, P. punjabense could be less competitive [65]. This could partly explain the small number of strains collected in fields.
To date, no specific molecular tools are available to detect P. punjabense [37,38]. They showed a negative or a weak signal around the expected value at 434 bp and additional nonspecific bands with the generic PCR Y1-Y2 assay widely used to detect Pectobacterium sp. [35]. In order to clarify this, the alignment of the pectate lyase sequences of the strains of the panel revealed six mismatches between the nucleotides of the Y1 primer and the corresponding pectate lyase region in P. punjabense strains ( Figure S3); this probably explains the negative PCR result observed with the Y1-Y2 primers. This information could also partly explain why P. punjabense strains remained poorly identified, sometimes for many years. In case of negative amplification obtained with the generic Y1-Y2 primers, some isolates have possibly been assimilated to contaminants. The Taqman 2b developed in this study has shown both its specificity and its sensitivity towards the P. punjabense species and will be helpful for rapid identification, seed testing, surveys and studies on epidemiology and disease management.
Overall, this study provides a first comprehensive look at P. punjabense in Europe and provides also information on its diversity. The data have shown that it has been present for almost 25 years, but the low frequency of isolation combined with the lack of specific identification tools have prevented its recognition. The Taqman 2b assay developed in this study will enable further studies in the potato ecosystem, including the capacity to cause blackleg and the host range.