Genomic Characterization of Extended-Spectrum Cephalosporin-Resistant Salmonella enterica in the Colombian Poultry Chain

Salmonella enterica serovars have been isolated from Colombian broilers and broiler meat. The aim of this study was to investigate the diversity of ESBL/pAmpC genes in extended-spectrum cephalosporin resistant Salmonella enterica and the phylogeny of ESBL/pAmpC-carrying Salmonella using Whole Genome Sequencing (WGS). A total of 260 cefotaxime resistant Salmonella isolates, obtained between 2008 and 2013 from broiler farms, slaughterhouses and retail, were included. Isolates were screened by PCR for ESBL/pAmpC genes. Gene and plasmid subtyping and strain Multi Locus Sequence Typing was performed in silico for a selection of fully sequenced isolates. Core-genome-based analyses were performed per ST encountered. blaCMY−2−like was carried in 168 isolates, 52 carried blaCTX−M−2 group, 7 blaSHV, 5 a combination of blaCMY−2−like-blaSHV and 3 a combination of blaCMY−2−like-blaCTX−M−2 group. In 25 isolates no ESBL/pAmpC genes that were screened for were found. WGS characterization of 36 selected strains showed plasmid-encoded blaCMY−2 in 21, blaCTX−M−165 in 11 and blaSHV−12 in 7 strains. These genes were mostly carried on IncI1/ST12, IncQ1, and IncI1/ST231 plasmids, respectively. Finally, 17 strains belonged to S. Heidelberg ST15, 16 to S. Paratyphi B variant Java ST28, 1 to S. Enteritidis ST11, 1 to S. Kentucky ST152 and 1 to S. Albany ST292. Phylogenetic comparisons with publicly available genomes showed separate clustering of Colombian S. Heidelberg and S. Paratyphi B var. Java. In conclusion, resistance to extended-spectrum cephalosporins in Salmonella from Colombian poultry is mainly encoded by blaCMY−2 and blaCTX−M−165 genes. These genes are mostly associated with IncI1/ST12 and IncQ1 plasmids, respectively. Evolutionary divergence is observed between Colombian S. Heidelberg and S. Paratyphi B var. Java and those from other countries.


