WGS-Based Lineage and Antimicrobial Resistance Pattern of Salmonella Typhimurium Isolated during 2000–2017 in Peru

Abstract

Salmonella Typhimurium is associated with foodborne diseases worldwide, including in Peru, and its emerging antibiotic resistance (AMR) is now a global public health problem. Therefore, country-specific monitoring of the AMR emergence is vital to control this pathogen, and in these aspects, whole genome sequence (WGS)—based approaches are better than gene-based analyses. Here, we performed the antimicrobial susceptibility test for ten widely used antibiotics and WGS-based various analyses of 90 S. Typhimurium isolates (human, animal, and environment) from 14 cities of Peru isolated from 2000 to 2017 to understand the lineage and antimicrobial resistance pattern of this pathogen in Peru. Our results suggest that the Peruvian isolates are of Typhimurium serovar and predominantly belong to sequence type ST19. Genomic diversity analyses indicate an open pan-genome, and at least ten lineages are circulating in Peru. A total of 48.8% and 31.0% of isolates are phenotypically and genotypically resistant to at least one antibiotic, while 12.0% are multi-drug resistant (MDR). Genotype–phenotype correlations for ten tested drugs show >80% accuracy, and >90% specificity. Sensitivity above 90% was only achieved for ciprofloxacin and ceftazidime. Two lineages exhibit the majority of the MDR isolates. A total of 63 different AMR genes are detected, of which 30 are found in 17 different plasmids. Transmissible plasmids such as lncI-gamma/k, IncI1-I(Alpha), Col(pHAD28), IncFIB, IncHI2, and lncI2 that carry AMR genes associated with third-generation antibiotics are also identified. Finally, three new non-synonymous single nucleotide variations (SNVs) for nalidixic acid and eight new SNVs for nitrofurantoin resistance are predicted using genome-wide association studies, comparative genomics, and functional annotation. Our analysis provides for the first time the WGS-based details of the circulating S. Typhimurium lineages and their antimicrobial resistance pattern in Peru.

1. Introduction

Salmonella Typhimurium, and other non-typhoidal Salmonella (NTS), are responsible for foodborne illnesses worldwide [1]. NTS can cause gastrointestinal disease, progressing to systemic infections in some patients. The World Health Organization (WHO) recently declared Salmonella a high-priority pathogen due to increased resistance to first-line antibiotics, fluoroquinolones, and third-generation cephalosporins [2]. In Peru, 234 foodborne outbreaks were reported between 2014 and 2018. In 2019, 1,204,136 diarrheal cases were reported mainly due to contaminated drinking water and food [3,4]. S. enterica serovars, mainly Infantis, Enteritidis, and Typhimurium, were associated with gastroenteritis in Lima hospitals during 2008–2013 and 2015–2017 [5,6]. Additionally, multi-drug resistant (MDR) Salmonella Infantis isolates resistant to first-line antibiotics, third-generation cephalosporins, and ciprofloxacin antibiotics are highly prevalent in Peru [7,8]. Unfortunately, limited studies are using whole-genome sequencing (WGS) strategies to understand diversity, alignment with antimicrobial resistance (AMR) phenotype and the AMR prevalence of circulating Salmonella.

S. Typhimurium has a broad host range and emerging dominant MDR phenotypes. For example, the MDR DT104 group disseminated rapidly globally [9], and ST313, an MDR group, is responsible for invasive diseases in Africa [10]. Apart from chromosomal mobile AMR genes mobile genetic elements, such as plasmids and pathogenic islands, are essential in expanding AMR distribution among the population and are often associated with hospital-acquired infections and foodborne outbreaks [8,9,10,11,12]. MDR Salmonella isolated in various countries largely contain genes for β-lactam, tetracycline, aminoglycoside, and quinolone resistance on plasmids [13,14].

The use of WGS in molecular epidemiology and AMR surveillance has several advantages as compared to conventional PCR, other molecular methods, or phenotypic approaches [15,16,17]. Since S. Typhimurium exhibits a broad and diverse host range, pathogenicity, and risk to human health [18,19,20]; WGS-based comparative and phylogenetic analysis as part of AMR surveillance is a promising approach to rapidly predict resistances that is much faster than phenotypic methods. Reports indicate that WGS-based approaches to predicting antimicrobial susceptibility with a good correlation between genotype and phenotype, detecting and tracing outbreaks, and determining the complement of AMR determinants are essential resources in the appropriate selection of antibiotic treatment in S. Typhimurium infections [8,21,22,23,24]. Furthermore, this approach allows the discovery of new AMR genes, or alleles of known AMR genes, as reported in various pathogenic bacteria [12,25,26]. Comparative genomics and genome-wide association studies (GWAS) have identified these potential causal variants associated with virulence and with AMR in multiple organisms [26,27].

In this study, we used whole genome sequences of S. Typhimurium isolates from a Peruvian surveillance study to determine its Peru-specific lineage and antimicrobial resistance pattern. More specifically, we classified the S. Typhimurium serovars using serotyping, multi-locus sequence typing, and average nucleotide identity (ANI) analysis approaches. Further, the nucleotide diversity, phylogenetic, pan-genome, and population structure analyses were carried out to determine the genome diversity of the circulating Peruvian S. Typhimurium isolates. Finally, identification and characterisation of AMR genes, GWAS-based prediction of new non-synonymous single nucleotide variation (SNVs) for AMR, prediction of AMR effect of newly identified SNVs, and genotype–phenotype correlations were performed for antimicrobial resistome profiling of these isolates. The overview of the methods applied, and the objectives of this study are represented in Figure 1.

2. Materials and Methods

2.1. Samples and Collection Sites

