Characterization and whole genome sequencing of closely related multidrug-resistant Salmonella enterica serovar Heidelberg isolates from imported poultry meat in the Netherlands

Multidrug-resistant Salmonella enterica serovar Heidelberg isolates are frequently recovered in the Netherlands from poultry meat imported from South America. Our aim was to retrospectively assess the characteristics of the antimicrobial determinants, gene content and the clonal relatedness of 122 unique S. Heidelberg isolates from chicken meat from Brazil (n = 119) and Argentina (n = 3) that were imported between 2010 and 2015. These isolates were subjected to antimicrobial susceptibility testing, PCR and Illumina HiSeq2500 whole genome sequencing. Draft genomes were assembled to assess the gene content, and the phylogenetic relationships between isolates were determined using single nucleotide polymorphisms. Ciprofloxacin-resistance was identified in 98.4% of the isolates and 83.7% isolates showed resistance to the extended-spectrum cephalosporins cefotaxime and ceftazidime (83.6% and 82.8% respectively). Of the latter, 97.1% exhibited an AmpC phenotype and contained blaCMY-2, whereas the remaining three isolates contained an extended spectrum beta-lactamase. Of the 99 extended-spectrum cephalosporins-resistant isolates harboring CMY-2 plasmids, 56.6% contained the incompatibility group I1 replicon. Phylogenetic cluster analysis showed that all isolates from Brazil clustered together, with 49% occurring in clusters larger than 5 isolates that revealed intra-cluster similarities based on geographical location and/or resistance profiles. The remaining isolates were classified in smaller clusters or as singletons, highlighting the large diversity of S. Heidelberg in the poultry chain in Brazil that was revealed by this study. Considering the potential public health risk associated with multidrug-resistant S. Heidelberg in imported poultry, collaborative whole genome sequencing-based surveillance is needed to monitor the spread, pathogenic properties and epidemiological distribution of these isolates.


