Comparative genomic analysis of Escherichia coli isolates from cases of bovine clinical mastitis identifies nine specific pathotype marker genes.

Escherichia coli is a major causative agent of environmental bovine mastitis and this disease causes significant economic losses for the dairy industry. There is still debate in the literature as to whether mammary pathogenic E. coli (MPEC) is indeed a unique E. coli pathotype, or whether this infection is merely an opportunistic infection caused by any E. coli isolate being displaced from the bovine gastrointestinal tract to the environment and, then, into the udder. In this study, we conducted a thorough genomic analysis of 113 novel MPEC isolates from clinical mastitis cases and 100 bovine commensal E. coli isolates. A phylogenomic analysis indicated that MPEC and commensal E. coli isolates formed clades based on common sequence types and O antigens, but did not cluster based on mammary pathogenicity. A comparative genomic analysis of MPEC and commensal isolates led to the identification of nine genes that were part of either the core or the soft-core MPEC genome, but were not found in any bovine commensal isolates. These apparent MPEC marker genes were genes involved with nutrient intake and metabolism [adeQ, adenine permease; nifJ, pyruvate-flavodoxin oxidoreductase; and yhjX, putative major facilitator superfamily (MFS)-type transporter], included fitness and virulence factors commonly seen in uropathogenic E. coli (pqqL, zinc metallopeptidase, and fdeC, intimin-like adhesin, respectively), and putative proteins [yfiE, uncharacterized helix-turn-helix-type transcriptional activator; ygjI, putative inner membrane transporter; and ygjJ, putative periplasmic protein]. Further characterization of these highly conserved MPEC genes may be critical to understanding the pathobiology of MPEC.

Data Summary

Sequencing data and genome assemblies are available from GenBank/ENA/DDBJ as BioProject PRJNA612640, under the accession numbers JAASLI000000000–JAASQG000000000.

Mammary pathogenic (MPEC) is a common cause of mastitis in dairy cattle. It is still controversial as to whether MPEC is a unique pathotype, since a core set of virulence factors that are unique to MPEC have not yet been defined. Our comparative genomics study of MPEC and bovine commensal identified nine unique MPEC genes. The nine genes are associated with nutrient intake, metabolism and fitness; in addition, we have identified that a few virulence factors common to uropathogenic are found in MPEC, but are absent from commensal bovine . These genes may also be highly conserved in the genomes of MPEC because of the absence of genomic islands in MPEC genomes. This study represents a significant step towards further understanding the pathobiology of MPEC, as well as designing MPEC targeted diagnostics and treatments.

Introduction

Bovine mastitis – inflammation of bovine udder usually caused by a bacterial infection – is a costly disease in the dairy industry [1], and results in annual losses of $665 million (CAD) [£386 million, 1 CAD=£0.58] for the Canadian dairy industry [2], $2 billion (USD) [£1.4 billion, 1 USD=£0.71] for the American dairy industry [3] and £168 million for the British dairy industry [4]. The aetiological agents of bovine mastitis can be categorized as either contagious or environmental pathogens. Contagious bovine mastitis is commonly caused by , , and , which are transmitted from infected to uninfected cows via milking equipment, direct contact or vectors like farm workers. Modern dairy farm practices, including early mastitis prevention programmes, were focused on controlling contagious mastitis, and now, as a result, environmental mastitis is the most common form of this disease [5]. Environmental mastitis pathogens originate from the farm environment, such as pasture, stable or bedding material. The bovine gastrointestinal tract is a common source for environmental pathogens [5]. is the most common aetiological agent of environmental mastitis [6, 7].

is a genetically and phenotypically diverse bacterial species. The range of diversity is particularly apparent in terms of host–bacteria relationships where it can be a mutualist, commensal, pathogen or occasional symbiont in the gastrointestinal tract of a variety of host species [8]. In humans, pathogenic strains are broadly categorized as either diarrhoeagenic or extraintestinal pathogenic (ExPEC). ExPEC typically reside asymptomatically within the intestine, but cause severe infection when allowed to colonize extraintestinal niches [9]. Within each broad group, there are several sub-groups of strains that share virulence factors and share similar clinical manifestations, which are known as pathotypes [10]. Uropathogenic (UPEC), which is the aetiological agent of about 90 % of human urinary tract infections [11], has been relatively recently recognized as a distinct ExPEC pathotype [10]. This infection was once thought to be an opportunistic infection caused solely by displacement of any intestinal into the urinary tract [10], but now it is known that only a distinct subset of , originating from the gastrointestinal tract, result in UPEC infections [12]. Four main UPEC phylogroups (A, B1, B2 and D) have been identified based on the presence of UPEC-specific virulence [13]. Most virulent UPEC strains are from the B2 lineage [14]. Many pathogenicity-associated islands (PAIs) are associated with UPEC, and these islands can carry important virulence factors, specifically: P fimbriae, type I fimbriae, haemolysins, iron-acquisition proteins, bacteriocins and the malX gene, which is associated with the phosphotransferase system enzyme II that uses glucose and maltose as the main substrates [15-17].

