# Population genomics and antimicrobial resistance in Corynebacterium diphtheriae

Melanie Hennart, Leonardo G. Panunzi, Carla Rodrigues, Quentin Gaday, Sarah L. Baines, Marina Barros-Pinkelnig, Annick Carmi-Leroy, Melody Dazas, Anne Marie Wehenkel, Xavier Didelot, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Melanie Hennart, Leonardo G. Panunzi, Carla Rodrigues, Quentin Gaday, Sarah L. Baines, et al. Population genomics and antimicrobial resistance in Corynebacterium diphtheriae. Genome Medicine, 2020, 12 (1), pp.107. 10.1186/s13073-020-00805-7. pasteur-03199447

HAL Id: pasteur-03199447

https://pasteur.hal.science/pasteur-03199447

Submitted on 15 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Population genomics and antimicrobial
resistance in *Corynebacterium diphtheriae*

Melanie Hennart¹,², Leonardo G. Panunzi¹,³, Carla Rodrigues¹, Quentin Gaday⁴, Sarah L. Baines⁵,
Marina Barros-Pinkelnig¹, Annick Carmi-Leroy¹,⁶, Melody Dazas¹,⁶, Anne Marie Wehenkel⁴, Xavier Didelot⁷,
Julie Toubiana¹,⁶,⁸, Edgar Badell¹,⁶ and Sylvain Brisse¹,⁶*

## Abstract

**Background:** *Corynebacterium diphtheriae*, the agent of diphtheria, is a genetically diverse bacterial species. Although antimicrobial resistance has emerged against several drugs including first-line penicillin, the genomic determinants and population dynamics of resistance are largely unknown for this neglected human pathogen.

**Methods:** Here, we analyzed the associations of antimicrobial susceptibility phenotypes, diphtheria toxin production, and genomic features in *C. diphtheriae*. We used 247 strains collected over several decades in multiple world regions, including the 163 clinical isolates collected prospectively from 2008 to 2017 in France mainland and overseas territories.

**Results:** Phylogenetic analysis revealed multiple deep-branching sublineages, grouped into a Mitis lineage strongly associated with diphtheria toxin production and a largely toxin gene-negative Gravis lineage with few toxin-producing isolates including the 1990s ex-Soviet Union outbreak strain. The distribution of susceptibility phenotypes allowed proposing ecological cutoffs for most of the 19 agents tested, thereby defining acquired antimicrobial resistance. Penicillin resistance was found in 17.2% of prospective isolates. Seventeen (10.4%) prospective isolates were multidrug-resistant (≥ 3 antimicrobial categories), including four isolates resistant to penicillin and macrolides. Homologous recombination was frequent (r/m = 5), and horizontal gene transfer contributed to the emergence of antimicrobial resistance in multiple sublineages. Genome-wide association mapping uncovered genetic factors of resistance, including an accessory penicillin-binding protein (PBP2m) located in diverse genomic contexts. Gene *pbp2m* is widespread in other *Corynebacterium* species, and its expression in *C. glutamicum* demonstrated its effect against several beta-lactams. A novel 73-kb *C. diphtheriae* multiresistance plasmid was discovered.

**Conclusions:** This work uncovers the dynamics of antimicrobial resistance in *C. diphtheriae* in the context of phylogenetic structure, biovar, and diphtheria toxin production and provides a blueprint to analyze re-emerging diphtheria.

**Keywords:** *C. diphtheriae*, Antibiotic resistance, Genome-wide association study, Phylogeny, Mobile genetic element, Biovar

<footnote>* Correspondence: sbrisse@pasteur.fr
¹Institut Pasteur, Biodiversity and Epidemiology of Bacterial Pathogens, Paris, France
⁶Institut Pasteur, National Reference Center for Corynebacteria of the Diphtheriae Complex, Paris, France
Full list of author information is available at the end of the article</footnote>
## Background

Diphtheria, if untreated, is one of the most severe bacterial infections of humans. It typically affects the upper respiratory tract causing pseudomembrane formation, sometimes leading to suffocation and death. The infection can be complicated by systemic symptoms, caused by the diphtheria toxin. Other forms of disease are skin and invasive infections, including endocarditis [1, 2].

The agent of diphtheria is *Corynebacterium diphtheriae*, a member of the phylum Actinomycetes [3, 4]. The diphtheria toxin, encoded by the *tox* gene, is carried by lysogenized corynephages within the chromosome of some *C. diphtheriae* strains [5, 6]. Concern exists about the possibility of lysogenic conversion of previously non-toxigenic strains during colonization, infection, or transmission chains [7]. However, knowledge on the microevolutionary dynamics between *tox*-positive and *tox*-negative strains is limited. The high genetic diversity of *C. diphtheriae* strains underlies their variable colonization, adhesion, and pathogenicity properties [8–10]. Although three main biovars (Mitis, Gravis, and Belfanti) are distinguished since the 1950s, their phylogenetic relationships are poorly defined [11–13].

Diphtheria used to be one of the deadliest infections in young children but has been largely controlled by vaccination with the highly effective toxoid vaccine [14]. Even so, thousands of cases of diphtheria are still reported annually [15], and large outbreaks can quickly follow the disruption of public health systems [14, 16–18]. In countries with high vaccination coverage, diphtheria cases are associated with travel and migration from endemic regions [19–21]. As diphtheria vaccination is performed using an inactivated form of diphtheria toxin, it is not considered to prevent asymptomatic colonization and silent transmission of the pathogen, which still circulates and is the object of intense epidemiological surveillance [4]. However, vaccine preparations may include other antigens, and the impact of vaccination on *C. diphtheriae* evolution deserves further studies [22].

Clinical management of infections with toxigenic isolates includes treatment with diphtheria antitoxin (DAT), which can prevent or reduce the systemic effects of the toxin [4]. Nevertheless, antimicrobial treatment is critical in the clinical management of both *tox*-positive and *tox*-negative infections, as it contributes to the elimination of the bacteria within the patient and limits transmission to novel individuals [23]. With DAT production being threatened [24], antimicrobial treatment might become even more critical in diphtheria therapy.

Penicillin is the first-line therapeutics to treat diphtheria, with erythromycin being recommended in case of allergy [25]. Both antimicrobial agents are generally effective for the treatment of diphtheria [23, 26]. However, reduced susceptibility or full resistance of *C. diphtheriae* to

penicillin, erythromycin, and other antimicrobial agents, sometimes combined in multidrug-resistant strains, has been reported from multiple world regions [23, 26–35].

Antimicrobial resistance genes have been described in *C. diphtheriae*, including the erythromycin resistance gene *ermX* on plasmid pNG2 [36] and genes *dfrA16*, *qacH*, and *sul1* carried on a class 1 integron, mobilized by IS6100 [37]. However, the prevalence and phylogenetic distribution of resistance genes in *C. diphtheriae* clinical isolates are unknown. Six chromosomal penicillin-binding proteins (PBPs) have been reported in *C. diphtheriae* [38], but so far, no association between *pbp* or other genetic variation and penicillin resistance has been described. Understanding the genetic basis of antimicrobial resistance in *C. diphtheriae* would improve our ability to diagnose and track its spread.

The aims of this study were (i) to characterize antimicrobial resistance phenotypes in a large collection of *C. diphtheriae* strains with diverse geographical and temporal origins and to uncover genomic determinants of resistance and (ii) to analyze the population structure of *C. diphtheriae* and define the associations between antimicrobial resistance, diphtheria toxin production, biovars, and phylogenetic sublineages.

## Methods

### Corynebacterium diphtheriae isolates and strains

A collection of 247 *C. diphtheriae* isolates were included (Additional file 1: Table S1), corresponding to three subsets. First, we included 163 clinical isolates (recent clinical isolates subset, Additional file 1: Table S1) collected prospectively between 2008 and 2017 by the French National Reference Center for Corynebacteria of the *Corynebacterium diphtheriae* complex (NRC-CCD). These isolates represented all isolates received at the NRC-CCD that corresponded to the *C. diphtheriae* species (*C. belfantii*, *C. ulcerans*, *C. pseudotuberculosis*, or other corynebacteria were excluded). They were collected from cutaneous (*n* = 136), respiratory (*n* = 23), and other (bones, blood; *n* = 4) infections. Of these, 74 were from Mainland France and 89 from French overseas territories, including Mayotte (*n* = 50), New Caledonia (*n* = 19), La Réunion Island (*n* = 11), French Guiana (*n* = 4), French Polynesia (*n* = 3), and Guadeloupe (*n* = 1); one isolate received from Institut Pasteur in Madagascar was also included (Fig. 1). Four isolates from New Caledonia collected between 2002 and 2006 (02-0322, 02-0338, 03-1641, and 06-1569) were included in a previous study [39]; the trimethoprim- and sulfamethoxazole-resistant isolate FRC0024 was previously shown to harbor an integron with gene *drfA16* [37].

Second, we included 15 clinical isolates collected in France between 1981 and 1991, 11 of which had been deposited in the Collection de l'Institut Pasteur (CIP;
<figure><img src="image_4.png" /></figure>

A. Strains per year and resistance phenotypes origins

B. Geographic origins

<caption>**Fig. 1** Temporal and geographical distribution of strains studied. a The number of strains per year, 1945–2017. The 2008–2017 clinical isolates are represented in red- or yellow-shaded patterns (see key), whereas the older clinical isolates are in gray and the reference strains in blue. The inset shows the pie charts with the frequency of resistance phenotypes among strains from the four most represented geographic origins; the remaining ones are pooled in the right-most pie chart. b Geographic origins of strains from the three subsets</caption>
historical clinical isolates subset in Additional file 1:
Table S1).