Introduction
A set of control samples was included to act as internal controls (Table 1). First, controls were included which consisted of isolates from three distinct shipments of poultry, referred to as "batches", which arrived as a single shipment on the same date from the same processing plant. Two of these batches were sampled eight days apart in 2014 and one batch was sampled in 2015. Five bacterial isolates were recovered from each individual batch, making a total of 15 isolates. To prevent calculation bias, as the isolates in each separate batch were expected to be near identical S. Heidelberg variants, only one randomly selected isolate from each batch was included in the counts and statistics mentioned throughout the study. The second part of the control samples consisted of individual isolates that were sequenced multiple times. One S. Heidelberg isolate from Brazil (2013) was taken from frozen stock and sequenced twice (indicated as Duplicate A). Next, one S. Heidelberg isolate from Argentina (2015) was subcultured and a single colony was selected. The genome sequence of the same DNA extract was determined twice to identify potential WGS artefacts (indicated as Duplicate B). Finally, one S. Heidelberg isolate from Brazil (2012) was selected to be included in each of the three sequencing runs performed for this study (indicated as Triplicate).
In this study, we considered 114 isolates from Brazil and two isolates from Argentina. Additionally, one isolate from each batch (A, B and C; n = 3) was added, as well as the three control isolates that were sequenced multiple times (Duplicate A, Duplicate B and the Triplicate). In total, we analysed 122 unique isolates that were recovered from imported poultry. Any statistics reported throughout this paper will refer only to these 122 unique isolates, unless explicitly stated otherwise.
In total, 134 S. Heidelberg isolates from Brazil and Argentina (including the batches; n = 15) were analysed by Whole Genome Sequencing, generating 138 SNP profiles (including batch A, B and C, Duplicate A and B, and the Triplicate) that were included in a phylogenetic tree. All isolates included in this study were sequenced to a mean read depth of at least 41. To assess the performance of the WGS-based SNP analysis the read data of 19 S. Heidelberg isolates from various publicly available sources (see S1 Table in   five isolates were randomly selected and included in this study. In addition, the genome data of two S. Heidelberg isolates from Colombia [7] collected from fecal content and chicken meat (2013, 2012, respectively) and two isolates from Thailand (source and collection date unknown) were included in this study (ebi.ac.uk 100K Pathogen Genome Project). Two Salmonella enterica serovar Senftenberg isolates were included in the phylogenetic analysis to serve as an outgroup.

Bacterial culturing, antimicrobial susceptibility testing, and plasmid analysis
Cryotubes containing S. Heidelberg isolates were taken from the -80˚C freezer for subculturing. One bead from each cryotube was transferred on tryptic soy agar (TSA) and incubated overnight for 22h at 37˚C. Next, a single colony was subcultured to TSA and incubated overnight for 22h at 37˚C. Of each bacterial isolate a 0.5 McFarland suspension was prepared containing 5 ml sterilized water. Next, 55 μl suspension was added to 11 ml

DNA library preparation, genome sequencing, and assembly
The complete DNA content of each bacterial isolate was extracted from overnight cultures using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) and sequenced using the HiSeq 2500 platform (Illumina, San Diego, USA) according to the manufacturer's protocol. SNP profiles were generated by mapping the Illumina reads against S. Heidelberg reference genome SL476 (NC_011083, excluding plasmids) based on the Genome Analysis ToolKit (GATK) best practises. In short, Illumina adapters and low quality reads were removed using Trimmomatic (v0.35). The resulting reads were mapped to S. Heidelberg reference genome SL476 using BWA-mem (v0.7.12-5) and PCR duplicates were removed via Picard Tools (v1.113-2). Next, indel realignment was performed, genomic VCF files were generated using the GATK (v3.60) in ERC mode, and SNP variants were called using GATK GenotypeGVCFs. SNPs with a low coverage (read depth <10) and conflicting reads (<90% of reads agree with the called genotype) were not included. The filtered SNPs were then concatenated into SNP profiles.
To determine the gene content, trimmed reads were assembled using ABySS (v1.9.01), and BLAST (v2.2.31) was used to extract genes of interest from the assembly. Blast hits with a nucleotide identity below 80% were rejected. The MLST scheme from pubmlst.org was used to predict the sequence types. The downloaded databases of ResFinder (2018-09-10) and Plas-midFinder (2018-11-06), available at https://cge.cbs.dtu.dk/services/, were used with in-house developed tools to determine the antibiotic resistance determinants and plasmid types, respectively. PointFinder (commit bd50f0b) and the PointFinder database (commit d1413d2) were used to identify point mutations in chromosomal genes that confer antibiotic resistance [24].

Phylogenetic cluster analysis
Randomized Axelerated Maximum Likelihood (RAxML v8.2.4.1) was used for maximum likelihood phylogenetic tree estimation [25]. In short, the Gamma model for rate heterogeneity was used to allow for varying rates of evolution in different regions of the profile. The Lewis ascertainment bias correction was used to compensate for the fact that the SNP profile only includes variable positions. The bootstrap convergence of the phylogenetic tree was calculated and the extended majority rule consensus was applied to determine the required number of bootstrap iterations [26].