The existence of a distinct mammary pathogenic (MPEC) pathotype has been proposed [9], but defining virulence factors of this group have not yet been identified [18]. The lack of a set of virulence genes common to all MPEC isolates, despite several attempts to identify them [18-22], has led to a proposed model for this disease where the mere introduction of any gastrointestinal-originating into the mammary gland and the resultant inflammatory response can result in clinical mastitis [23, 24]. In this model, the severity of clinical mastitis is primarily dependent on host factors. However, this model fails to explain several aspects of clinical mastitis. For example, not all strains can cause clinical mastitis in experimental models of the disease [25], and mastitis strains are much less genetically diverse than bovine commensal [18, 22]. The Fec system appears to be much more common in MPEC isolates than in other isolates derived from dairy cow environments; in addition the Fec system is overexpressed when MPEC strains are grown in milk, and Fec knockouts are unable to induce clinical mastitis [26, 27]. Thus, the complex aetiology of mastitis caused by is not fully understood.

In this study, we advance previous work by performing a detailed genomic analysis of 113 novel MPEC isolates. To identify the genetic traits that differentiate MPEC isolates from other bovine isolates, we performed a comparative genomic analysis in which MPEC isolates were compared to 100 isolates from dairy cattle habitats that were not associated with disease.

Methods

MPEC isolates and genomes of bovine commensal

MPEC isolates (n=113) were obtained in 2019 from the Mastitis Pathogen Culture Collection, which is maintained and curated by the Canadian Bovine Mastitis Research Network [28]. Each isolate was obtained from milk samples originating from 113 different cows from 57 herds (Alberta=9, Ontario=17, Quebec=17 and Atlantic provinces=14) experiencing clinical mastitis either on the day of diagnosis (n=100) or on subsequent post-clinical mastitis follow-up sampling (within 14 days, n=7; between 14–28 days, n=6) between 2007 and 2008 [28]. As previously described, MPEC isolates were isolated on bi-plates containing Columbia agar with 5 % sheep blood and MacConkey agar, and biochemical tests were performed to confirm the isolates were (lactose and indole positive, oxidase and citrate negative) [29]. Bovine metadata, including herd number and location, cow ID, quarter position, sampling data, mastitis severity score [30], days in milk (DIM) at sampling and cow’s parity, are summarized in Table S1 (available with the online version of this article).

The whole genomes of 100 bovine isolates not associated with bovine disease were obtained from the National Center for Biotechnology Information (NCBI) database. These genomes were from isolates from bovine faeces, skin, cow sheds and milking areas as described in previous studies, and came from a variety of international locations excluding Canada [18, 26, 31, 32] (Table S2). The sequenced reads for the bovine commensal genomes were assembled using Platanus v1.2.2, Newbler v2.3 [31], CLC Genomics Workbench v.6.5.2 [26], and SPAdes v3.1.1 [18] and v3.5.0 [32] (Table S2).

Whole-genome sequencing, assembly and annotation

Each MPEC isolate was streaked on tryptic soy agar (TSA) (Becton Dickinson) and incubated overnight at 37 °C. A single colony was picked and incubated in tryptic soy broth (TSB) (Becton Dickinson) overnight at 37 °C at 200 r.p.m. DNA was extracted from each isolate with a culture that had >1×108 cells ml−1 (OD600 > 0.8) at the time of extraction using DNAzol reagent (Invitrogen), following the manufacturer’s instructions. DNA was further purified using the Qiagen DNeasy PowerClean Pro Cleanup kit (Qiagen), as per the manufacturer’s instructions. DNA from isolates that did not produce high-quality DNA via this method was re-extracted using the Maxwell RSC instrument and the recommended Blood DNA kit (Promega), according to the manufacturer’s instructions. A DNA concentration between 10 and 100 ng μl−1, with corresponding purity measurements of A260/A280 >1.8 and A260/A230 between 1.8 and 2.2 based on Nanodrop measurements (Thermofisher), was achieved prior to each sequencing library preparation.