INTRODUCTION
Salmonella Enteritidis and S. Typhimurium have been reported as the most frequent serovars causing salmonellosis in humans worldwide (Hendriksen et al., 2011). According to data collected in the European Union (EU) in 2015 and the United States (USA) in 2013, S. Enteritidis accounted for 46 and 15% of Salmonella infections and S. Typhimurium for 16 and 13%, respectively (Centers for Disease Control and Prevention (CDC), 2016; European Food Safety Authority (EFSA) and European Centre for Disease Prevention and Control (ECDC), 2016). Likewise, S. Typhimurium and S. Enteritidis are the most isolated serovars from human cases in Colombia. Their overall prevalence among human isolates between 2005 and 2011 was 32 and 28%, respectively (Rodríguez et al., 2017).
Among foods of animal origin, poultry products (e.g., eggs) have been primarily associated with S. Enteritidis infections in humans while S. Typhimurium infections are associated with a wider range of sources including pork, beef and poultry products (Mughini-Gras et al., 2014;Antunes et al., 2016). Nevertheless, in broilers and chicken meat the prevalence of serotypes other than S. Enteritidis and S. Typhimurium have been on the rise over the last 2 decades (van Pelt et al., 2003;Foley et al., 2011;Wagenaar et al., 2013). Baseline studies of Salmonella  Canada, 2016). In Colombia, baseline studies performed by the Colombian integrated program for antimicrobial resistance surveillance (Coipars), demonstrated S. Paratyphi B var. Java, and S. Heidelberg to be the most prevalent serovars in broiler farms and meat at retail. Serovar distribution at farm level was 76 and 23%, respectively (Donado-Godoy et al., 2012b), and at the retail level 45 and 19%, respectively (Donado-Godoy et al., 2014). A similar distribution of S. Paratyphi B var. Java and S. Heidelberg, has been reported in chicken meat at retail in Guatemala (35 and 16% respectively) (Jarquin et al., 2015) and broilers at slaughter in Venezuela (62 and 31% respectively) (Boscán-Duque et al., 2007).
Isolates of S. Paratyphi B var. Java and S. Heidelberg are often multi drug resistant (Denny et al., 2007;Dutil et al., 2010;Antunes et al., 2016;Liakopoulos et al., 2016) and have been associated with carriage of plasmid-mediated extended spectrum β-lactamases (ESBLs) and plasmid associated AmpC β-Lactamases in poultry (pAmpC) (Antunes et al., 2016). For instance, S. Paratyphi B var. Java isolates from poultry in Europe have been found carrying the ESBL gene bla CTX−M−1 and the pAmpC gene bla CMY−2 (Doublet et al., 2014;Mevius et al., 2015;Veldman et al., 2016). In turn, S. Heidelberg has been associated to bla CMY−2 in North America (Andrysiak et al., 2008;Folster et al., 2012) and poultry products imported from South America into Europe (Liakopoulos et al., 2016). In South America S. Heidelberg has been associated with bla CTX−M−2 (Antunes et al., 2016). To date, no data of the genetic determinants of Extended Spectrum Cephalosporin (ESC) resistance in Salmonella from broilers and chicken meat in Colombia are available. These data are highly relevant to understand the epidemic spread of ESBL/pAmpC-producing Salmonella in poultry resulting in frequent occurrence in chicken meat in multiple countries.
Serotyping has traditionally been used for the epidemiological investigation of Salmonella, but does not provide information about the evolutionary relatedness of strains. Sequence based methodologies such as Multi Locus Sequence Typing (MLST) and Whole Genome Sequencing (WGS) have been proposed as a replacement of serotyping to identify evolutionary and epidemiological relatedness (Achtman et al., 2012;Ashton et al., 2016;Nadon et al., 2017). Additionally, information on the genetic basis of antimicrobial resistance (AMR) and plasmids harboring AMR genes can be readily obtained from WGS assemblies (Zankari et al., 2012;Carattoli et al., 2014). Altogether, the objectives of this study were to investigate the diversity of ESBL/pAmpC genes and encoding plasmids found in ESCresistant S. enterica from broilers and broiler meat in Colombia and to determine the genetic relatedness with Salmonella strains from other countries using MLST and core-genome alignments.

