Genomic Characterization of Antibiotic Resistant Escherichia coli Isolated From Domestic Chickens in Pakistan

Poultry husbandry is important for the economic health of Pakistan, but the Pakistani poultry industry is negatively impacted by infections from Escherichia coli. We performed Illumina whole genome sequencing on 92 E. coli isolates obtained from the livers of deceased chickens originating in five Pakistani geographical regions. Our analysis indicates that the isolates are predominantly from the B1 and A clade and harbor a diverse number of antibiotic resistance and virulence genes, with no linkage between phylogeny and antibiotic resistance gene presence but some association between phylogeny and virulence gene and SNP presence for the B1 and E phylogroups. The colistin resistance gene mcr-1 and the quinolone resistance gene qnrS1 were both found in 13/92 isolates. Alarmingly, 82/92 of the E. coli strains characterized in this study are multidrug resistant with 100% (92/92) resistance to lincomycin, 81.5% (75/92) to streptomycin, 79.3% (73/92) to ampicillin and 66.3% (61/92) to ciprofloxacin. These results provide a high-resolution analysis of poultry-associated E. coli isolates in an area with a high endemic burden of antibiotic resistance. Surveillance of antibiotic resistance in poultry associated E. coli isolates is an important pillar of the One Health concept to integrate analysis of potential pathogens in human, animal, and environmental niches.

Introduction

The poultry industry is an important component of Pakistan’s gross domestic product (Hussain et al., 2015). However, the Pakistani poultry industry faces several environmental conditions that threaten continued economic output, livestock health, and human health, including pathogen Escherichia coli (Hussain et al., 2015; Manges, 2016). E. coli is particularly relevant to human health, as E. coli can cause a diverse array of infections, exist as a gut commensal, and is often antibiotic resistant through horizontally acquired antibiotic resistance genes (ARGs) (Croxen and Finlay, 2009). In poultry, these infections often manifest as colibacillosis, which can lead to many health abnormalities, often resulting in chicken death (Dho-Moulin and Fairbrother, 1999). Numerous studies have investigated pathogenicity of avian associated E. coli isolates and identified factors such as iss (serum survival), ibeA (cell invasion), sitA (iron acquisition), and iroN (iron acquisition), which are associated with increased virulence in both chickens and humans (Mellata et al., 2003; Sabri et al., 2006; Cieza et al., 2015; Sarowska et al., 2019). Unlike the strong associations identified between human uropathogenic E. coli and the B2 E. coli phylogroup, previous analysis of European E. coli from birds has not identified genetic signatures strongly associated with any single phylogroup (Cordoni et al., 2016).

Given the economic importance of the domestic poultry industry to Pakistan and the frequent identification of ARGs emanating from the central Asian region (CDDEP, 2017), a comparative phenotypic and genotypic analysis of multidrug-resistant E. coli from poultry infections throughout Pakistan is warranted. Total consumption of antimicrobials by chickens in Asia is expected to increase by 129% between 2010 and 2030, driven in large part by increasing chicken consumption in India and Pakistan concomitant with population growth and poor antibiotic stewardship (Hussain et al., 2015; Van Boeckel et al., 2015). Although policies have recently been introduced to encourage restricted antibiotic use in food animals in Pakistan, challenges to regulation and surveillance remain such that accurate estimates of antibiotic usage in poultry rearing are elusive (Rahman and Mohsin, 2019). The types of antibiotics used might be inferred from the multi-drug resistance profiles identified in E. coli isolates by previous studies, alarmingly with high prevalence to antibiotics of clinical importance (Idrees et al., 2011; Akhtar et al., 2016; Azam et al., 2019). However, these studies have been limited by the geographic range of their isolates and lack of whole genome sequencing and analysis.

Because antibiotic resistance transmission rates in poultry-associated bacteria may be high in Pakistan due to poultry rearing conditions, constant antibiotic selection, and horizontal gene transfer, it is important to closely monitor the landscape of antibiotic resistance in the Pakistan poultry industry as a whole using genomic methods. To address the knowledge gap on ARG burden in E. coli from chickens in Pakistan, we collected E. coli isolates from the livers of chicken in five regions of Pakistan and performed Illumina whole-genome sequencing. We then analyzed the population structure, identified relevant gene presence, and assayed for phenotypic antibiotic resistance in these isolates, with the primary goals of identifying any associations between virulence or resistance determinants and phylogroup, or between phylogroup and geographic region.

Materials and Methods

E. coli Cohort

A total of 1,219 liver samples from culled layer and broiler chickens that had poor birth growth and reduced appetite but not otherwise symptomatic for colibacillosis were collected from the National Reference Laboratory for Poultry Diseases in Pakistan via federal and provincial sentinel surveillance laboratories under a national surveillance program from 2015 to 2017. As the National Reference Laboratory for Poultry Disease does not handle live animals for experimentation, ethics approval was not sought or obtained for this study. The chickens originated in the Pakistani provinces of Balochistan, Sindh, Punjab, Khyber Pakhtunkhwa, and the federally administered Islamabad Capital Territory. Whole liver samples were cultured in nutrient broth (Sigma-Aldrich, St. Louis, MO, United States) at 37°C for 24 h. A loopful of the overnight culture was then plated onto eosin-methylene blue agar and grown at 37°C for 24 h. Suspected E. coli isolates were identified using the API 20E assay (bioMérieux, Durham, NC, United States), and glycerol stocks were generated. Frozen cultures were sent to Washington University in St. Louis for further analysis. Isolates were plated onto blood agar and E. coli identification was confirmed with the MALDI-TOF VITEK MS IVD v2.3.3 (bioMérieux, Durham, NC, United States) mass spectrometry system (Richter et al., 2013).