DNA was further quantified using the Quant-iT dsDNA assay kit prior to library preparation (Thermofisher). DNA library preparation was performed using a Nextera DNA Flex Library Prep kit (Illumina) optimized for short-read sequencing by the Illumina MiSeq system, as per the manufacturer’s instructions. The tagmentation step was optimized to 15 min to achieve a DNA target length of 500–600 bp, this was followed by a clean-up step. Tagmented DNA was amplified using Nextera DNA CD indexes via PCR, followed by a clean-up step and concentration check. A pooled library was made combining all samples into one 1.5 ml tube, and a final quantification step was performed to ensure a final concentration of 1.6 ng μl−1 (4 nM). After library pool denaturation was performed by adding 5μL of 0.2N sodium hydroxide, a final concentration of 12 pM was obtained and a PhiX control was added to a concentration of 20 pM. The library and PhiX control were loaded into a MiSeq v3 reagent kit, and 600 cycles (300 forward and 300 reverse) of sequencing was conducted using a MiSeq benchtop sequencer (Illumina).

Sequence reads were de novo assembled using the software pipeline ProkaryoteAssembly version 0.1.6 (https://github.com/bfssi-forest-dussault/ProkaryoteAssembly). This pipeline includes quality control and trimming of low-quality sequences (Q value <20) using BBDuk (BBMap v38.79), error-correction using Tadpole (BBMap), assembly using Skesa v2.4, alignment of error-corrected reads against draft assembly BBMap and polishing of assembly using Pilon v1.23 [33-35]. After assembly, contigs shorter than 1 kbp were discarded, and the coverage and contigs were quantified using Qualimap [36]. Prokka was used to annotate the assembled contigs of genomes of MPEC and bovine commensal [37]. The pipeline includes annotation of protein-encoding genes by identifying coordinates of candidate genes from ISfinder, UniProt, Pfam and TIGRFAMs [38-42].

Pan-genome analysis

Roary was used to construct a pan-genome for MPEC and bovine commensal isolates to allow for a direct comparison between the two groups of genomes [43]. The predicted functional proteins encoded in the pan-genome of MPEC and the commensal sets were identified by Clusters of Orthologous Groups (COGs) on eggNOG-mapper (E value >1×10−10) [44, 45]. The core genome alignment file was used as input for iq-tree, which can use the ModelFinder Plus algorithm, selecting the best performing substitution model and building a tree with it [46, 47]. Specifically, the GTR+F+R10 model was used on iq-tree to build a phylogenomic tree of MPEC and the commensal genomes. To visualize the tree, interactive Tree Of Life (iTOL) v4 (https://itol.embl.de) was used [48].

Core and soft-core genes from MPEC and commensal isolates were identified using Roary, and compared using Venny v2.1, determining the genes of MPEC isolates that are either accessory genes (shell or cloud genes) or unique genes (not shared with any commensal isolates) in the commensal isolates by Venn diagram [49]. The genes in the unique group of MPEC genome in the diagram that are copies of the same gene in the group of the commensal genome were discarded after local blast using the pan-genome of the commensal group as reference by BioEdit v7.2 [50]. The unique genes that are annotated as ‘hypothetical protein’ were searched on NCBI blastx (https://blast.ncbi.nlm.nih.gov/Blast.cgi) using their nucleotide sequences as query to find closely related protein (identity and coverage >98 % and E value <1×10−10) [51]. One hundred additional MPEC genomes from previous studies were obtained to be searched to determine whether the identified MPEC marker genes were applicable beyond the 113 isolates investigated in this study, using the command line version of blast [20, 26].

Identification of sequence type (ST), O and H antigens, plasmid replicons and genomic islands (GIs)

STs of each isolate were identified using the tool mlst (https://github.com/tseemann/mlst), which incorporates data from the PubMLST database [52]. To identify the distribution of O and H serotypes, ABRicate v1.0 (https://github.com/tseemann/abricate) was used with the EcOH database for O and H serotypes [53]. Minimum coverage and identity settings for the screening were set to 90 %. ABRicate was used to identify replicons of plasmids, using the PlasmidFinder v2.1 database, and plasmid multilocus sequence typing (pMLST) was performed on the most prevalent replicon of plasmids in MPEC and the bovine commensal genomes [54]. Putative GIs were predicted using IslandViewer 4 using IslandPath-DIMOB and SIGI-HMM as island prediction methods. The previously closed ECC-1470 genome was used as a reference strain (accession no. NZ_CP010344.1) (>8 kbp as cut-off) [55]. The identified predicted GIs were screened if the unique genes of MPEC, Fec operon genes and the ail gene were present.

Results

Quality of sequenced genomes of MPEC and bovine commensal

The assembly of each draft genome for MPEC isolates was evaluated; the coverage and number of contigs are reported in Table S1. The range of coverage for individual genomes was between 22× and 360×, and the number of contigs ranged from 28 to 149. The genomes of bovine commensal were selected from those available in the NCBI database based on isolation from dairy cattle environments including cowsheds, faeces, skin, gastrointestinal tracts or from the milking room, having coverage between 20× and 90×, and having less than 419 contigs [18, 19, 26, 31, 32, 56, 57].

Absence of major clusters of MPEC by origin, herds and provinces

A phylogenomic tree that illustrates the phylogenomic relatedness of the 113 MPEC and 100 bovine commensal isolates examined in this study was created by comparing core-genomeSNPs across the entire genome of each isolate (Fig. 1). There was no significant clustering of MPEC or commensal isolates based on origin, herds or provinces, and MPEC isolates were not phylogenetically differentiated from commensal isolates. There was a large range in the diversity of isolates in this study, which included 102 different STs, 88 different O antigens and 38 different H antigens. The STs, O antigens and H antigens of each genome are indicated in Tables S1 and S2.

Fig. 1.

Phylogenomic tree of clinical mastitis-related MPEC and bovine commensal isolates by core-genomeSNPs. The phylogenomic tree was reconstructed using iq-tree based on the core genomes of MPEC and commensal genomes. The tree was visualized using iTOL v4 and each genome was annotated with STs by multilocus sequence typing (n=102), O antigens (n=88) and H antigens (n=38). The scale bar is subistutions per site. The nd is not determined.

Comparative genomic analysis between clinical mastitis-related MPEC and bovine commensal isolates

The pan-genomes of MPEC and the commensal isolates were constructed using Roary after assembly and annotation of each individual genome. A total of 17 532 and 20 042 genes were identified in the pan-genome of MPEC and the commensal isolates, respectively. The MPEC pan-genome included 3391 core and soft-core genes (a core gene is defined as a gene that is shared by 99–100 % of genomes, and a soft-core gene is a gene that is found in 95–99 % genomes), 1638 shell genes (a shell gene is a gene that is shared by 15–95 % of included genomes) and 12 503 cloud genes (a cloud gene is a gene that is shared by 0–15 % of genomes). The commensal pan-genome included 3538 core and soft-core genes, 1539 shell genes and 14 965 cloud genes.

The pan-genome of MPEC and the commensal genomes were compared via functional classification by COGs. There was no significant difference between the COGs of MPEC and the commensal genomes (P=0.85; P >0.05) (Fig. 2).

Fig. 2.

COGs of pan-genes of clinical mastitis-related MPEC and bovine commensal . The groups were identified using eggNOG-mapper with E value >1×10−10. The COGs are related to information storage and processing (group B, J, K, L), cellular processes and signalling (group D, V, T, M, N, W, U, O), metabolism (group C, G, E, F, H, I, P, Q) and uncharacterized functions (S).

To identify core and soft-core genes that are unique to MPEC relative to the other bovine-associated isolates, a gene-by-gene pairwise comparison between MPEC and other commensal isolates core and soft-core genes was performed (Fig. 3). This analysis identified 91 genes that were both unique to and widely conserved in MPEC isolates. Each of these 91 genes was individually compared to each commensal genome using command-line blast to identify any that were identical to a commensal cloud gene. Hypothetical genes were also manually removed from this set. Refining the pool of these 91 genes left 22 potential MPEC marker genes. Of these 22 genes, 13 were identified as part of the shell genes in the commensal genomes (Table 1) and 9 genes were unique to only MPEC isolates (Table 2).

Fig. 3.

The number of core and soft-core genes in clinical mastitis-related MPEC and bovine commensal genomes illustrated by a Venn diagram. The core genes of each genome set were extracted from the pan-genome result by Roary. Local blast against each set of genomes was conducted to distinguish hypothetical protein genes that were identical but with the same gene name. Then, the names of core and soft-core genes with annotation from each pan-genome result were used with Venny v2.1 to generate a Venn diagram showing the genes that are only for MPEC (n=91), only for bovine commensal (n=233) or for both sets of genomes (n=3288). A total of 91 genes were identified as being unique to MPEC. However, when each of these genes were manually annotated and compared to commensal genomes, this number was reduced to nine marker genes.

Table 1.

List of core and soft-core genes of clinical mastitis-related MPEC identified in the commensal genomes as shell genes

Gene	Putative annotation	Relative abundance in MPEC isolates (n/113)	Relative abundance in commensal E. coli isolates (n/100)
ail*	Putative phage portal protein	111/113	22/100
appX	Putative cytochrome bd-II ubiquinol oxidase subunit	109/113	91/100
dedA	Uncharacterized protein	113/113	62/100
fecA	Fe(3+) dicitrate transport protein	110/113	27/100
fecC	Fe(3+) dicitrate transport system permease protein	108/113	26/100
fecI	Putative RNA polymerase sigma factor	110/113	27/100
fecR	Regulator of iron dicitrate transporter	110/113	27/100
folK	2-Amino-4-hydroxy-6-hydroxymethyldihydropteridine pyrophosphokinase	109/113	93/100
ghoT	Toxic component of a type V toxin–antitoxin (TA) system	110/113	93/100
higA	Antitoxin	109/113	83/100
yjiK	Putative protein	108/113	90/100
yqeI*	Transcriptional regulatory protein, C terminal protein	108/113	92/100
ybfB*	Uncharacterized MFS-type transporter	113/113	89/100

*Genes that were initially annotated as hypothetical proteins. The identical genes were identified manually by blast on NCBI and UniProt.

Table 2.

Relative abundance of clinical mastitis-related MPEC isolates from this study and previous studies

Gene	Putative annotation	Relative abundance in clinical mastitis-related MPEC isolates from this study (n/113)	Relative abundance in clinical mastitis-related MPEC isolates from previous studies (n/100)
				adeQ	Adenine permease	Core (112/113)	100/100
yfiE	HTH-type transcriptional activator	Core (113/113)	100/100
nifJ	Pyruvate-flavodoxin oxidoreductase	Core (113/113)	100/100
ygjI	Putative inner membrane transporter	Core (112/113)	97/100
yhjX	Putative MFS-type transporter	Core (112/113)	100/100
ygjJ*	Putative periplasmic protein	Core (113/113)	100/100
pqqL*	Zinc metallopeptidase	Core (112/113)	100/100
fdeC*	Intimin-like adhesin	Soft-core (110/113)	95/100
Group_69†	–	Soft-core (108/113)	96/100

*Genes that were initially annotated as hypothetical proteins. The identical genes were identified manually by blast on NCBI and UniProt.

†Pseudogene.

Relative abundance in commensal isolates (n/100)

–

To verify whether the nine unique MPEC marker genes identified in this study were indeed good markers for MPEC, an additional 100 clinical mastitis-related MPEC genomes sequenced by previous studies were downloaded and marker genes were identified using a local blast (Tables 2, S1 and S2) ygjI, fdeC and group_69 genes were identified from 97, 95 and 96 genomes, respectively, while the other marker genes were present in 100 % of additional MPEC genomes.

Analysis of mobile genetic elements in MPEC

The replicons of plasmids in MPEC and the commensal genomes were screened using ABRicate (https://github.com/tseemann/abricate) based on the PlasmidFinder database [54]. In the MPEC isolates from this study, 25 different types of plasmid replicons were identified, and 31 types of replicons were identified in the commensal isolates (Tables S1–S3). The plasmid IncF was the most prevalent replicon of plasmid type in both MPEC (n=79) and the commensal (n=76) isolates. IncF replicon sequence typing was conducted to identify the difference between IncF-type plasmid replicons in the 79 MPEC and 76 commensal genomes. The most prevalent replicon of IncF type plasmid was IncFIB (AP001918), which was identified in 61 out of 79 MPEC genomes and 67 out of 76 commensal genomes (Table S4). Twenty MPEC genomes had novel alleles of FIA (similar replicon with >97 % identity: 67) and FII (similar replicon with >97 % identity: 64) replicons, while only one commensal bovine genome had these novel alleles. These MPEC genomes with two novel alleles of FIA and FII contained replicons of IncFIB(AP001918) and IncFIC(FII), except one genome.

IslandPath-DIMOB and SIGI-HMM prediction methods were used within the IslandViewer 4 platform to predict the presence of GIs in MPEC genomes using an alignment-based strategy and the closed MPEC genome ECC-1470. Each MPEC isolate contained between 14 and 35 predicted GIs. No MPEC unique genes, except fdeC, were located in the predicted GIs in MPEC genomes. Ten MPEC genomes contained fdeC in the predicted GI that contains ykgOMR (50S ribosomal protein L36, L31 type B, putative membrane protein), ecpABCDER ( common pilus), paoABCD (aldehyde oxidoreductase). The presence of Fec operon genes was also identified in predicted GIs from MPEC genomes. All of the Fec operon genes (fecABCDEIR) were identified in the predicted GIs in 65 out of 110 MPEC genomes that contained the operon. The rest of the genomes contained either partial Fec operons or did not contain the Fec operon in the predicted GIs (Tables S5 and S6). One hundred and one MPEC genomes contained an ail gene on a predicted GI and these GIs commonly contained the genes related to environment adaptation: ydfO, ydfR, gnsA (putative proteins); cspB, cspG, cspJ (cold shock-like proteins); rrrD (lysozyme); hokC (toxic compound of a type I toxin–antitoxin system); relE, relB (type I toxin–antitoxin system); and flxA (phage or prophage related protein).

Discussion

In this study, 113 novel clinical mastitis-related MPEC genomes were characterized, and a comparative genomics approach was used to identify marker genes that could potentially differentiate MPEC isolates from bovine commensal . Nine MPEC marker genes were ultimately identified. These marker genes are involved in a variety of cellular processes including the uptake of nutrients, metabolism, transcriptional regulation and virulence. The adeQ gene encodes adenine permease, which may be involved with uptake of adenine. However, it may not be essential for the MPEC pathotype, since both MPEC and commensal genomes have an isozyme adeP that encodes a second adenine permease with higher affinity to adenine than AdeQ [58]. Two of the potential marker genes, nifJ and yhjX, are induced by pyruvate and involved in metabolism. The nifJ gene encodes pyruvate flavodoxin oxidoreductase, which catalyses the oxidation of pyruvate to acetyl-coenzyme A, followed by reduction of flavodoxin (NifF) providing an electron to dinitrogenase reductase, which then provides an electron to dinitrogenase [59, 60]. Dinitrogenase is a nitrogen-fixing protein that reduces N2 to form ammonia, and MPEC may utilize this pathway to obtain nitrogen from milk, which contains bounded nitrogen in the form of casein and whey, non-protein nitrogen and urea [61].

Although its function is not fully characterized, the yhjX gene encodes a major facilitator superfamily (MFS) type transporter and is targeted by PyrSR, which is reported to be induced after pyruvate uptake during the exponential growth phase [62]. There could be potential interplay between YhjX and YjiY, which is a pyruvate/H+ symporter regulated by BtsSR resulting in pyruvate uptake followed by expression of YhjX [63]. As yhjX and yjiY are core genes of MPEC, there might be further regulatory processes in response to pyruvate uptake for their survival in the bovine mammary gland compared to non-MPEC bovine isolates.

Two MPEC marker genes identified in this study, pqqL and fdeC, have previously been identified in UPEC isolates, and can contribute to the fitness and virulence of UPEC [64, 65]. The pqqL gene likely encodes a zinc metallopeptidase, and is reported to act with yddA and yddB to form an ABC transporter ATPase and an outer membrane β-barrel protein, respectively, as a locus (yddABpqqL) in the UPEC genome [64]. Even though the effect of the yddABpqqL locus on fitness and growth of UPEC is not fully characterized, it was reported that expression is highly upregulated under iron-limiting conditions, which are similar in urine and milk [66, 67]. Unlike the UPEC genome, which contains pqqL and yddA as core genes, the MPEC genomes from this study contained pqqL and yddA as a core and a soft-core gene, respectively; however, the yddB gene was not present in any MPEC genomes. However, yddB has a high degree of sequence similarity to another outer membrane β-barrel protein, ferrienterobactin fepA; therefore, it might not be necessary for MPEC to possess yddB while containing fepA and other iron-uptake systems. Another common gene in UPEC also observed in MPEC was fdeC, an adhesin that shares similarity to intimin and other intimin-like adhesins such as eaeH of ETEC (94 % of similarity) [65, 68, 69]. It has been shown that fdeC is expressed by UPEC when bound to the plasma membrane of human bladder and urethral epithelial cells in vitro, and that it is associated with an aggressive UPEC phenotype [65]. The fdeC gene has also been found in human gastrointestinal isolates, from healthy individuals, indicating the presence of this gene is not necessarily associated with pathogenic [70, 71]. The enterohaemorrhagic (EHEC) isolate, N39 (also known as EC673), with FdeC from bovine faeces of Australian calf was also characterized and its expression level of FdeC is significantly higher at >39 °C [72]. This indicates that MPEC with FdeC can potentially originate from the bovine rectum where there is a consistent temperature above 37 ℃, and as the temperature of bovine udder with mastitis is above 38℃, fdeC could be upregulated during clinical mastitis [72, 73].

The functions of proteins encoded by other unique genes, such as yfiE, ygjI, ygjJ and group_69 genes, are not characterized yet. YfiE is reported to be an uncharacterized helix-turn-helix (HTH)-type transcriptional activator that is predicted to be involved with a repressor for the metabolism of cofactors, vitamins and amino acids based on computational analysis [74]. This might play a role in regulating the uptake of nutrients and utilization of metabolism along with nifJ and yhjX, which may be involved with metabolism and uptake of nutrients. ygjI is localized in the ebg operon, β-galactosidase genes, and encodes a putative transporter localized in the inner membrane, and ygiJ encodes a periplasmic protein of unknown function [75, 76]. The other unique gene, group_69, was identical to a pseudo gene from strains (isolated from stray dog and fox) and shares 58 % coverage with ShET-2 gene in tblastx, indicating this gene exists as a pseudo gene in MPEC genome [77].

The Fec operon, which was identified by previous studies to have a higher prevalence in MPEC than in non-MPEC , was identified in our current study as part of the MPEC soft-core genome and in the shell genome of commensal . This result agrees with that of Leimbach et al. [18], who found that fecIRABCDE genes were present in at least 50 % of commensal genomes [18]. It was also reported that from other mammals, fish, frogs, turtles, snakes and lizards, crocodiles, birds and lakes, especially the ones that belong to phylogenomic group A, contained Fec operons [78]. Considering the presence of the operon not only in the dairy environment but also in other animal hosts and the environment, the Fec operon alone is not a good marker of MPEC. The Fec operon also could not be strictly essential for MPEC as other genes such as efeUOB, which encodes an iron (Fe2+) uptake system, were contained as core genes in the isolates [79].

The ail gene, which encodes a putative phage portal protein, was a soft-core gene in MPEC, but in the shell genome of commensal . Multiple copies of the ail gene were identified in the predicted GIs with other phage genes such as nohA, a prophage DNA-packing protein gene, and tfaE, prophage tail fibre assembly protein gene. The gene is also found to be identical to the one in complete genome of isolated from bovine clinical mastitis (accession no. CP009166.1), indicating that this might be the potential coliphage specifically targeting MPEC [80]. The predicted GIs that contained ail also commonly had genes that are related to stress tolerance such as cold shock-like protein (cspB, cspG, cspJ), toxic component of type I toxin–antitoxin system (hokC), type II toxin–antitoxin system (relB, relE) and lysozyme (rrrD) genes. However, it is unclear whether these genes are crucial for inducing mastitis and adaptation in the mammary gland, as these were shell genes in both MPEC and bovine commensal genomes. Therefore, interactions with these genes in GIs and other prevalent genes, including nine unique genes in MPEC that contribute to the pathogenicity, needs to be characterized.

This study has identified the unique genes in the MPEC genome that differentiate it from the other bovine commensal . While there were no significant differences in COGs by pan-genes and phylogenomic relationship, nine unique genes were conserved as core and soft-core genes in the MPEC genome. In the future, the presence of these genes in can possibly be used to make advances in the diagnosis and therapeutics for MPEC.

Click here for additional data file.