Antimicrobial resistance
Antimicrobial susceptibility testing of the 122 unique S. Heidelberg isolates showed that all isolates from Brazil (n = 119) and two of the three isolates from Argentina were multidrug-resistant to four or more classes of antimicrobial agents ( Table 2). Ciprofloxacin resistance (MIC of >0.064 mg/L) was identified in 98.4% (120 of 122) of the isolates and extended spectrum cephalosporin (ESC) resistance was identified in the majority of the isolates; 83.6% (102/122) were resistant to cefotaxime (MIC of >0.5 mg/L) and 82.8% (101/122) showed resistance to ceftazidime (MIC of >2 mg/L). Of the cefotaxime-resistant isolates, 99.0% (101/102) also showed resistance to ceftazidime.
Assessment of the antimicrobial resistance phenotypes of these isolates showed that 97.1% (99 of 102) ESC-resistant isolates exhibited an AmpC phenotype. Genome analysis revealed the presence of bla CMY-2 in all AmpC S. Heidelberg isolates and no discrepant results were identified following PCR confirmation. Of the three remaining ESC-resistant isolates that contained no bla CMY-2 , one isolate harbored bla CTX-M-2 while another isolate harbored bla CTX-M-8 . In the third isolate, association of bla CTX-M-8 and bla TEM-1 genes was revealed. In 20 of the 122 isolates that were susceptible to cefotaxime and ceftazidime, no β-lactamases were identified with the exception of one Argentinian isolate that harbored bla TEM-1 , which confers resistance to several penicillins. This isolate was the only isolate of these 20 isolates with resistance to ampicillin. Overall, 99 of 122 (81.1%) of the multidrug-resistant S. Heidelberg isolates from Brazil and Argentina harbored bla CMY-2 and exhibited an AmpC β-lactamase phenotype.

Identification of antimicrobial resistance genes and plasmids
A comparison between the phenotypic antimicrobial resistance data and the resistance gene WGS data showed that in addition to the bla genes, 16 different plasmid-associated antibiotic resistance genes were identified ( Table 2). Of the antibiotics included in this study, the identified resistance genes were expected to produce resistance to sulfamethoxazole, tetracycline, chloramphenicol, gentamycin, ciprofloxacin and trimethoprim. Our results confirmed that all resistance phenotypes correlated with the genotypes identified for these antibiotic agents. Ciprofloxacin resistance is often induced by mutations in gyrA and parC and in Salmonella these mutations are also related to nalidixic acid resistance. All of the 122 S. Heidelberg isolates had a chromosomal parC mutation (T57S). Of these, 118 (96.7%) isolates carried a chromosomal gyrA mutation (S83F) as well. Interestingly, three of the 118 (2.5%) isolates had also a plasmid that carried QnrB19. Two (1.6%) ciprofloxacin-resistant isolated lacked the gyrA mutation but harbored the QnrB19 gene. In the two (1.6%) ciprofloxacin-susceptible isolates, both QnrB19 and gyrA mutations were absent. QnrB19 is known to confer a low level of resistance to ciprofloxacin while also facilitating mutations in gyrA [27]. In the colistin-resistant isolate no mcr genes (1 to 5) were identified and the MIC of >4 mg/L is just one step above the cut-off (MIC of >2 mg/L).
After de novo assembly of the sequencing data, the presence of several plasmid replicons was revealed by using the PlasmidFinder database (Table 3). Co-carriage of plasmids was common in this isolate collection. Large multidrug-resistant plasmids of the IncX1 and IncA/C incompatibility groups were most frequently identified and were present in 120 of 122 (98.4%) and 118 of 122 (96.7%) of the isolates, respectively. IncI1 was present in 65 of 122 (53.3%) of the isolates. Of the 99 ESC-resistant isolates harboring bla CMY-2 , IncI1 was identified in 56 (56.6%) isolates. To verify whether the plasmid replicons of these incompatibility groups identified were present, transformation experiments were conducted on an arbitrarily selected subset of 45 isolates. All identified IncI1 replicons (n = 32) were confirmed and no discrepant results were identified. Ninety-seven (98.0%) of these 99 isolates carried IncX1 and IncA/C plasmids of which 54/97 (54.5%) co-carried IncI1. Notably, the four different plasmid profiles that occur in 91.0% of the isolates (111/122) exhibited 11 different antimicrobial phenotypes. Additional genomic analysis of the 65 isolates containing an IncI1 plasmid replicon revealed that the primary plasmid multilocus sequence type (pMLST) was ST12 (n = 56; 86.2%). The remaining isolates exhibited ST178 (n = 4; 6.2%), ST113 (n = 3; 4.6%), ST26 (n = 1; 1.5%) and one of the isolates from Argentina showed a combination of known plasmid alleles that was not part of any of the known types, suggesting this plasmid is a novel sequence type.