Illumina Whole-Genome Sequencing

Frozen stocks of E. coli isolates were plated onto blood agar using four-quadrant streaking and ∼10 morphologically similar colonies from the fourth quadrant were used as input for the QIAamp BiOstic Bacteremia DNA Kit (Qiagen, Germantown, MD, United States). A total of 0.5 ng of genomic DNA per isolate was used to create sequencing libraries with the Nextera Kit (Illumina, San Diego, CA, United States) (Baym et al., 2015). The libraries were pooled together at equimolar concentrations and sequenced on a NextSeq 500 to obtain 25–183X coverage of each genome with 2 bp × 150 bp reads. The reads were demultiplexed by barcode and Illumina adaptors and contaminating sequences were removed with Trimmomatic v.38 (Bolger et al., 2014) and Deconseq v.4.3 (Schmieder and Edwards, 2011), respectively. The processed reads were used to construct de novo assemblies of each genome with SPAdes v3.13.0 (Bankevich et al., 2012). The assembly metrics of the scaffolds.fasta files were assessed with QUAST v4.5 (Gurevich et al., 2013) and open reading frames identified with Prokka v1.12 (Seemann, 2014). A total of 92 genomes with less than 300 contigs were chosen for downstream genomic and phenotypic analysis (Supplementary Table S1).

In silico Analysis

To obtain phylogroup information for each E. coli genome, we gathered 11 publicly available genomes from known E. coli phylogroups and identified open reading frames using Prokka (Supplementary Table S2) (Schreiber et al., 2017; Hutton et al., 2018). The gff files from Prokka for the phylogroup reference strains and the genomes sequenced in this study were used as input for Roary v3.12.0 to construct a core-genome alignment of the 2,755 core-genes with PRANK v1.0 (Loytynoja, 2014; Page et al., 2015). The core-genome alignment file was converted into an approximate maximum likelihood tree with FastTree v2.1.10 and the resulting newick file was uploaded to iToL¹ (Letunic and Bork, 2007; Price et al., 2010). In parallel, we identified in silico antibiotic resistance determinants for acquired antimicrobial genes using ResFinder v4.0 and for E. coli single nucleotide polymorphisms with PointFinder v4.0 (Kleinheinz et al., 2014; Zankari et al., 2017). Additionally, we identified known virulence genes with VirulenceFinder v1.5 and Enterobacteriaceae plasmid replicons with PlasmidFinder v4.0 (Carattoli et al., 2014; Kleinheinz et al., 2014). Hypergeometric tests were used to determine the significant enrichment of isolate groups within phylogenetic clades or dendrogram clusters, with Bonferroni correction for multiple hypothesis testing.

Antibiotic Susceptibility Testing

Phenotypic antibiotic resistance assessment was performed using a protocol similar to those previously described for various Enterobacteriaceae on a variety of antibiotics relevant to human and veterinary use (Potter et al., 2018a, b). Briefly, the phenotypic antibiotic resistance of E. coli isolates was assayed through growth on Mueller-Hinton agar in the presence of antibiotic-laden Kirby-Bauer disks in accordance with Clinical Laboratory and Standards Institute requirements for human derived isolates and also for veterinary derived isolates (Supplementary Tables S3, S4) (CLSI, 2018). The interpretation of zone of clearance was used to create a heatmap with hierarchical clustering for each isolate in pheatmap (R studio). ComASP^TM (Liofilchem) was used exactly in concordance with manufacturer’s instructions to quantify colistin resistance in isolates with identified in silico resistance determinants. Disk diffusion values for ceftriaxone/ceftiofur and ciprofloxacin/norfloxacin were plotted as XY coordinates in Prism v8.0 and linear regression analysis was performed using default conditions to quantify the R² value.

Results

Of the 92 genomes that were used for full genomic and phenotypic analysis, 41.3% (38/92) originated in Punjab and 26% (24/92) came from Islamabad. The remainder of the isolates originated from Khyber Pakhtunkhwa (14%; 13/92), Balochistan (12%; 11/92), and Sindh (7%; 7/92) (Figure 1). Phylogenetic reconstruction of the similarity between the E. coli within our cohort and with known phylogroup strains indicates that there is not a clear association between phylogroup and geographic region, although we note isolates from Balochistan fell exclusively in the B1 Clade (Figure 2). Overall, 53.2% (49/92) of the cohort are in the B1 Clade. The rest of the isolates are in Clade A (22/92), B2 (11/92), E (9/92), and F (1/92). The most abundant sequence types were ST115 (16/92) in B1 and ST117 (9/92) in phylogroup B2 (Figure 2A).

FIGURE 1

Figure 1. Escherichia coli cohort originates from five regions in Pakistan. Map depicting the four provinces of Pakistan and the Islamabad Capital Territory from which isolates were obtained. Number of isolates are shown adjacent to the provinces, and colors correspond to isolate prevalence.

FIGURE 2

Figure 2. Escherichia coli genomes are predominantly in B1 and A phylogroup with a mosaic of antibiotic resistance determinants and virulence genes. (A) Population structure of the E. coli cohort depicting the phylogroup, MLST, and region obtained. Presence absence of antibiotic resistance determinants (B), virulence genes (C), and plasmid replicons (D) identified using ResFinder, PointFinder, VirulenceFinder, and PlasmidFinder on the cohort.

