High prevalence of mcr-1-encoded colistin resistance in commensal Escherichia coli from broiler chicken in Bangladesh

Colistin is a last-resort antimicrobial used for the treatment of human infections caused by multidrug-resistant Gram-negative bacteria. However, colistin is still widely used in intensive poultry production in Bangladesh. We aimed to investigate the dynamics and genetic diversity of colistin-resistant commensal Escherichia coli from broiler chickens. A total of 1200 E. coli strains were characterized from 20 broiler farms at three-time points along the production period. All strains were screened for mcr-1 to mcr-5 genes by a multiplex PCR, and their genetic diversity was measured by repetitive extragenic palindromic (REP)-PCR fingerprinting. Genomic diversity and characterization were performed by whole genome sequencing (WGS). Twenty-five percent of the commensal E. coli strains harbored mcr-1 genes. Frequency of mcr-1 gene detection correlated positively (odds ratio 1.71; 95% CI 0.96–3.06; p = 0.068) with the use of colistin in poultry flocks. REP-PCR profiles and WGS analysis showed diverse E. coli population carrying multiple antimicrobial resistance genes. Phylogenetic comparison of mcr-1-bearing strains recovered from this study with a global strain collection revealed wide phylogenetic relationship. This study identified a high prevalence of mcr-1 gene among genetically diverse E. coli populations from broiler chickens in Bangladesh suggesting a massive horizontal spread of mcr-1 rather than by clonal expansion.


Scientific Reports
| (2020) 10:18637 | https://doi.org/10.1038/s41598-020-75608-2 www.nature.com/scientificreports/ In Bangladesh, small-scale broiler farms are the major source of poultry meat. Colistin sulphate is used massively for the treatment and prevention of diseases in this production system 9 . Owing to the absence of strict regulation for AMU in Bangladesh, farmers can acquire colistin without prescription from a registered veterinarian 9 . To date, there is no systematic investigation on the colistin resistance commensal E. coli population in animals in this production system where imprudent use of antimicrobial is very frequent. Moreover, the overall genetic diversity of the commensal E. coli population in broiler chicken over the entire production stage has not been investigated before.
In this study, the genetic diversity of commensal E. coli from broiler chicken in Bangladesh was thoroughly investigated throughout the production period, and the prevalence of plasmid-encoded colistin resistance genes (mcr) was determined. The objectives of the current study were to investigate the genetic diversity of commensal E. coli in poultry, to explore the distribution and genetic background of colistin-resistant strains among the commensal E. coli population, and to characterize the phylogenetic relationship of the colistin-resistant commensal E. coli strains in the global population structure of E. coli.

