In silico Analysis of Virulence Associated Genes in Genomes of Escherichia Coli Strains Causing Colibacillosis in Poultry.

INTRODUCTION: Colibacillosis - the most common disease of poultry, is caused mainly by avian pathogenic Escherichia coli (APEC). However, thus far, no pattern to the molecular basis of the pathogenicity of these bacteria has been established beyond dispute. In this study, genomes of APEC were investigated to ascribe importance and explore the distribution of 16 genes recognised as their virulence factors.
MATERIAL AND METHODS: A total of 14 pathogenic for poultry E. coli strains were isolated, and their DNA was sequenced, assembled de novo, and annotated. Amino acid sequences from these bacteria and an additional 16 freely available APEC amino acid sequences were analysed with the DIFFIND tool to define their virulence factors.
RESULTS: The DIFFIND tool enabled quick, reliable, and convenient assessment of the differences between compared amino acid sequences from bacterial genomes. The presence of 16 protein sequences indicated as pathogenicity factors in poultry resulted in the generation of a heatmap which categorises genomes in terms of the existence and similarity of the analysed protein sequences.
CONCLUSION: The proposed method of detection of virulence factors using the capabilities of the DIFFIND tool may be useful in the analysis of similarities of E. coli and other sequences deriving from bacteria. Phylogenetic analysis resulted in reliable segregation of 30 APEC strains into five main clusters containing various virulence associated genes (VAGs).

Introduction

Extraintestinal pathogenic Escherichia coli (ExPEC) strains cause various infectious diseases in poultry, among which avian colibacillosis is the most common and regarded as the most significant in terms of economic losses worldwide. The significance of avian colibacillosis is deserved mainly because it markedly increases mortality and morbidity of birds in all ages. Moreover, it can indirectly affect human health as ExPEC strains share many common features, regardless of origin (1, 2, 16, 19, 24). However, in spite of the ubiquity of these bacteria, their virulence mechanisms are still not defined or understood sufficiently. Currently, diagnosis of colibacillosis is mainly achieved with serotyping, but designation of a pathogenic strain to a certain serogroup is not reliable. Inconsistency exists because its virulence and serotypes differ according to geographic regions and because there are many avian pathogenic Escherichia coli(APEC) strains that do not belong to any known serogroup (1, 3). Therefore, molecular methods are widely used to identify the pathogens. A lot of genes have been already defined as virulence factors in infected poultry, among which toxins, adhesins, or iron-related genes can be distinguished (12, 14). It is known that adhesins are important in extra-intestinal E. coli pathogenesis because once the bacterium binds to host cells, it can progress to its cell infection stage. Genes responsible for iron acquisition mechanisms are also of importance here because the extra-intestinal environment is deficient in this metal and therefore the bacteria need to develop mechanisms to survive under such unfavourable conditions. Additionally, bacteria are protected with genes coding protectins, serum resistance proteins, and outer membrane proteins (1, 17, 28). Despite such a plethora of already described virulence factors, the exact pattern of colibacillosis intensity in poultry remains diagnostically insoluble. There are no clear predictors for the pathogenicity of the strains. They are very diversified according to geographic location and, to make it even more complicated, AFEC strains (avian faecal E. coli) isolated from healthy birds can also contain certain genes coding virulence mechanisms (31). Moreover, the complexity of the molecular basis of pathogenesis makes the task of defining the exact APEC pathotype extremely difficult (4, 19).

The present study, aiming at sequencing and using whole genomes of bacterial strains isolated from poultry with symptoms of colibacillosis, might provide key information for further identification of virulence-associated genes (VAGs) or structural variants involved in the APEC pathotype such as: cvaC, iroN, iucC, iutA, sitA, traT, tsh, feoB, ireA, irp2, hlyD, hlyF, fliC, chuA, iss, and ompT. The availability of APEC genome sequences permits the comparison of genome contents within APEC strains, and helps to understand the evolutionary processes involved in shaping the phenotypes of different E. coli pathotypes (5, 6, 11, 18, 19, 26, 29, 30, 32). The use of next generation sequencing data may result in a better understanding of the evolution and formation of colibacillosis mechanisms and can lead to simplification of molecular fingerprint methods and acceleration of diagnosis in the future.

Material and Methods

Bacterial strains

Previously unreported Escherichia coli strains were isolated from poultry from different geographic regions of Poland. All 14 strains, which were obtained from infected birds (Gallus gallus domesticus and Meleagris gallopavo) with symptoms of colibacillosis, were subjected to detailed analysis. Sampling and isolation procedures were performed on clinically diagnosed individuals: breeder hens (n = 6), broiler breeders (n = 2), commercial layers (n = 2), and turkeys (n = 4). Strains were isolated from yolk sac, marrow, liver, trachea, spleen, and lungs (Table 1). The isolation procedure was performed according to a horizontal method for the detection and enumeration of presumptive Escherichia coli(9). Afterwards, the isolates were cultured in LB medium and grown overnight at 37°C.