To identify a genotypic basis for phenotypic antibiotic resistance, following Illumina whole-genome sequencing, we applied ResFinder to identify acquired ARGs and PointFinder to locate relevant SNP resistance determinants in the assembled genomes. Consistent with previous reports on genomic analysis of E. coli isolates, we identified a mosaic of antibiotic resistant determinants and virulence genes within our cohort with no clear association between ARG composition, phylogroup, and region source, other than the identification of all isolates from Balochistan as members of the B1 phylogroup (p < 0.01, hypergeometric test) (Figure 2A). In total, we identified 49 unique ARGs and 5 previously validated antibiotic resistance-conferring SNPs (Figure 2B) (Zankari et al., 2017). The median number of ARGs per isolate was 7 and the median prevalence for each ARG was 8. We found that 17/49 of the ARGs are predicted to have activity against aminoglycosides, 11/49 against folate-synthesis inhibitors, 6/49 against β-lactams, 5/49 against amphenicols, and 5/49 against lincosamides. In addition, we identified 3 tet ARGs (76 isolates), and the quinolone resistance gene qnrS1 (13 isolates). PointFinder identified amino acid changes in GyrA (D87G 4/92 isolates, D87N 44/92 isolates, D87Y 3/92 isolates, and S83L 53/92 isolates), ParC (E84A 1/92 isolates, S80I 51/92 isolates, S80R 3/92 isolates), ParE (I464F 2/92 isolates), and PmrB (V161G 1/92 isolates). The aminoglycoside ARGs include representatives of the aac, aadA, and aph families. Notably, no known carbapenem resistance genes were identified within this E. coli cohort. The most prevalent β-lactamase was bla_TEM–1B (found in 66/92 isolates). The class A β-lactamases bla_CTX–M–15 and bla_CTX–M–1 were found in 1/92 and 12/92 isolates. The only non-class A β-lactamase present was bla_CMY–2, which was found in a single isolate. The most conserved gene, mdf(A), conferring multidrug resistance, was found in 92/92 of the E. coli isolates, indicating that it is a core-gene within this cohort (Edgar and Bibi, 1997). mdf(A) was the sole ARG in 7/92 of the isolates. Despite this conservation, the other lincosamide ARGs, lnu(F) (1/92), mef(B) (2/92), mph(A) (7/92), and mph(B) (3/92) were at or below the median ARG prevalence. The colistin ARG mcr-1 was found in 13/92 of the E. coli isolates. The prevalence of mcr-1 was much higher than SNPs in genes associated with colistin resistance as only 1/92 isolates contained the V161G amino acid change in PmrB. When hierarchically clustered on ARG presence (k = 5), no segregation of isolates by region or phylogroup is observed, except for an enrichment of B2 isolates (p < 0.01) in a cluster characterized by a high ARG load that includes bla_TEM–1B, tet(A), and mdf(A), as do most of the isolates, but also aadA2, aadA1, sul3, and cmlA1 (Figure 3A). When hierarchically clustered on antibiotic resistance-conferring SNPs, the isolates are segregated into two major clades, one with a preponderance of SNPs in gyrA and parC, and a clade of isolates most of which carry no SNPs (Figure 3C). B2 isolates are enriched in the latter clade (p < 0.001).

FIGURE 3

Figure 3. Escherichia coli phylogroups segregate better by virulence gene presence than ARG or SNP presence. Heatmaps depicting isolates as rows and ARGs (A), virulence genes (B), or SNPs (C) as columns. Rows and columns are hierarchically clustered by Euclidian distance. Region, phylogroup, or expected phenotypic resistance conferred by ARGs are portrayed as metadata, as indicated.

We applied VirulenceFinder on the sequenced E. coli cohort to annotate genes putatively involved in poultry infections (Figure 2C). We identified 21 virulence genes and found that the median number of virulence genes per isolate was four. Only 2/92 of the isolates had no known virulence genes identified. Consistent with their previously identified roles in E. coli virulence, the serum survival gene iss was found in 78.2% (72/92) of the isolates, the iron acquisition gene iroN was in 60.9% (56/92) of the cohort, and the long polar fimbriae gene lpfA was in 59.8% (55/92) of the cohort. When hierarchically clustered on virulence gene presence (k = 5), we again observed no segregation of the isolates by geographic region (Figure 3B). However, the B2 and E isolates were each enriched in their own clusters (p < 1E-9 and p < 0.01, respectively) characterized by different virulence gene profiles. While the larger B1 phylogroup was more distributed across the clusters, there was one cluster that was exclusively comprised of B1 isolates (p < 0.001). Together these data indicate that the phylogroups segregate better by virulence gene presence than ARG presence; the lack of segregation by geographic region suggests that chicken-borne E. coli strains are readily transmitted across Pakistan.

We applied PlasmidFinder using the Enterobacteriaceae database on our E. coli cohort to identify known plasmid replicons (Figure 2D). 90/92 isolates had plasmid replicons identified, with a maximum of 7 in EC_10, EC_44, and EC_67, and a median of 4.26 different replicons were identified, among which members of the IncF plasmid replicon family were the most prevalent (12/26). IncFIB(AP001918) and IncFII were the most prevalent among the isolates with 65/92 and 41/92 identified within the cohort, respectively.