Third, to increase the genetic diversity and geographic range of the sample, the 65 available reference strains of ribotypes that belong to *C. diphtheriae* were included [40]. These reference strains represent an international collection of isolates collected over several decades and originating from multiple world regions including the Americas, Europe, Asia, Africa, and Oceania. Our sub-cultures of these strains were controlled for *tox* gene presence, toxin production, and biovar, leading to modifications of published characteristics in some instances (Additional file 1: Table S1; Fig. 1). Finally, four reference strains were included: strain NCTC13129, which is used as genomic sequence reference [38]; strain NCTC10648, which is used as the *tox*-positive and toxigenic reference strain in PCR and Elek tests, respectively; strain NCTC11397<sup>T</sup>, which is the taxonomic type strain of the *C. diphtheriae* species; and the vaccine production strain PW8, which corresponds to CIP A102 [41]. This third subset is referred to as “ribotype and reference strains” subset (Additional file 1: Table S1).

## Bacterial cultures, identification, and biovar

Bacteria were cultivated on Trypto-Casein-Soy (TCS) agar for 24 h at 35–37 °C. Bacterial identification was performed at the NRC-CCD as described previously [42] by multiplex polymerase chain reaction (PCR) combining a *dtxR* gene fragment specific for *C. diphtheriae* and a multiplex PCR that targets a fragment of the *pld* gene specific for *C. pseudotuberculosis*, the gene *rpoB* (amplified in all species of the *C. diphtheriae* complex) and a fragment of 16S rRNA gene specific for *C. pseudotuberculosis* and *C. ulcerans*, respectively. Isolates collected since 2014 were confirmed as *C. diphtheriae* by matrix-assisted laser desorption-ionization time-of-flight mass spectrometry (MALDI-TOF MS) using Bruker technology. In order to exclude strains initially identified as *C. diphtheriae* but now classified as *C. belfantii* [42] or *C. rouxii* [43], genome-wide average nucleotide identity (ANI) was used as described previously [42]. Strains were characterized biochemically for pyrazinamidase, urease, and nitrate reductase and for utilization of maltose and trehalose using API Coryne strips (BioMérieux, Marcy l'Etoile, France) and the Rosco Diagnostica reagents (Eurobio, Les Ulis, France). The Hiss serum water test was used for glycogen fermentation. The biovar of isolates was determined based on the combination of nitrate reductase (positive in Mitis and Gravis, negative in Belfanti) and glycogen fermentation (positive in Gravis only). The rare biovar Intermedius was not identified, as its distinction from other biovars is based on colony morphology, which is considered subjective, or on lipophily, which was not tested.

## Determination of the presence of the *tox* gene

Determination of the diphtheria toxin gene (*tox* gene) presence was achieved by a conventional *tox* PCR assay [44], while its phenotypic production was assessed by the modified Elek test [45]. We also confirmed *tox* PCR results by BLASTN (query: *tox* gene sequence from strain NCTC13129, RefSeq accession number: DIP_RS12515) analysis of the genomic assemblies.

## Antimicrobial susceptibility testing

Phenotypic susceptibility was tested for the following agents: penicillin G (10 IU), amoxicillin, oxacillin, cefotaxime, imipenem, erythromycin, azithromycin, clarithromycin, spiramycin, pristinamycin, kanamycin, gentamicin, rifampicin, tetracycline, ciprofloxacin, clindamycin, sulfonamide, trimethoprim, and trimethoprim + sulfamethoxazole. The 19 antimicrobial agents tested (Table S1) corresponded to seven classes, as described hereafter: β-lactams: penicillin G (PEN), amoxicillin (AMX), oxacillin (OXA), cefotaxime (CFT), imipenem (IMP); macrolides: azithromycin (AZM), clarithromycin (CLR), erythromycin (ERT), and spiramycin (SPR); lincosamides: clindamycin (CLD); streptogramins: pristinamycin (PRT); aminoglycosides: gentamicin (GEN) and kanamycin (KAN); folate pathway inhibitors: sulfonamide (SUL), trimethoprim (TMP), and trimethoprim + sulfamethoxazole (cotrimoxazole, TMP-STX); ansamycins: rifampicin (RIF); tetracyclines: tetracycline (TET); and fluoroquinolones: ciprofloxacin (CIP).

Antimicrobial susceptibility was determined using the disk diffusion method with impregnated paper disks (Bio-Rad, Marnes-la-Coquette, France) on Mueller-Hinton agar plates supplemented with 5% of sheep blood and 20 mg/L β-NAD, as recommended. Minimum inhibitory concentrations (MICs) were determined using E-test strips (BioMerieux, Marcy l'Etoile, France). The control strain used is *Streptococcus pneumoniae* ATCC 49619. The zone diameter (ZD) data were interpreted into S, I, and R categories in the following way. First, we used the CA-SFM/EUCAST V.1.0 (January 2019) document (https://www.sfm-microbiologie.org/wp-content/uploads/2019/02/CASFM2019_V1.0.pdf), which contains interpretative criteria for *Corynebacterium* spp. only for CIP, GEN, CLD, TET, RIF, and TMP-STX. Second, for the other agents, we used the interpretative criteria published in Table III of the CA-SFM 2013 recommendations (https://resapath.anses.fr/resapath_uploadfiles/files/Documents/2013_CASFM.pdf). Note that for RIF, we used the 2013 breakpoints, as they fitted better with the observed distribution of ZD values. Clarithromycin breakpoints were taken from those for erythromycin, as recommended. The breakpoint for oxacillin was derived from the one used for *Staphylococcus* spp. ZD interpretation breakpoints are given in Additional file 1: Table S2.
Penicillin susceptibility was initially determined using 10-UI (6 µg) disks (resistance breakpoint, 18 mm), but CA-SFM/EUCAST recommendations were changed in 2014 to use 1-UI disks, while the resistance breakpoint was increased from 18 to 29 mm. As all *C. diphtheriae* strains end up in the resistant category following this recommendation, E-test strips were used to define the penicillin MIC since 2014 (Additional file 1: Table S2: isolates starting from FRC0259); the EUCAST breakpoint of 0.125 g/L was used as a cutoff. Penicillin E-test was also performed systematically for strains tested as resistant before 2014 (using 10-UI disks), as well as for some susceptible isolates (Additional file 1: Table S1).

Multidrug-resistant *C. diphtheriae* (MDR-DIP) were defined as strains resistant to three or more antimicrobial agent categories (defined in Additional file 1: Table S2), excluding intrinsic resistance to fosfomycin [46]. Note that we used ecological cutoffs rather than currently proposed clinical breakpoints (see Additional file 1: Table S2 for a comparison of both types of breakpoints).

## Whole-genome sequencing by Illumina and Oxford Nanopore Technologies

DNA was extracted from broth cultures, by making use of DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany). However, a lysis step was added to the extraction protocol described by the manufacturer as previously described [47]: a 1-µL loopful of bacterial colonies was emulsified in 180 µL of lysis buffer containing 20 mM Tris-HCl, pH 8, 2 mM EDTA, 1.2% Triton X-100, and 20 mg/mL lysozyme, in a DNase/RNase-free 1.5-mL Eppendorf tube and incubated in a heating block at 37 °C for 1 h, with mixing every 20 min. After extraction, DNA concentration was measured with the Qubit 3.0 Fluorometer (Invitrogen), employing the Qubit dsDNA BR Assay Kit (Invitrogen). Besides, the DNA quality was verified using a D-One spectrophotometer (Nanodrop). Multiplexed paired-end libraries (2 × 150 bp) were prepared using the Nextera XT DNA kit (Illumina, San Diego, CA, USA) and eventually sequenced with an Illumina NextSeq-500 instrument at a minimum of 50× coverage depth. Trimming and clipping were performed using AlienTrimmer v0.4.0 [48]. Redundant or over-represented reads were reduced using the khmer software package v1.3 [49]. Finally, sequencing errors were corrected using Musket v1.1 [50]. A de novo assembly was performed for each strain using SPAdes v3.12.0 [51]. The genomic sequences of the four reference strains were retrieved from public repositories (Additional file 1: Table S1).

Additionally, the multidrug-resistant isolate FRC0402 was subjected to long-read sequencing using Oxford Nanopore Technologies (ONT). Genomic DNA was extracted using the phenol-chloroform protocol combined

with Phase Lock Gel tubes (Qiagen GmbH). Libraries were prepared using a 1D ligation sequencing kit (SQK-LSK-108) without fragmentation and sequenced using a MinION FLO-MIN-106 flow cell. Finally, ONT and Illumina short reads were combined to generate a hybrid assembly using Unicycler v0.4.4 (normal assembly mode, default parameters).

## Phylogeny, recombination, and genomic sequence analyses

We built a core genome multiple sequence alignment (cg-MSA) from the assembled genome sequences. For this, the genome sequences were annotated using PROKKA v1.14.2 [52] with default parameters, resulting in GFF files. Roary v3.6 [53] was used to define protein-coding gene clusters, with a threshold set at 70% amino acid identity. Core genes were defined as being present in 95% of genomes and were concatenated into a cg-MSA by Roary. ClonalFrameML v1.11 [54] was used to build a phylogenetic tree based on the cg-MSA, which quantifies and accounts for the effects of recombination events. PhyML v20131022 [55] was used to build an initial tree. To evaluate bootstrap support, IQtree version 2 [56] with best-fit model GTR+F+R10 was used.

We used Kleborate v1.0.0-beta (https://github.com/katholt/Kleborate), with the --resistance option, to identify (identity > 80% and coverage > 90%) known resistance genes in *C. diphtheriae* genomic sequences, based on the August 1, 2019, update of the ARG-Annot database. We used BLASTN (identity > 80% and coverage > 95%) to search for the presence of the *tox* gene and of genes associated with biovar Gravis (DIP351, DIP354, and DIP357) and for nitrate utilization (*narKGHJI*).

MLST genotypes were defined using the international MLST scheme for *C. diphtheriae* and *C. ulcerans* [12].

## Genome-wide association studies

The software treeWAS [57] was used to find genome-wide associations between either antimicrobial resistance phenotypes or biovar on the one hand and genetic variants (both core-genome SNPs and accessory genome gene presence/absence) on the other hand. Core-genome SNPs were derived either from a mapping approach (Samtools v1.9 and GATK v3.4-0), which comprises intergenic regions, or from the alignment of core coding sequences found using Roary. We ran treeWAS v1.1 with default parameters, using as input the previously computed ClonalFrameML phylogenetic tree and distribution of homoplasies, in order to account for both the population structure and effect of recombination. For this analysis, susceptibility phenotypes were classified into resistant or susceptible categories based on zone diameter phenotypes using the CA-SFM/EUCAST 2019 cutoffs (Fig. 3, Additional file 3: Table S3). The
seven chromosomal PBP coding sequences (including gene with locus tag RS14485 in RefSeq NC_002935.2, or DIP0637 in the original GenBank file, renamed by us as PBP4b) were extracted from the genomic sequences and translated into amino acid (AA) sequences, which were also analyzed for association with penicillin resistance.

## Cloning and transformation experiments

For ectopic expression in *C. glutamicum*, the *pbp2m* gene was amplified from *C. diphtheriae* strain FRC0402 and put under the control of the inducible *PgntK* promoter on the shuttle vector pTGR5 [58] (Additional file 4: Fig. S1). *pbp2m* was assembled in this plasmid by Gibson assembly using the primers PBPdi_Fw (CAA AGA AAG GAT AAG ACC ATA TGA TGA CTA AGC ACA ATC GTT TCC GTC), PBPdi_Rv (TAC CTT AAG CGG CCG CTT TAT TGA ATT CCA GAG AAT TTC TGA ACA TCC G), pTGRdi_Fw (TAA AGC GGC CGC TTA AGG TAC C), and pTGRdi_Rv (ATG GTC TTA TCC TTT CTT TGG TGG CG).

*Escherichia coli* CopyCutter EPI400 (Lucigen) was used for the cloning of the *pbp2m* gene and was grown in Luria-Bertani (LB) broth or agar plates at 37 °C supplemented with 50 µg/mL kanamycin. The pTGR5_pbp2m plasmid was sequenced and electroporated into *Corynebacterium glutamicum* ATCC 13032. Positive colonies were grown in brain heart infusion (BHI) at 30 °C and 120 rpm supplemented with 25 µg/mL kanamycin and 1% (w/v) gluconate when required for ectopic expression of Pbp2m.

## Mapping of SNPs on *C. diphtheriae* PBP sequences

Functional annotation of the different sequences was performed with InterPro [59]. Conserved transpeptidation motifs SxxK, SxN, and KTG were identified and mapped on the PBP sequences from *Corynebacteriales* based on the results of multiple sequence alignments performed with Clustal Omega [60]. When uncertainty between a transmembrane domain and a signal peptide existed, a decision was made based on the previous characterization of the homologous PBP in other *Corynebacteriales* in the literature.

## Results

### Provenance and microbiological characteristics of *C. diphtheriae* isolates

We studied 247 *C. diphtheriae* strains of diverse geographic and temporal origins (Fig. 1). This collection included 163 isolates prospectively collected between 2008 and 2017 from French mainland and overseas territories, 15 older (1981–2006) French clinical isolates, 65 ribotype reference strains [40], and 4 other reference strains. All isolates were confirmed as *C. diphtheriae* (excluding *C. belfantii* and *C. rouxii*) based on an average

nucleotide identity (ANI) value higher than 96% with the *C. diphtheriae* type strain NCTC11397T.

Approximately one third (n = 78, 32%) of isolates were tox-positive (as defined by the detection of the *tox* gene by PCR), whereas the remaining 169 isolates (68%) were tox-negative. The proportions of *tox*-positive isolates were 42%, 34%, and 2% among reference strains, 2008–2017 clinical isolates, and older clinical isolates, respectively (Additional file 4: Fig. S2). Of the 78 *tox*-positive isolates, 17 (21.8%) were negative for toxin production and thus correspond to non-toxigenic toxin gene-bearing (NTTB) isolates. Six of the NTTB isolates had a stop codon within the *tox* gene sequence (Additional file 1: Table S1, Additional file 5: Table S4). However, for the 11 remaining strains, we found no explanation for the observed lack of toxin production.

Upon biotyping, 154 (62.3%) isolates belonged to biovar Mitis, 87 to biovar Gravis (35.2%), and 6 (2.4%) to biovar Belfanti (Additional file 1: Table S1). Biovar proportions were similar among the three datasets. Mitis isolates were more frequently tox-positive than Gravis isolates (56/154 versus 18/87, chi-squared test, *p* value 0.01; Additional file 4: Fig. S2). Among tox-positive isolates, NTTB were more frequent among Mitis isolates (13/56, 23.2%) than among Gravis isolates (1/18, 5.6%) although this difference was not statistically significant (*p* value 0.09). Three out of four tox-positive Belfanti isolates were NTTB.

## Phylogenetic structure of *C. diphtheriae* and distribution of the toxin gene

To infer a phylogenetic tree, we first aimed to detect and remove homologous recombination events among *C. diphtheriae* genomic sequences. ClonalFrameML inferred a relative rate of recombination to mutation (R/theta) of 0.86, with an average length of recombination segments (delta) of 287 bp. The mean genetic distance between donor and recipient of recombination (nu) was estimated by ClonalFrameML to be 0.02, resulting in a relative impact of recombination to mutation in genomic diversification (r/m = R/theta × delta × nu) of 5.01 [61].

The recombination-corrected phylogeny was similar to the uncorrected phylogeny but with shorter branches, as expected (Additional file 4: Fig. S3). The phylogeny, rooted with *C. belfantii* and *C. rouxii* (Additional file 4: Fig. S4), was star-like, with a multitude of sublineages branching off deeply. The deepest branching sublineages corresponded to two ribotype reference strains of biovar Mitis: CIP107521 (ribotype Dagestan) and CIP107534 (ribotype Kaliningrad). Remarkably, isolates of biovars Mitis and Gravis were mostly distributed in two distinct branches of the tree. We therefore named the two major branches, lineage Mitis (156 strains, of which 86% were of biovar Mitis) and lineage Gravis (91 strains, of which
77% were of biovar Gravis). The Gravis lineage branched off from within the Mitis lineage (Fig. 2; Additional file 4: Fig. S5). However, although most shallow branches were strongly supported, bootstrap support values of the deep nodes of the phylogeny were generally weak, implying that the position of the Gravis branch within the Mitis diversity was uncertain, as were the relative positions of the early-branching Gravis and Mitis sublineages (Additional file 4: Fig. S5, Fig. S6). Reference strains PW8 and NCTC11297T belonged to the Mitis lineage, whereas NCTC13129 (from the ex-Soviet Union 1990's outbreak) and NCTC10648 belonged to the Gravis lineage. The Belfanti isolates were scattered in three distinct sublineages within the Mitis lineage and one within the Gravis lineage.

The isolates carrying the *tox* gene belonged mostly to the Mitis lineage (68 of 78, 87.2%), in which they were distributed in multiple sublineages. In the Mitis lineage, 69 (44.2%) were *tox*-positive. In contrast, within the Gravis lineage, only 10 (11%) isolates were *tox*-positive, and they corresponded to the earliest-branching Gravis sublineages with only one exception. Interestingly, this exception corresponded to the large ex-Soviet Union outbreak in the 1990s (Fig. 2; Additional file 4: Fig. S5). This phylogenetic pattern is consistent with an evolutionary scenario where Mitis is the ancestral biovar of *C. diphtheriae* and where Gravis evolved from the Mitis lineage as an initially *tox*-positive sublineage, with subsequent loss of the toxin gene. In this scenario, the ex-Soviet Union outbreak sublineage would have re-acquired the *tox* gene. All NTTB isolates belonged to the Mitis lineage except strain CIPA99 (ribotype Rhone, biovar Belfanti; Fig. 2), and they were distributed in multiple sublineages, showing convergent evolution towards the loss of toxin production.

## Genetic events linked to biovar status

Biovar Mitis and Gravis are distinguished by the ability to utilize glycogen (positive in Gravis, negative in Mitis). The *spuA* gene, which codes for a putative alpha-1,6-glycosidase, was reported as being specific for biovar Gravis isolates [62]. Our genome-wide association study (GWAS) of accessory genes with the biovar phenotype revealed a strong association of a cluster of genes that includes *spuA* (DIP357; Additional file 4: Fig. S7) with biovar Gravis isolates, providing statistical support to the discovery of Santos et al. [62]. This association was stronger within the Gravis lineage; in contrast within the Mitis lineage, few (5 out of 17) of the biovar Gravis isolates possessed *spuA* (Fig. 2). GWAS analysis of core SNPs further demonstrated that a SNP (at position 324,487, Additional file 4: Fig. S7) downstream of the *spuA* cluster insertion point was also associated with biovar, suggesting homologous

recombination among core genes as a mechanism for the *spuA* cluster insertion event.

The nitrate reductase activity differentiates Mitis and Gravis isolates, which are positive, from Belfanti isolates, which are nitrate-negative. We found that the nitrate reduction *narKGHJI* gene cluster [62] was disrupted in three of the six isolates assigned to the biovar Belfanti: strains FRC0480 and FRC0481 had a G to A mutation at position 675 of the *narG* gene, leading to a stop codon, whereas in strain CIPA99, approximately 100 nucleotides were inserted at position 446 in *narG*. No molecular explanation was found for the lack of nitrate reductase ability of the three other Belfanti strains when scrutinizing the *narKGHJI* gene cluster and adjacent mobility cofactor biosynthesis genes [63].

## Antimicrobial susceptibility variation

Susceptibility to 19 antimicrobial agents was determined for the 247 clinical isolates and reference strains (Additional file 3: Table S3). For each agent, the distribution of zone diameter (ZD) values (Fig. 3) revealed a predominant mode located towards the right end of the distribution. This mode likely corresponds to the natural susceptibility distribution within the *C. diphtheriae* population and was used to define tentative epidemiological cutoffs (ECOFFs, also called ecological cutoff [64]). The proposed ECOFFs and their comparison with clinical breakpoints are presented in Additional file 2: Table S2. For each antimicrobial agent except cefotaxime, this approach led to the identification of strains that had ZD values smaller than those in the main mode, thus potentially corresponding to acquired resistance (Fig. 3).

Penicillin was exceptional in that the predominant mode (centered around 36 mm) was less neatly defined, due to the partial overlap with a second mode of smaller diameter values centered around 24 mm. This second mode corresponds mostly to the “intermediate” interpretative category ($18 \le ZD < 29$ mm) but also overlaps with the “resistant” category ($< 18$ mm). The distribution of ZD values for tetracycline also showed a clear second mode. For multiple other agents (amoxicillin, oxacillin, imipenem, kanamycin, rifampicin, ciprofloxacin, clindamycin and more evidently sulfonamide, trimethoprim, and the trimethoprim-sulfamethoxazole combination), some strains had the minimal diameter (6 mm, corresponding to growth at the disk contact). For trimethoprim, we observed both a second mode centered around 14 mm and a group of even more resistant strains with growth at disk contact.

Antimicrobial resistance levels were similarly distributed between *tox*-positive and *tox*-negative isolates (Additional file 4: Fig. S8) as well as between the two main phylogenetic lineages or biovars (Additional file 4: Fig. S9; Additional file 1: Table S1). In particular, the
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Phylogenetic tree of *C. diphtheriae*. The tree was obtained using ClonalFrameML and was rooted using *C. rouxii* and *C. belfantii* isolates (see supplementary material). The main lineages Mitis and Gravis are labeled, and their branches are drawn using purple and green, respectively. The first (internal) circle around the tree corresponds to the three strain subsets (red, recent clinical isolates; blue, reference strains; gray, older clinical isolates). The second circle (stars) gives the toxigenic status. The third circle corresponds to biovars Mitis (purple), Gravis (green), and Belfanti (yellow). The next three circles indicate the presence of the spuA-associated gene cluster; DIP357 = spuA gene. The positions of reference strains PW8, NCTC13129, and NCTC10648 are indicated. The scale bar gives the number of nucleotide substitutions per site</figcaption></figure>
<figure><img src="image_6.png" /></figure>

proportion of multidrug-resistant strains was not statisti-
cally different between biovars Mitis (n = 14) and Gravis
(n = 4; chi-square test *p* value 0.2) or between tox-posi-
tive and tox-negative isolates (6 versus 12, *p* value 0.87).

Resistance rates were 17.2%, 2.5%, and 2.5% for peni-
cillin, amoxicillin, and erythromycin, respectively, among
the prospectively collected 2008–2017 clinical isolates
(Fig. 4). Reference strains were generally susceptible to
most agents, including penicillin, but were partially re-
sistant to tetracycline (18%) and sulfonamide (35%). The
resistance profile distribution showed that approximately
half (121/247) of the strains had a fully susceptible
phenotype, whereas 18 isolates (7.3%) were multidrug-
resistant (Table S1). Notably, four isolates were resistant
at the same time to penicillin and macrolides (Fig. 4
inset), and two of them (FRC0402 and FRC0466) add-
itionally had a reduced susceptibility to amoxicillin. Two
of these isolates (FRC0475 and FRC0478) were collected
from a foot arch wound and in respiratory carriage in
the same patient (French mainland, with recent travel
from New Caledonia). The two others came from a

patient living in La Réunion Island (FRC0402) and from
a patient living in Paris, who had recently traveled to
Tunisia (FRC0466). Of the 18 multidrug-resistant iso-
lates, 1 was a 1994 isolate from India and 17 belonged to
the 2008–1017 prospective subset of clinical isolates
(representing 10.4% of these). Isolates from La Réunion
Island and mainland France showed resistance to mul-
tiple antimicrobial agents more often than isolates from
other geographic origins (Fig. 1a, inset).

Genomic associations with antimicrobial resistance phenotypes

We first searched for the presence in the genomic se-
quences, of previously described antimicrobial resistance
genes (ARGs). This approach led to the detection of 12
ARGs (Additional file 6: Table S5). We identified three
tetracycline resistance genes (tetW, tet33, and tetO), four
aminoglycoside resistance genes [aph (3')-Ia, aph (3'')-Ib,
aph (6)-Id, and aadA1 = ant (3')-Ia], and also ermX,
dfrA16, dfrA15b, dfrA1, and sul1 genes. We observed a
strong correlation between the presence of ARGs and the
<figure><img src="image_7.png" /></figure>

expected resistance phenotypes (Additional file 4: Fig. S10), particularly for ermX (macrolide resistance; phi coefficient = 0.66 with erythromycin resistance), sul1 (sulfonamide resistance; phi = 0.80), and aph (3')-Ia (kanamycin resistance; phi = 0.84). In strains FRC0137 and FRC0375, this latter gene was linked to aph (3')-Ib (strA), aph (6)-Id (strB), and ermX on a Tn5432-like genomic region with an IS1249 insertion sequence [65]. The phylogenetic distribution of ARGs (Fig. 5; Additional file 4: Fig. S6) revealed their presence in multiple unrelated sublineages, consistent with independent acquisitions by horizontal gene transfer. Gene ermX was present either in proximity to gene pbp2m (see below) or in a fragmented insertion-sequence rich accessory region. Gene dfrA16 was associated with sul1 on a reported [37] class 1 integron (see the “Discovery of a multidrug-resistant conjugative plasmid carrying the gene pbp2m” section). Tetracycline resistance was associated either with the ribosomal protection protein genes tet(O) or tet(W) or with the efflux pump gene tet33. These three genes were present in distinct strain subsets and appear to contribute independently to tetracycline resistance in C. diphtheriae (Fig. 5); they were mostly associated with insertion sequences but not with other ARGs.

Next, in order to identify novel genetic determinants potentially associated with antimicrobial resistance in C. diphtheriae, a GWAS approach was followed, based on either core genome SNPs or accessory gene presence/absence. SNPs that were strongly associated with ciprofloxacin, trimethoprim, and rifampicin resistance were identified within genes for gyrase subunit A, dihydrofolate reductase, and RNA polymerase subunit B, respectively (Additional file 6: Table S5; Additional file 4: Fig. S11), consistent with known mechanisms and validating our approach. SNPs were also found to be associated with penicillin, kanamycin, and tetracycline resistance (Additional file 6: Table S5), but based on the annotations of the genes where they are located, the functional implications of these SNPs are unclear in these cases. No association was found for penicillin resistance within the core PBP genes using the genome-wide approach. However, using a concatenation of the amino acid sequences of the seven identified putative PBPs of C. diphtheriae, we identified amino acid positions that were statistically associated with penicillin resistance (Additional file 4: Fig. S12). The identified positions were mapped onto the predicted functional domains of the different PBPs (Additional file 4: Fig. S13) revealing several mutation hotspots. A number of significant SNPs were
<figure><img src="image_8.png" /><figcaption>**Fig. 5** (See legend on next page.)</figcaption></figure>
(See figure on previous page.)

**Fig. 5** Phylogenetic distribution of antimicrobial resistance phenotypes and genes. The phylogenetic tree is the same as in Fig. 1; the Mitis branch is in purple, the Gravis branch in green. To the right of the tree, each bloc indicates first the phenotype (red, resistant; see key) and relevant corresponding genotypes (orange, gene or mutation presence)

found within the transpeptidase (TP) domains of the differ-
ent PBPs, but none of the mutations mapped to the con-
served transpeptidation motif. Other mutations were found
outside the TP domains, for instance, in the transglycosy-
lase and PASTA domains of PBP1b or in the dimerization
domain of PBP2b.

GWAS analysis of accessory genes demonstrated signifi-
cant associations with phenotypic resistance. Associated
genes included those mentioned above for erythromycin,
tetracycline, kanamycin, sulfonamide, and trimethoprim
(Additional file 6: Table S5), showing that these genes are
the main mechanisms of resistance to these antimicrobial
agents in *C. diphtheriae* populations. In addition, an
accessory penicillin-binding protein gene (which we name
*pbp2m*) was strongly associated with penicillin resistance.
This gene is described in more detail below.

Discovery of a penicillin-binding protein (PBP2m)
associated with penicillin resistance

The PBP gene *pbp2m* was observed in 11 isolates, 8 of
which were penicillin resistant with minimum inhibitory
concentrations (MIC) ranging from 0.19 to 1.5 mg/L.
Two of these isolates were also resistant to amoxicillin
and one was in addition resistant to oxacillin

(Additional file 2: Table S2). The phylogenetic distribu-
tion of *pbp2m*-positive strains was compatible with mul-
tiple independent acquisitions of the gene through
horizontal gene transfer (Fig. 5). Sequence analysis
showed that the newly identified PBP2m is almost iden-
tical (3 differences out of 593 amino acids, 99.5%) to
PBP2c from *C. jeikeium*, a class B PBP with an N-
terminal signal peptide followed by a lipobox domain
and the C-terminal transpeptidase domain (Add-
itional file 4: Fig. S13). The *C. jeikeium* PBP2c is a low-
affinity PBP and was associated with beta-lactam resist-
ance in *C. jeikeium* [66].

To demonstrate the role of PBP2m in penicillin resist-
ance, its gene was PCR amplified from FRC0402 and
cloned into the pTGR5 plasmid (Additional file 4: Fig.
S1). Transformation of the plasmid into C. glutamicum
strain ATCC 13032 raised the MIC for penicillin from
0.125 to 1.5 mg/L, and the MICs of the other beta-
lactams amoxicillin, oxacillin, and cefotaxime also in-
creased importantly (Fig. 6; Additional file 7: Table S6).
In contrast, MICs of non-beta-lactam agents were not
changed. Imipenem was less effective against the trans-
formant based on disk diffusion but not based on E-test.
Transformation with the empty plasmid used as control

<figure><img src="image_9.png" /></figure>

**Fig. 6** Phenotypic effect of *pbp2m* expression. Compared susceptibility phenotypes for *C. glutamicum* transformants with plasmid pTGR5 containing, or not, the *pbp2m* gene. Left, the shift in zone diameter size; right, the shift in the minimum inhibitory concentration (MIC). Diamonds are positioned on the scales, at positions corresponding to the difference of zone diameters (without *pbp2m* – with *pbp2m*) or the log2 of MIC ratios (with *pbp2m*/without *pbp2m*). Red, penicillins or cephalosporins; blue, other agents. Tmp-Stx, trimethoprim-sulfamethoxazole
did not affect the MIC of any agent. These results show
that PBP2m confers resistance to a broad range of beta-
lactams.

Discovery of a multidrug-resistant conjugative plasmid
carrying the gene pbp2m

Strain FRC0402, a *tox*-negative isolate from La Réunion
Island, stood out as being resistant to 12 agents (Fig. 4
inset). In addition to *pbp2m*, this isolate carried genes
*sul1*, *ermX*, and *dfrA16* and a *tetA* family *tet(Z)*-like
(71%) tetracycline efflux gene. To define the genomic

context of resistance genes, a complete genome se-
quence was obtained. The assembly revealed a chromo-
some of 2,397,465 bp and a circular plasmid of 73,763 bp
(Fig. 7), which we propose to name pLRPD (for large re-
sistance plasmid of *C. diphtheriae*).

The *pbp2m* gene was located on the large plasmid in a
region comprising three other genes, a *blaB* beta-
lactamase family gene, a LysR family regulator gene, and
the *ermX* gene, flanked by two insertion sequences
(IS1628) of the IS6 family (Fig. 7). With disparate direct
repeat sequences, it remains unclear if this region

<figure><img src="image_10.png" /></figure>
represents a single transposable unit or a mosaic of gene acquisition events in pLRPD. A nearly identical PBP was observed in 78% of publicly available *C. jeikeium* genomes, in 57% of *C. striatum* genomes and in multiple other *Corynebacterium* genomes (Additional file 8: Table S7). However, the genetic context of PBP2m was highly variable in *C. diphtheriae* and among other *Corynebacterium* species (Additional file 4: Fig. S14). A putative transposable PBP-containing unit (PCU) comprising genes *pbp2m*, *blaB*, and *lysR*, commonly flanked by IS3503 (IS256 family) with a fragment identified in pLRPD (Fig. 7), appeared to be highly conserved and was associated variably with *ermX* in *C. diphtheriae* and with a helicase in *C. diphtheriae* and other *Corynebacterium* species. The PCU was sometimes found in 2 or 3 tandem copies and was chromosomally located in most genomes.

Further elements carried by pLRPD included an integer carrying genes *dfrA16*, *qacL*, and *sul1*, as well as elements of a putative conjugation apparatus gene cluster (Fig. 7). Our conjugation experiments aiming to demonstrate the transfer of pLRPD into recipient *C. diphtheriae* isolates failed. This plasmid was not found in other *C. diphtheriae* strains.

## Discussion

Strains of *C. diphtheriae* that are resistant to antimicrobial therapy may compromise the management of diphtheria cases and the control of pathogen transmission. Here, we aimed to define the genomic determinants of resistance to penicillin and other antimicrobial agents in *C. diphtheriae* and to analyze the relationships of resistance with diphtheria toxin production, biochemical variants, and phylogenetic sublineages. To this aim, we characterized phenotypically and genotypically a large sample of *C. diphtheriae* isolates from diverse geographic and temporal origins. We confirmed that the species is made of multiple phylogenetic sublineages [9, 13, 67] and showed that homologous recombination contributes five times more to their diversification than mutation, consistent with previous evidence of recombination in *C. diphtheriae* populations [12, 13].

Historically, *C. diphtheriae* isolates have been classified into three main biovars, but the links between biovars and phylogenetic structure have remained obscure. Whereas previous work concluded on the absence of an association [11], our phylogenetic analyses reveal that Gravis and Mitis, the two main biovars of *C. diphtheriae*, are associated strongly with two phylogenetic lineages. Lineage Gravis appears to have acquired ancestrally a gene cluster comprising the extracellular glycogen debranching enzyme gene *spuA* [62]. Although the Gravis phenotype was associated with *spuA* within lineage Gravis, other genomic determinants of glycogen

utilization remain to be discovered within the Mitis lineage. Whereas most of our biovar Belfanti isolates were excluded from this work because they belonged to *C. belfantii* or *C. rouxii*, a few Belfanti isolates did belong to *C. diphtheriae*. For three of them, mutations leading to a non-functional nitrate reductase were uncovered, explaining the phenotypic switch to biovar Belfanti. Our results show that biotyping is subject to parallel evolution, perhaps due to a yet-unknown selective pressure, and has limited species-level identification (as Belfanti is found in the three species *C. diphtheriae*, *C. belfantii*, and *C. rouxii*) or epidemiological typing value.

The most important factor of *C. diphtheriae* pathogenicity is the diphtheria toxin. Despite early realization that it is encoded on a prophage [5], few studies have investigated the phylogenetic distribution of the *tox* gene in *C. diphtheriae* [8, 9]. Here, we show that *tox*-positive strains mainly belong to the Mitis lineage and to early-diverging branches of the Gravis lineage. The distribution of *tox*-positive isolates into multiple Mitis sublineages is strongly indicative of independent acquisitions of the toxin gene. Alternately, this pattern might result from initial acquisition of the *tox* gene, followed by secondary loss in multiple sublineages. The phylogenetic pattern is also consistent with an ancestral presence of the *tox*-bearing phage in the Gravis lineage, with subsequent loss of the *tox* gene in the branch leading to the ancestor of most Gravis isolates. Future work should investigate the dynamics of the lysogenic corynephages and molecular determinants of their sublineage distribution. One important open question is the likelihood of *tox*-negative strains acquiring the *tox* gene during colonization, infection, or short-term epidemiological timeframes [7].

Remarkably, except for early-diverging sublineages, only one Gravis sublineage was found to carry the *tox* gene. This sublineage happens to correspond to the largest outbreak in recent times, which occurred in Newly Independent States of the ex-Soviet Union in the 1990s [13, 14, 38]. This sublineage, which comprises the ST8 reference strain NCTC13129, is genetically distant from other *tox*-positive lineages, which belong to the Mitis lineage. Hence, its antigenic structures or other pathogenicity properties may have diverged from those of more common *tox*-positive isolates, which might have contributed to its exceptional transmission in the 1990s, in addition to the decline in vaccine coverage [14]. Of note, biovar Gravis was named to reflect a perceived higher severity of infection compared to diphtheria cases caused by biovar Mitis isolates [68, 69]. Recently, it was shown that most diphtheria vaccines contain, besides the anatoxin, multiple other *C. diphtheriae* immunogens [22]. The impact of vaccination on the evolution of *C.*
*diphtheriae* populations, and possible variations of cross-protection as a function of strain diversity, is currently undefined. This work provides a framework onto which future studies can build to address this important question.

Although antimicrobial-resistant *C. diphtheriae* strains have been reported on numerous occasions [23, 32], knowledge on antimicrobial resistance in *C. diphtheriae* is largely fragmented and suffers from the lack of harmonization. Breakpoints used to define resistance vary according to world region and have changed over time within single countries [70]. The lack of consensus on the definition of resistance restricts our ability to define the magnitude of the problem and its global significance.

We aimed to define biologically meaningful cutoffs [64] based on susceptibility phenotype distributions, taking advantage of our large and diverse sample. Our data allowed us to propose tentative ecological cutoffs for *C. diphtheriae*. This approach should in the future be extended to MIC values [71] and should use larger and more diverse strain collections. Nevertheless, our analyses suggest that non-susceptibility to at least one antimicrobial agent was acquired by half of *C. diphtheriae* strains, regardless of lineage, biovar, or toxigenic status. This study further suggests that acquired resistance to penicillin, the first-line therapy against diphtheria, is far from being rare, affecting > 15% of *C. diphtheriae* isolates collected in the last decade in France and its overseas territories. Besides, 10.4% of these isolates were multidrug-resistant. The high prevalence of resistance to penicillin, tetracycline, and trimethoprim/sulfamethoxazole found here is consistent with susceptibility surveys of recent *C. diphtheriae* isolates in Algeria [29], Indonesia [31, 72], and India [73]. Many high-income countries such as France have chosen to use amoxicillin as the first choice for antibiotic therapy [74], as this molecule remains highly active. Still, widespread penicillin resistance is concerning, since diphtheria mainly occurs in resource-poor settings where penicillin G is largely used. In contrast, resistance to erythromycin and other macrolides remains rare. Our results call for concerted research into the magnitude of the antimicrobial resistance threat in *C. diphtheriae*.

Knowledge of the genetic mechanisms of antimicrobial resistance is critical for defining appropriate treatments, refining diagnostics, and conducting epidemiological studies of antimicrobial resistance. Resistance genes to several antimicrobial classes have been described in *C. diphtheriae* [36, 37, 75], while additional genes described in other *Corynebacterium* species [76] might also be present in *C. diphtheriae*. Here, we defined the prevalence and phylogenetic distribution of previously reported and newly identified resistance determinants in

*C. diphtheriae*. We demonstrate the co-occurrence of resistance phenotypes and genes, suggesting a causative link in multiple instances. We further show that resistance genes have been acquired independently in multiple sublineages, demonstrating a dynamic resistome in *C. diphtheriae*. In addition, we demonstrate an association between alterations in chromosomally encoded targets and phenotypic resistance for quinolone, trimethoprim, and rifampicin. Fluoroquinolone resistance was previously linked to mutations in the *gyrA* gene in *C. amycolatum* [77], *C. striatum* [78], and *C. belfantii* [79] but seemingly never for *C. diphtheriae*. Finally, we demonstrate the co-occurrence within some strains of multiple resistance determinants and uncover a previously undescribed large resistance plasmid in *C. diphtheriae*. The mechanism of genetic transfer of this plasmid remains to be investigated. This work provides the first overview of the *C. diphtheriae* resistome and will facilitate further studies into the evolutionary emergence of multiresistant *C. diphtheriae* strains.

Here, we discovered an accessory PBP (PBP2m) and experimentally showed it to confer resistance to penicillin and other beta-lactam antimicrobial agents. Its distribution in multiple sublineages, and its presence in other *Corynebacterium* species, clearly demonstrates its horizontal transfer, and we revealed a multiplicity of genomic contexts in which it is found within *Corynebacterium*. PBP2m is a putative low-affinity PBP, which would explain why it is less affected by beta-lactam antibiotics. Very recently, Forde et al. [80] have reported high-level penicillin resistance (MIC > 256 g/L) in *C. diphtheriae* strain BQ11 and have attributed this phenotype to a penicillin-binding protein located on transposon Tn3503. We found that the *pbp2m*-containing unit is 99.92% similar to Tn3503, which additionally harbors a helicase (Additional file 4: Fig. S14). Besides, PBP2m is 100% identical to the PBP from BQ11 and therefore also possesses the V361A alteration (using *C. jeikeium* as reference) that was hypothesized to confer high-level penicillin resistance [80]. Further studies on the expression of *pbp2m* in its variable genetic backgrounds are required to understand the observed difference in penicillin resistance levels. The location of *pbp2m* on a mobile element and its association with beta-lactam and carbapenem resistance represents a public health threat.

The seven chromosomal PBPs of *C. diphtheriae* (including the newly annotated PBP4b) were investigated to identify amino acid sequence polymorphisms associated with penicillin resistance. Although several alterations were significantly associated, none was directly implicated with the catalytic residues of the transpeptidase or transglycosylase domains. The association with resistance may not be directly linked to these catalytic residues but could be due to secondary sites that are
thought to interfere with beta-lactam ligand binding. While biochemical and structural studies are necessary to understand how these mutations affect penicillin susceptibility, we postulate that some changes in these domains could lead to allosteric effects ultimately resulting in beta-lactam resistance, as described for *Staphylococcus aureus* PBP2a [81, 82] or *Streptococcus pneumoniae* PBP2x [83]. Other SNPs might simply have been hitchhiking due to their physical linkage with functionally important SNPs [84].

## Conclusions

Our phylogenetic analyses showed that *C. diphtheriae* is highly diverse and comprises multiple sublineages that are grouped into two major lineages, characterized by distinctive associations with biovars Mitis and Gravis and with the diphtheria toxin gene. Antimicrobial resistance remains rare in France and its overseas territories, but multidrug-resistant strains were found and have a wide phylogenetic distribution. Analyses of the genetic bases of antimicrobial resistance defined the contribution of previously described mechanisms of resistance and uncovered a penicillin-binding protein, PBP2m, associated at the population level with penicillin resistance and shown experimentally to confer resistance to several beta-lactam antimicrobial agents. PBP2m is encoded in diverse genetic contexts, including on a large putatively conjugative resistance plasmid. This study reveals how biovars, diphtheria toxin gene, and antimicrobial resistance are associated with the phylogenetic diversity of *C. diphtheriae* strains, and provides novel data into the evolutionary dynamics of these medically important features. This framework will facilitate future studies investigating the emergence of sublineages of *C. diphtheriae* involved in diphtheria outbreaks, antimicrobial resistant infections, and the patterns and drivers of global spread of *C. diphtheriae* sublineages.

## Supplementary Information

Supplementary information accompanies this paper at https://doi.org/10.1186/s13073-020-00805-7.

**Additional file 1: Table S1.** Strains used in this study and their characteristics.

**Additional file 2: Table S2.** Antimicrobial agents and zone diameters interpretative breakpoints used.

**Additional file 3: Table S3.** Susceptibility data and interpretative categories based on current EUCAST recommendations (SIR) and proposed ecological cutoffs (ECOFF).

**Additional file 4.** Contains Supplementary Figures S1 to S14.

**Additional file 5: Table S4.** Molecular explanations for non-toxigenic tox gene-bearing (NTTB) strains.

**Additional file 6: Table S5.** Associations with resistance phenotypes based on treeWAS results for accessory genes, SNPs in core genes and whole-genome SNPs; or based on Kleborate or direct analysis of potential target sequences.

**Additional file 7: Table S6.** Susceptibility phenotypes of transformed *C. glutamicum* ATCC 13032 strains.

**Additional file 8: Table S7.** Detection of the pbp2m gene in other *Corynebacterium* genomes.

## Acknowledgements

We thank Vincent Enouf and the P2M core facility of Institut Pasteur for genomic sequencing. We are indebted to Collection de l'Institut Pasteur for providing reference strains of ribotypes, which were deposited following their described by Grimont and colleagues in 2004 [40]. We thank Valerie Bouchez for help with the Oxford Nanopore Technologies sequencing.

## Authors' contributions

S.B. conceived, designed, and coordinated the study. A.C-L., E.B., M.B., and M.D. performed the microbiological cultures of the isolates and their biochemical and molecular characterizations. M.H., L.G.P., C.R., S.L.B., M.B.-P., M.D., and S.B. analyzed the genomic and phenotypic data. J.T. reviewed the clinical source data of the isolates. Q.G. and A.-M.W. performed the pbp2m cloning experiments and analyzed the PBP sequences. X.D. provided help with the phylogenetic, recombination, and GWAS analyses. S.B. wrote the initial version of the manuscript. All authors provided input to the manuscript. All authors read and approved the final manuscript.

## Funding

MH was supported financially by a PhD grant from the European Joint Programme One Health, which has received funding from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 773830. LGP was supported financially by a grant from the Institut Français de Bioinformatique, the national infrastructure for services in bioinformatics created in the frame of the French Government's Investissement d'Avenir program. QG was funded by MTCI PhD school (ED 563). CR was supported by a Pasteur-Roux fellowship from Institut Pasteur. SLB received support from a Victorian Fellowship, provided by VESKI and funded by the State Government of Victoria, Australia.
The National Reference Center for Corynebacteria of the *Diphtheriae* complex receives support from Institut Pasteur and Public Health France (Santé Publique France, Saint Maurice, France). This work was supported financially by the French government's Investissement d'Avenir program Laboratoire d'Excellence "Integrative Biology of Emerging Infectious Diseases" (ANR-10-LABX-62-IBEID). AW and QG receive support from Institut Pasteur and the CNRS (France).

## Availability of data and materials

The genomic sequencing data generated in this study were deposited in the European Nucleotide Archive (ENA) database and are accessible through the BioProject PRJEB22103 (https://www.ebi.ac.uk/ena/browser/view/PRJEB22103) [85].

## Ethics approval and consent to participate

To conduct the research, we used bacterial strains, which are not considered human samples. Accordingly, this research was not considered human research and is out of the scope of the decree n° 2016-1537 of November 16, 2016, implementing law n° 2012-300 of March 5, 2012, on research involving human subjects. Therefore, no ethics approval was needed and no informed consent was required.

All French bacteriological samples were collected in the frame of the French national diphtheria surveillance and are collected, coded, shipped, managed, and analyzed according to the French National Reference Center protocols that received approval from the French competent body CNIL (Commission Nationale Informatique et Libertés; approval number: 1474671). Other strains were obtained from culture collections.

The research conformed to the principles of the Helsinki Declaration.

## Consent for publication

Not applicable

## Competing interests

The authors declare that they have no competing interests.
Author details

¹Institut Pasteur, Biodiversity and Epidemiology of Bacterial Pathogens, Paris, France. ²Collège doctoral, Sorbonne Université, F-75005 Paris, France. ³Institut Français de Bioinformatique, CNRS UMS 3601, Evry, France. ⁴Unité de Microbiologie Structurale, Institut Pasteur, CNRS UMR 3528, Université de Paris, F-75015 Paris, France. ⁵Doherty Applied Microbial Genomics, Department of Microbiology & Immunology, The University of Melbourne at The Peter Doherty Institute for Infection & Immunity, Melbourne, Victoria, Australia. ⁶Institut Pasteur, National Reference Center for Corynebacteria of the Diphtheriae Complex, Paris, France. ⁷School of Life Sciences and Department of Statistics, University of Warwick, Coventry, UK. ⁸Department of General Pediatrics and Pediatric Infectious Diseases, Hôpital Necker-Enfants Malades, APHP, Université de Paris, Paris, France.

Received: 29 May 2020 Accepted: 11 November 2020

Published online: 27 November 2020

References

1. Patey O, et al. Clinical and molecular study of Corynebacterium diphtheriae systemic infections in France. Coryne Study Group. J Clin Microbiol. 1997;35: 441-5.

2. Hadfield TL, McEvoy P, Polotsky Y, Tzinserling VA, Yakovlev AA. The pathology of diphtheria. J Infect Dis. 2000;181(Suppl 1):S116-20.

3. Burkovski, A. *Corynebacterium diphtheriae* and related toxigenic species: genomics, pathogenicity and applications. (Springer, 2014).

4. Sharma NC, et al. Diphtheria. Nat Rev Dis Primer. 2019;5:81.

5. Freeman VJ. Studies on the virulence of bacteriophage-infected strains of Corynebacterium diphtheriae. J Bacteriol. 1951;61:675-88.

6. Parveen S, Bishai WR, Murphy JR. Corynebacterium diphtheriae: diphtheria toxin, the tox Operon, and its regulation by Fe2+ activation of apo-DtxR. Microbiol Spectr. 2019;7(4):GPP3-0063-2019.

7. Pappenheimer AM, Murphy JR. Studies on the molecular epidemiology of diphtheria. Lancet Lond. Engl. 1983;2:923-6.

8. Trost E, et al. Pangenomic study of Corynebacterium diphtheriae that provides insights into the genomic diversity of pathogenic isolates from cases of classical diphtheria, endocarditis, and pneumonia. J Bacteriol. 2012; 194:3199-215.

9. Sangal V, Hoskisson PA. Evolution, epidemiology and diversity of Corynebacterium diphtheriae: new perspectives on an old foe. Infect. Genet. Evol. J. Mol. Epidemiol. Evol. Genet. Infect. Dis. 2016;43:364-70.

10. Timms VJ, Nguyen T, Crighton T, Yuen M, Sintchenko V. Genome-wide comparison of Corynebacterium diphtheriae isolates from Australia identifies differences in the Pan-genomes between respiratory and cutaneous strains. BMC Genomics. 2018;19:869.

11. Sangal V, et al. A lack of genetic basis for biovar differentiation in clinically important Corynebacterium diphtheriae from whole genome sequencing. Infect Genet Evol J Mol Epidemiol Evol Genet Infect Dis. 2014;21:54-7.

12. Bolt F, et al. Multilocus sequence typing identifies evidence for recombination and two distinct lineages of Corynebacterium diphtheriae. J Clin Microbiol. 2010;48:4177-85.

13. Grosse-Kock S, et al. Genomic analysis of endemic clones of toxigenic and non-toxigenic Corynebacterium diphtheriae in Belarus during and after the major epidemic in 1990s. BMC Genomics. 2017;18:873.

14. Dittmann S, et al. Successful control of epidemic diphtheria in the states of the former Union of Soviet Socialist Republics: lessons learned. J Infect Dis. 2000;181(Suppl 1):S10-22.

15. WHO. Surveillance and burden of diphtheria. http://www.who.int/ immunization/monitoring_surveillance/burden/diphtheria/en/. (2018).

16. Galazka A. The changing epidemiology of diphtheria in the vaccine era. J Infect Dis. 2000;181(Suppl 1):S2-9.

17. Rahman MR, Islam K. Massive diphtheria outbreak among Rohingya refugees: lessons learnt. J Travel Med. 2019;26(1):1-3.

18. Dureab F, Müller O, Jahn A. Resurgence of diphtheria in Yemen due to population movement. J. Travel Med. 2018;25(1):1-2.

19. Zakikhany K, Efstratiou A. Diphtheria in Europe: current problems and new challenges. Future Microbiol. 2012;7:595-607.

20. Meinel DM, et al. Outbreak investigation for toxigenic Corynebacterium diphtheriae wound infections in refugees from Northeast Africa and Syria in Switzerland and Germany by whole genome sequencing. Clin. Microbiol. Infect. Off. Publ. Eur. Soc. Clin. Microbiol. Infect. Dis. 2016;22:1003.e1-1003.e8.

21. Scheifer C, et al. Re-emergence of Corynebacterium diphtheriae. Med Mal Infect. 2019;49:463-6.

22. Möller J, et al. Proteomics of diphtheria toxoid vaccines reveals multiple proteins that are immunogenic and may contribute to protection of humans against Corynebacterium diphtheriae. Vaccine. 2019;37:3061-70.

23. Zasada, A. A. Antimicrobial susceptibility and treatment. in *Corynebacterium diphtheriae* and related toxigenic species genomics, pathogenicity and applications. Burkovski A., editor. 239-246 (New York: Springer, 2014).

24. Kupferschmidt K. Life-saving diphtheria drug is running out. Science. 2017; 355:118-9.

25. PAHO/WHO, P. A. H. O. / W. H. O. Epidemiological update: diphtheria. (2018).

26. Kneen R, et al. Penicillin vs erythromycin in the treatment of diphtheria. Clin Infect Dis Off Publ Infect Dis Soc Am. 1998;27:845-50.

27. Pereira GA, et al. Antimicrobial resistance among Brazilian Corynebacterium diphtheriae strains. Mem Inst Oswaldo Cruz. 2008;103:507-10.

28. von Hunolstein C, Scopetti F, Efstratiou A, Engler K. Penicillin tolerance amongst non-toxigenic Corynebacterium diphtheriae isolated from cases of pharyngitis. J Antimicrob Chemother. 2002;50:125-8.

29. Benamrouche N, et al. Microbiological and molecular characterization of Corynebacterium diphtheriae isolated in Algeria between 1992 and 2015. Clin. Microbiol. Infect. Off. Publ. Eur. Soc. Clin. Microbiol. Infect. Dis. 2016;22: 1005.e1-1005.e7.

30. Paveenkittiporn W, Sripakdee S, Koobkratok O, Sangkitporn S, Kerdsin A. Molecular epidemiology and antimicrobial susceptibility of outbreak-associated Corynebacterium diphtheriae in Thailand, 2012. Infect. Genet Evol J Mol Epidemiol Evol Genet Infect Dis. 2019;75:104007.

31. Husada D, et al. First-line antibiotic susceptibility pattern of toxigenic Corynebacterium diphtheriae in Indonesia. BMC Infect Dis. 2019;19:1049.

32. Bernard K, Pacheco AL. In vitro activity of 22 antimicrobial agents against Corynebacterium and microbacterium species referred to the Canadian National Microbiology Laboratory. Clin Microbiol Newsl. 2015;37:187-98.

33. Schiller J, Groman N, Coyle M. Plasmids in Corynebacterium diphtheriae and diphtheroids mediating erythromycin resistance. Antimicrob Agents Chemother. 1980;18:814-21.

34. Jellard CH, Lipinski AE. Corynebacterium diphtheriae resistant to erythromycin and lincomycin. Lancet Lond Engl. 1973;1:156.

35. Patey O, et al. Antibiotic susceptibilities of 38 non-toxigenic strains of Corynebacterium diphtheriae. J Antimicrob Chemother. 1995;36:1108-10.

36. Tauch A, Bischoff N, Brune I, Kalinowski J. Insights into the genetic organization of the Corynebacterium diphtheriae erythromycin resistance plasmid pNG2 deduced from its complete nucleotide sequence. Plasmid. 2003;49:63-74.

37. Barraud O, Badell E, Denis F, Guiso N, Ploy M-C. Antimicrobial drug resistance in Corynebacterium diphtheriae mitis. Emerg Infect Dis. 2011;17: 2078-80.

38. Cerdeño-Tárraga AM, et al. The complete genome sequence and analysis of Corynebacterium diphtheriae NCTC13129. Nucleic Acids Res. 2003;31:6516-23.

39. Farfour E, et al. Characterization and comparison of invasive Corynebacterium diphtheriae isolates from France and Poland. J Clin Microbiol. 2012;50:173-5.

40. Grimont PAD, et al. International nomenclature for Corynebacterium diphtheriae ribotypes. Res Microbiol. 2004;155:162-6.

41. Park WH, Williams AW. The production of diphtheria toxin. J Exp Med. 1896; 1:164-85.

42. Dazas M, Badell E, Carmi-Leroy A, Criscuolo A, Brisse S. Taxonomic status of Corynebacterium diphtheriae biovar Belfanti and proposal of Corynebacterium belfantii sp. nov. Int J Syst Evol Microbiol. 2018;68: 3826-31.

43. Badell, E. et al. *Corynebacterium rouxii* sp. nov., a novel member of the diphtheriae species complex. Res. Microbiol. 2020. https://doi.org/10.1016/j. resmic.2020.02.003.

44. Hauser D, Popoff MR, Kiredjian M, Boquet P, Bimet F. Polymerase chain reaction assay for diagnosis of potentially toxigenic Corynebacterium diphtheriae strains: correlation with ADP-ribosylation activity assay. J Clin Microbiol. 1993;31:2720-3.

45. Engler KH, Glushkevich T, Mazurova IK, George RC, Efstratiou A. A modified Elek test for detection of toxigenic corynebacteria in the diagnostic laboratory. J Clin Microbiol. 1997;35:495-8.

46. Magiorakos AP, et al. Multidrug-resistant, extensively drug-resistant and pandrug-resistant bacteria: an international expert proposal for interim
standard definitions for acquired resistance. Clin Microbiol Infect. 2012;18:268-81.

47. Badell E, et al. Improved quadruplex real-time PCR assay for the diagnosis of diphtheria. J Med Microbiol. 2019;68:1455-65.

48. Criscuolo, A. & Brisse, S. AlienTrimmer: a tool to quickly and accurately trim off multiple short contaminant sequences from high-throughput sequencing reads. Genomics. 2013. https://doi.org/10.1016/j.ygeno.2013.07.011.

49. Crusoe MR, et al. The khmer software package: enabling efficient nucleotide sequence analysis. F1000Research. 2015;4:900.

50. Liu Y, Schröder J, Schmidt B. Musket: a multistage k-mer spectrum-based error corrector for Illumina sequence data. Bioinforma Oxf Engl. 2013;29:308-15.

51. Bankevich A, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455-77.

52. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinforma. Oxf. Engl. 2014;30:2068-9.

53. Page AJ, et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinforma. Oxf. Engl. 2015;31:3691-3.

54. Didelot X, Wilson DJ. ClonalFrameML: efficient inference of recombination in whole bacterial genomes. PLoS Comput Biol. 2015;11:e1004041.

55. Guindon S, et al. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59:307-21.

56. Minh BQ, et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37:1530-4.

57. Collins C, Didelot X. A phylogenetic method to perform genome-wide association studies in microbes that accounts for population structure and recombination. PLoS Comput Biol. 2018;14:e1005958.

58. Sogues A, et al. Essential dynamic interdependence of FtsZ and SepF for Z-ring and septum formation in Corynebacterium glutamicum. Nat Commun. 2020;11:1641.

59. Mitchell AL, et al. InterPro in 2019: improving coverage, classification and access to protein sequence annotations. Nucleic Acids Res. 2019;47:D351-60.

60. Madeira F, et al. The EMBL-EBI search and sequence analysis tools APIs in 2019. Nucleic Acids Res. 2019;47:W636-41.

61. Ansari MA, Didelot X. Inference of the properties of the recombination process from whole bacterial genomes. Genetics. 2014;196:253-65.

62. Santos AS, et al. Searching whole genome sequences for biochemical identification features of emerging and reemerging pathogenic Corynebacterium species. Funct Integr Genomics. 2018;18:593-610.

63. Tagini F, et al. Distinct genomic features characterize two clades of Corynebacterium diphtheriae: proposal of Corynebacterium diphtheriae subsp. diphtheriae subsp. nov. and Corynebacterium diphtheriae subsp.lausannense subsp. nov. Front. Microbiol. 2018;9:1743.

64. Kahlmeter G, et al. European harmonization of MIC breakpoints for antimicrobial susceptibility testing of bacteria. J Antimicrob Chemother. 2003;52:145-8.

65. Tauch A, Kassing F, Kalinowski J, Pühler A. The Corynebacterium xerosis composite transposon Tn5432 consists of two identical insertion sequences, designated IS1249, flanking the erythromycin resistance gene ermCX. Plasmid. 1995;34:119-31.

66. Lavollay M, et al. The beta-lactam-sensitive D,D-carboxypeptidase activity of Pbp4 controls the L,D and D,D transpeptidation pathways in Corynebacterium jeikeium. Mol. Microbiol. 2009;74:650-61.

67. Seth-Smith HMB, Egli A. Whole genome sequencing for surveillance of diphtheria in low incidence settings. Front Public Health. 2019;7:235.

68. Barksdale L. Corynebacterium diphtheriae and its relatives. Bacteriol Rev. 1970;34:378-422.

69. McLeod, J. W. The types mitis, intermedius and gravis of corynebacterium diphtheriae: a review of observations during the past ten years. Bacteriol Rev 7, 1-41 (1943).

70. Zou J, et al. Phenotypic and genotypic correlates of penicillin susceptibility in nontoxigenic Corynebacterium diphtheriae, British Columbia, Canada, 2015-2018. Emerg Infect Dis. 2020;26:97-103.

71. Marosevic DV, et al. Antimicrobial susceptibility of Corynebacterium diphtheriae and Corynebacterium ulcerans in Germany 2011-17. J Antimicrob Chemother. 2020. https://doi.org/10.1093/jac/dkaa280.

72. Sariadji K, Puspandari SN, Sembiring M. Antibiotic susceptibility pattern of Corynebacterium diphtheriae isolated from outbreaks in Indonesia 2010-2015. Indones Biomed J. 2018;10:51-5.

73. Parande MV, Roy S, Mantur BG, Parande AM, Shinde RS. Resurgence of diphtheria in rural areas of North Karnataka, India. Indian J Med Microbiol. 2017;35:247-51.

74. Anonymous. Conduite à tenir lors de l'apparition d'un cas de diphtérie. https://www.hcsp.fr/explore.cgi/hcspr20110304_conduitediphterie.pdf. (2011).

75. Tauch A, Götker S, Pühler A, Kalinowski J, Thierbach G. The 27.8-kb R-plasmid pTET3 from Corynebacterium glutamicum encodes the aminoglycoside adenyltransferase gene cassette aadA9 and the regulated tetracycline efflux system Tet 33 flanked by active copies of the widespread insertion sequence IS6100. Plasmid. 2002;48:117-29.

76. Tauch A, Kassing F, Kalinowski J, Pühler A. The erythromycin resistance gene of the Corynebacterium xerosis R-plasmid pTP10 also carrying chloramphenicol, kanamycin, and tetracycline resistances is capable of transposition in Corynebacterium glutamicum. Plasmid. 1995;33:168-79.

77. Yoon S, Kim H, Lee Y, Kim S. Bacteremia caused by Corynebacterium amycolatum with a novel mutation in gyrA gene that confers high-level quinolone resistance. Korean J Lab Med. 2011;31:47-8.

78. Sierra JM, Martinez-Martinez L, Vázquez F, Giralt E, Vila J. Relationship between mutations in the gyrA gene and quinolone resistance in clinical isolates of Corynebacterium striatum and Corynebacterium amycolatum. Antimicrob Agents Chemother. 2005;49:1714-9.

79. Pivot D, et al. Carriage of a single strain of non-toxicigenic Corynebacterium diphtheriae biovar Belfanti (Corynebacterium belfantii) in four patients with cystic fibrosis. J Clin Microbiol. 2019. https://doi.org/10.1128/JCM.00042-19.

80. Forde BM, et al. Fatal respiratory diphtheria caused by β-lactam-resistant Corynebacterium diphtheriae. Clin Infect Dis Off Publ Infect Dis Soc Am. 2020. https://doi.org/10.1093/cid/ciaa1147.

81. Fishovitz J, et al. Disruption of allosteric response as an unprecedented mechanism of resistance to antibiotics. J Am Chem Soc. 2014;136:9814-7.

82. Otero LH, et al. How allosteric control of Staphylococcus aureus penicillin binding protein 2a enables methicillin resistance and physiological function. Proc Natl Acad Sci U S A. 2013;110:16808-13.

83. Bernardo-García N, et al. Allostery, recognition of nascent peptidoglycan, and cross-linking of the cell wall by the essential penicillin-binding protein 2x of Streptococcus pneumoniae. ACS Chem Biol. 2018;13:694-702.

84. Chewapreecha C, et al. Comprehensive identification of single nucleotide polymorphisms associated with beta-lactam resistance within pneumococcal mosaic genes. PLoS Genet. 2014;10:e1004547.

85. Brisse, S., Badell, E., Hennart, M. & Criscuolo, A. Genomic content and phylogenetic diversity of C. diphtheriae and related species were compared to clarify taxonomic relationships and study evolutionary trends. European Nucleotide Archive https://www.ebi.ac.uk/ena/browser/view/PRJEB22103 (2020).

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_11.png" /></figure>