We examined 90 pathogenic S. Typhimurium isolated from humans (n = 78), animals (n = 3) and the environment (n = 9). The isolates were collected between 2000 and 2017 from 14 cities in Peru, as shown in the global map (https://glenjasper.github.io/leaflet-salmonella).

2.2. Antibiotic Susceptibility Test

The disk diffusion method [28] was used to determine AMR following Committee for Clinical Laboratory Standards (CLSI) guidelines for ten commonly used antibiotics: ampicillin (AM), chloramphenicol (C), ciprofloxacin (CIP), trimethoprim-sulfamethoxazole (SXT), cefotaxime (CTX), nalidixic acid (NA), amoxicillin-clavulanate (AMC), nitrofurantoin (N), tetracycline (TE), and ceftazidime (CAZ). MDR was defined as resistance to at least three antibiotics. The intermediate or reduced susceptibility phenotype was considered a susceptible isolate (Supplementary Table S1).

2.3. Whole Genome Sequencing

We used Peruvian Salmonella Typhimurium genome (90 isolates) sequence raw reads generated under the 10K Salmonella Project (BioProject: PRJEB35182). After assembly and annotation, we submitted these genomes to the BioProject (PRJNA635403) as “Peruvian Salmonella spp. Genome sequencing and assembly”. The genome sequencing was performed within the 10K Salmonella Project as described by Perez-Sepulveda et al., 2021 [29]. Additionally, a global set of 50 genomes were obtained from the Genbank database and included in the further analysis to estimate the diversity and ancestry of Peruvian isolates globally (Supplementary Table S2).

2.4. Genome Assembly, Annotation, and Plasmid Detection

The raw fastq sequences from the Illumina 150 bp paired-end were checked for quality using FastQC v0.11.8 [30]. De novo genome assembly was performed using Unicycler v0.4.5 [31], and the quality of the assembly was evaluated using Quast v3.2 [32] and BUSCO v4.0.6 [33]. Contigs < 200 bp were removed, and the assembled genomes were submitted under BioProject: PRJNA635403. For further analysis, the annotation was performed using PROKKA v1.11 [34], and only the genome sequences with >20× depth coverage were considered. We considered 20× to be the minimum coverage to include the maximum number of our isolates (n = 90 out of a total n = 109) for core genome analysis, variant determination, and monitoring of infectious outbreaks in other previous studies [35,36,37]. Plasmids were predicted and reconstructed from the assembled genomes using MOB-suite [38] and classified as conjugative, mobilisable, and non-mobilisable plasmids. PlasmidFinder [39] and ABRicate (https://github.com/tseemann/abricate; accessed on 4 April 2022) were also used to crosscheck the MOB-suite results.

2.5. In Silico Serotyping, MLST, and ANI Analysis

The serovar and sequence types were predicted using SeqSero [40] and MLST (https://github.com/tseemann/mlst; accessed on 5 April 2022) [41], respectively. Pairwise average nucleotide identity (ANI) values were calculated using FastANI v.1.1 [42] to determine the degree of genomic relatedness. Results were visualised using the ggplot2 V.3.3.5 (https://ggplot2.tidyverse.org/reference/index.html; accessed on 5 April 2022). The respective tools’ assembled genomes and default parameters were used for these analyses.

2.6. Pan-Genome and Phylogenetic Analysis

Core and accessory genes were identified using Roary v3.12 [43] with default settings. The R package micropan v.2.1 [44] was used to model the openness of the pan-genome using Heaps’ law as described by Tettelin et al., [45] with the number of permutations set to 1000. Values α ≤ 1 representant an open pan-genome, where adding new genomes will increase the pan-genome substantially. For this analysis, we used our Peruvian samples and 50 S. Typhimurium isolates from 29 countries distributed on six continents (Asia, Europe, Africa, Australia, and North and South America). We used S. Typhimurium (LT2, GenBank: NC_003197) as the reference. Core genome-SNPs were predicted using Snp-sites [46] from the Roary v3.12 [43] core genome output. Phylogenetic analysis was based on core genome genes of the Peruvian and 140 global samples. The RaxML v8.2.12 [47], maximum likelihood method, GTR + Gamma model, and 1000 bootstrap replicates were applied to create the Phylogenetic trees that were visualised using the ggtree package in R (R Development Core Team, 2016).

2.7. Population Structure and Diversity Analyses

The population structure for two sample sets was determined separately using core gene SNPs: (i) Peruvian and (ii) Peruvian + 50 global samples. The Bayesian analysis of population structure (BAPS) [48] model was defined using RhierBAPS v1.0.1 [49]. The nucleotide variation analysis within the Peruvian and global S. Typhimurium populations was calculated using the pairwise similarity (inverse diversity calculation) and median pairwise similarities using the core genome SNPs and the MEGA-X tool [50]. Ggplot2 package in R (R Development Core Team, 2016) was used to visualise the results. These analyses allow estimating the diversity and ancestry of Peruvian isolates from global representative isolates.

2.8. Identification of Known Antimicrobial Resistance Genes and Single Nucleotide Variations

For each genome, the AMR-associated genes were identified in the chromosome and plasmids using ABRicate (https://github.com/tseemann/abricate; accessed on 4 April 2022) using CARD [51] and ARG-ANNOT [52] databases. More than 85.0% of sequence coverage and identity were considered the lower limit. Additionally, AMR Finder [53] and ResFinder [54], contain a large number of Salmonella spp. sequences, and were included to corroborate the results. The AMR genes were classified according to resistance mechanisms and drug class using the CARD database and manual curation [51].

To identify the SNPs in the AMR genes from the chromosomal DNA, first, we extracted the AMR gene sequences using Roary [43], followed by analysis with MAFFT v7.307 [55] for sequence alignment. Finally, SNP-sites [46] was used for the variant calling. Subsequently, we used PointFinder [56] to identify alleles and associate them with AMR.

2.9. Genotype–Phenotype Correlation

We performed the genotype–phenotype correlation (GPC) analysis using disc diffusion antimicrobial susceptibility test results in combination with the WGS-informed AMR analysis. These assays were compared with the occurrence of known AMR genes and alleles associated with the resistance to the respective drug. False-positive, false-negative, sensitivity, specificity, and accuracy of the GPC were calculated as described previously [21]. Finally, phenotypically resistant isolates that did not contain known genes for this resistance were identified as candidates for use in finding new variants that may confer resistance.

2.10. Genome-Wide Association Study to Identify New Single Nucleotide Variations Associated with Resistance Phenotype

In isolates where genotype–phenotype association did not match, we used the genome-wide association study (GWAS) approach to identify possible genes and alleles for specific AMR. A similar approach, as described by Bandoy and Weimer [57], was used for the GWAS analysis with a chi-square test in R (R Development Core Team, 2016) to identify mutant alleles that conferred a specific phenotype. In this process, the variant calling was performed for all the isolates using LT2 isolate as the reference. Only the phenotype observed for the ten tested drugs was considered. The Snippy v3.2 (https://github.com/tseemann/snippy; accessed on 16 April 2022) was used for variant calling, and the SNPs and indels were filtered considering a minimum sequencing depth > 20x and minor allele frequency less than 0.02 were removed. A significant association was considered when p < 0.05, and the values were automatically corrected for multiple testing using the Bonferroni method [58]. The results were visualised with the Manhattan plot in R (R Development Core Team, 2016) using the “qqman” package. To minimise the false positive association in the GWAS analysis, we applied the population structure as a covariable in Firth’s logistic regression analysis using the “logistf” v1.24 package in R (https://github.com/vicbp1/Genetic-Arquitecture-of-Zika; accessed on 16 April 2022). The population structure was predicted using principal component analysis (PCA) and PLINK v1.90b6.9 [59], considering the first six principal components (PC1–PC6) as continuous covariables.

2.11. Association of New Single Nucleotide Variations and Drug Resistance

The GWAS represent a powerful approach to identifying new genetic variants in isolates that demonstrated phenotypic drug resistance but did not contain determinants associated with known antimicrobials. We examined the potential for additional gene variants that may explain the specific drug resistance in those isolates. When found, we adopted the strategy described by Ferla et al. [57,60]. Additionally, we modelled the 3D structure of the protein harbouring the variation using the corresponding amino acid sequence and Swiss Model homology server [61]. The 3D structure was then used to predict the stability of these proteins for the new variations using the Dynamute2 tool [62]. Finally, the stability/instability property of the new allele was correlated with the observed drug resistance phenotype.

3. Results

3.1. Genomic Characterisation Shows Peruvian Salmonella Samples Belong to Typhimurium Serovar and Sequence Type ST19

Ninety WGS samples passed the quality metrics and were used for further analysis (Supplementary Table S1). The SeqSero and MLST analysis confirmed that 83 belong to serovar Typhimurium and 85 belong to sequence type ST19, respectively (Supplementary Table S1). For the global samples, the majority were serovar Typhimurium and sequence type ST19. Use of ANI analysis found the Peruvian genomes were >99.6% identical, while the global samples (n = 140) were >99.4% identical (Supplementary Figure S1). These analyses confirm that both samples belong to serovar S. Typhimurium and use all the isolates in subsequent analyses.

3.2. Pan-Genome Structure and Nucleotide Diversity of Peruvian Samples

The pan-genome analysis examined the content of gene diversity, and the core genome was used as input data to analyse the phylogeny, population structure, and nucleotide variation within these isolates. The pan-genome of 140 S. Typhimurium was open (α = 0.53) and contained 11,168 orthologous genes with 3455 core genes and 7713 accessory genes (softcore = 495, shell = 1000, and cloud = 6218) that represent the 30.9% and 69.1%, respectively (Figure 2A,B). Specifically for the Peruvian isolates (n = 90), the core and accessory genome constitutes 42.6% and 57.4%, respectively, where there are 7462 orthologous genes, 3181 core genes and 4281 accessory genes (softcore = 889, shell = 781, and cloud = 2611) (Figure 2C,D). We identified values > 85.0% and 70.0% of pairwise nucleotides dissimilarity (Supplementary Figure S2) between the overall and within-population nucleotide diversity to be 0.057 and 0.135 for global and Peruvian samples, respectively, indicating that our Peruvian isolates are more genetically diverse and present an open pan-genome with α value minor to 1 (α = 0.73).

3.3. Ten Lineages of S. Typhimurium Classified under Two Major Clades Are Circulating in Peru

The un-rooted phylogenetic tree based on the core genome clustered similarly for both the global and Peruvian samples. The Peruvian isolates were analysed with the other global samples and contained ten distinct sub-clades (Supplementary Figure S3). It is also observed that the lineages of Peruvian isolates are distributed across the continents (Figure 3). When analysed exclusively for the 90 Peruvian isolates, the phylogenetic tree shows a similar number of sub-clades (Figure 3). Both analyses also found that the Peruvian samples are mostly distributed in two major clades. When we determined the population structure as sequence clusters (SC, assigned as lineages) by a two-level hierarchical Bayesian approach (BAPS) using the core gene SNPs, we observed at least 19 sequence clusters for the global samples and 10 sequence clusters for the 90 Peruvian samples within this global sample pool (Supplementary Figure S3). The first clade (n = 47) includes the sequences cluster SC13 (n = 24), SC14 (n = 9), SC15 (n = 7), SC16 (n = 4), and SC9 (n = 3); while the second clade (n = 43) includes SC2 (n = 6), SC3 (n = 4), SC8 (n = 6), SC17 (n = 26), and SC18 (n = 1) (Figure 3).

While we further analysed the sequence clusters only for the 90 Peruvian samples, we identified four additional sub-populations without much phylogenetic difference, where SC2 consists of three sub-sequence clusters and SC3 and SC8 had two sub-sequence clusters each (Figure 3). Therefore, ten sequence clusters or lineages are circulating in Peru. The major clonal group CG-I (SC3) consists of 24 isolates, distributed in seven different Peruvian cities, mostly isolated during 2005–2006, and closely related to a Switzerland sample. The second crucial clonal group, CG-II (SC17), consists of 26 isolates, is found in eight different Peruvian cities, isolated during 2000–2001, and is closely related to a Chilean sample (Figure 3).

3.4. The Phenotypic Profile Shows 48.9% Resistance to at least One Drug in Isolates Circulating in Peru for Ten Drugs

Out of the 90 Peruvian samples, 44 isolates (48.9%) contained resistance to at least one drug from the tested drugs using the disc diffusion assay. Out of these 44 isolates, the highest phenotypic resistance profile was for NA (n = 18, 40.9%; or 20.0% considering total 90 isolates), TE (n = 17, 38.6% or 18.9% considering total 90 isolates), AM (n = 14, 31.8%; or 15.6% considering total 90 isolates), N (n = 13, 29.6%; or 14.4% considering total 90 isolates), and 11 samples (25.0%; or 12.0% considering total 90 isolates) show MDR (Supplementary Table S1, Figure 4). These MDR isolates mostly belong to SC8 (n = 5), followed by SC9 (n = 2), SC15 (n = 2), SC17 (n = 1). In sample collection site-based analysis, Lima shows most (n = 10) of the MDR isolates and only one from Huánuco (Supplementary Table S1).

Figure 3. The phylogenetic tree was generated using a maximum likelihood method with 1000 replicates of bootstrap using GTR + GAMMA to estimate the evolutionary distance between Peruvian isolates (n = 90). The phylogenetic tree was clustered into at least two large clades and separated into nine sub-clades. Each sub-clade corresponds to a population group, except for SC18. Two emerging clades (CG-I and CG-II) are also found. Additionally, we identified subgroups based on the prediction of only the Peruvian population structure, but they did not show the phylogenetic distinction.

Figure 3. The phylogenetic tree was generated using a maximum likelihood method with 1000 replicates of bootstrap using GTR + GAMMA to estimate the evolutionary distance between Peruvian isolates (n = 90). The phylogenetic tree was clustered into at least two large clades and separated into nine sub-clades. Each sub-clade corresponds to a population group, except for SC18. Two emerging clades (CG-I and CG-II) are also found. Additionally, we identified subgroups based on the prediction of only the Peruvian population structure, but they did not show the phylogenetic distinction.

Figure 4. Individual antibiotic resistance profile of Peruvian isolates and their genotype–phenotype correlations for our ten tested antibiotics. Antibiotic abbreviation: ampicillin (AM), chloramphenicol (C), ciprofloxacin (CIP), trimethoprim-sulfamethoxazole (SXT), cefotaxime (CTX), nalidixic acid (NA), amoxicillin-clavulanate (AMC), nitrofurantoin (N), tetracycline (TE), and ceftazidime (CAZ). AMR profile abbreviation: Gen, genotypic; Phen, phenotypic; Gen-Phen, genotypic-phenotypic.

Figure 4. Individual antibiotic resistance profile of Peruvian isolates and their genotype–phenotype correlations for our ten tested antibiotics. Antibiotic abbreviation: ampicillin (AM), chloramphenicol (C), ciprofloxacin (CIP), trimethoprim-sulfamethoxazole (SXT), cefotaxime (CTX), nalidixic acid (NA), amoxicillin-clavulanate (AMC), nitrofurantoin (N), tetracycline (TE), and ceftazidime (CAZ). AMR profile abbreviation: Gen, genotypic; Phen, phenotypic; Gen-Phen, genotypic-phenotypic.

3.5. Genotypic Profile Showed 31.0% of Drug Resistance Isolates Are Circulating in Peru

The genotypic profile for 90 isolates showed the presence of a total of 63 (chromosomal + plasmid) different AMR genes as per drug class and resistance mechanisms. We found 31 were chromosomal, 32 were acquired in the 90 samples (Supplementary Table S3). However, since we have only the DNA sequence data, we over-ruled the expression-based drug resistance mechanism and have considered only the genes that confer drug resistance if it is present or when alleles were present. Based on these criteria, we found 19 isolates containing mobile AMR genes, two isolates harboured known single nucleotide variations in chromosomal AMR genes that conferred drug resistance, and seven isolates contained both the mobile AMR gene as well as variations in chromosomal AMR genes (Supplementary Table S3). Therefore, 28 isolates (31.0%) showed resistance to at least one drug gene/variation, including an untested drug.

While we considered the use of 10 antibiotics and DNA-based criteria, we found 26 isolates that have at least one AMR gene/variation in the genome. Out of these 26 isolates, 17 isolates contained mobile AMR genes, two isolates had known variations in chromosomal AMR genes, and seven isolates had both the mobile AMR gene and variations in chromosomal AMR genes (Supplementary Table S3). In this analysis, the resistance genotypic profile was observed: nalidixic acid was (n = 17, 65.4% or 18.9% considering total 90 isolates), tetracycline (n = 13, 50.0% or 14.4% considering total 90 isolates), and ampicillin (n = 10, 11.0%) (Figure 4). These MDR isolates mostly belong to SC8 (n = 5) and SC9 (n = 3) (Supplementary Table S3). The details of the mobile AMR gene and variations in chromosomal AMR genes of our samples are given in

Table 1. The genotype to phenotype concordance using accuracy, sensitivity, and specificity for Peruvian isolates for ten tested drugs.

Antibiotic	New Variations in Chromosomal Genes	Known Variations in Chromosomal Genes	Mobile Resistance Genes	Susceptible Phenotype		Resistant Phenotype		Accuracy	Sensitivity	Specificity
				Resistant Genotype	Susceptible Genotype	Resistant Genotype	Susceptible Genotype
				FP	TN	TP	FN	(TP + TN)/TOTAL	TP/(TP + FN)	TN/(TN + FP)
Tetracycline (T)	-	-	tetA (n = 15), tetB (n = 2), tetD (n = 3), tetR (n = 4)	1	72	13	5	94.4%	72.2%	98.6%
Ampicillin (AM)	ompA (n = 90), ampH (n = 46), golS (n = 90), mdsA (n = 87), mdsB (n = 90), mdsC (n = 69)	-	blaTEM-176 (n = 1), blaTEM-181 (n = 9), blaSHV-134 (n = 1), blaTEM-181 (n = 1), blaCTX-M-15 (n = 1)	5	71	5	9	84.4%	35.7%	93.4%
Amoxicillin-clavulanate (AMC)		-	blaCTX-M-15 (n = 1), blaSHV-12 (n = 1), blaSHV-134 (n = 1)	2	86	0	2	95.6%	0.0%	97.7%
Cefotaxime (CTX)		-		0	85	2	3	96.7%	40.0%	100.0%
Ceftazidime (CAZ)		-		1	88	1	0	98.9%	100.0%	98.9%
Chloramphenicol (C)	golS (n = 90), mdsA (n = 87), mdsB (n = 90), mdsC (n = 69)	-	florR (n = 5)	2	81	3	4	93.3%	42.9%	97.6%
Ciprofloxacin (CIP)	mdtK (n = 90), crp (n = 90), emrA (n = 90), emrB (n = 87) [p.G509D, G510D (n = 1)], emrR (n = 90), gyrB (n = 89) p.S347P (n = 1), parC (n = 89) [p.A554T (n = 3), p.R360H, R365H (n = 1), p.T571S (n = 1)]	gyrA (n = 90) p.S83Y (n = 6), p.S83F (n = 1), p.D87G (n = 1), p.D87Y (n = 1)	qnrB5 (n = 8), qnrE2 (n = 3), qnrB19 (n = 2)	3	85	2	0	96.7%	100.0%	96.6%
Nalidixic acid (NA)				3	69	14	4	92.2%	77.8%	95.8%
Trimethoprim-sulfamethoxazole (STX)	-	-	dfr1 (n = 5), dfr12 (n = 2), dfr14 (n = 1), sul2 (n = 4), sul3 (n = 5)	2	84	3	3	96.7%	50.0%	97.7%
Nitrofurantoin (N)	nfsA (n = 90), nfsB (n = 88), ribE (n = 89)	-	-	0	77	0	13	85.6%	0.0%	100.0%
Aminoglycoside	aac(6′)-Iaa (n = 90), kdpE (n = 74) [p.A115G, A115E (n = 2), p.L39P (n = 3), p.R81H, S100R (n = 1)]	-	aac(3) IIe (n = 6), aph(3′)-Ia (n = 4), aac(3)-Iie (n = 6), aac(6′)-Ian (n = 2), ant(3″)-Iia (n = 1), aph(3″)-Ib (n = 7), aph(6)-Id (n = 6)
Aminocoumarin and Aminoglycoside	baeR (n = 90), cpxA (n = 86)	-	-
aminocoumarin	mdtB (n = 84) [p.T69A (n = 5), p.S157A (n = 5), p.N199H (n = 1), p.R512H (n = 3), p.Q315R (n = 1), p.R590H (n = 5)], mdtC (n = 90) [p.T81S (n = 1), p.N113K (n = 1), p.N133D (n = 1)]	-	-
Multiclass	acrB2 (n = 90) [p.L332F (n = 3), p.V482A, A491T (n = 1)], sdiA (n = 90), tolC (n = 89), H-NS (n = 90), marA (n = 90), acrB1 (n = 90) [p.T599P (n = 27), p.R418H (n = 4), p.L845F (n = 1)], mdtC, marA	-	kpnH (n = 3)
Bacitracin	bacA (n = 89)	-	-
Nitroimidazole	msbA (n = 90)	-	-
Microcin	yojI (n = 90) [p.H431Y (n = 6), p.A366D (n = 2)]	-	-
Lincosamide	-	-	linG (n = 3)
Bleomycin	-	-	BLMT (n = 1)
Fosfomycin	-	-	fosA3 (n = 1)
Colistin	-	-	mcr-1 (n = 3)

and Table S3.

3.6. Ciprofloxacin and Ceftazidime Resistance Shows the Best Genotype–Phenotype Correlation, but Nitrofurantoin Does Not

While we compared our genotype-based drug resistance isolates to the phenotype data for ten tested antibiotics, we observed that 22 isolates were phenotypically resistant; they did not contain a known AMR plasmid or mobile gene or any mutation in the chromosomal AMR gene. On the other hand, out of the 28 genotypically resistant isolates, four isolates did not demonstrate phenotypically correlation to NA (n = 2), multi-drug resistance (n = 2) and additional drugs not testing (n = 2). Therefore, we considered a total of 50 isolates (44 phenotypically and 4 genotypically resistant to 10 tested drugs, and 2 genotypically resistant to untested antibiotics) for our further analysis (Table 2).

While we considered the genotype and phenotype data for ten drugs, 22 (n = 44–22) isolates had no correlation between the genotype and phenotype. In these 22 phenotypically drug-resistant isolates, ten isolates showed resistance exclusively to N, 3 to AM, 2 to CTX, and one each to NA, TE, AMC, AM + N, AM+ CTX+ C, AM + C+ TE, and AM+ SXT (Table 2, Supplementary Table S3). According to our calculation, as described in the method, the genotype–phenotype correlations for the ten tested antibiotics, the accuracy was 84.4% to 98.9%, specificity is between 93.4% and 100.0%, and the sensitivity reached up to 100.0% only for ciprofloxacin and ceftazidime resistance (Table 1). Ciprofloxacin and ceftazidime show >96.0% accuracy, specificity, sensitivity, and the lowest values were for beta-lactam resistance, followed by chloramphenicol and trimethoprim-sulfamethoxazole. In addition to that, we did not obtain a good sensitivity for nitrofurantoin as most of the phenotypical nitrofurantoin resistance isolates do not have any known nitrofurantoin resistance marker (Table 1).

3.7. Seventeen Different Plasmids Carrying 30 AMR Genes Were Identified in Peruvian Isolates

We identified a total of 47 different plasmids in the 90 Peruvian isolates of which 30 plasmids did not carry any AMR gene, while 17 contained at least one AMR gene. A total of 30 AMR genes were found in the 17 plasmids from 28 isolates. A maximum occurrence of eight AMR genes was found in one plasmid (lncHI2A family, isolate-FD01846422) (Supplementary Table S4). These isolates mainly belong to SC9, SC2 and SC8 lineages. Among the 90 isolates, 83 isolates had the lncFIB virulence plasmid, and this lncFIB plasmid was observed in all the sequencing clusters except the isolates that belong to SC18 and one isolate under SC2 (Supplementary Table S1, Figure 5). The details of the isolates and their corresponding plasmid AMR genes are given in Table 3 and Table S4.

In addition, we detected plasmid-mediated resistance to third-generation cephalosporin, lncI-gamma/k1_P7 plasmid (n = 1/SC8) carrying blaCTX-M-15 gene and IncI1-I(Alpha)_P14 plasmid (n = 1/SC3), blaSHV-12 and blaSHV-134 genes (Supplementary Table S4). Eleven isolates (n = 5/SC8, n = 3/SC2, n = 2/SC16, and n = 1/SC13) showed presence of Col(pHAD28) plasmid carrying fluoroquinolone-resistant either qnrB19 or qnrB5 gene. The fluoroquinolone-resistant qnrE2 gene was found in three isolates belonging to the SC9 group that contained the IncHI2_P1 plasmid (Supplementary Table S4). Another critically important antibiotic-resistant gene, fosA3, for fosfomycin was identified in IncFIB(pN55391)_P12 (n = 1/SC15), and mcr-1 gene for colistin resistance was found in IncHI2_P1 (n = 1/SC9), lncI2_P3 (n = 2/SC9) (Supplementary Table S4).

Table 2. Profile of 50 Peruvian drug resistance isolates (44 phenotypically and 4 genotypically resistant to 10 tested drugs, and 2 genotypically resistant to untested antibiotics).

SC	Isolate	Resistance Genes Profile	Resistance Plasmids Profile	Resistance Phenotypic Profile	Institution/Hospital	City	Year
SC9	FD01846446	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1	IncHI2_P1, IncI2_P3, P6, P8, IncFIB(S)_P11	-	Hospital Emergencias Pediátricas	Lima	2017
	FD01846422	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1, sul2	IncHI2_P1, P6	NA, C, SXT, TE	Hospital Emergencias Pediátricas	Lima	2017
	FD01872670	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1, sul2, tetR	IncHI2_P1, IncI2_P3, P6, P8, IncFIB(S)_P11	NA, C, SXT, TE	Hospital Emergencias Pediátricas	Lima	2017
SC8	FD01852492	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3″)-Ib, dfrA1, floR, qacL, sul3, qnrB5, ANT(3″)-Iia, blaTEM-176	P4, IncI-gamma/K1_P7, P6, Col(pHAD28)_P13	NA, N	LRR Apurimac	Apurimac	2015
	FD01852476	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, tetA, qnrB19, blaCTX-M-15	P4, IncI-gamma/K1_P7, P10, Col(pHAD28)_P13, P16	NA, TE, AM, CTX, CAZ	INSN	Lima	2015
	FD01852484	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, dfrA1, floR, qacL, sul3	P4, P6, IncI-gamma/K1_P7, Col(pHAD28)_P13	NA, C, SXT	Hospital Emergencias Pediátricas	Lima	2015
	FD01852637	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, tetA, qnrB5	P9, Col(pHAD28)_P13	NA, CIP, TE, AM	Hospital Emergencias Pediátricas	Lima	2015
	FD01852500	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, tetA	P4, P15, P16	NA, TE, AM	INSN	Lima	2016
	FD01852499	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, qnrB19	P4, Col(pHAD28)_P13	NA, CIP, AM, AMC	Hospital Emergencias Pediátricas	Lima	2016
SC2	FD01852647	aph(3″)-Ib, aph(6)-Id, sul2, tetA, tetR	P2	TE	INEN	Lima	2015
	FD01852582	aph(3″)-Ib, aph(6)-Id, sul2, tetA, tetR	P2, IncN_P17	TE	Hospital Emergencias Pediátricas	Lima	2013
	FD01851394	aph(3″)-Ib, aph(6)-Id, sul2, tetA	P2	TE	Hospital Emergencias Pediátricas	Lima	2013
	FD01852653	qnrB5, tetA, tetR	Col(pHAD28)_P13, IncN_P17	NA, TE	DISA Lima Ciudad	Lima	2012
	FD01852523	qnrB5, tetA	IncFIB(S)_P11, Col(pHAD28)_P13	NA, TE	INEN	Lima	2016
	FD01852539	qnrB5, tetA	IncFIB(S)_P11, Col(pHAD28)_P13	NA, TE, AM, C, SXT	Hospital Emergencias Pediátricas	Lima	2016
SC18	FD01851477	aph(3′)-Ia, dfrA14	IncFIB(pN55391)_P12	NA	DIRESA Trujillo	La Libertad	2008
SC16	FD01852549	qnrB5, gyrA p.S83F	Col(pHAD28)_P13	NA	UNMSM	Lima	2017
SC16	FD01852587	gyrA p.D87Y	-	NA	INEN	Lima	2014
SC14	FD01852461	gyrA p.D87Y	-	NA	INEN	Lima	2015
SC16	FD01852535	qnrB5	Col(pHAD28)_P13	-	Cusco	Cusco	2016
SC13	FD01852545	qnrB5	Col(pHAD28)_P13	-	INEN	Lima	2017
SC15	FD01852865	sul2, tetB, aac(6′)-Ian	IncC_P5	AM, C, TE	Hospital Dos de Mayo	Lima	2010
	FD01852857	sul2, tetB, aac(6′)-Ian	IncC_P5	STX, NA, TE, N	DIRESA Huanuco	Huanuco	2010
	FD01852558	fosA3	IncFIB(pN55391)_P12	NOT TESTED	DIRESA Callao	Callao	2013
	FD01852530	KpnH	-	SXT, NA	INEN	Lima	2016
SC13	FD01851538	KpnH	-	N	INEN	Lima	2016
SC3	FD01852748	blaSHV-12, blaSHV-134	IncI1-I(Alpha)_P14	AM, CTX	INEN	Lima	2011
	FD01852600	aph(3′)-Iia, BLMT	IncI1-I(Alpha)_P14	NOT TESTED	INEN	Lima	2014
SC17	FD01851425	-	-	NA	CENAN/INS	Lima	2007
SC16	FD01852548	-	-	TE	Hospital Emergencias Pediátricas	Lima	2017
SC13	FD01851503	-	-	N	DIRESA Trujillo	La Libertad	2005
	FD01848616	-	-	N	LIMA CIUDAD	Lima	2005
	FD01851500	-	-	N	Puno	Puno	2005
	FD01851509	-	-	N	Direccion de Salud I	Callao	2005
	FD01851516	-	-	N	LRR Chiclayo	Lambayeque	2005
	FD01851519	-	-	N	Huaraz	Huaraz	2005
	FD01851527	-	-	N	Direccion de Salud I	Callao	2006
	FD01851529	-	-	N	Huaraz	Huaraz	2006
	FD01851530	-	-	N	INEN	Lima	2006
	FD01851538	-	-	N	LIMA	Lima	2006
SC17	FD01848677	-	-	CTX	Cajamarca	Cajamarca	2000
	FD01848679	-	-	CTX	LIMA ESTE	Lima	2000
	FD01851388	-	-	AM	INEN	Lima	2010
SC15	FD01851386	-	-	AM	Hospital Emergencias Pediátricas	Lima	2010
SC13	FD01851541	-	-	AM	Hospital Santa Maria Del Socorro	Ica	2006
SC14	FD01851320	-	-	AMC	CENAN/INS	Lima	2009
SC13	FD01851499	-	-	AM, N	LRR Chiclayo	Lambayeque	2005
SC13	FD01848615	-	-	AM, SXT	Hospital San Bartolomé	Lima	2005
SC17	FD01852858	-	-	AM, C, TE	DISA Lima Ciudad	Lima	2010
SC17	FD01848690	-	-	AM, C, CTX	Hospital San Bartolomé	Lima	2001

Antibiotic abbreviation: ampicillin (AM), chloramphenicol (C), ciprofloxacin (CIP), trimethoprim-sulfamethoxazole (SXT), cefotaxime (CTX), nalidixic acid (NA), amoxicillin-clavulanate (AMC), nitrofurantoin (N), tetracycline (TE), and ceftazidime (CAZ).

Table 2. Profile of 50 Peruvian drug resistance isolates (44 phenotypically and 4 genotypically resistant to 10 tested drugs, and 2 genotypically resistant to untested antibiotics).

SC	Isolate	Resistance Genes Profile	Resistance Plasmids Profile	Resistance Phenotypic Profile	Institution/Hospital	City	Year
SC9	FD01846446	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1	IncHI2_P1, IncI2_P3, P6, P8, IncFIB(S)_P11	-	Hospital Emergencias Pediátricas	Lima	2017
	FD01846422	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1, sul2	IncHI2_P1, P6	NA, C, SXT, TE	Hospital Emergencias Pediátricas	Lima	2017
	FD01872670	tetA, tetD, sul3, linG, floR, dfrA12, blaTEM-181, qnrE2, mcr-1, aph(6)-Id, aph(3″)-Ib, aadA2, dfrA1, sul2, tetR	IncHI2_P1, IncI2_P3, P6, P8, IncFIB(S)_P11	NA, C, SXT, TE	Hospital Emergencias Pediátricas	Lima	2017
SC8	FD01852492	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3″)-Ib, dfrA1, floR, qacL, sul3, qnrB5, ANT(3″)-Iia, blaTEM-176	P4, IncI-gamma/K1_P7, P6, Col(pHAD28)_P13	NA, N	LRR Apurimac	Apurimac	2015
	FD01852476	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, tetA, qnrB19, blaCTX-M-15	P4, IncI-gamma/K1_P7, P10, Col(pHAD28)_P13, P16	NA, TE, AM, CTX, CAZ	INSN	Lima	2015
	FD01852484	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, dfrA1, floR, qacL, sul3	P4, P6, IncI-gamma/K1_P7, Col(pHAD28)_P13	NA, C, SXT	Hospital Emergencias Pediátricas	Lima	2015
	FD01852637	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, tetA, qnrB5	P9, Col(pHAD28)_P13	NA, CIP, TE, AM	Hospital Emergencias Pediátricas	Lima	2015
	FD01852500	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, aph(3′)-Ia, tetA	P4, P15, P16	NA, TE, AM	INSN	Lima	2016
	FD01852499	aac(3)-Iie, blaTEM-181, gyrA p.S83Y, qnrB19	P4, Col(pHAD28)_P13	NA, CIP, AM, AMC	Hospital Emergencias Pediátricas	Lima	2016
SC2	FD01852647	aph(3″)-Ib, aph(6)-Id, sul2, tetA, tetR	P2	TE	INEN	Lima	2015
	FD01852582	aph(3″)-Ib, aph(6)-Id, sul2, tetA, tetR	P2, IncN_P17	TE	Hospital Emergencias Pediátricas	Lima	2013
	FD01851394	aph(3″)-Ib, aph(6)-Id, sul2, tetA	P2	TE	Hospital Emergencias Pediátricas	Lima	2013
	FD01852653	qnrB5, tetA, tetR	Col(pHAD28)_P13, IncN_P17	NA, TE	DISA Lima Ciudad	Lima	2012
	FD01852523	qnrB5, tetA	IncFIB(S)_P11, Col(pHAD28)_P13	NA, TE	INEN	Lima	2016
	FD01852539	qnrB5, tetA	IncFIB(S)_P11, Col(pHAD28)_P13	NA, TE, AM, C, SXT	Hospital Emergencias Pediátricas	Lima	2016
SC18	FD01851477	aph(3′)-Ia, dfrA14	IncFIB(pN55391)_P12	NA	DIRESA Trujillo	La Libertad	2008
SC16	FD01852549	qnrB5, gyrA p.S83F	Col(pHAD28)_P13	NA	UNMSM	Lima	2017
SC16	FD01852587	gyrA p.D87Y	-	NA	INEN	Lima	2014
SC14	FD01852461	gyrA p.D87Y	-	NA	INEN	Lima	2015
SC16	FD01852535	qnrB5	Col(pHAD28)_P13	-	Cusco	Cusco	2016
SC13	FD01852545	qnrB5	Col(pHAD28)_P13	-	INEN	Lima	2017
SC15	FD01852865	sul2, tetB, aac(6′)-Ian	IncC_P5	AM, C, TE	Hospital Dos de Mayo	Lima	2010
	FD01852857	sul2, tetB, aac(6′)-Ian	IncC_P5	STX, NA, TE, N	DIRESA Huanuco	Huanuco	2010
	FD01852558	fosA3	IncFIB(pN55391)_P12	NOT TESTED	DIRESA Callao	Callao	2013
	FD01852530	KpnH	-	SXT, NA	INEN	Lima	2016
SC13	FD01851538	KpnH	-	N	INEN	Lima	2016
SC3	FD01852748	blaSHV-12, blaSHV-134	IncI1-I(Alpha)_P14	AM, CTX	INEN	Lima	2011
	FD01852600	aph(3′)-Iia, BLMT	IncI1-I(Alpha)_P14	NOT TESTED	INEN	Lima	2014
SC17	FD01851425	-	-	NA	CENAN/INS	Lima	2007
SC16	FD01852548	-	-	TE	Hospital Emergencias Pediátricas	Lima	2017
SC13	FD01851503	-	-	N	DIRESA Trujillo	La Libertad	2005
	FD01848616	-	-	N	LIMA CIUDAD	Lima	2005
	FD01851500	-	-	N	Puno	Puno	2005
	FD01851509	-	-	N	Direccion de Salud I	Callao	2005
	FD01851516	-	-	N	LRR Chiclayo	Lambayeque	2005
	FD01851519	-	-	N	Huaraz	Huaraz	2005
	FD01851527	-	-	N	Direccion de Salud I	Callao	2006
	FD01851529	-	-	N	Huaraz	Huaraz	2006
	FD01851530	-	-	N	INEN	Lima	2006
	FD01851538	-	-	N	LIMA	Lima	2006
SC17	FD01848677	-	-	CTX	Cajamarca	Cajamarca	2000
	FD01848679	-	-	CTX	LIMA ESTE	Lima	2000
	FD01851388	-	-	AM	INEN	Lima	2010
SC15	FD01851386	-	-	AM	Hospital Emergencias Pediátricas	Lima	2010
SC13	FD01851541	-	-	AM	Hospital Santa Maria Del Socorro	Ica	2006
SC14	FD01851320	-	-	AMC	CENAN/INS	Lima	2009
SC13	FD01851499	-	-	AM, N	LRR Chiclayo	Lambayeque	2005
SC13	FD01848615	-	-	AM, SXT	Hospital San Bartolomé	Lima	2005
SC17	FD01852858	-	-	AM, C, TE	DISA Lima Ciudad	Lima	2010
SC17	FD01848690	-	-	AM, C, CTX	Hospital San Bartolomé	Lima	2001

Antibiotic abbreviation: ampicillin (AM), chloramphenicol (C), ciprofloxacin (CIP), trimethoprim-sulfamethoxazole (SXT), cefotaxime (CTX), nalidixic acid (NA), amoxicillin-clavulanate (AMC), nitrofurantoin (N), tetracycline (TE), and ceftazidime (CAZ).

Figure 5. Heatmap showing presence and absence of antimicrobial resistance (due to chromosomal AMR gene variations, mobile genes, and plasmids), their corresponding phenotype, and sequence cluster classification in 50 Salmonella Typhimurium isolates.

Figure 5. Heatmap showing presence and absence of antimicrobial resistance (due to chromosomal AMR gene variations, mobile genes, and plasmids), their corresponding phenotype, and sequence cluster classification in 50 Salmonella Typhimurium isolates.

3.8. Hospital Emergencias Pediatrica Shows the Presence of Most of the MDR Isolates in Peru

Out of the 26 genotypically AMR isolates, 21 were from Lima, and 1 was from Apurimac, La libertad, Callao, Huanuco, and Cusco. Of these twenty-one AMR isolates from Lima, ten are MDR isolates, with seven found at Hospital Emergencias Pediatrica, two found at Instituto Nacional de Salud del niño, and one is present in LRR Apurimac. Out of the phenotypically resistant 44 isolates, 12 were MDR (Lima n = 11, Huanuco n = 1). From these twelve MDR isolates, six isolates were from Lima, Hospital Emergencias Pediatrica, two were from Instituto Nacional de Salud del niño and another four MDR isolates were distributed in four different hospitals in Lima. AM, NA and TE resistance were the most prevalent in these MDR isolates. Three isolates showed phenotypic resistance to maximum five antibiotics (FD01852539: AM + C + SXT + NA + TE; FD01852476: AM + CTX + NA + TA + CAZ; and FD01852499: AM + CIP + NA + AMC + TE) and were prevalent at Hospital Emergencias Pediatrica and Instituto Nacional de Salud del niño. Importantly, we also found that Salmonella acquired MDR genes after 2015 (Table 2 and Table S3).

3.9. Probable NA Resistance New Single Nucleotide Variations from Core AMR Gene Analysis

Based on the core genome analysis, we identified 30 non-reported non-synonymous single nucleotide variations (SNVs) in 11 AMR genes in our samples. Among these, four aminoglycoside resistance genes (acrB1, kdpE, mdtB, mdtC) had 17 SNVs and the microcin resistance gene (yojI) showed two SNVs. However, we were unable to proceed with these variations as we did not perform the phenotypic tests for these drugs (Supplementary Tables S3 and S5). The other 11 new SNVs detected in 6 AMR genes (emrB, parC, gyrB, acrB2, sdiA, nfsA) are associated with fluoroquinolone, multiclass, and nitrofurantoin resistance were selected for further analysis. The isolates that showed a specific drug resistance phenotype but had no known marker for that phenotype but had a new variation were selected. Following this strategy, finally, a total of three SNVs in two genes parC (A554T, T571S for NA), and sdiA (K104Q for NA), were selected for structure-based functional annotation (Supplementary Table S5).

3.10. Probable Nitrofurantoin Resistance Eight New Genes and Their Variations from GWAS Analysis

Out of the 44 phenotypically resistant isolates, 13 isolates showed N resistance but had no known AMR markers for N. Therefore, in the GWAS analysis, we focused on these 13 isolates to predict the potential causal SNVs associated with N resistance. We observed eight non-synonymous substitution variations in eight different genes in eleven isolates that were significantly associated with N resistance (p < 0.0000526). It is also important to note that these 11 isolates belong to SC13, and the SC13 group had 24 isolates. Therefore, any SC13 isolate containing all eight variations has a 46.0% chance of showing the N-resistant phenotype. The SNVs were MdtH L15P, hypothetical protein G6V, QseC1 L19V, PpnN R116C, BioH A236T, WecA B284I, PurA A103G, and Tsr2 D161G. The details of these genes and SNVs are provided in Table 3. The first logistic regression-based negative results (false positive) of N resistance association are shown in Supplement Table S6.

Table 3. Potential nitrofurantoin resistance eight genes and their SNVs identified with GWAS.

Locus	Gene	Product	Effect	p-Value
NC_003197.2:1252216	mdtH	Multi-drug resistance protein MdtH	missense_variant c.44T > C p. Leu15Pro	1.85112239616509e-06
NC_003197.2:2585776		hypothetical protein	missense_variant c.17G > T p. Gly6Val	3.9922233235481e-06
NC_003197.2:2933724	qseC_1	Sensor protein QseC	missense_variant c.55A > C p.Ile19Leu	1.85112239616509e-06
NC_003197.2:3119194	ppnN	Pyrimidine/purine nucleotide 5′-monophosphate nucleosidase	missense_variant c.346C > T p.Arg116Cys	1.85112239616509e-06
NC_003197.2:3667964	bioH	Pimeloyl-[acyl-carrier protein] methyl ester esterase	missense_variant c.706G > A p.Ala236Thr	1.85112239616509e-06
NC_003197.2:4127935	wecA	Undecaprenyl-phosphate alpha-N-acetylglucosaminyl 1-phosphate transferase	missense_variant c.850G > A p.Val284Ile	1.85112239616509e-06
NC_003197.2:4609418	purA	Adenylosuccinate synthetase	missense_variant c.308C > G p.Ala103Gly	1.85112239616509e-06
NC_003197.2:4790581	tsr_2	Methyl-accepting chemotaxis protein I	missense_variant c.482A > G p.Asp161Gly	1.85112239616509e-06

Table 3. Potential nitrofurantoin resistance eight genes and their SNVs identified with GWAS.

Locus	Gene	Product	Effect	p-Value
NC_003197.2:1252216	mdtH	Multi-drug resistance protein MdtH	missense_variant c.44T > C p. Leu15Pro	1.85112239616509e-06
NC_003197.2:2585776		hypothetical protein	missense_variant c.17G > T p. Gly6Val	3.9922233235481e-06
NC_003197.2:2933724	qseC_1	Sensor protein QseC	missense_variant c.55A > C p.Ile19Leu	1.85112239616509e-06
NC_003197.2:3119194	ppnN	Pyrimidine/purine nucleotide 5′-monophosphate nucleosidase	missense_variant c.346C > T p.Arg116Cys	1.85112239616509e-06
NC_003197.2:3667964	bioH	Pimeloyl-[acyl-carrier protein] methyl ester esterase	missense_variant c.706G > A p.Ala236Thr	1.85112239616509e-06
NC_003197.2:4127935	wecA	Undecaprenyl-phosphate alpha-N-acetylglucosaminyl 1-phosphate transferase	missense_variant c.850G > A p.Val284Ile	1.85112239616509e-06
NC_003197.2:4609418	purA	Adenylosuccinate synthetase	missense_variant c.308C > G p.Ala103Gly	1.85112239616509e-06
NC_003197.2:4790581	tsr_2	Methyl-accepting chemotaxis protein I	missense_variant c.482A > G p.Asp161Gly	1.85112239616509e-06

3.11. Functional Annotation of Nalidixic Acid and Nitrofurantoin Resistance Possible New Variations

We tested the effect of the newly identified three SNVs from two core AMR genes (for NA resistance) and eight new SNVs from eight genes (N resistance) from our GWAS analysis following the methods we described.

We determined the 3D structure of all these proteins, except the hypothetical protein and tsr_2. Therefore, we succeed in checking the effect of twelve new SNVs in nine genes (Supplementary material S1). Out of the four known fluoroquinolone resistance SNVs in gyrA, we found three SNVs that are destabilising (

Table 4. Most of the non-synonymous variations for the NA and N resistance are predicted to be destabilising based on ΔΔG (Kcal/mol) calculation.

Analysis	Gene	Non-Synonymous Variations	Isolates	Prediction Stability ΔΔG (Kcal/mol)	Stability	Antibiotic Resistance
Variants in core AMR genes	parC	A554T	3	−1.82	Destabilising	Fluoroquinolone (NA, CIP)
		R360H, R365H	1	−1.86	Destabilising
		T571S, A554T	1	−1.01	Destabilising
	gyrA	D87G	1	−0.84	Destabilising
		D87Y	1	0.18	Stabilising
		S83F	1	−0.85	Destabilising
		S83Y	6	−0.91	Destabilising
	sdiA	K104Q	1	−0.18	Destabilising	Multiclass
	nfsA	E99K	7	−0.61	Destabilising	Nitrofurantoin
New SNPs using GWAS	mdtH	L15P	24	−1.34	Destabilising
	hypothetical protein	G6V	24	-	-
	qseC_1	I19L	24	−0.09	Destabilising
	ppnN	R116C	24	0.24	Stabilising
	bioH	A236T	24	−1.04	Destabilising
	wecA	V284I	24	−0.47	Destabilising
	purA	A103G	24	−1.19	Destabilising
	tsr_2	N161G	24	-	-

), indicating that many of the variations that have destabilising effects may be associated with drug resistance. Only the new SNVs for parC p.T571S, A554T, and sdiA p.K104Q were present in strains phenotypically resistant to NA and N, respectively, with no other known resistance gene. While we tested the effect of the identified unknown SNVs in other genes, except ppnN p.R116C, all 11 new SNVs were found to destabilise their corresponding protein. Therefore, considering the destabilising effect of the gyrA mutations, we may conclude that our identified new SNVs may be associated with NA and N resistance (Table 4).

4. Discussion

Our study is the first using WGS analysis and antimicrobial susceptibility test (AST) of S. Typhimurium isolates (n = 90) between 2000 and 2017 from different cities in Peru to examine the diversity, WGS-AST correspondence, resistome profile, and emerging lineages. Our study concludes that a considerable nucleotide, gene, and phylogenetic group diversity circulates in Peru.

Peruvian samples belong predominant to ST19, with a diverse accessory content, constituted by plasmids and segregated by phylogenetic groups. We reported that 94.4% of S. Typhimurium isolates in Peru belong to sequence type ST19. This was expected because ST19 is among the most predominant ST associated with gastroenteritis cases worldwide, and it is the most ancestral and diverse phylogenetically ST for the serovar Typhimurium [11,19,20].

The population contained an open pan-genome with many repertoire genes for global and Peruvian isolates comparable to previous reports [19,63,64]. However, the core gene content was smaller than previously reported (3672 genes) [19] Fu et al. (3846 genes) [64]. Likewise, we defined the accessory genome as constituted of diverse plasmids without predicting other mobile elements. Other studies highlight that the composition of the accessory genome for S. Typhimurium is mainly based on the diversity contribution of prophage genes, up 23.4%, followed by other elements such as plasmids and mobile islands, up 13.3% [19,64]. In addition, we did not corroborate that phylogenetic groups segregate the accessory genome by PCA analysis (analysis not shown).

Initially, a low accumulation nucleotide of 400–600 SNPs was reported for the serovar [20]. Although, other studies showed a large accumulation of SNPs by isolate (1232 SNPs) [64] and a total of 62,884 SNPs for the African population [10]. These studies are comparable to the vast number of polymorphisms (3045 SNPs) within genomes circulating in Peru with a moderate median pi value (0.135) at the intrapopulation diversity level described by Pons, 1996 [65].

The phylogenetic topology showed two high-order clades and population structures identified at least 10 lineages supported by subclades with depth branches to multiple branches. This topology has already been discussed previously in studies; this includes an alpha clade basal with livestock samples with a diversity of terminal branches corresponding to clonal expansions. A distinct beta clade is characterised by multiple lineages from the vast host (including wild avian) that are deeply rooted [19,20]. Likewise, we report the emergence of SC13 and SC17 lineages by their relative widespread in the country and the prevalent MDR phenotype of SC9 and SC8 lineages. It has already been shown that subclades are under different anthropogenic selection pressures. Antibiotic use might provide the selection pressure driving the emergence of sub-lineages containing AMR genes [10,11,19,66]. MDR strains harbour variable numbers of resistant plasmids reported mainly from livestock samples and outbreaks in hospitals and foods [19,20].

S. Typhimurium presents considerable resistance to first-line antibiotics, susceptibility to cephalosporins and ciprofloxacin, and low sensitivity value of phenotypic prediction due to the unknown resistance mechanism. This study explored resistance trends, transmission, profile genotypic and phenotypic, correlation, and discovery of new genetic bases for resistance phenotypes with phylogeny. The study reports a high number of strains (48.8%) with phenotypic resistance to at least one antibiotic compared to the variable prevalence (26.0%/n = 95, 61.6%/n = 193, 37.3%/n = 3491) for NTS clinical and food samples from Peru, US, and England, respectively [21,22,67]. We found considerable resistance to first-line antibiotics (nitrofurantoin, tetracycline, nalidixic acid, and ampicillin); only 12.0% were MDR strains. The persistence of resistance to first-line antibiotics over the years in this study was expected due to the continued use of treatment in Salmonella [24].

Previous studies of NTS and S. Typhimurium strains show a significant number (24.3%-43.0%) of MDR strains to first-line antibiotics except for nitrofurantoin, and with the addition of sulphonamide, streptomycin, and chloramphenicol (in some cases) in American Latin, USA, and England [21,22,24,66,68]. This minor prevalence of MDR strains compared to other studies of NTS samples should include different serovars that present high and diverse MDR. In recent years, S. Infantis has been the most predominant serovar detected in clinical samples associated with high resistance to first-line antibiotics, third-generation cephalosporins [7,8], and ciprofloxacin [8,68,69]. Although 85.0% of our samples were clinical isolates, S. Typhimurium, compared to other serovars, still presents a reduced resistance to priority antibiotics. On the contrary, S. Typhimurium is a relevant pathogen in guinea pigs with distinct and moderate AMR profiles that include colistin and enrofloxacin [70], erythromycin and nitrofurantoin resistance [71]. In addition, Salmonella isolates from chicken meat show high quinolone resistance (enrofloxacin and NA). Currently, the resistance spectra of the MDR strains of Salmonella serovars have been emerging in farm animals [24]. Because quinolones, chloramphenicol, aminoglycosides, and nitrofurantoin are exhaustively used as treatment and/or prophylactic in farm animals [68,70,71,72,73,74], constituting the leading resource of transmission on the emergence of resistant strains is reported.

WGS showed the content of AMR chromosomal genes with an impact on the resistance based on their expression. These AMR genes encoded most efflux pumps reported by Seribelli [11]. Complementary, we identified 32 AMR acquired genes (Table 1, this profile AMR genes were similar and minor to previous works [21,22,23,66,73]. These studies detecting additionally other EBSLs, ribosomal protection mechanisms, PMQR genes, phosphotransferase, efflux pump, and an acetyltransferase that confer resistance to spectrum extended beta-lactam, tetracycline, fluoroquinolone, macrolides, and phenicol, respectively. Genotypic AMR profile to priority antibiotics shows the nalidixic acid resistance associated with the presence of PMQR genes (qnrB5, qnrE2, qnrB19) or mutations (D87G, D87Y, S83F, S83Y) in the gyrA gene, and presence combined of both markers are associated with ciprofloxacin resistance [21]. Interestingly, the considerably reduced ciprofloxacin susceptibility did not exhibit these markers. This is because the efflux pump’s overexpression was related to intermediate resistance [75], and the disc diffusion test did not adequately detect reduced susceptibility to fluoroquinolones [76]. In only two strains, extended-spectrum beta-lactam resistance was associated with three EBSLs genes (blaCTX-M-15, blaSHV-12, blaSHV-134). Additional WGS allowed the detection of genes associated with additional drug resistance as aminoglycoside, aminocoumarin, bacitracin, nitroimidazole, microcin, lincosamide, bleomycin, fosfomycin, colistin, and multiclass.

A worrying trend is an increase in resistance to treatment antibiotics (extended-spectrum cephalosporins and ciprofloxacin) for Salmonella. Despite them, our study shows susceptibility to ciprofloxacin and third-generation cephalosporins, reported in other works [21,22,24,66,67,68,73] without co-resistance to both antimicrobial classes. Whereby the use of third-generation cephalosporins for treatment in S. Typhimurium infections would be recommended.

WGS strategy allows monitoring and complementing the prediction of phenotypic AMR profiles [23,73]. We found good accuracy values (prediction of true positives and negatives) and specificity (the absence of known markers predicting true negatives). Comparable successful correlation values were reported, 99.0% [73], 97.8% [21], 95.4% [22], 89.9% [23], and 85.4% [69]. Nonetheless, we found the best sensitivity values to predict resistance to only ceftazidime and ciprofloxacin and the highest discrepancy values to predict resistance to beta-lactam, chloramphenicol, SXT, and nitrofurantoin. Previous work also found discord in the prediction of resistance to beta-lactams [66,73], sulfamethoxazole [23], ciprofloxacin [73] and tetracycline [66]. Due to the high number of contradictions between genotype–phenotype compared with previous reports cannot conclude that this would be an alternative method of predicting antimicrobial susceptibility.

Mismatch categories with a lower sensitivity, where an isolate is genotypically predicted to be susceptible but exhibits phenotypic resistance, highlights limitations based on sequence quality [77], partial assemblies, lack of updated AMR database, and rely on the prediction based only on the genome [12,78,79]. Continuous findings should be carried out to identify novel resistance mechanisms and be incorporated into the reference databases to maintain a high level of prediction sensitivity. Other inconsistent results support the need for combined AST strategies, such as the microdilution test, an efficient method with quantitative results [28]. This study also highlights the performance of the routine antibiotic susceptibility test and works with a balanced number of phenotypic samples and population representations [80].

The diversity of plasmid families carrying AMR genes in isolates belong to SC8 and SC9 lineages and was detected in two hospital centres considered the focus of transmission of antimicrobial resistance. WGS also allows identifying AMR genes commonly present on plasmids, primary transmission resources, and related to the emergency lineages. Here, we reported at least 17 resistant plasmid families, including conjugative and mobilisable plasmids, with the potential threat of spreading AMR genes in NTS [14,22]. These resistance plasmids belong to F, ColE, I1, C, HI1, HI2, and N families, such as a previous work of Salmonella isolates from food animals in the USA predicted 212 resistance plasmids [22]. In Peru, a conservative virulence plasmid (pSV) in S. Typhimurium [67] and MDR Mega plasmid in S. Infantis [8] have been described. In Peru, S. Typhimurium isolates are prevalent that contains the virulence plasmid (pSV) [67,81]. However, the absence of plasmids in the phylogenetic group could be due to the competence with other AMR plasmids or other genetic and environmental factors that modulate the residence [81].

Interestingly, the genotypic profile is linked to family plasmid in some lineages. For example, two only strains contain replicon plasmid carrying ESBLs genes and unnamed replicon plasmid P4 carrying beta-lactamases genes in lineage SC8. Lineages S8 and SC2 harbour PMQR genes in Col(pHAD28) plasmid. Col(pHAD28) plasmids related to fluoroquinolone resistance were reported [82]. MCR-1-carrying lncI2 and IncHI2 belong to the SC9 lineage. These dominant mobilisable plasmids show colistin resistance [83,84]. Likewise, only a strain carries fosA3-carrying IncFIB plasmid. Previous studies have reported antibiotic’s last line resistance as the fosA3 gene in AMR plasmids [85] and IncFIB plasmids [86,87]. Intriguingly, SC9 isolates harbour the IncHI2 plasmid that carries many AMR genes. This plasmid is a dominant mobilisable detected among MDR Salmonella, playing a role in the acquisition of ARGs, and has been reported recently encoding ESBLs [10,88,89,90]. Thus, active surveillance is needed to minimise the global spread of lncIA-I(Alpha), IncI-gamma/K1, lncFIB, lncI2, and IncHI2 resistant plasmids link an SC9, SC8, and SC2 lineage. MDR strains were found mainly from the Hospital de Emergencias Pediátricas, and Instituto Nacional de Salud del Niño. The previous report shows the presence of MDR isolates from serovar Infantis in the Hospital de Emergencias Pediátricas [91]. Therefore, antimicrobial screening routines should be implemented to mitigate the spread of these healthcare-associated MDR strains in Peruvian hospitals.

GWAS analysis allows us to identify new non-synonymous mutations that can potentially improve resistance fitness; however other resistance confirmation strategies are necessary. We reported new SNVs in parC p.R360H, R365H gene, and sdiA p.K104Q that show destabilising effect protein with a possible impact on the protein function. Mutations in gyrA-parC genes decrease the binding affinity of quinolones with DNA-topoisomerase enzymes [92,93]. Efflux pumps are encoded in chromosomes and play intrinsic roles in multi-drug resistant Gram-negative bacteria [25]. sdiA gene acts as a positive regulator of the constitutive expression of the AcrAB–TolC pump system [94]. Nucleotide variations in this regulator [25] act in high-level fluoroquinolone resistance and other antimicrobials [94,95]. Genetic variation targets could result in overexpression of these proteins. For instance, SNVs in acrR regulator or multi-drug pump AcrAB were associated with high-level fluoroquinolone resistance [94].

We found that considerable nitrofurantoin resistance could not be associated with known AMR markers such as those previously reported [96]. This would happen because the ARM databases have few AMR genes since the mechanisms of action of nitrofurans are poorly studied [94]. Thereby, GWAS analysis identified eight non-synonymous substitutions potentially associated with resistance to nitrofurantoin; six show the destabilising effect protein. Genetic variation is within multi-drug efflux pump [97]; regulator purine homeostasis and biosynthesis [98]; biotin ring assembly [99]; pathway LPS O-antigen biosynthesis [100], and chemotactic-signal transducers added. However, we did not find that the mutational effect on these functional mechanisms could confer nitrofurantoin resistance. Likewise, these variants were present in the lineage SC13. To avoid spurious associations, GWAS analysis using the population stratification covariate was performed subsequently without associating any mutation to nitrofurantoin resistance. Hence, we suggest that the maintenance of these combined mutations in phylogenetic group strains would be fixed randomly and consequently could confer a resistance advantage.

5. Conclusions

Our work based on the WGS analysis has allowed us to understand the dynamics and determinants of antimicrobial resistance distinguished by the population diversity for S. Typhimurium. However, due to low sensitivity values from genotype–phenotype resistance correlation, it is still necessary to evaluate the use of WGS to predict AMR susceptibility. We recommend the third-generation cephalosporin antibiotic as a potential treatment against infection by S. Typhimurium. We can reinforce that WGS constitutes a complement but not an alternative to traditional methods to infer antimicrobial susceptibility, as a powerful tool that allows genome-based epidemiological study, monitoring AMR genes, virulence, plasmid typing, outbreaks, understanding of resistance mechanism, and transmission patterns. Pathogen genomic surveillance should be expanded globally and continuously monitored for better treatment of Salmonellosis and control strategies against the disseminating AMR. Future work should be based on better discordant prediction due to the absence of AMR genes on resistance phenotypes, adding new resistance mechanisms and improving the database’s reliability, replicating assay, large-scale samples, and supplementary technical methods.