To assess the effect of ARG burden on phenotypic antibiotic resistance, we performed antimicrobial susceptibility testing using the Kirby-Bauer Disk Diffusion method and the Clinical Laboratory and Standards Institute (CLSI) interpretative criteria from the M100 (Edition 29) and VET01 (Edition 5) for Enterobacteriaceae on a variety of antibiotics relevant for human and veterinary use. Using the definition of multidrug-resistant (MDR) as non-susceptibility to at least one agent in three or more antimicrobial classes, and extensively drug resistant (XDR) as susceptibility to at least one agent in only one or two classes assayed, we found that 82/92 are MDR but only 1/92 are XDR (Figure 4) (Magiorakos et al., 2012). Consistent with the presence of mdf(A) in all the genomes, all isolates were resistant to the lincosamide antibiotic lincomycin. Hierarchal clustering on the antibiotics with CLSI interpretive criteria (using 1 for resistant, 0 for intermediate, and −1 for susceptible) with a cluster cutoff just below the second node of the dendrogram resulted in seven clusters (k = 7) (Figure 4). The first cluster (from left to right) was characterized by susceptibility to β-lactams. There was widespread resistance to quinolones and the aminoglycoside streptomycin across all clusters, while the two rightmost clusters were characterized by additional resistance to the β-lactams cefazolin, ceftriaxone, and ceftiofur. There was no significant association between any cluster and any region or phylogroup, unlike the associations detected between clusters formed on genetic features and phylogroup (Figure 3). Given problems using disk diffusion testing for colistin resistance, we performed a broth minimum inhibitory concentration (MIC) assay using the ComASP^TM colistin test on all mcr-1 positive isolates, the pmrB SNP isolate 55, quality control strain E. coli ATCC 25922, and the mcr-1 positive E. coli AR Bank #0350 from the CDC and FDA Antibiotic Resistance Isolate Bank (Table 1). All mcr-1 positive isolates had MIC values of 4 or 8 μg/mL but the pmrB SNP isolate 55 had an MIC of 0.25 μg/mL. This cohort demonstrated 100% in vitro susceptibility to meropenem, imipenem, cefotetan, piperacillin-tazobactam, and amikacin. We found strong association between phenotypic resistance within the 3^rd generation cephalosporin ceftriaxone and ceftiofur (R² = 0.9004) and quinolones ciprofloxacin and norfloxacin (R² = 0.8897) (Supplementary Figure S1). One isolate (EC_44) was discordant for the 3^rd generation cephalosporins and tested as resistant to ceftriaxone but intermediate to the veterinary antibiotic ceftiofur and contained bla_CMY–2 and bla_TEM–1B. While there are not CLSI interpretative criteria for norfloxacin, one isolate (EC_72) without any identified quinolone resistance determinants tested as ciprofloxacin susceptible but had a comparatively low disk diffusion radius to norfloxacin.

FIGURE 4

Figure 4. Escherichia coli isolates are predominantly (82/92) MDR. Heatmap depicting each isolate as a column and each row as an antibiotic. Columns are hierarchically clustered by Euclidian distance. Region and phylogroup are portrayed as metadata.

TABLE 1

Table 1. Colistin susceptibility testing results.

Discussion

Antibiotics are widely used in the poultry industry of developing nations for the prevention of disease and growth promotion (CDDEP, 2017). Despite growing concern over the selective pressure this contributes to antimicrobial resistance and the potential for this to expand the reservoir for the dissemination of MDR pathogens, global annual consumption of antimicrobials by food animals is expected to increase by 67% between 2010 and 2030 (Van Boeckel et al., 2015). To our knowledge, this study represents the first use of next-generation sequencing and bioinformatics techniques in combination with phenotypic susceptibility testing to assess the extent of antibiotic resistance in poultry E. coli isolates from multiple provinces of Pakistan.

Given that chicken infections by E. coli can be explained by atypical environmental conditions leading to increased stress-related immunosuppression, some groups argue that the concept of avian pathogenic E. coli has no strong basis (Hoerr, 2010; Shini et al., 2010; Collingwood et al., 2014). This is consistent with our analysis that did not find any common genetic signature associated with chicken-borne E. coli isolates (Zhu Ge et al., 2014; Paudel et al., 2016). To understand if the chicken-borne E. coli isolates within our cohort have an identifiable population structure, we applied Roary and FastTree to the core-genome of the E. coli cohort and a selection of reference phylogroup strains to construct a phylogenetic tree. Our results indicate that a majority of the isolates belong to the B1 or A clades. It is possible that a greater extent of genetic signatures related to phylogroup may be identified with the more stringent criteria of sequencing isolates derived from active chicken infections, and not just birds with failure to grow. These findings are similar to those describing E. coli from bovine mastitis in Ireland, which also predominantly yielded B1 and A clade genomes (Keane, 2016). An analysis of avian associated E. coli isolates from several countries in Europe found variation within each country but that A1 and B2 phylogroups were the most prevalent, indicating that on a continental scale geographic location may affect E. coli background (Cordoni et al., 2016). A broader analysis of European avian associated E. coli isolates came to similar results with B2 and A as the dominant phylogroups (Mora et al., 2013). Since there was not a control chicken gut arm of our study, it is possible that the B1 and A clade isolates we found are also found as commensals in a greater proportion of the chickens. We found that the phylogroups are represented evenly in the different regions of Pakistan, except for the finding that all Balochistan isolates were B1. This may be simply due to the low sample size from that region (11 isolates), but the fact that we did observe multiple phylogroups in the similarly sized isolate sets from Sindh and Khyber Pakhtunkhwa, and that most of Balochistan is geographically separated from the rest of Pakistan by the Sulaiman and Brahui ranges (southern offshoots of the Hindu Kush Himalayan Region), suggest the possibility that phylogenetic segregation of chicken-associated E. coli at the scale of Pakistan’s geography may be detected in a larger collection of isolates. The B2 isolates repeatedly segregated when the isolates were hierarchically clustered by ARG, virulence gene, or SNP presence, as well as phenotypic resistance, reflecting the longer branch lengths of these isolates in the phylogeny in Figure 2A compared to most of the cohort. The phylogroups segregated better by virulence gene presence than ARG presence, possibly reflecting the increased mobility of ARGs. We found that IncFIB plasmid replicons were the most prevalent in our cohort. An analysis of E. coli in the United States found that IncFIB plasmids were significantly more present in E. coli suspected as avian pathogens or from retail poultry than compared to human and avian commensals (Johnson et al., 2007, 2008).