Isolates of Salmonella enterica
The isolates included in this study originated from different cross-sectional baseline studies conducted between 2008 and 2013 at different stages during the development of Coipars. In these studies, non-clinical samples were obtained from three different levels of poultry production in Colombia: Broiler farms (Donado-Godoy et al., 2012b), broilers at slaughter (Donado-Godoy et al., 2015b) and broiler meat at retail (Donado-Godoy et al., 2012a, 2015a. The methodology for sampling, random isolation (i.e., without antimicrobials during enrichment) of Salmonella and antimicrobial susceptibility testing with the BD Phoenix automated system, was previously described in detail for the studies in broiler farms, broilers at slaughter and broiler meat at retail mentioned above. Previous results of the prevalence of Salmonella and resistance to ESC at the different levels of poultry production, are summarized in Table 1.
For the present study, all available S. enterica isolates (n = 673) were considered. The 673 isolates belonged to 578 production flocks, 28 flocks from broiler farms, 140 from slaughterhouses and 410 from retail (Table 1). Next, ESC-resistant isolates . For all flocks with multiple ESC-resistant isolates, each first isolate was selected to make sure we included epidemiologically unrelated isolates only. As a result, a total of 260 isolates were selected for ESCresistance characterization. The selected isolates from farms originated from drag swabs of fecal material from the floor of broiler houses and fresh feces directly from the chicken. These isolates originated from 19 broiler farms (n = 19) (Donado-Godoy et al., 2012b). Isolates at slaughter, originated from cecal content and carcass rinse from 32 slaughterhouses (n = 84) (Donado-Godoy et al., 2015b). At retail, isolates originated from carcass rinse from 143 retail suppliers (Donado-Godoy et al., 2012a and chicken thighs meat from 8 retail suppliers (n = 157) (Donado-Godoy et al., 2015a). The isolates originated from 18 out of 32 departments (i.e. provinces) of Colombia. Together, these 18 departments are responsible for more than 90% of the chicken population in the country. A map of Colombia with the location of these departments is available in Supplementary Figure 1.
Information related to the origin of the samples and the prevalence of Salmonella is shown in Table 1.

Selection of Strains and WGS
Based on the diversity of gene families found after PCR screening, a representative selection of isolates was made and further subjected to WGS. Due to its high prevalence, a random selection was made using the Random function in Microsoft R Excel to assign random numbers to isolates positive for bla CMY−2−like and bla CTX−M−2 group. The size of selection was fixed to the square root of the number of resulting positive strains for these genes. This resulted in 13 bla CMY−2−like -and 8 bla CTX−M−2 group-carrying isolates. In addition and due to their low prevalence, all positive isolates for bla SHV family and combinations of bla CMY−2−like -bla SHV family and bla CMY−2−likebla CTX−M−2 group were included. An overview of the selection of isolates is shown in Table 1.
Identification of plasmids associated to the ESBL/pAmpC genes was based on co-localization on the same contig as resulted from ResFinder and PlasmidFinder analysis. Since this was not possible for all genomes, transformation of plasmids harboring each of the ESBL/pAmpC variants identified with ResFinder was performed together with selective culturing as described before (Castellanos et al., 2017). This was done to obtain the plasmid types identified with PlasmidFinder in transformed Escherichia coli DH10B harboring the different ESBL/pAmpC variants identified. Afterwards, the sequences of the transformed plasmids were used as reference to map against the genomes of the selected strains. To this purpose, transformants were sequenced with Illumina MiSeq and NextSeq sequencing as described above. Chromosomal contigs in transformants were detected and removed by mapping against the E. coli DH10B genome sequence (GenBank accession number: CP000948.1) using BLAST. The remaining contigs were considered to be part of the ESBL/pAmpC carrying plasmid and were used as a reference. Next, the contigs of the initially sequenced selected strains were aligned to the obtained reference plasmid sequences using MUSCLE (Edgar, 2004). Alignments were made on selected strains according to the ESBL/pAmpC gene variants they harbored. Resulting aligned contigs were selected and considered as newly inferred ESBL/pAmpC-carrying plasmids.

In silico MLST and Serotype Prediction
To determine the population structure of the selected ESCresistant Salmonella, 7-gene MLST at the strain level was performed in silico with MLST 1.8 (Larsen et al., 2012). Serotype was predicted using the Salmonella In Silico Typing Resource (SISTR) (Yoshida et al., 2016). Whole genome phylogenetic analyses were performed for S. Heidelberg ST15 and S. Paratyphi B var. Java ST28. Given the limited number of isolates (n = 1), this was not performed for S. Enteritidis ST11, S. Kentucky ST152 or S. Albany ST292.

Collection of Publicly Available Genomes for Phylogenetic Comparisons
Genome sequences of S. Paratyphi B var. Java ST28, were downloaded from Enterobase 1 (last accessed: 12-Sept-2017) for comparison. Likewise, sequences of S. Heidelberg ST15 were obtained from Enterobase 1 (last accessed: 11-Jan-2017). Only genomes with data available for their year of isolation, country of origin and source were considered for S. Heidelberg ST15. For both S. Paratyphi B var. Java ST28 and S. Heidelberg ST15, the genomes were collected disregarding their susceptibility to 3 rd generation cephalosporin's or to other antimicrobials. Additionally, the quality of genomes obtained in this study and the downloaded genomes was assessed with CheckM (Parks et al., 2015). Only genomes with >98% completeness score, when compared against the set of genomic markers for S. enterica, were included. MLST designation was amended using a custom BLAST-based tool 2 .

Core Genome Phylogenetic Analysis
Whole genome analysis was performed by a core-genome alignment using Parsnp v1.2 (Treangen et al., 2014). Recombination regions in the core genome alignment were detected and filtered using Gubbins (Croucher et al., 2015). Phylogenetic maximum likelihood (mid-point rooted) trees were constructed with the recombination-filtered core genomes alignments using FastTree2 3 (Price et al., 2010) and visualized with FigTree 4 .

Data Availability
The obtained genome sequences of the Salmonella strains selected for WGS ( Table 2) and those of the transformed E. coli DH10B strain with the ESBL/pAmpC-carrying plasmids used as reference have been deposited in the short-read archive of the ENA under Project Number: PRJEB23610.

PCR Screening of ESBL/pAmpC Genes
After PCR screening of the 260 isolates, 235 isolates were positive for the genes screened for. 168 were positive for bla CMY−2−like , 52 for bla CTX−M−2 group, 7 for bla SHV family, 5 harbored a combination of bla CMY−2−like -bla SHV family and 3 a combination of bla CMY−2−like -bla CTX−M−2 group. In 48 isolates, bla TEM was co-located with bla CMY−2−like, bla CTX−M−2 group or bla SHV . In 25 isolates no ESBL/pAmpC genes that were screened for were encountered.
The distribution of positive samples for Salmonella according to their source, year and ESBL/pAmpC genes is shown in Table 1.

Selection of Strains for WGS
After PCR characterization, a random selection of bla CMY−2−likepositive isolates (n = 13) and bla CTX−M−2 group (n = 8) together with all positive isolates for bla SHV family (n = 7), bla CMY−2−likebla SHV family (n = 5) and bla CMY−2−like -bla CTX−M−2 group (n = 3) were included for WGS. In total, 36 isolates were selected and subjected to WGS. The list of selected isolates and the results of the characterization of ESBL/pAmpC genes, plasmid types, serotypes, and strain MLST based on WGS is shown in Table 2.

In silico Characterization of ESBL/pAmpC-Carrying Plasmids
Co-localization of ESBL/pAmpC and plasmid replicon genes in the same contig was observed for 17 out of 36 strains selected for WGS. For the remaining strains, co-localization was determined by analyzing selected contigs harboring ESBL/pAmpC and plasmid replicon genes using MUSCLE. Genes conferring resistance to other antimicrobials only co-localized with bla CTX−M−165 -harboring contigs, not with other ESBL/pAmpC genes (Supplementary Table 1). In detail, one transformant was obtained carrying either a bla CMY−2 -, bla CTX−M−165or bla SHV−12 -harboring plasmid. The plasmids isolated from strains UGBOG4 (bla CMY−2 ), UGBAR1160 (bla CTX−M−165 ) and UGBOG327 (bla SHV−12 ) (Supplementary Table 1) were used as a reference to map against the genomes of all selected strains. The bla SHV−129 gene present in five selected strains was not transferable by transformation from any of the strains, suggesting chromosomal localization of this gene. After characterizing the plasmid contigs of the selected strains, 18 of 21 bla CMY−2 -carrying plasmids were found to belong to IncI1/ST12, 1 was non-typeable based on the PCR Based Replicon Typing (PBRT) scheme used in PlasmidFinder (Carattoli et al., 2014). Two were designated IncI1, but the plasmid from strain UGVAL515 lacked the pilL allele and the plasmid from FBOG7 lacked the sogS allele. Nonetheless, these two plasmids remained single-allele variants of IncI1/ST12. For bla CTX−M−165 , 9 of 11 plasmids harbored the IncQ1 plasmidreplicon and 2 remained non-typeable. Four of seven bla SHV−12 plasmids belonged to IncI1/ST231 and three remained nontypeable ( Table 2 and Supplementary Table 1).

Strain MLST, Serotype Characterization, and Core Genome Phylogeny
After using 7-gene MLST and the Salmonella In Silico Typing Resource (SISTR), 17 strains belonged to S. Heidelberg ST15, 16 to S. Paratyphi B var. Java ST28, 1 to S. Enteritidis ST11, 1 to S. Kentucky ST152 and 1 to S. Albany ST292 ( Table 2). Further, whole genome analysis was performed for ST28 and ST15 isolates. For the phylogenetic analysis, additional genomes for ST28 (n = 60) and ST15 (n = 1221) were selected from Enterobase.warwick.ac.uk disregarding their characteristics of susceptibility to 3rd generation cephalosporins and used to construct the phylogenetic maximum likelihood trees. All Colombian genomes belonging to ST28 and ST15 formed a single cluster in the phylogenetic analysis. Phylogenetic trees for ST28 and ST15 with data regarding the source, year and country of the strains, are shown in Figures 1 and 2 respectively. No clustering related to the presence or absence of an ESBL/pAmpC gene is observed in Figure 1, suggesting the observed clustering is related to the geographical origin of S. Paratyphi B var. Java ST28 strains and not to the presence of an AMR gene. Likewise, in Figure 2, a cluster of S. Heidelberg ST15 strains originating from Colombian poultry is observed. Furthermore, the Colombian strains from ST28 and ST15 were disseminated in multiple departments within the country. A map of Colombia with the location of origin of ST28 and ST15 isolates selected for WGS is available in Supplementary Figure 2. An additional table with the metadata of strains selected for the construction of the ST28 and ST15 phylogenies is available as Supplementary Table 2.

DISCUSSION
In summary, bla CMY−2 , bla CTX−M−165 , bla SHV−12, and bla SHV−129 are described as the most prevalent ESBL/pAmpC genes conferring resistance to ESC in S. enterica isolated from the Colombian poultry chain between 2008 and 2013. According to the objectives of the present study, the collection of isolates served to reflect maximum diversity of ESBL/pAmpC genes in different years, departments and levels of the poultry production chain in Colombia.
The finding of bla CMY−2 as the main cause of ESC resistance is comparable to a previous report of ESBL/pAmpC-producing E. coli in Colombian poultry. In that study, it was encountered a prominent association of this gene with IncI1/ST12 plasmids (Castellanos et al., 2017). Those results suggested occurrence of horizontal gene transfer of this plasmid lineage between heterogeneous E. coli STs. The results from the present study, suggest that transfer of bla CMY−2 -IncI1/ST12 plasmids between E. coli and Salmonella is also likely to occur, and could be considered a driver of the frequent occurrence of this resistance gene along the poultry chain. This particular gene-plasmid association has been described in E. coli from Brazil (da Silva et al., 2017), Salmonella from the USA (Folster et al., 2010) and E. coli and Salmonella from Europe (Accogli et al., 2013;Smith et al., 2015), suggesting an epidemiological link between the presence of bla CMY−2 -IncI1/ST12 and poultry from different countries. Nevertheless, WGS of these plasmids is necessary to assess the level of genetic relatedness among them and estimate the potential transmission of these plasmids between E. coli and Salmonella.
After PCR screening, bla CTX−M−2 group was found to be the most prevalent ESBL gene among Colombian isolates. After subtyping a selection of these isolates with WGS, these genes were found to be bla CTX−M−165 . To date, this variant has been solely reported in an isolate of Klebsiella pneumoniae from a urine sample in Chile and reported in 2016 (Accession number: KP727572) without further epidemiological records. Noticeably, all bla CTX−M−165 -positive strains, alone or in combination with other bla genes, were identified in S. Heidelberg ST15. In our collection of isolates comprising the years 2008 to 2013, bla CTX−M−165 is only detected from the year 2010 onwards ( Table 1). Not taking into account the potential bias that could occur by having isolates comprising a period of time no longer than 5 years, this finding may suggest a recent introduction of this gene in Colombian poultry, and until 2013, is limited to S. Heidelberg ST15. However, analysis of recently collected isolates of S. enterica and other Enterobacteriaceae is necessary to confirm this hypothesis.
After electroporation experiments, bla SHV−129 remained nontransferable and no plasmid markers were identified in its harboring contigs, which ranged in size between 2100 and 8594 bp. Therefore, it is likely that this gene is chromosomally located. Furthermore, the gene was found in two different serovars, S. Paratyphi B var. Java ST28, and S. Heidelberg ST15. It can be hypothesized that its transfer could be associated to an integrative or transposable element. Initial screening of transposases and Insertion Sequences (IS) using BLAST FIGURE 1 | Core-genome phylogeny of S. Java ST28 using the genomes obtained in the present study and those available in Enterobase 1 (last accessed: 12-Sept-2017).
FIGURE 2 | Core-genome phylogeny of S. Heidelberg ST15 using genomes obtained in this study and those with complete metadata selected from Enterobase 1 (last accessed: 11-Jan-2017). (Siguier, 2006) on the contigs harboring these genes detected several IS families flanking the sequences inside the contigs. Nevertheless, given the restricted size of the contigs no definite association with a unique IS element (e.g., AMR-associated IS) was possible. In such a case, complementing the short-read WGS data with additional data obtained through long-read sequencing is necessary to confirm the chromosomal location of this gene and its association to a particular mobile genetic element. This approach could also be used for further characterization of ESBL/AmpC-plasmids that were non-typeable according to the PBRT scheme used in PlasmidFinder, which could have also been affected by the limitations of genome and plasmid assembly of short-reads.
As mentioned previously, ESC and non-ESC S. Paratyphi B var. Java and S. Heidelberg were reported as the most prevalent serovars in the Colombian poultry chain (Donado-Godoy et al., 2012b. In the present study, investigation of resistance to cefotaxime showed a total of 235 (41%) resistant isolates from which, 17 (61%) originated from farms, 82 (59%) from slaughterhouses and 136 (33%) from retail (Table 1). Noteworthy, the prevalence of resistance diminishes from one level of production to the other. From previous studies, S. Paratyphi B var. Java and S. Heidelberg were the most prevalent serovars encountered at the farm level (Donado-Godoy et al., 2012b) and a larger diversity of serovars was found at retail, with more than 10 different serovars isolated repeatedly (Donado-Godoy et al., 2014). As observed in the present study, strains belonging to S. Paratyphi B var. Java and S. Heidelberg had a higher prevalence of ESC-resistance in comparison to the other serovars. These results indicate that the higher prevalence of S. Paratyphi B var. Java and S. Heidelberg, accounted in large part for the higher prevalence of ESC-resistance at the farm level and the presence of different serovars resulted in the reduction of resistance along the production chain, which is reflected in the lower prevalence of ESC-resistance at retail.
As anticipated, the analysis of Salmonella strains using MLST in addition to the resolution provided by WGS data has proven to be very useful in showing clustering of Colombian strains belonging to S. Paratyphi B var. Java ST28. According to the phylogenetic analysis including ESBL/pAmpC-positive and-negative strains (Figure 1), the clustering seems to occur independently of ESC-susceptibility and may be related to the geographical origin of the strains. Whether the cluster of Colombian S. Paratyphi B var. Java ST28 represents a particular separate lineage circulating in the country, or is present elsewhere in Latin America is a question that requires further investigation. At the time of publication of the present study no genomes of ST28 from other Latin American countries were publicly available and the comparisons within the region were limited to the strains we sequenced and analyzed.
In conclusion, resistance to ESC in S. enterica from Colombian poultry is mainly caused by bla CMY−2 and bla CTX−M−165 genes. These genes are mostly associated with IncI1/ST12 and IncQ1 plasmids, respectively. The resolution provided by WGS was appropriate to assess the evolutionary divergence of strains from Colombian poultry belonging to S. Paratyphi B var. Java ST28. Dissemination of ESBL/pAmpC genes in Salmonella is mainly due to the carriage of plasmids encoding these genes in strains belonging to S. Paratyphi B var. Java ST28 and S. Heidelberg ST15.