MLST analysis
Genome data analysis to establish the multilocus sequence type (MLST) revealed that all but one of the isolates from Brazil and Argentina were indistinguishable, as all of these isolates were assigned sequence type 15 (n = 121). The only exception was a single isolate from Brazil that was more distantly related to the other Brazilian isolates, and was assigned ST4632 (n = 1) due to a single SNP in thrA1. This outlier was confirmed by resequencing, and all other control samples were assigned MLST type 15 (Table 4).

Cluster analysis of the control set
To assess the validity of the high resolution cluster classification based on SNP analysis, we first investigated the clustering of the control samples. The control samples, derived from the sampling of three distinct shipments of poultry meat, assigned batch A, B, and C (Table 4), were grouped together with a bootstrap support values of 95, 96, and 100 and with a maximal SNP distance between isolates of 1, 2, and 2 SNPs, respectively (Fig 1).
Next, the classification of the control samples that were sequenced multiple times (Table 4) was assessed as illustrated in Figs 1 and 2.
The SNP profiles of the control sample duplicate A from Brazil, which was sequenced twice, clustered together with a bootstrap support value of 100, and 0 SNPs. The SNP profiles from duplicate B, derived from the same DNA extraction step, showed zero SNPs difference and were paired together with a bootstrap support value of 100. The three SNP profiles that were derived from the single isolate indicated as triplicate, were placed together with a bootstrap support value of 94, and a maximal distance of one SNP. We therefore decided to set the cut-off value for identical isolates at 5 SNPs, since we are not able to reliably differentiate between isolates that differ less than five SNPs, using our pipeline. Next we set out to establish a terminology to describe clusters of more distantly related isolates. Note that this terminology   Characterization and WGS-based typing of multidrug-resistant S. Heidelberg was only used when the reliability of a cluster had been established, specifically if the bootstrap support value was at least 80. In order of increasing distance, we described isolates as identical (0-5 SNPs), highly similar (6-10 SNPs) or similar (11-30 SNPs). The arbitrary aspect of these cut-off values mean that they are well suited for the analysis of this sample set using our pipeline and parameters. Therefore, these cut-off values should not be taken as a general guideline, but rather as a situational convention [28].