Of greatest concern from our analysis of the ARG content within the cohort was the presence of the mobilizable colistin ARG, mcr-1, in 14% (13/92) of the isolates, all of which displayed phenotypic resistance to colistin. Initially discovered in E. coli obtained from swine farms in China, it has since been identified in a wider cohort of bacteria including in a cohort of avian associated E. coli isolates from Egypt and China, although at a much lower prevalence of 1% (12/1220 isolates) (Liu et al., 2016; Barbieri et al., 2017). A wider analysis of E. coli from chickens in 13 Chinese provinces found that mcr-1 was detected in 4% (58/1136) of the genomes (Yang et al., 2017), while an analysis of 100 E. coli chicken isolates from Faisalabad, Pakistan found an mcr-1 prevalence of 8% (Lv et al., 2018). An E. coli isolate from a poultry farm in Tunisia was found to harbor both bla_CMY–2 and mcr-1 and exhibited resistance to several other antibiotic classes (Maamar et al., 2018). bla_CMY–2 was identified once in our cohort and that isolate did not have mcr-1 but it did have the pmrB V161G mutation which has been previously shown to confer colistin resistance, but did not in our study (Delannoy et al., 2017). Similar to previous analyses of E. coli isolates from poultry, iss was the most common virulence gene identified in our cohort (Keane, 2016). Deletion of iss from an isolate significantly perturbed E. coli growth in serum (Huja et al., 2015).

We found that 82/92 isolates are MDR using the definition of MDR as non-susceptibility to at least one agent in three or more antimicrobial classes (Magiorakos et al., 2012). This result is similar to analysis of avian associated E. coli isolates from Nepal, which determined that 94% (47/50) of their cohort was MDR (Subedi et al., 2018); and from Hebei, China which had 100% (87/87) MDR (Li et al., 2015). One analysis of avian associated E. coli from the Punjab region in Pakistan found almost universal resistance to ampicillin (98.6%) while we found 79.3% (73/92) (Azam et al., 2019). In Egypt, one analysis of 116 avian associated E. coli isolates found 100% (116/116) resistance to ampicillin (Awad et al., 2016). A study of retail poultry products from the United States found that E. coli originating from turkey products had an ampicillin resistance rate (62%) higher than E. coli originating from chicken products (20%) (Davis et al., 2018). Thankfully a number of clinically relevant antibiotics (meropenem, imipenem, cefotetan, piperacillin-tazobactam, and amikacin) had no phenotypic resistance observed. This result is consistent with a previous report of environmental E. coli isolates from Japan, described by the authors as harboring a virulence gene profile similar to isolates associated with avian infections, which were 100% susceptible to carbapenems and aminoglycosides (Hayashi et al., 2019).

A limitation of this study was that with short Illumina reads, we were unable to unequivocally implicate ARGs or virulence genes as present on mobilizable plasmids. As we did not sequence isolates originating from the intestinal contents of these chickens, we are not able to discuss the similarity and differences between the liver-borne isolates and those from commensal sites. Additionally, we do not have access to a chicken model of infection to demonstrate links between comparative pathogenicity of the strains and their virulence gene mosaic. Given the limitation of clinical breakpoints some isolates that possess genetic determinants of resistance may test phenotypically susceptible. Since we tested multiple antibiotics in the same antibiotic classes there is not always a clear relation between phenotypic resistance and genotypic presence of a resistant determinant. As we used disk diffusion to assess antimicrobial susceptibility instead of quantitative broth microdilution (except in cases of mcr-1 positive isolates) we are not able to report MIC values. In summation, we assembled a cohort of chicken associated E. coli isolates obtained from multiple provinces in Pakistan. Genomic analysis of these isolates identified that most of the cohort are in the B1 and A clades and harbor a mosaic of ARGs and virulence genes which may complicate treatment in a human infection.

Data Availability Statement

All genomes sequenced in this study have been uploaded to the NCBI WGS database associated with BioProject PRJNA522294.

Author Contributions

MR and RP sequenced the isolates, performed the in silico analysis, and wrote the manuscript. AF performed the in silico analysis and helped to write the manuscript. MW carried out the culture work and susceptibility testing. AR, AA, NS, and MA collected and stored the samples. AD helped to write the manuscript and generate the figures. C-AB, NA, and GD devised the study.

Funding

This work was supported in part by an International Research Support Initiative Program (IRSIP) award of the Higher Education Commission (HEC) of Pakistan to MR, and an Awards through the National Institute of Allergy and Infectious Diseases, and the Eunice Kennedy Shriver National Institute of Child Health & Human Development, of the National Institutes of Health under award numbers R01AI123394 and R01HD092414 to GD. RP received support from an NIGMS training grant through award T32 GM007067 (PI: James Skeath) and the Bayer/Monsanto Excellence Fund graduate fellowship. AF received support from the Chancellor’s Graduate Research Fellowship Program at Washington University in St. Louis. AD received support from the Institutional Program Unifying Population and Laboratory-Based Sciences Burroughs Welcome Fund grant to Washington University. The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding agencies.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank The Edison Family Center for Genome Sciences & Systems Biology staff Eric Martin, Brian Koebbe, Jessica Hoisington-López, and MariaLynn Crosby for their technical support and sequencing expertise.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.03052/full#supplementary-material

FIGURE S1 | Strong concordance between disk diffusion relationships. Scatter plots depicting the disk diffusion relationship for the 3^rd generation cephalosporins ceftriaxone and ceftiofur (A) and quinolones ciprofloxacin and norfloxacin (B).

TABLE S1 | Genomes analyzed in this study.