Table 1

Summary of the whole genome sequencing and genome assembly results of 14 APEC strains from Poland

Isolate	Animal	Source of material	Scaffolds (n)	Assembly length (Mb)	N50 (bp)	GC (%)	No. of CDS	No. of tRNA	No. of ncRNA	Accession no.	SPAdes version
001PP2015	Gallus gallus domesticus	Yolk-sac	56	4.84	236150	50.67	4744	73	13	MBPT00000000	3.8.0
002PP2015	Gallus gallus domesticus	Marrow	95	5.16	237499	50.37	5121	81	6	MBPU00000000	3.8.0
004PP2015	Meleagris gallopavo	Liver	128	5.39	209164	50.36	5395	86	5	MBPW00000000	3.8.0
007PP2015	Meleagris gallopavo	Trachea	161	5.27	142232	50.54	5266	76	11	LWSE00000000	3.7.1
009PP2015	Meleagris gallopavo	Liver	95	5.12	172385	50.56	5063	75	14	LWUK00000000	3.7.1
011PP2015	Meleagris gallopavo	Liver	99	4.88	157543	50.61	4807	74	11	LXWO00000000	3.7.1
012PP2015	Gallus gallus domesticus	Spleen	112	5.1	224574	50.71	5038	79	10	LXWP00000000	3.7.1
015PP2015	Gallus gallus domesticus	Liver	68	5.08	203619	50.51	5035	82	12	LXYD00000000	3.7.1
016PP2015	Gallus gallus domesticus	Liver	106	5.17	173467	50.56	5124	73	13	LXYC00000000	3.7.1
017PP2015	Gallus gallus domesticus	Lungs	140	5.08	255774	50.51	5046	80	7	LYWM00000000	3.7.1
018PP2015	Gallus gallus domesticus	Liver	110	5.14	268770	50.52	5165	75	6	MBPX00000000	3.7.1
019PP2015	Gallus gallus domesticus	Spleen	72	4.85	204306	50.52	4772	77	13	LYWL00000000	3.7.1
020PP2016	Gallus gallus domesticus	Spleen	111	5.06	127463	50.61	5022	73	11	MBPY00000000	3.7.1
022PP2016	Gallus gallus domesticus	Liver	81	4.73	196459	50.59	4661	76	13	MAYE00000000	3.7.1

MP-PCR differentiation

Genomic DNA was extracted and purified from the bacteria using the Wizard Genomic DNA Purification Kit (Promega, USA) according to the manufacturer’s instructions. Afterwards, the strains were differentiated with the PCR MP (PCR melting profile) method (15). In the first step, 500–1,000 ng of genomic DNA was digested with Hind III restriction enzyme, then, after thermal inactivation of the enzyme, a Hind III adapter was ligated and ligation products were thermally purified. The PCR was performed in the following cycling parameters: 72°C for 7 min, 85°C for 1.5 min, 24 cycles of 85°C for 60 s, 72°C for 120 s, 72°C for 135 s, and a final cycle of 72°C for 5 min. Then, PCR products were subjected to electrophoresis in 1.5% agarose gel and TAE buffer. Data were analysed using the BioNumerics package (Version 6.01, Applied Maths, Belgium) based on images of PCR MP electrophoretic band patterns obtained for the analysed strains. Dendrograms were generated with BioNumerics software using the Dice similarity coefficient and clustering by the unweighted pair group method with arithmetic mean (UPGMA), with 1% tolerance for differences in the band position (Fig. 1).

Fig. 1

Differentiation of 14 potentially pathogenic E. coli strains using MP-PCR method (Hind III) with 1% tolerance for differences in the band position (UPGMA)

Whole-genome shotgun sequencing

Genomic DNA obtained from 14 potentially pathogenic strains was used to prepare Illumina paired-end libraries with the Nextera XT DNA Library Preparation Kit (Illumina, USA) protocol. Whole-genome shotgun sequencing was performed in-house on the Illumina NextSeq 500 platform (Illumina, USA) at a read length of 2 × 150 bp. The draft genomes were further assembled de novo by SPAdes versions 3.7.1 and 3.8.0 (St. Petersburg State University, Russia) (27), with necessary manual editing and evaluation accomplished by FA_TOOL (21). Annotation was carried out automatically using the NCBI pipeline (PGAAP). Nucleotide sequence accession numbers and the results of de novo assemblies are presented in Table 1. The draft genome sequences of the E. coli strains which are pathogenic for poultry, namely: 001PP2015, 002PP2015, 004PP2015, 007PP2015, 009PP2015, 011PP2015, 012PP2015, 015PP2015, 016PP2015, 017PP2015, 018PP2015, 019PP2015, 020PP2016, and 022PP2016, were deposited at DDBJ/EMBL/GenBank (Table 1).