Results
Farms, sampling and E. coli isolates. A total of 20 commercial broiler chicken farms at Chattogram division in Bangladesh were investigated in this study. Farms were located across two administrative districts of the division with a minimum and maximum distance between farms of 0.5 km and 80 km, respectively. All the birds were raised in small-scale intensive systems of rearing. The farm size varies from 800 to 4000 birds per farm (median size 1060 birds/farm). Most farms used antimicrobials for treatment and prophylaxis purpose. No farms used antimicrobials for the purpose of growth promotor. A detail list of different classes of antimicrobials used during the production is given in Supplementary Table 1.
One pooled faecal sample were collected in each of the three sampling times at day1, day15, and day28 of the production which comprises 60 pooled samples from 20 farms.
Quantification of E. coli from faecal specimens shown that the average count of E. coli was significantly (p < 0.0001) higher at day1 (8.7 ± 0.43 log 10 CFU/gm faeces) compared with day15 (7.6 ± 0.74 log 10 CFU/gm faeces) and day28 (7.4 ± 0.55 log 10 CFU/gm faeces) (Fig. 1a). A total of 20 confirmed E. coli isolates from each faecal sample, comprising an overall total collection of 1200 isolates, were characterized in this study.
A high prevalence of mcr-1 gene was detected in the commensal E. coli population. We identified 305 (25%, 95% CI, 23-28%) mcr-1-positive E. coli out of the 1200 isolates. No strains were positive for the other genes (mcr-2 -mcr- 5) investigated. The prevalence of the mcr-1 gene in the E. coli population was significantly higher (p < 0.0001) at day15 (45%) compared with day28 (25%) and day1 (6%) in the production (Fig. 1b). The distribution of mcr-1-bearing E. coli isolates across the farms at three sampling times was highly diverse (Fig. 1c). Colistin-resistant E. coli could be detected in all farms at least one sampling time. Seven variables were tested in the univariable analysis (Table 1). Of them, two parameters: "antimicrobial use" and "colistin use" were identified as eligible to be included for the multivariable analysis. In the final model "colistin use" was found to Core and pangenome comparison. The degree of genomic flexibility of the 32 mcr-1-positive commensal E. coli strains was assessed by comparing the core and pangenome structure of the isolates. The overall pangenome consisted of 13,132 genes which is four times larger than the core genome of the same strains (2737 genes (20.84%) (Fig. 2a). The accessory gene pool was highly variable among the isolates (N = 10,395). Of the genes identified by Roary, 50.88% (n = 6681), 25.66% (n = 3370), and 2.62% (n = 344) were found in less than 15% of the isolates (referred to as cloud genes), between 15 and 95% of the isolates (referred as shell genes), and between 95 and 99% of the isolates (referred as soft core genes), respectively (Fig. 2a). We compared the genome of mcr-1-positive E. coli isolates based on the presence or absence a gene which is depicted in Fig. 2b  www.nature.com/scientificreports/ comparison based on the concatenated core gene alignment also showed large genomic diversity among the mcr-1-positive commensal E. coli isolates (Fig. 2c).
Virulence genes harboured by commensal E. coli. In Total, 13 virulence genes were detected in the genomes of the 32 commensal E. coli (Supplementary Table 4). The most frequent virulence determinants were astA (EAST-1 heat-stable toxin) and iss (Increased Serum Survival) which were found in 50% and 44% of the isolates, respectively. A major proportion of the isolates (53%) harboured multiple virulence genes. Notably, seven of the 32 isolates harboured at least 4-6 virulence determinants, but none of the strains carried combinations of virulence genes known to be characteristic for pathogenic subtypes.
In silico MLSTs, serotypes, and plasmid genotyping. We identified 16 different sequence types (STs) among the 32 E. coli isolates. The most frequent ST was ST43 (n = 6) followed by ST4965 (n = 5). We could not infer the ST types of three isolates (listed as unknown types in Supplementary Table 4). In silico serotyping identified 11 complete serotypes. The most frequent serotype was O13:H30 (n = 6) followed by O6:H10 (n = 4). We could not identify the O antigen specific serotype of 11 of the isolates. However, based on the H type, three isolates were O-:H10 and two were O-:H21 and seven were singletons. We identified 31 plasmid replicons in the 32 isolates. All the isolates harboured multiple plasmid replicons ranging from four to 11. The most abundant replicon type was ColRNAI which was found in 24 isolates followed by IncFIB (AP001918) in 23 isolates. Other dominant plasmid replicon types were IncHI2 (n = 19), IncHI2A (n = 19), IncN (n = 16), IncX1 (n = 15), and IncI2 (n = 13).
High prevalence of antimicrobial resistance determinants. We identified a total of 42 different antimicrobial resistance determinants in the 32 strains (Fig. 3a). The strain selection was based on the presence of mcr-1 and in accordance, this gene was detected in all 32 strains. The most abundant antimicrobial resistance gene was mdf(A) (31/32) followed by tet(A) (23/32) and bla TEM-1B (23/32) which confer resistance to a broad spectrum of antimicrobials, tetracycline, and beta-lactam antibiotics, respectively. All 32 strains carried multiple acquired antimicrobial resistance genes with a minimum of five to a maximum of 20 genes (Fig. 3a). According to the database used, more than 50% of the strains (n = 18) harboured at least 13 resistance genes. Twenty-seven of the isolates exhibited at least one known point mutation and 15 isolates exhibited at least two point mutations in any of the three-resistance determinants (gyrA, parC, and parE) which may confer resistance to nalidixic acid and ciprofloxacin (Fig. 3b). The insertion sequence ISApl1 was found in 30 of the 32 mcr-1-bearing E. coli strains (Supplementary Table 4). The other two had no ISApl1 insertion sequence in their genomes. A total of eight, fifteen, and six of the isolates contained one, two, and three ISApl1 insertion sequences, respectively.
Presence of chromosomal point mutations in genes previously associated with colistin resistance. Since MIC varied considerably between mcr-1-positive isolates, we carried out a search for point muta- Phylogenetic position of mcr-1-positive E. coli at the global population structure. The WGS of mcr-1-positive strains identified in the current study were compared to a global collection originated from Cambodia (6 human isolates), Canada (one human isolate), China (two environment, 16 chicken, 33 human, and 15 pig isolates), Denmark (one human and five chicken isolates), Japan (three cattle isolates), the Netherlands (three chicken isolates), Singapore (five human isolates), Thailand (three isolates from unknown source), the USA (one human isolate), and Vietnam (32 isolates from poultry farms and farmers) (Supplementary File 1). A total of 93,130 core SNPs was detected among the 158 isolates. The pairwise comparison of SNP between isolates varied from 1 to 46,621 SNPs. A total of six and five E. coli isolates from this study, respectively, formed two separate clades. However, most of the strains of this study were spread across the phylogenetic tree (Fig. 4).

Discussion
Colistin sulphate is widely used for treatment and for the prevention of diseases in poultry in Bangladesh 9 . Poultry farmers in Bangladesh usually administer colistin or other antimicrobials through the water, and often in the absence of disease symptoms. Once they start using any antimicrobials treatment, they continuously keep it in water for a certain period of time ranging from three to seven days. Little is known on how this imprudent use of antimicrobial over the decades have affected the genetic diversity of commensal gut bacteria in the poultry.
Our study detected a diverse population of commensal E. coli in broiler chicken according to genetic typing. The genetic diversity of E. coli increased with the age of the birds. Previous studies [20][21][22] have also shown the association of age of host (e.g., pigs) with the diversity of commensal E. coli. The increased diversity in the older age group of chicken might be due to the establishment and persistence of adapted strains in the gut over time.
Apart from the wide genetic diversity, we detected a high prevalence of the mcr-1 gene (25%) in the commensal E. coli population. Although the study design and the host species were different, a similar alarming high frequency of mcr-1 was detected in E. coli from human and livestock sources in Thailand 23 A242T  A360V  D283G  H2R  M287I  S138N  T119P  T246I  V351I  V88I  Y358N   S001  S002  S003  S004  S005  S024  S025  S026  S027  S028  S029  S030  S031  S032  S033  S034  S035  S036  S037  S038  S039  S040  S041  S042  S043  S044  S045  S046  S047  S048  S049  S050 pmrA pmrB www.nature.com/scientificreports/ higher than previously reported studies from Germany (3.8%) 25 , China (1% or less) and other countries [26][27][28] . Also, this study shows high MIC against colistin (≥ 8 mg/L) in more than 50% of the mcr-1 carrying E. coli which is very similar with another study in Bangladesh 29 , and the overall MIC is higher than other studies in some European countries 25,30,31 . It suggests that the commensal E. coli population in livestock in countries with a high and unregulated use of colistin constitutes a hotspot for selection for carriage of transferrable resistance. This is worrisome and calls for immediate action to reduce the preventive use of antimicrobials in general and colistin in particular. The study showed that the occurrence of the mcr-1 gene increased significantly as the birds grew older, most likely simply reflecting the selection pressure. In accordance with this, colistin use is positively associated with the isolation frequency of the mcr-1 gene in this study. This association should be cautiously interpreted as the selection and spread of AMR gene is a complex biological process. Although the selection of mcr-1 gene is primarily influenced under the pressure of colistin, other antimicrobials can also co-select mcr-1 gene 32 . Additionally, many factors, such as habitat, co-selection by metals, and biocides, can influence the selection of antimicrobial resistance genes 33 . There is no veterinary monitoring of antimicrobial usage for the studied farms, therefore, we mainly relied on the data provided by the farmers which should also be considered during interpreting this association. The presence of mcr-1 gene in the E. coli population in the day-old broiler chicks raises questions about the origin of the colistin resistance strains in such early age of birds. Vertical transmission of colistin resistance E. coli from broiler breeders to their offspring could be one explanation for the presence of mcr-1 gene in the day-old chicks as indicated by previous studies 34, 35 . The pangenome analysis of 32 mcr-1-bearing commensal E. coli strains showed a highly flexible accessory genome. The E. coli core genome of this study is four times smaller than the pangenome. Previous studies have also shown a smaller size core genome of E. coli [36][37][38] . However, the core genome size is a comparative measurement because the core would be reduced when more genomes are added to the comparison 39   www.nature.com/scientificreports/ Almost all the isolates (91%) harboured at least one and several isolates harboured multiple virulence genes. The most frequent virulence gene was EAST-1 heat-stable toxin (astA) and the increased serum survival (iss) was the second most frequent virulence marker among the commensal strains. The astA gene is commonly found in enterotoxigenic E. coli (ETEC) 40 . ETEC is one of the leading causes of severe diarrhoea in young children in developing countries like Bangladesh 41 . Furthermore, increased serum survival gene iss is a well-known virulence marker of extraintestinal pathogenic E. coli (ExPEC) in avian species 42 , however, the presence of this gene was not found in combination with other virulence markers of avian pathogenic E. coli (APEC). Additionally, the dominant ST type among the isolates was ST43 (n = 6). This type is often found in human clinical infections as carbapenemase-producing E. coli 43,44 . The presence of virulence markers in the commensal E. coli strains isolated from poultry that have genetic similarity with diverse pathogenic clones indicates their potential of transferring virulence genes to pathogenic clones of E. coli in humans and animals.
Our study revealed that all the mcr-1-positive E. coli isolates carry multiple ARGs. The presence of tetracycline resistance genes in more than 70% of the mcr-1-bearing isolates indicated an alarming spread of tetracycline resistance genes among the commensal population. Although tetracycline is less commonly used in broiler production in Bangladesh 9 , its increased resistance in commensal E. coli could be an indication of resistance selection by a bystander effect of other antimicrobials 45 . A previous study from Bangladesh has identified tetracycline resistance genes in multidrug resistance E. coli strains isolated from children 46 . Earlier studies identified a varied prevalence of tetracycline resistance genes in E. coli isolated from both humans and animals 47 . One of the strains harboured extended-spectrum β-lactamase-producing gene bla CTX-M-55 . Co-harbouring of mcr-1 and bla CTX-M-55 gene was previously reported in the pyelonephritis case in a three-year child in France, caused by E. coli 48 . A combination of a mutation in the gyrA gene and the parC gene was detected in several isolates suggesting resistance to fluroquinolones 49,50 . The presence of single mutations in the gyrA gene altering serine 83 to leucine (S83L) and the double mutations in the gyrA gene altering aspartic acid (D87N) and serine (S83L) is known to confer high resistance capability of E. coli strains to fluroquinolones 51 .
All the 32 isolates harboured multiple plasmid replicons. The abundant plasmid types such as ColRNAI, IncFIB, IncHI2, and IncI2 were previously shown as mcr-1-carrying plasmid replicon 52,53 . The mcr-1 gene mobilizing transposon component ISApl1 was detected in all the isolates except two. Evidence shows that ISApl1 plays a pivotal role in the spread of mcr-1 gene 54 . In some cases, one or both copies of ISApl1 can be lost, however, a single copy of upstream ISApl1 is capable of mobilizing mcr-1 genes 55 .
Phylogenetic comparison of the strains identified in the present study with an international collection of colistin-resistant E. coli showed that the mcr-1-bearing E. coli strains of this study were highly diverse. Only two small clusters of the isolates from this study were observed; these encompassed strains of ST43 (n = 6) and ST4965 (n = 5), respectively. All other isolates were found to be diffusely distributed across the tree according to the SNP phylogeny. Similarly, a high level of phylogenetic diversity of mcr-1-bearing E. coli was reported previously 53 . The scattering of mcr-1-bearing E. coli into genetically diverse strains suggests that the transfer of the gene between strains is a frequent event in the gut of broilers in the poultry farms in Bangladesh.
This study has some limitations. First, the study was conducted in a short time frame which did not allow us to observe if there is any seasonal variation on the AMR and genetic diversity of E. coli. Second, the study was conducted in a certain geographical region in Bangladesh. The addition of more farms from different regions of the country would bring a more detailed picture of the E. coli diversity. Third, only broiler chickens were investigated in this study, whereas layer and backyard chickens are also a major part of poultry production in Bangladesh. Finally, although the study identified the mcr-1 gene in E. coli from day-old chicks, however, it was not possible to clarify whether the source of the resistant-strains was a vertical transmission from breeder flocks or not.
In conclusion, this first population level study on commensal E. coli and colistin resistance in broiler chicken encompassing different stages of production has revealed a high genetic diversity of commensal E. coli in broiler chicken in Bangladesh. The study also revealed a high frequency of mcr-1 carriage among the commensal E. coli population and showed that the carriage was positively associated with colistin use in chicken production. The presence of mcr-1 gene in a diverse E. coli population suggests a massive horizontal spread of mcr-1 gene rather than clonal expansion. These results call for immediate action from the policy makers to stop imprudent use or to actively regulate the rational use of colistin along with other critically important antimicrobials in poultry production in Bangladesh.

Materials and methods
Study farms and sampling. A longitudinal study was conducted between August to September 2018.
Commercial broiler chicken farms in the Chattogram division of Bangladesh was chosen for this study. Broiler farm owners were invited to participate in the study. Farms were selected based on the starting date of a new flock during the study period, minimum flock size of 500 birds, and geographically well-representing across the division. Pooled faecal samples from each farm were collected longitudinally at three sampling times: at day1, day15, and day28 of the production. The pooled sample consists of five randomly picked fresh faecal droplets collected from four corners and the centre of the flock. A trained veterinarian collected the samples and brought them to the laboratory on the same day with a proper transporting system in a cooling box. A questionnaire survey was also conducted for epidemiological data on each sampling time. The parameters investigated includes strains of birds, use of antimicrobials, sources of food and water, and farm biosecurity practices. All farm owners were agreed to participate in the study and informed consent was taken before sampling from farms and questionnaire survey on farming data. The animal ethical committee of Chattogram Veterinary and Animal Sciences University (CVASU), Bangladesh approved the study protocol and collection of animal related data by a questionnaire.

PCR for mcr genes and minimum inhibitory concentration (MIC) determination. All isolates
were screened by a multiplex PCR to detect five mcr genes (mcr-1 to mcr-5) as described previously 58  Repetitive extragenic palindromic-PCR (REP-PCR) fingerprinting. The genetic diversity of the E.
coli isolates was determined by REP-PCR. The primers used for REP-PCR fingerprinting were Rep1R-I (5′-III ICG ICG ICA TCI GGC-3′) and Rep2-I (5′-ICG ICT TAT CIG GCC TAC-3′). The PCR was performed according to a previously described protocol 60 . E. coli ATCC 25922 and sterile MilliQ water were used as positive and negative controls, respectively. As an external reference control, GeneRuler 100 bp plus DNA ladder (Thermo Fisher Scientific, USA) was used to standardize the fingerprint profiles. GelJ, an image analyzing program, was used to analyze REP-PCR DNA fingerprints data 61 . The normalization of every gel images was done by 100 bp plus DNA ladder (Thermo Fisher Scientific, USA) as an external reference. Unique REP-type was assigned using the Dice similarity method with more than 90% band similarity and 4% tolerance. The REP-PCR-based genetic diversity of commensal E. coli strains as well as the diversity of colistin resistant E. coli strains were determined by the Shannon diversity index (H´). The following formula 62 was used to calculate the diversity-S denotes the number of unique genotype and p i is the number of isolates sharing the same genotype [i] over the total number of isolates.
Whole genome sequencing (WGS) and analysis. A total of 32 mcr-1-bearing E. coli isolates among the dominant REP-types (REP-type with at least three isolates) were selected for WGS. The genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). The sequencing library was prepared according to the Illumina protocol and paired-end next generation sequencing was performed on the Illumina MiSeq platform (Illumina, San Diego, USA).
Quality checking and de-novo assembly of sequencing reads. The quality of the sequencing reads was tested using fastqc (Galaxy Version 0.72 + galaxy1) 63 in the Galaxy platform 64 . We used "Galaxy Europe Instance" for all other downstream analyses related to the Galaxy platform. The raw reads that passed quality control were denovo assembled into a draft genome using a hybrid assembly pipeline Unicycler (Galaxy Version 0.4.7.0) 65 . The pipeline uses both Illumina reads and long reads to produce complete and accurate assemblies. We used default assembly parameters. Genome assembly quality was assessed using Quast (Galaxy Version 5.0.2 + galaxy0) 66 .
Genome annotation and pangenome analysis. The de-novo assembled genomes were annotated using the prokaryotic genome annotation pipeline (PROKKA, Galaxy Version 1.13) 67 . Core and accessory genome comparison of the 32 isolates isolated from of this current study was performed using the Roary pangenome pipeline (Galaxy Version 3.10.2) 68 . Roary generates a core gene alignment from gff3 files produced by PROKKA annotation. A 99% identity cut-off was used to define the "core" gene among the isolates. A concatenated core gene alignment of each of the core genes of all the isolates was generated. The combined core gene alignment was used to construct a maximum-likelihood phylogenetic tree using RAxML (Galaxy Version 1.0.0) 69 . We used the default 1000 fast bootstraps on the best likelihood tree constructed with the General Time Reversible (GTR) substitution model and with a Gamma rate of correction heterogeneity. The gene presence/absence file generated by the Roary pangenome annotation pipeline and the core gene phylogenetic tree were visualized using a web based interactive visualization tool Phandango 70 . 3) and the sequencing reads of interest. Multiple Snippy outputs were combined into a "core SNPs" alignment using Snippy Core (Galaxy Version 4.3.6). The "core site" represents a common genomic position present in all the genomes. The "core SNPs" alignment was used to build a high-resolution phylogeny (ignoring possible recombination). We reconstructed a maximum likelihood phylogenetic tree using FastTree 73 (Galaxy Version 2.1.10) with GTR + CAT Nucleotide evolution model. The phylogenetic tree was visualized using an interactive webtool iTOL 74 .
In silico typing. Acquired antimicrobial resistance genes and chromosomal mutations were detected using Resfinder version 3.2 on Center for Genomic Epidemiology (CGE) web server 75 . CGE webserver was used for further molecular typing of the isolates. Multi-locus sequence typing (MLST), serotyping, virulence determination, plasmid replicon identification and typing were performed using MLST version 2.0, SerotypeFinder version 2.0, VirulenceFinder version 2.0, PlasmidFinder version 2.0 and pMLST version 2.0, respectively [76][77][78][79] . The presence of the insertion sequence ISApl1 belonged to the IS30 family of transposons was identified using the ISfinder online tool 80 .
Statistical analysis. Epidemiological data were analyzed in R 3.5.1 81 . Univariable analysis followed by multivariable logistic regression analysis was performed to identify possible risk factors associated with the prevalence of the mcr-1 gene in the E. coli population. First, we used the univariable logistic regression analysis to identify potential risk factors to be included in the multivariable analysis. Variables with a p-value of less than 0.1 in the univariable analysis were selected for multivariable analysis. Due to the hierarchical data structure, we used three level logistic regression model. The individual observation (1200 E. coli isolates) were nested within three sampling times (day1, day15, and day28) at level-2 which were in turn nested within twenty farms at level-3. The R package lme4 82 was used for logistic regression analysis.
Ethical approval. The study protocol and a questionnaire to collect animal related data were developed in accordance with relevant guidelines and regulations in Bangladesh which was approved by the animal ethical committee of Chattogram Veterinary and Animal Sciences University (CVASU), Bangladesh (Approval No. CVASU/Dir(R&E)EC/-2019/39(2/6)). Informed consent was taken from all participatory farm owners before sampling and data collection from farms. No animals were handled or harmed in this study as only fecal droppings were collected from the farm floor.

Data availability
The genome sequencing data were submitted to the European Nucleotide Archive (ENA) under the project accession number PRJEB34000. All other data generated or analysed in this study are included in the manuscript and in the Supplementary Materials. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.