TABLE S2 | Publicly available genomes used in this study.

TABLE S3 | Antimicrobial susceptibility testing information.

TABLE S4 | Full antimicrobial susceptibility testing results.

Footnotes

1. ^ https://itol.embl.de/

References

Akhtar, F., Rabbani, M., Muhammad, K., Younus, M., Ghafoor, A., Sheikh, A. A., et al. (2016). Comparative antibiotic resistance profile of the multidrug resistant E. coli isoalted from commercial and backyard poultry. J. Anim. Plant Sci. 26, 1628–1632.

Google Scholar

Awad, A., Arafat, N., and Elhadidy, M. (2016). Genetic elements associated with antimicrobial resistance among avian pathogenic Escherichia coli. Ann. Clin. Microbiol. Antimicrob. 15, 59. doi: 10.1186/s12941-016-0174-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Azam, M., Mohsin, M., Sajjad-Ur-Rahman, and Saleemi, M. K. (2019). Virulence-associated genes and antimicrobial resistance among avian pathogenic Escherichia coli from colibacillosis affected broilers in Pakistan. Trop. Anim. Health Prod. 51, 1259–1265. doi: 10.1007/s11250-019-01823-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bankevich, A., Nurk, S., Antipov, D., Gurevich, A. A., Dvorkin, M., Kulikov, A. S., et al. (2012). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477. doi: 10.1089/cmb.2012.0021

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbieri, N. L., Nielsen, D. W., Wannemuehler, Y., Cavender, T., Hussein, A., Yan, S. G., et al. (2017). mcr-1 identified in avian pathogenic Escherichia coli (APEC). PLoS One 12:e0172997. doi: 10.1371/journal.pone.0172997

PubMed Abstract | CrossRef Full Text | Google Scholar

Baym, M., Kryazhimskiy, S., Lieberman, T. D., Chung, H., Desai, M. M., and Kishony, R. (2015). Inexpensive multiplexed library preparation for megabase-sized genomes. PLoS One 10:e0128036. doi: 10.1371/journal.pone.0128036

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Carattoli, A., Zankari, E., García-Fernández, A., Voldby Larsen, M., Lund, O., Villa, L., et al. (2014). In silico detection and typing of plasmids using PlasmidFinder and plasmid multilocus sequence typing. Antimicrob. Agents Chemother. 58, 3895–3903. doi: 10.1128/AAC.02412-14

PubMed Abstract | CrossRef Full Text | Google Scholar

CDDEP (2017). Situation Analysis Report on Antimicrobial Resistance in Pakistan: Findings and Recommendations for Antibiotics Use and Resistance. Washington, DC: Global Antibiotic Resistance Partnership.

Google Scholar

Cieza, R. J., Hu, J., Ross, B. N., Sbrana, E., and Torres, A. G. (2015). The IbeA invasin of adherent-invasive Escherichia coli mediates interaction with intestinal epithelia and macrophages. Infect. Immun. 83, 1904–1918. doi: 10.1128/IAI.03003-14

PubMed Abstract | CrossRef Full Text | Google Scholar

CLSI (2018). Performance Standards for Antimicrobial Disk Susceptibility Tests. Wayne, PA: Clinical and Laboratory Standards Institute.

Google Scholar

Collingwood, C., Kemmett, K., Williams, N., and Wigley, P. (2014). Is the concept of avian pathogenic Escherichia coli as a single pathotype fundamentally flawed? Front. Vet. Sci. 1:5. doi: 10.3389/fvets.2014.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cordoni, G., Woodward, M. J., Wu, H., Alanazi, M., Wallis, T., and La Ragione, R. M. (2016). Comparative genomics of European avian pathogenic E. Coli (APEC). BMC Genomics 17:960. doi: 10.1186/s12864-016-3289-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Croxen, M. A., and Finlay, B. B. (2009). Molecular mechanisms of Escherichia coli pathogenicity. Nat. Rev. Microbiol. 8, 26–38. doi: 10.1038/nrmicro2265

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, G. S., Waits, K., Nordstrom, L., Grande, H., Weaver, B., Papp, K., et al. (2018). Antibiotic-resistant Escherichia coli from retail poultry meat with different antibiotic use claims. BMC Microbiol. 18:174. doi: 10.1186/s12866-018-1322-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Delannoy, S., Le Devendec, L., Jouy, E., Fach, P., Drider, D., and Kempf, I. (2017). Characterization of colistin-resistant Escherichia coli isolated from diseased pigs in france. Front. Microbiol. 8:2278. doi: 10.3389/fmicb.2017.02278

PubMed Abstract | CrossRef Full Text | Google Scholar

Dho-Moulin, M., and Fairbrother, J. M. (1999). Avian pathogenic Escherichia coli (APEC). Vet. Res. 30, 299–316.

Google Scholar

Edgar, R., and Bibi, E. (1997). MdfA, an Escherichia coli multidrug resistance protein with an extraordinarily broad spectrum of drug recognition. J. Bacteriol. 179, 2274–2280. doi: 10.1128/jb.179.7.2274-2280.1997

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurevich, A., Saveliev, V., Vyahhi, N., and Tesler, G. (2013). QUAST: quality assessment tool for genome assemblies. Bioinformatics 29, 1072–1075. doi: 10.1093/bioinformatics/btt086

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayashi, W., Tanaka, H., Taniguchi, Y., Iimura, M., Soga, E., Kubo, R., et al. (2019). Acquisition of mcr-1 and cocarriage of virulence genes in avian pathogenic Escherichia coli isolates from municipal wastewater influents in Japan. Appl. Environ. Microbiol. 85:e01661-19. doi: 10.1128/AEM.01661-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoerr, F. J. (2010). Clinical aspects of immunosuppression in poultry. Avian Dis. 54, 2–15.