Cluster analysis
The RAxML tree in S1 Fig shows that all S. Heidelberg isolates formed a clade, with an average SNP distance of 40,076 SNPs from the out-group that comprised Salmonella enterica subsp. enterica serovar Senftenberg. An average SNP distance of 27 and 29 SNPs was found between the isolates from Brazil and between the isolates from Argentina, respectively, whereas 131 SNPs were identified between the isolates of both countries. An average number of 7350, 100, 34, 85, 69 SNPs were found between the Brazilian S. Heidelberg isolates and isolates from Thailand, Colombia and Canadian outbreak I, II and III, respectively. Fig 2 also shows that our method is able to differentiate between the SNP profiles from isolates from various geographical areas and outbreaks by grouping them with good bootstrap support. All SNP profiles from Brazil cluster exclusively together, as do the SNP profiles from Argentina, Colombia, Thailand and the three Canadian outbreaks.
Eighty-four of the 122 isolates (68.9%) were assigned to at least one cluster, while the remaining 38 isolates (31.1%) were singletons. Note that isolates can be a member of multiple clusters, e.g. when a highly similar cluster contains two or more identical isolates. The isolates from Brazil and Argentina were classified in a total of 14 clusters (not counting the batches) with identical isolates (0 to 5 SNPs between isolates), ranging in size from 2 to 7 isolates. The period between sampling in these clusters varied from the same day to isolates that were sampled 20 months apart. A total of 38 (31.1%) isolates were assigned to 11 clusters with highly similar isolates (a SNP difference between 6 and 10), ranging from 2 to 9 isolates. The period between sampling in these clusters varied from the same day to isolates that were sampled approximately 23 months apart. When we considered similar clusters, with a SNP difference of 11 to 30 SNPs between isolates, a total of 66 (54.1%) isolates were assigned to ten clusters with similar isolates, ranging from 2 to 13 isolates. The period between sampling in these clusters varied from approximately three months to almost 5.5 years between the isolates from Argentina. In the following section we only describe six clusters (I to VI) that contain at least 5 similar isolates (11-30 SNPs, see Table 5).
The 119 unique isolates from chicken meat from Brazil originated from two adjacent states, Paraná (n = 79) and Santa Catarina (n = 40). The isolates from chicken meat that originated from Paraná came from six different poultry processing plants, located in Cafelândia, Cascavel, Marechal Cândido Rondon, Palotina, Rolândia and Toledo. The isolates from chicken meat that originated from Santa Catarina came from three different poultry processing plants, located in Chapecó, Concórdia and Nova Veneza.
In the clusters with similar isolates, clusters I, III, IV, V and VI, a total of 1/7 (14.3%), 13/13 (100%), 7/7 (100%), 5/6 (83.3%) and 3/12 (25.0%) of the isolates originated from Paraná, respectively. These clusters all contained isolates from chicken meat that originated from a poultry processing plant in the city of Toledo in the Paraná state. The remaining isolates from Paraná were classified in two small clusters with identical isolates, four clusters with highly similar isolates, and as 34 (43.0%) singletons. These clusters and singletons were unique and thus not part of a larger cluster.
We found that the isolates from chicken meat that originated from Santa Catarina (n = 40) showed a strong tendency to group together (Figs 1 and 2). In the clusters with similar isolates, clusters I, II, V and VI, a total of 6/7 (85.7%), 13/13 (100%), 1/6 (16.7%) and 9/12 (75.0%) of the isolates from chicken meat originated from Santa Catarina, respectively. Interestingly, the subcluster with highly similar isolates that was part of cluster VI, exclusively contained isolates from Santa Catarina. Cluster I contained three isolates from chicken meat from Chapecó and three isolates from Concórdia that all originated from the Santa Catarina state. Cluster II exclusively contained isolates from a poultry processing plant in Chapecó. Cluster V only contained one isolate from Concórdia and cluster VI contained nine unique isolates from two poultry processing plants in the city of Chapecó. The remaining isolates from chicken meat that originated from Santa Catarina (n = 11) were classified in one cluster with identical Location poultry processing plant (n) Chapecó (3) Chapecó (13) Toledo (13) Toledo (7) Concórdia (1) Chapecó (9) Concórdia (3) Toledo (5) Toledo (3) Toledo (1) Approx. sampling period in months The three primary antimicrobial resistance phenotypes (Table 2) also grouped together in the tree as illustrated in Figs 1 and 2. Clusters I to V shared the most common ESC resistant antimicrobial resistance phenotype AMP, CTX, CAZ, CIP, NAL, SMX, TET. All the isolates in cluster VI contained the phenotype CIP, NAL, SMX, TET with the exception of one ESC resistant isolate that exhibited an AmpC β-lactamase phenotype, harbored bla CTX-M-8 and displayed the AMP, CTX, CAZ, CIP, CST, NAL, SMX, TET phenotype. None of the other isolates contained bla CMY-2 or were ESC resistant. The phenotype CIP, NAL, SMX, TET was only identified in two other isolates. One isolate belonged to a cluster with three similar isolates and one isolates was classified in a highly similar cluster of only two isolates. Isolates with antimicrobial resistance phenotype CIP, GEN, NAL, SMX, TET clearly grouped together, even though only three of the six isolates were classified to a small cluster with identical isolates. Six different plasmid profiles were identified in clusters I to VI, ranging from two to five per cluster. Notably, seven of the 13 isolates in cluster II, belong to a subcluster with identical isolates that contained plasmid replicons IncX, IncA/C, ColpVC and IncI1. The sampling period between the isolates in clusters I to VI varied from approximately 2.5 to 27.5 months. The three Argentinian S. Heidelberg isolates grouped together despite a delay of almost 5 years and 7 months between the first and last isolate sampled. A maximal SNP distance of 63 SNPs between the isolates indicates they are more distantly related.
The Canadian outbreak isolates were classified in three separate clusters. The maximal SNP distance between the isolates of each separate outbreak cluster was one SNP. The isolates from Thailand and Colombia were grouped in pairs with a maximal SNP distance of 7 and 43 between isolates, respectively. Interestingly, the Colombian isolates were grouped next to the isolates of Canadian outbreak III and a maximal SNP distance of 75 SNPs was found between the Colombian and the Canadian outbreak III isolates. This makes these clusters distantly related according to the terminology used for the phylogenetic classification.