Virulence genotyping

Based on already published data (12, 14, 26, 28), genes representing the various groups such as toxins, iron acquisition genes, protectins, or adhesins, known as virulence factors of ExPEC strains, were chosen for in silico analysis: cvaC(UniProt accession no. P22522), iroN(A0A0H2XJD7), iucC(B7NCN3), iutA(A0A0H2XK75), sitA(Q1RCU3), traT(B7LI89), tsh(A0A0H2XJN4), feoB(Q8FCT7), ireA(B7MN59), irp2(Q1RAG0), hlyD(A0A0H2YYP7), hlyF(Q849X5), fliC (Q1RAL7), chuA(Q1R5A7), iss(P19592) and ompT(Q3L7I1). Amino acid sequences of APEC strains, annotated in PGAAP and deposited in NCBI database (Table 1) and original sequences of sequenced ExPEC genomes also annotated in the NCBI annotation system were tested for the presence of the mentioned virulence factors.

DIFFIND

A total of 30 genomes of APEC strains were chosen for in silico analyses of virulence associated with selected genes (14 genomes of APEC strains from Poland presented for the first time in this study and 16 previously reported genomes from USA, Denmark, Brazil, China, and Italy (GenBank accessions nos. ACN 001 CP007442.1; APEC IMT5155 CP005930.1; APECO1 CP000468.1; APEC O1 NC008563.1; APEC O2 CP006834.1; APEC O2 LSZR00000000.1; APEC O08 AOGM00000000.1; APEC O18 CP006830.1; APEC O78 CP004009.1; APEC RS76 CP013048.1; APEC S17 AOGN00000000.1; APEC SCI-07 GCA00252805.2; APEC SEPT362 AOGL00000000.1; APEC ST410 LBBC00000000.1; E. coli 44 LXWV00000000.1; and E. coli 51 LYPJ00000000.1) (5, 6, 11, 18, 19, 26, 29, 30, 32). For in silico examination, amino acid sequences of selected genes were used. Mapping to the constructed reference sequence built from amino acid sequences of selected VAGs was performed by using the newly designed tool DIFFIND (-c 0.7 -g 1 -s2 0.8) (20) based on the Luigi framework involving the Python packages Pathos (23) and Scipy (13) and the novel modules pycdhit (22) and FA_TOOL (21). It allowed estimation of the degree of similarity between compared sequences by using farthest point algorithm calculation completed by generation of a heatmap plotted with the Matplotlib module (8). The analysis was facilitated by visualisation of the matching results of the sequence comparison with minimal sequence similarity ℑ70%. Branches visualised on the dendrogram were separated if the distance value between two consecutive nodes exceeded 0.7 × maximum distance between all nodes. This calculation was considered an independent cluster reflected into the various colours of each branch (Fig. 2).

Fig. 2

In silico examination of prevalence of chosen virulence factors within 30 analysed APEC strains using DIFFIND tool. 001 − 001PP2015 (MBPT00000000), 002 − 002PP2015 (MBPU00000000), 003 − 004PP2015 (MBPW00000000), 004 − 007PP2015 (LWSE00000000), 005 − 009PP2015 (LWUK00000000), 006 − 011PP2015 (LXWO00000000), 007 − 012PP2015 (LXWP00000000), 008 − 015PP2015 (LXYD00000000), 009 − 016PP2015 (LXYC00000000), 010 − 017PP2015 (LYWM00000000), 011 − 018PP2015 (MBPX00000000), 012 − 019PP2015 (LYWL00000000), 013 − 020PP2016 (MBPY00000000), 014 − 022PP2016 (MAYE00000000), 015 − ACN_001 (CP007442), 016 − APEC_IMT5155 (CP005930), 017 − APEC_O1 (CP000468), 018 − APEC_O1 (NC008563), 019 − APEC_O2 (CP006834), 020 − APEC_O2 (LSZR00000000), 021 − APEC_O08 (AOGM00000000), 022 − APEC_O18 (CP006830), 023 − APEC_O78 (CP004009), 024 − APEC_RS76 (CP013048), 025 − APEC_S17 (AOGN00000000), 026 − APEC_SCI-07 (GCA00252805), 027 − APEC_SEPT362 (AOGL00000000), 028 − APEC_ST410 (LBBC00000000), 029 − E_coli_44 (LXWV00000000), and 030 − E_coli_51 (LYPJ0000000000)

Results

The results of differentiation with the MP-PCR method are presented in Fig. 1 in the form of dendrograms (UPGMA) produced in BioNumerics software. Estimation of chosen virulence factors in the reported genomes in comparison to known APEC sequences was based on the presence or absence of 16 VAGs, which resulted in a formation of clusters with clearly visible strain differentiation. During cluster analysis, strains were distributed into five main clusters: yellow (no. 1), violet (no. 2), blue (no. 3), red (no. 4), and green (no. 5), differentiated by the selected gene content (Fig 2.). Results show that the strains under comparison contained at least two genes described as virulence indicators for ExPEC strains. The virulence associated genes with the highest frequencies (≥70%) in all analysed genomes, both original and previously reported, were: feoB 100% (30/30), ompT 93.3% (28/30), iss 86.6% (26/30), sitA 80.0% (24/30), iutA 70.0% (21/30), and iucC 70% (21/30). The VAGs with the lowest presence in examined genomes were fliC13.3% (4/30) and ireA 40 % (12/30). Phylogenetic typing demonstrated that the analysed strains contained various numbers of selected virulence associated genes. Clusters 1 and 2 (Fig. 2) included the highest contents of toxins, adhesins, protectins, iron acquisition, and serum resistance genes ranging from 12 to 15 and from 10 to 12, respectively. Clusters 3, 4, and 5 contained VAGs in the range from 2 to 7 for cluster 3 and from 5 to 9 for clusters 4 and 5.

Discussion

It is well known that the mechanisms of colibacillosis, a common and significant bacterial disease of poultry, have so far been only poorly explained (1, 16). Many extraintestinal E. coli share common virulence genes. It is reported that APEC and human UPEC (uropathogenic E. coli) share many virulence traits (28). Moreover, many human and animal-source ExPEC isolates exhibit similar VAGs which may suggest that they represent zoonotic pathogens (10). What is more, previously described virulence factors associated with this disease do not always serve to distinguish APEC from commensal strains as they can be found in both groups (16, 31). The genomes announced in this paper provide original information concerning sequences of extraintestinal pathogenic Escherichia coli strains (ExPEC) isolated from colibacillosis lesions. The methods chosen in this study, based on sequencing of the whole genomes and assembling them de novo, seem to be reasonable choices, enhancing knowledge about the DNA sequences of pathogenic strains and simultaneously recognising chosen genes coding for virulence factors. Moreover, it is an innovative approach to the problem as previous reports mainly used PCR-based analysis (4, 12, 17, 24, 28). There are already many known and described virulence factors (1, 3, 12, 14, 16, 25) and PCR-based tests aimed at defining if a chosen strain is a representative of APEC. For example, Ewers et al. (4) designed a multiplex PCR in which the presence of a minimum of five out of eight genes determined the virulence of a strain, while Johnson et al. (12) reported five genes (iutA, hlyF, iss, iroN, and ompT) to be the minimal predictors of pathogenic strains, and Silveira et al. (31) compared APEC and AFEC strains and revealed that a bigger set of genes might be needed to differentiate pathogenic from non-pathogenic bacteria well. Therefore, in this study, a set of 16 factors reported previously (7, 12, 14, 26, 28) as virulent for poultry were chosen and in silico analysis of 14 original APEC genomes with 16 strains regarded as APEC was performed and the results deposited in the NCBI database. Our results characterised the analysed strains by the varying content of selected virulence associated genes ranging from 2 to 15 VAGs. It is worth noticing that the substantial number of genes involved in the iron acquisition mechanisms are essential for expression of pathogenicity of strains, which is in accordance with predictions as the extra-intestinal environment is iron deficient, and it is clear to see that they are multiply present in all clusters. However, if virulent strains are to be compared between themselves, genes coding serum resistance or adhesins such as ompT and iss could be chosen to differentiate them. In future work more traits should be subjected to such analysis to distinguish strains with high- and low-pathogenic profiles. However, even such a small set of chosen virulence factors indicate that strains causing colibacillosis can be highly different in their genotypes, even if they are isolated from proximate geographic regions (Poland). Efforts to define the pathotype of APEC and hence better understand colibacillosis can contribute to a reduction in spread of foodborne agents from contaminated poultry to people and an increase in food safety and animal well-being. The DIFFIND tool provides a reliable method for fast estimation of amino acid sequence similarities in at least two representatives of different species and may be useful for detection of virulence factors.