Google Scholar

Huja, S., Oren, Y., Trost, E., Brzuszkiewicz, E., Biran, D., Blom, J., et al. (2015). Genomic avenue to avian colisepticemia. MBio 6:e01681–14. doi: 10.1128/mBio.01681-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, J., Rabbani, I., Aslam, S., and Ahmad, H. A. (2015). An overview of poultry industry in Pakistan. Worlds Poul. Sci. J. 71, 689–700. doi: 10.1017/s0043933915002366

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutton, T. A., Innes, G. K., Harel, J., Garneau, P., Cucchiara, A., Schifferli, D. M., et al. (2018). Phylogroup and virulence gene association with clinical characteristics of Escherichia coli urinary tract infections from dogs and cats. J. Vet. Diagn. Invest. 30, 64–70. doi: 10.1177/1040638717729395

PubMed Abstract | CrossRef Full Text | Google Scholar

Idrees, M., Shah, M. A., Michael, S., Qamar, R., and Bokhari, H. (2011). Antimicrobial resistant Escherichia coli strains isolated from food animals in Pakistan. Pak. J. Zool. 43, 303–310.

Google Scholar

Johnson, T. J., Wannemuehler, Y., Johnson, S. J., Stell, A. L., Doetkott, C., Johnson, J. R., et al. (2008). Comparison of extraintestinal pathogenic Escherichia coli strains from human and avian sources reveals a mixed subset representing potential zoonotic pathogens. Appl. Environ. Microbiol. 74, 7043–7050. doi: 10.1128/AEM.01395-08

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, T. J., Wannemuehler, Y. M., Johnson, S. J., Logue, C. M., White, D. G., Doetkott, C., et al. (2007). Plasmid replicon typing of commensal and pathogenic Escherichia coli isolates. Appl. Environ. Microbiol. 73, 1976–1983. doi: 10.1128/aem.02171-06

PubMed Abstract | CrossRef Full Text | Google Scholar

Keane, O. M. (2016). Genetic diversity, the virulence gene profile and antimicrobial resistance of clinical mastitis-associated Escherichia coli. Res. Microbiol. 167, 678–684. doi: 10.1016/j.resmic.2016.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleinheinz, K. A., Joensen, K. G., and Larsen, M. V. (2014). Applying the ResFinder and VirulenceFinder web-services for easy identification of acquired antibiotic resistance and E. coli virulence genes in bacteriophage and prophage nucleotide sequences. Bacteriophage 4:e27943. doi: 10.4161/bact.27943

PubMed Abstract | CrossRef Full Text | Google Scholar

Letunic, I., and Bork, P. (2007). Interactive tree of life (iTOL): an online tool for phylogenetic tree display and annotation. Bioinformatics 23, 127–128. doi: 10.1093/bioinformatics/btl529

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Chen, L., Wu, X., and Huo, S. (2015). Molecular characterization of multidrug-resistant avian pathogenic Escherichia coli isolated from septicemic broilers. Poult. Sci. 94, 601–611. doi: 10.3382/ps/pev008

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y. Y., Wang, Y., Walsh, T. R., Yi, L. X., Zhang, R., Spencer, J., et al. (2016). Emergence of plasmid-mediated colistin resistance mechanism MCR-1 in animals and human beings in China: a microbiological and molecular biological study. Lancet Infect. Dis. 16, 161–168. doi: 10.1016/S1473-3099(15)00424-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Loytynoja, A. (2014). Phylogeny-aware alignment with PRANK. Methods Mol. Biol. 1079, 155–170. doi: 10.1007/978-1-62703-646-7_10

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, J., Mohsin, M., Lei, S., Srinivas, S., Wiqar, R. T., Lin, J., et al. (2018). Discovery of a mcr-1-bearing plasmid in commensal colistin-resistant Escherichia coli from healthy broilers in Faisalabad, Pakistan. Virulence 9, 994–999. doi: 10.1080/21505594.2018.1462060

PubMed Abstract | CrossRef Full Text | Google Scholar

Maamar, E., Alonso, C. A., Hamzaoui, Z., Dakhli, N., Abbassi, M. S., Ferjani, S., et al. (2018). Emergence of plasmid-mediated colistin-resistance in CMY-2-producing Escherichia coli of lineage ST2197 in a Tunisian poultry farm. Int. J. Food Microbiol. 269, 60–63. doi: 10.1016/j.ijfoodmicro.2018.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Magiorakos, A. P., Srinivasan, A., Carey, R. B., Carmeli, Y., Falagas, M. E., Giske, C. G., et al. (2012). Multidrug-resistant, extensively drug-resistant and pandrug-resistant bacteria: an international expert proposal for interim standard definitions for acquired resistance. Clin. Microbiol. Infect. 18, 268–281. doi: 10.1111/j.1469-0691.2011.03570.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Manges, A. R. (2016). Escherichia coli and urinary tract infections: the role of poultry-meat. Clin. Microbiol. Infect. 22, 122–129. doi: 10.1016/j.cmi.2015.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellata, M., Dho-Moulin, M., Dozois, C. M., Curtiss, R., Lehoux, B., and Fairbrother, J. M. (2003). Role of avian pathogenic Escherichia coli virulence factors in bacterial interaction with chicken heterophils and macrophages. Infect. Immun. 71, 494–503. doi: 10.1128/iai.71.1.494-503.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, A., Viso, S., Lopez, C., Alonso, M. P., Garcia-Garrote, F., Dabhi, G., et al. (2013). Poultry as reservoir for extraintestinal pathogenic Escherichia coli O45:K1:H7-B2-ST95 in humans. Vet. Microbiol. 167, 506–512. doi: 10.1016/j.vetmic.2013.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Page, A. J., Cummins, C. A., Hunt, M., Wong, V. K., Reuter, S., Holden, M. T., et al. (2015). Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics 31, 3691–3693. doi: 10.1093/bioinformatics/btv421