Discussion
In this retrospective study, we show the emergence of clonally related multidrug-resistant (MDR) S. Heidelberg isolates from poultry meat that was imported to the Netherlands between 2010 and 2015. Most of these isolates originated from ten different cities of two adjacent states from the Southern part of Brazil, Paraná (n = 79) and Santa Catarina (n = 40) and were all but one assigned ST15, a sequence type commonly identified among S. Heidelberg strains and often reported in Asia, Europe and the Americas in both food and human isolates [5,6,29]. Interestingly, the only exception was the most distant isolate relative to the other isolates from Brazil that was assigned ST4632 based on a single SNP in one of the MLST genes. Nearly all isolates exhibited fluoroquinolone resistance (98.4%), a worrying observation as these antibiotics are considered the first-line antimicrobial agents for the treatment of salmonellosis [30][31][32]. In 2016, the EFSA also reported a substantial number (64.7%) of ciprofloxacin-resistant isolates (breakpoint MIC of >0.064 μg/ml, EUCAST) from broiler meat in 16 of the 19 European Member States that participated. In contrast, the NARMS reported that in 2015 only 0.7% of the retail chicken meat samples in the United States were resistant to ciprofloxacin (breakpoint MIC of �0.12 μg/ml, CLSI) [33,34]. Even more disturbingly, 83.6% of the isolates were resistant to cefotaxime. Since ESCs are often administered in severe cases of salmonellosis in risk groups, or in cases of infection caused by fluoroquinolone-resistant strains, these results show the serious limitations concerning effective treatment in case of severe infection with these S. Heidelberg strains imported from Brazil [35]. In the isolates imported from Brazil, 81.5% displayed an AmpC phenotype and bla CMY-2 was the most frequently identified β-lactamase in these isolates. Although limited information is available on the prevalence of bla  in Salmonella in Brazil, Giuriatti et al. 2017 reported CMY-2-positive S. Heidelberg isolates that exhibited ESC resistance from Paraná state, Brazil, which is where most isolates in this study were imported from [10]. Moreover, bla CMY-2 was reported to be widely distributed in Escherichia coli isolates detected in Brazilian chicken carcasses [36]. Besides bla CMY-2 , several ESBL genes (e.g. TEM-, SHV-and CTX-M-type ESBLs) have been associated with the appearance of ESC resistance in Salmonella [37]. In our study, we identified only four isolates containing ESBL genes (bla CTX-M-2, bla TEM-1, and bla CTX-M-8 ) using genome analysis. Despite the limited number of ESBLs detected in this study, the CTX-M-2 and CTX-M-8 enzymes are nowadays highly prevalent in South America and contribute to ESC resistance in several Salmonella enterica serovars [8,38,39]. Until recently CTX-M-8 enzymes were rarely identified in S. Heidelberg, but a recent report also described the presence of bla CTX-M-8 in an S. Heidelberg isolate from Brazilian poultry meat [40]. Notably, no ESC resistance was identified in absence of a corresponding bla CMY-2 or ESBL gene in this study.
It appears these resistance determinants are widely disseminated among various Salmonella serovars in Brazil, including S. Heidelberg [8,38,40,41]. The identification of these clonally related multidrug-resistant S. Heidelberg variants originating from several Brazilian states suggests the successful expansion and persistence of these apparently biologically fit variants. The extensive use of antibiotics in the poultry industry for the treatment and prevention of disease, as well as growth promotion in broilers, can create favorable conditions for the acquisition of resistance determinants [42][43][44]. Commensal microbiota in the gastrointestinal tract serve as an especially dangerous reservoir for antimicrobial resistance as these bacteria can rapidly be altered by antibiotic pressure [45]. Interspecies horizontal transfer and/or inter-serovar exchange of resistance determinants can provide Salmonella with a selective advantage to facilitate spread and persistence. Although the logistics concerning the Brazilian poultry industry are not entirely clear, contamination with MDR S. Heidelberg could occur during (inter) national exchange of broiler stock and/or equipment between poultry farms or during the transport and storage of poultry carcasses [9]. The environmental persistence of MDR S. Heidelberg in poultry farms may allow for continued exposure and the often asymptomatic infection of live poultry [46,47]. Another reason that allowed these S. Heidelberg variants to survive might be inadequate hygienic procedures at the poultry processing plants. Especially problems with the pre-chilling procedure, which has been reported as a critical processing step to eradicate Salmonella spp., can contribute to the contamination of the final poultry products with viable Salmonella spp. [48]. In turn, this might also lead to the persistence of pathogens on processing equipment, a known risk for cross contamination [48,49].
Of the plasmids belonging to known incompatibility groups, transmissible plasmids IncX1, IncA/C and IncI1 were most common. These plasmids are often associated with persistence, as well as carrying virulence and antimicrobial resistance in poultry by disseminating ESBL genes [50,51]. Our results also strongly suggest that most of the bla CMY-2 identified were located on the plasmids of the incompatibility groups IncI1 and IncA/C, although in some cases a chromosomal association may exist, as transformation experiments and replicon typing were only performed on a representative subset of samples collected [29, [52][53][54].
To further assess the potential variation between these clonally related isolates, we performed high-resolution subtyping by using 138 SNP profiles (including the control samples) for phylogenetic cluster analysis. The correct classification of the control samples was a strong indication that the phylogenetic allocation of the filtered SNP profiles provided the means to discriminate between different S. Heidelberg variants while maintaining meaningful groups.
For the purpose of discussion, we adopted an arbitrary terminology for cluster classification, partially based on the differences of the SNP profiles between matching control samples, to assess three different levels of genotyping resolution. Interestingly, instead of large groupings of Brazilian S. Heidelberg isolates, many small (sub)clusters with less than five isolates were observed. Only when we considered clusters which contained isolates that differed between 11 to 30 SNPs, six substantial clusters (of at least five isolates) with good bootstrap support were identified, meaning that only 58 isolates were assigned to a substantial cluster. Thus, even after using the latter cluster classification, 31.1% (37/119) of the Brazilian S. Heidelberg isolates appeared as a singleton or were assigned to a small cluster (20.2%; 24/119). The limited number of substantial clusters strongly suggests that in a large clonally related population of S. Heidelberg strains, subtyping resulted in the identification of numerous S. Heidelberg variants that circulate in the states of Paraná and Santa Catarina. It appears that in Brazil several S. Heidelberg lineages, presumable of common ancestry, contained an epidemic strain that underwent clonal expansion and diversification. Resistance determinants might have been acquired by ancestral strains, while maintaining biological fitness, conferring antimicrobial resistance to the epidemic strains. Microevolution during the expansion of the epidemic clones could explain the high number of S. Heidelberg variants identified [55]. Moreover, it is likely that the substantial number of variants we presented in this study are an underestimation of the real variation of S. Heidelberg isolates that are circulating in the Brazilian poultry chain, taking into account the high number of singletons and the fact that a relatively small sample size was presented provided by random sampling of imports from the southern part of Brazil.
We were able to differentiate six clusters, containing five or more isolates, which showed similarities in geographical location and antimicrobial resistance phenotype. Assessment of the geographical location of the isolates showed that the isolates from different poultry processing plants in Santa Catarina had a tendency to cluster together, which is clearly illustrated by the separate clusters I, II and VI. These results suggest that, although not exclusively, a persistent S. Heidelberg clone was circulating on the poultry production farms and/or poultry processing plants located in this state. Although the resistance markers might not be ultimately conclusive, the classification of cluster VI was reinforced by the antimicrobial resistance phenotypic data. Eleven of the 12 (91.7%) isolates in cluster VI were susceptible for the ESCs tested and showed an MDR phenotype (CIP, NAL, SMX, TET) that diverged from the dominant MDR phenotype found in this study (CTX, CAZ, CIP, AMP, NAL, SMX, TET). The isolates in this cluster clearly contained an S. Heidelberg variant that is, with a maximal SNP distance of 54 SNPs, more distantly related to the isolates in clusters I to V. However, we are aware that the lack of profound epidemiological data concerning our isolate collection hindered in-depth epidemiological characterization of these high-resolution (sub)clusters.
The S. Heidelberg isolates from Argentina, Thailand, Colombia and Canada showed a strong geographical signature and could be clearly differentiated from the Brazilian isolates, despite the fact that almost all Brazilian isolates were classified as ST15. Each of the three individual Argentinian isolates presented a unique and divergent MDR phenotype compared to the Brazilian MDR phenotypes, as well as an average SNP distance of 90 SNPs, indicating that the Argentinian S. Heidelberg isolates had no direct link to the Brazilian isolates. The isolates from different Canadian outbreaks were assigned to three different clusters, which is consistent with the findings of Bekal et al. 2016, supporting the validity of the SNP analysis in this study [23]. Notably, the Colombian isolates that were collected from poultry fecal matter and chicken meat, appear to be closest related (an average of 65 SNPs) to the human clinical isolates from Canadian outbreak III. The link between the consumption of poultry products and human infection by Salmonella is well documented. Several studies reported that the presence of Salmonella in healthy poultry is considered as one of the major transmission risks for human infection by direct transmission or via eggs and poultry meat [2,56,57].
In this study, we demonstrated the high discriminative power of WGS-based SNP analysis to distinguish clonally related S. Heidelberg isolates from imported chicken meat that were multidrug-resistant against medically important antibiotics, including ESC's. We revealed a high genetic diversity among the S. Heidelberg variants that circulate in Brazil, while maintaining meaningful groups of S. Heidelberg variants that were outbreak related and originated from other geographical locations over a substantial period of time. The phylogenetic relationships suggest that there is evidence that implies the variants from Brazil originated from a common ancestor. Yet, the different plasmid profiles also suggest that the ESC resistance can be attributed to the wide dissemination of several plasmids that are transferred by horizontal transmission between identified and still unidentified S. Heidelberg variants as well as interspecies plasmid transfer. Although we are confident that our method is applicable for the efficient monitoring of local outbreaks and the tracing of variants in a contaminated production environment, we are aware that the current data set is limited in size. Studies with a larger sample size and extensive epidemiological data, collected over a longer time-span, are required to determine suitability of the method for large scale analysis. Nevertheless, our data provides additional insight into the spread of MDR S. Heidelberg variants and the reservoir of resistance determinants in the Brazilian poultry industry. These findings are serious reason for concern and should be addressed. The MDR character of these S. Heidelberg variants pose a potential health hazard, and a collaborative WGS-based surveillance by food safety authorities is needed to provide a deeper understanding of the pathogenicity, antimicrobial resistance and epidemiological spread of these isolates.  Table. Public read data that was included in this study. The first column is the sample name used in this study, the remaining columns are data taken from the SRA archive. Heymans.