PubMed Abstract | CrossRef Full Text | Google Scholar

Paudel, S., Stessl, B., Hess, C., Zloch, A., and Hess, M. (2016). High genetic diversity among extraintestinal Escherichia coli isolates in pullets and layers revealed by a longitudinal study. BMC Vet. Res. 12:221. doi: 10.1186/s12917-016-0859-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Potter, R. F., Lainhart, W., Twentyman, J., Wallace, M. A., Wang, B., Burnham, C. A., et al. (2018a). Population structure, antibiotic resistance, and uropathogenicity of Klebsiella variicola. MBio 9:e02481–18. doi: 10.1128/mBio.02481-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Potter, R. F., Wallace, M. A., McMullen, A. R., Prusa, J., Stallings, C. L., Burnham, C. A. D., et al. (2018b). blaIMP-27 on transferable plasmids in Proteus mirabilis and Providencia rettgeri. Clin. Microbiol. Infect. 24, 1019.e5–1019.e8. doi: 10.1016/j.cmi.2018.02.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, M. N., Dehal, P. S., and Arkin, A. P. (2010). FastTree 2–approximately maximum-likelihood trees for large alignments. PLoS One 5:e9490. doi: 10.1371/journal.pone.0009490

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahman, S., and Mohsin, M. (2019). The under reported issue of antibiotic-resistance in food-producing animals in Pakistan. Pak. Vet. J. 5:e9490.

Google Scholar

Richter, S. S., Sercia, L., Branda, J. A., Burnham, C. A. D., Bythrow, M., Ferraro, M. J., et al. (2013). Identification of Enterobacteriaceae by matrix-assisted laser desorption/ionization time-of-flight mass spectrometry using the VITEK MS system. Eur. J. Clin. Microbiol. Infect. Dis. 32, 1571–1578.

Google Scholar

Sabri, M., Léveillé, S., and Dozois, C. M. (2006). A SitABCD homologue from an avian pathogenic Escherichia coli strain mediates transport of iron and manganese and resistance to hydrogen peroxide. Microbiology 152, 745–758. doi: 10.1099/mic.0.28682-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarowska, J., Futoma-Koloch, B., Jama-Kmiecik, A., Frej-Madrzak, M., Ksiazczyk, M., Bugla-Ploskonska, G., et al. (2019). Virulence factors, prevalence and potential transmission of extraintestinal pathogenic Escherichia coli isolated from different sources: recent reports. Gut Pathog. 11, 10. doi: 10.1186/s13099-019-0290-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmieder, R., and Edwards, R. (2011). Quality control and preprocessing of metagenomic datasets. Bioinformatics 27, 863–864. doi: 10.1093/bioinformatics/btr026

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreiber, H. L. IV, Conover, M. S., Chou, W. C., Hibbing, M. E., Manson, A. L., Dodson, K. W., et al. (2017). Bacterial virulence phenotypes of Escherichia coli and host susceptibility determine risk for urinary tract infections. Sci. Transl. Med. 9:eaaf1283.

Google Scholar

Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069. doi: 10.1093/bioinformatics/btu153

PubMed Abstract | CrossRef Full Text | Google Scholar

Shini, S., Huff, G. R., Shini, A., and Kaiser, P. (2010). Understanding stress-induced immunosuppression: exploration of cytokine and chemokine gene profiles in chicken peripheral leukocytes 1. Poult. Sci. 89, 841–851. doi: 10.3382/ps.2009-00483

PubMed Abstract | CrossRef Full Text | Google Scholar

Subedi, M., Luitel, H., Devkota, B., Bhattarai, R. K., Phuyal, S., Panthi, P., et al. (2018). Antibiotic resistance pattern and virulence genes content in avian pathogenic Escherichia coli (APEC) from broiler chickens in Chitwan, Nepal. BMC Vet. Res. 14:113. doi: 10.1186/s12917-018-1442-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Boeckel, T. P., Brower, C., Gilbert, M., Grenfell, B. T., Levin, S. A., Robinson, T. P., et al. (2015). Global trends in antimicrobial use in food animals. Proc. Natl. Acad. Sci. U.S.A. 112, 5649–5654. doi: 10.1073/pnas.1503141112

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y. Q., Li, Y. X., Song, T., Yang, Y. X., Jiang, W., Zhang, A. Y., et al. (2017). Colistin resistance gene mcr-1 and its variant in Escherichia coli isolates from chickens in China. Antimicrob. Agents Chemother. 61:e01204–16

Google Scholar

Zankari, E., Allesoe, R., Joensen, K. G., Cavaco, L. M., Lund, O., and Aarestrup, F. M. (2017). PointFinder: a novel web tool for WGS-based detection of antimicrobial resistance associated with chromosomal point mutations in bacterial pathogens. J. Antimicrob. Chemother. 72, 2764–2768. doi: 10.1093/jac/dkx217

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu Ge, X., Jiang, J., Pan, Z., Hu, L., Wang, S., Wang, H., et al. (2014). Comparative genomic analysis shows that avian pathogenic Escherichia coli isolate IMT5155 (O2:K1:H5; ST complex 95, ST140) shares close relationship with ST95 APEC O1:K1 and human ExPEC O18:K1 strains. PLoS One 9:e112048. doi: 10.1371/journal.pone.0112048

PubMed Abstract | CrossRef Full Text | Google Scholar