Comparative Genome Analysis of Enterococcus cecorum Reveals Intercontinental Spread of a Lineage of Clinical Poultry Isolates

To obtain insight into the gene repertoire of E. cecorum, we performed whole-genome sequencing of 118 isolates, including 100 clinical isolates collected from 16 broiler farms in western France between 2007 and 2017, 6 clinical isolates of human origin, and 8 clinical and 4 nonclinical isolates from other studies that had been previously sequenced (see Table S1 in the supplemental material). A total of 118 genomes sequenced by Illumina technology had sequencing coverage of >160× and an N₅₀ of between 62 kbp and 276 kbp (Table S2A). Hybrid assemblies of Illumina and Nanopore data for 14 genomes allowed the reconstruction of 10 complete genomes (CIRMBP-1212, CIRMBP-1228, CIRMBP-1246, CIRMBP-1261, CIRMBP-1274, CIRMBP-1281, CIRMBP-1283, CIRMBP-1287, CIRMBP-1292, and CIRMBP-1302) and improvements of the genome assembly completeness of 4 others (Table S2A). The estimated average length of the genomes is ~2.4 Mb and varies between ~2.05 and ~2.8 Mbp. Each genome had an average of 2,345 predicted protein coding DNA sequences (CDSs). A total of 277,011 CDSs were annotated. Comparison of the chromosomal architecture of the 11 complete E. cecorum strains using the NCTC12421 strain as a reference revealed that strain CIRMBP-1261 has a large chromosomal inversion of ~1.8 Mbp between the second and sixth rRNA operons (Fig. S1). Strain CIRMBP-1287 also has a chromosomal inversion of 280 kbp from genes DQL78_RS05120 to DQL78_RS06585, involving an insertion sequence (IS) of the IS3 family.

Thirty further nonredundant E. cecorum genomes available at the time of this work were included. These comprised genomes of 9 clinical and 21 nonclinical isolates from Belgium, Germany, and the United States (Table S2B). Comparative genomics analysis of the 351,733 CDSs from the 148 genomes identified 8,523 gene clusters in the pangenome, composed of a strict core genome (present in all E. cecorum genomes) of 1,207 CDSs, an accessory genome of 4,664 CDSs, and a unique genome of 2,652 (31.1%) CDSs (Fig. 1). The pangenome curve displayed an asymptotic trend after the 140-genome iteration, indicating the stabilization of the pangenome within this data set. Consistently, the core genome stabilized after the 125-genome iteration (Fig. 1A). These trends confirm that the genome data set used here provides a comprehensive overview of the gene repertoire of the E. cecorum species. The distribution of the gene clusters revealed that ~47% of the members of the pangenome (n = 3,979, including the unique genes) are present in one to three isolates (Fig. 1B), indicating that the genetic diversity is partly attributable to gene acquisition. This hypothesis is further supported by the high proportion (42%) of genes of unknown function.

The E. cecorum neighbor-joining (BioNJ) phylogenetic tree was constructed using the concatenated sequence of the core genes (Fig. 2). The isolates clustered into five distinct phylogenetic clades (A to E), supported by a bootstrap value of 100% and the higher maximum pairwise genetic distance between clades (0.026 to 0.041) than within clades (<0.021) (Fig. S2). Clade E was further divided into 13 well-supported subclades (E1 to E13) with intrasubclade pairwise genetic distances of <0.012, indicating the greater clonality of these isolates. Clades A to D contain 2 to 16 genomes, only 35% of which were avian clinical isolates. While clades B and D contain mainly genomes of nonclinical isolates from the United States, clade A contains two European clinical human isolates, and clade C contains both nonclinical and clinical poultry isolates. In contrast, clade E comprises 117 genomes, 95% of which belong to avian clinical isolates from the United States, Belgium, Germany, Poland, and France, suggesting the widespread distribution of this clade. Of note, the type strain NCTC12421, isolated from the cecal contents of a dead chicken from a farm in Belgium (1, 31), is part of subclade E3, and subclade E12 contains only avian isolates from the United States. Although the number of isolates per subclade is limited, almost all isolates from subclades E6, E10, and E13 were isolated after 2009, while those from subclades E4 and E11 were isolated before 2014 and after 2015, respectively, indicating that the dominant subclasses have varied over the years. Due to sequence data availability, single nucleotide polymorphism (SNP) analysis was possible only for genomes sequenced in this study, thus excluding genomes from clades B and D and subclade E12. Of the 65,226 SNPs of the core genome, 2,443, 168, and 62 were specific to isolates of clades A (n = 2), C (n = 9), and E (n = 107), respectively. Of the 13 nonsynonymous clade E-specific SNPs, 8 are nonconservative, and 2 are predicted to be nonneutral in a phage shock protein (PspC [DQL78_RS04285 in NCTC12421]) and an ATP-binding cassette domain-containing protein (DQL78_RS04370 in NCTC12421).

Altogether, these data revealed that clinical isolates of E. cecorum from poultry from different countries are grouped phylogenetically and belong mainly to clade E, confirming their clonality and suggesting either the dissemination of well-adapted clones or convergent selection/adaptation to poultry genetics and breeding methods.

To further explore the genetic basis for the epidemic success of the E. cecorum clade E isolates, we searched for genes significantly associated with clade membership. The majority of unique genes (71.9%) are carried by isolates from clades A, B, C, and D, indicating higher genetic diversity in these clades. However, this trend may result from the enrichment of clade E isolates in the collection. Consistently, we identified 97 accessory genes significantly associated with clade E isolates (Table S3A). The most abundant group of the 83 clade E-enriched genes after the unclassified hypothetical protein-encoding genes (n = 17) are predicted carbohydrate metabolism and transport genes (n = 13) and cell wall/membrane/envelope biogenesis genes (n = 12). A high proportion of the clade E-enriched genes are clustered loci, including the capsular polysaccharide (CPS) biosynthesis locus (CIRMBP1228_00568 to CIRMBP1228_00581 in strain CIRMBP-1228), a large gene cluster comprising biotin biosynthesis genes (CIRMBP1228_01807 to CIRMBP1228_01837), and one putative operon of carbohydrate metabolism (CIRMBP1228_02727 to CIRMBP1228_02733). Alignment of the gene products of the representative capsule biosynthesis loci from the complete genomes identified two closely related loci represented by the genomes of CIRMBP-1228 and CIRMBP-1212 (Fig. S3), which account for 66.7 and 14.5% of the clade E isolates, respectively. The remaining 14 genes significantly associated with the origin of the isolates (clinical and nonclinical poultry isolates or clinical human isolates) are predominantly absent in clade E isolates. They include genes encoding a pseudouridine synthase (CIRMBP1294_00547), a predicted nucleotide sugar dehydrogenase (CIRMBP1294_00386), a diacylglycerol kinase family lipid kinase (CIRMBP1294_00623), and two stress-related proteins (CIRMBP1294_00560 and CIRMBP1294_00758). In line with the large proportion of clade E clinical isolates in the collection, 30 clade E-enriched genes are also enriched in clinical poultry isolates compared to nonclinical poultry and clinical human isolates (Table S3B). These include genes of the CPS biosynthesis locus (CIRMBP1228_00573, CIRMBP1228_00574, and CIRMBP1228_00575), 27 genes of the biotin gene cluster, and an H protein gene of the glycine cleavage system generally involved in protein lipoylation. A total of 65 genes are significantly associated with the origin of the isolates. In addition to those enriched in clade E isolates, other genes enriched in avian clinical isolates encode a phosphoenolpyruvate:carbohydrate phosphotransferase system (PTS), a transketolase (CIRMBP1228_00604 to CIRMBP1228_00607), as well as proteins of unknown function. Although no specific gene signature for avian clinical isolates was found, we identified 6 accessory genes (CIRMBP1228_00573, CIRMBP1228_00586, CIRMBP1228_00757, CIRMBP1228_01816, CIRMBP1228_02735, and CIRMBP1283_01819) that allowed the identification of 94% of avian clinical isolates (Table S3C). With the exception of 7 avian clinical isolates from clade C or D, all other genomes of clinical isolates have an average of 4 of the 6 genes mentioned above (range, 2 to 6 genes), while nonclinical isolates have 1 gene at most. We also found 12 genes that may be enriched in clinical human isolates. These include 8 genes predicted to be involved in the import and utilization of ascorbate under anaerobic conditions (32, 33).

Despite the difficulty in identifying genes specific to the origin of the isolates, this work highlights six genes encoding a glycosyltransferase (CIRMBP1228_00573), two PTS EIIC components (CIRMBP1320_01424 and CIRMBP1228_02735), and three hypothetical proteins (CIRMBP1228_00586, CIRMBP1228_01816, and CIRMBP1228_00757) that can discriminate most of the isolates associated with poultry disease from those that are not.

Next, we searched the 148 sequenced genomes for the presence of acquired ARGs using ResFinder, BLAST, and the PLSDB database. A total of 18 ARGs were identified (Fig. 3). Tetracycline and macrolide-lincosamide-streptogramin B (MLS) resistance genes were detected in 95% and 75% of the isolates, respectively. In total, 70% of the genomes had at least one gene for resistance to these two families of antimicrobials, with tet(M) and erm(B) being the most prevalent. Genes for resistance to other antimicrobial classes such as aminoglycosides (10%), bacitracin (30%), and vancomycin (0.6%) were also detected. Only two isolates (CIRMBP-1244 and CIRMBP-1314) did not contain any ARG searched. Overall, 39 genomes carried genes conferring resistance to three antimicrobial families, and 7 genomes had at least four genes for resistance to different antimicrobial families (Fig. 3). The most common combinations of ARGs in multidrug-resistant isolates cover the tetracycline and MLS families in combination with aminoglycosides in clades C and D or with bacitracin in clades C and D and subclades E10, E11, and E12. Aminoglycoside resistance genes are carried by nonclinical genomes with enrichment in U.S. strains. The near-systematic presence of ARGs in the genomes of E. cecorum isolates suggests that the species may have undergone strong selection for antibiotics, tetracyclines and MLS in particular, but that resistance to these two antibiotic families did not provide a selective advantage to clinical isolates.

Several E. cecorum genome regions carrying ARGs conserved in Gram-positive bacteria were identified in the PLSDB database, but none were carried by a previously identified plasmid. This analysis highlighted different ARG clusters. They comprised the erm(B), vat(D), and msr(D) genes; the tet(L), tet(M), and bcr genes; or the tet(L), tet(M), ant(6), and bcr genes. Further characterization of the vancomycin resistance operon carried by strain CIRMBP-1294 revealed that the vanA operon was integrated into the chromosome 12 kbp away from the genes that confer resistance to narasin and erythromycin in between (Fig. S4). This chromosomal region was closely related to that of the E. faecium plasmid pVEF3 (34) identified in broiler isolates from Sweden (35) and was flanked by two transposase genes of the IS3 family. However, no conjugative or mobilization element was detected close to this ARG cluster. Although about 20 plasmid-related gene clusters were identified in the pangenome, a single putative 4.4-kbp plasmid was detected in the complete genome of strain CIRMBP-1292 and 6 other isolates (CIRMBP-1233, CIRMBP-1259, CIRMBP-1286, CIRMBP-1289, and CIRMBP-1304), but no ARG was detected on this small plasmid. Together, these results suggest that plasmids are rare in this bacterial collection and that ARGs are located on mobile genomic islands (GIs).

In Gram-positive bacteria, ARGs and related transposons are frequently integrated into complex MGEs forming integrative conjugative elements (ICEs) or integrative mobilizable elements (IMEs) that may be difficult to identify in draft genomes (36). To evaluate the contribution of these elements to the spread of ARGs and E. cecorum genome diversity, we predicted ICEs and IMEs in the 11 complete genomes and the 2 large contigs of the genome of CIRMBP-1320 using ICEscreen (37). ICEscreen ICE and IME detection relies on the presence of signature protein CDSs (integrase, coupling protein, relaxase, and VirB4) grouped on the genomes and previously shown to be a valuable clue to the presence of an integrative element (38, 39). ICEs were defined by the superfamily and family of their signature proteins (38). Each E. cecorum genome contained at least one ICE (Fig. 4 and Table S4A). In total, 33 mobile genetic elements, including 5 types of ICEs (belonging to the Tn916, Tn5252, and TnGBS1 superfamilies), 2 IMEs, 3 Tn917 elements, and 1 partial conjugative element, were identified. Their sizes ranged from ~5 kbp for Tn917 to ~103 kbp for the TnGBS1 ICE of the genome of CIRMBP-1302. The Tn916 ICE of the genomes of CIRMBP-1228 and CIRMBP-1281, integrated upstream of the 30S ribosomal protein S6-encoding gene, included another ICE of an undescribed family encoding a DDE transposase, a MobC-like relaxase, a VirD4-like coupling protein, and a type IV secretion system (T4SS) protein, VirB4. We also investigated the presence of other GIs in the complete genomes (Fig. 4 and Table S4B). A total of 42 GIs were identified, with sizes ranging from ~8 kbp to ~122 kbp for the complex GI that comprises an ICE (Tn916, ICEBs1, and ICESt3) in the genome of CIRMBP-1228. Each genome harbored from 3 ICE-related elements or genomic elements (strain CIRMBP-1212 of subclade E5) up to 12 (strain CIRMBP-1228 of subclade E4). They were mainly integrated into intergenic regions or the 3′ ends of genes, with no effect on the encoded sequences. As expected, all Tn917 elements and ICEs of the Tn916 family carried erm(B) and tet(M), respectively. Other ARGs, such as tet(L), tet(O), ant(6), and the bcr operon, were carried on Tn916-related ICEs or on GIs. Besides genes involved in GI transfer, a few had predicted functions related to biotin biosynthesis, restriction-modification enzymes, cadmium/arsenate resistance, toxin-antitoxin systems, redox enzymes, type 2 secretion systems, and flagella and chemotaxis. Analysis of the distributions of the ICE-related elements and the GIs in the other 136 E. cecorum genomes (Table S4B) revealed three highly dispersed elements: the Tn916-related ICE of strain CIRMBP-1292 in 113 genomes, the GI inserted near the rpsB gene and carrying the biotin biosynthesis genes found in 90 genomes, and the transposon Tn917 in 65 genomes (Table S4B). The GI inserted between dnaX and sufB of strains CIRMBP-1287 and CIRMBP-1274 was highly prevalent in subclades E6 and E10, respectively. The same is true for the cognate element of CIRMBP-1283 and CIRMBP-1302 that is detected in all isolates of subclade E11, to which CIRMBP-1283 belongs. Other elements are enriched in specific subclades, such as the GI near the 23S rRNA methyltransferase gene rlmD and ICEs of the TnGBS1 and Tn916 families of CIRMBP-1274 enriched in subclade E6, the ICE of the Tn916 family of CIRMBP-1212 enriched in subclade E13, and the GI inserted at the 3′ end of CIRMBP1228_01030 and the ICE of the TnGBS2 family of CIRMBP-1228 close to gene CIRMBP1228_01622 enriched in subclade E4. According to their gene contents, 15 elements appear to be more strain specific in strains CIRMBP-1320 (n = 4), CIRMBP-1281 (n = 4), CIRMBP-1292 (n = 3), and NCTC12421 (n = 4). Although the genomes used to identify ICEs and GIs are not representative of all clades, homologs were found in almost all clades except for clade B, which gathers only five nonclinical poultry isolates from the United States. However, these genomic elements were less conserved in isolates from the United States, suggesting that they were acquired independently. We conclude that the majority of the identified ICEs and GIs are enriched in some clades, while only a few are shared across clades.

As prophages are easier to detect than GIs and mobile genetic elements, prophages were searched for in all E. cecorum genomes sequenced, including the genome of NCTC12421. Totals of 103 complete prophages and 30 incomplete prophages were identified in 78 (66.1%) E. cecorum genomes (Table S5). The prophage size ranged between 31 and 59 kbp. The prophages were related to 19 different types that varied across the phylogenetic groups. Subclade E4 was enriched with prophages related to Faecalibacterium phage oengus (40), which was the most commonly found prophage (n = 31), followed by prophages related to the Siphoviridae temperate phage EFC-1 of Enterococcus faecalis (41) (n = 16). Prophages related to the Myoviridae temperate bacteriophage EJ-1 of Streptococcus pneumoniae (42) (n = 13) were found mostly in the genomes of subclade E10, while those related to the Siphoviridae temperate phage 50101 of Lactococcus (n = 12) were most prevalent in subclade E6. Conversely, clade C had the most diverse prophages, and subclades E10 and E12 contained the smallest number of prophages. Prophages related to Siphoviridae temperate phage 5093 of Streptococcus thermophilus (43) (n = 9) were identified in several phylogroups. A total of 13 different integration sites for temperate bacteriophages were identified, mainly in intergenic regions and at the 3′ ends of genes, without changing the reading frame (Table S5). However, the prophage related to PHAGE_Strept_EJ_1_NC_005294 integrated into the PreQ1 riboswitch is likely to change the expression of the downstream nucleoside hydrolase gene. Like ICEs and GIs, prophages are less prevalent in the ~700 kbp surrounding the origin of replication of the chromosome (Fig. 4), indicating that MGEs are not randomly distributed in the E. cecorum genome.

To obtain phenotypic insight into the 118 isolates, we independently examined phenotypes relevant to E. cecorum pathogenesis: biofilm robustness, adhesion to type II collagen, and growth in chicken serum (see Text S1 in the supplemental material). Hierarchical clustering of these phenotypes (Fig. 5) revealed 17, 11, and 9 groups of strains for biofilm robustness, type II collagen adhesion, and growth in chicken serum, respectively. The lack of concordance between these phenotypic groups and clades strongly suggests that none of the phenotypes discriminate between clades. In an attempt to rank isolates, the most robust biofilm-forming isolates were CIRMBP-1302 (subclade E13), CIRMBP-1277 (clade C), and CIRMBP-1228 (subclade E4), and the most fragile biofilm-forming isolates were CIRMBP-1204 (subclade E1), CIRMBP-1205 and CIRMBP-1211 (subclade E4), and CIRMBP-1225 (subclade E6). The most collagen-adherent isolates were CIRMBP-1274 (subclade E6) and CIRMBP-1277 (clade C). All of the E. cecorum isolates had a serum growth index below the index of the Escherichia coli control strain sensitive to the bactericidal effect (190.34 ± 123.88) of serum and in the range of that of the nonsensitive E. coli control strain (3.172 ± 4.2). The E. cecorum isolates with the lowest indices were strains CIRMBP-1318 (0.229 ± 0.21) and CIRMBP-1206 (0.467 ± 0.001), belonging to subclades E3 and E4, respectively. The E. cecorum strains with the highest indices were strains CIRMBP-1312 (3.85 ± 0.44) of subclade E12, CIRMBP-1197 (4.21 ± 2.2) of clade A, and CIRMBP-1298 (5.71 ± 2.9) of subclade E10. Taken together, these results indicate that the growth of E. cecorum is not affected by the presence of chicken serum. Besides strain CIRMBP-1277, which formed strong biofilms and had a high capacity to adhere to collagen, no relationship between the expression of the three phenotypes and the strains was observed. Moreover, no accessory genes were associated with the phenotypic groups or with binary transformed values, suggesting functional redundancy between genes and/or differential gene expression between isolates.

To evaluate the pathogenic potential of E. cecorum strains isolated from diseased broilers in French farms, we used the chicken embryo lethality assay (CELA) reported previously by Borst et al. using strain CIRMBP-1309 (original name, CE3) and strain CIRMBP-1311 (original name, SA2) as biological controls for a nonpathogenic and pathogenic poultry isolate, respectively (21). We tested the 11 E. cecorum strains with the most complete genomes (see Text S1 in the supplemental material). No embryo died in the control group. The two E. cecorum control strains behaved as expected: 53% of the embryos were still alive 6 days after inoculation with the commensal strain CIRMBP-1309, whereas 100% of the embryos died after 2 days when inoculated with strain CIRMBP-1311 isolated from spondylitis lesions (Fig. 6). Strain CIRMBP-1294, isolated from infected vertebrae, and strain CIRMBP-1320, isolated from human infection, induced <27% embryonic lethality, indicative of low virulence in the CELA. In contrast, strains CIRMBP-1228, CIRMBP-1274, CIRMBP-1292, CIRMBP-1302, and CIRMBP-1304, isolated mainly from infected vertebrae, showed the same virulence as that of strain CIRMBP-1311, with an overall average lethality of 85% at day 2 postinfection. They thus could be considered virulent isolates. The lethality of strains CIRMBP-1212, CIRMBP-1281, CIRMBP-1283, and CIRMBP-1287 was intermediate and higher than that observed for the least virulent strains CIRMBP-1294 (P values of between 0.026 and 0.064) and CIRMBP-1320 (P values of between 0.026 and 0.073). Close analysis of the gene contents of the 10 strains exhibiting virulence and the 3 strains considered nonvirulent in the CELA identified 33 genes, 18 of which were already found to be enriched in clade E isolates and 6 of which were enriched in avian clinical isolates. The predicted function of these genes indicates a bias in favor of carbohydrate transport and metabolism (Table S3D).

In line with the data obtained for biofilm robustness, binding to collagen, and growth in chicken serum, the clear gradient of virulence observed in the CELA for clinical isolates supports the hypothesis of the multifactorial nature of E. cecorum pathogenesis.

A total of 118 strains were collected from various laboratories and deposited at the International Center for Microbial Resources—Bacterial Pathogens (CIRM-BP) (https://www6.inrae.fr/cirm_eng/BRC-collection-and-catalogue/CIRM-BP) (see Table S1 in the supplemental material). The majority of them (n = 100) were isolated between 2007 and 2017 from diseased birds from 16 broiler farms located in Brittany, France, and were provided by Labofarm (Loudéac, France). Other poultry strains were isolated in Poland (n = 5), Belgium (n = 1), and the United States (n = 6). Six strains from human infections were isolated in France (n = 3), Belgium (n = 2), and Germany (n = 1). The sources and the original names of the strains from abroad are indicated in Table S1. Additional details are available in the supplemental material.

All genomic DNA was subjected to random shotgun library preparation using the TruSeq DNA PCR-free kit (Illumina). Ready-to-load libraries were sequenced on an Illumina MiSeq or HiSeq 3000 platform at GeT-PlaGe (Toulouse, France) and a HiSeq 2500 platform at Eurofins Genomics (Germany) using 150-bp paired-end chemistry. DNA of 14 isolates was also sequenced using Oxford Nanopore Technology platforms. The preparation of libraries and sequencing were performed at the GeT-PlaGe core facility (INRAE, Toulouse, France) or MetaGenoPolis (INRAE, Jouy-en-Josas, France) according to the manufacturer’s instructions. Additional details are available in the supplemental material.

The 104 genomes with only Illumina reads were assembled using RiboSeed v0.4.73 (91). RiboSeed uses a reference genome to resolve rRNA operons and globally improve the whole-genome assembly. Assemblies were performed using NCTC12421 (GenBank accession number NZ_LS483306) as the reference genome for rRNA operons, using SPAdes v3.13.0 (92) as the assembler in “careful” mode with k values of 21, 33, 55, 77, and 99. The 14 genomes with Illumina and Nanopore reads were assembled using Unicycler 0.4.4 (93), an assembly pipeline for bacterial genomes that uses SPAdes for short-read assembly and Miniasm and Racon for long-read assemblies and polishing. Unicycler was launched with default parameters. Genomes were compared using Mummer 3.2.3 (94) with parameter to compute all maximal matches (-maxmatch) and a minimal length of match of 100 nucleotides (-l parameter). Inversions in the genomes were checked by realigning long reads on the assembly using minimap 2.2.17 (95). The veracity of the assembly on the breaking points of the inversion was manually checked by visual inspection of the alignment of the long reads encompassing the inversion.

Genome annotation was conducted with Prokka v1.12 (96). First, coding DNA sequences were identified on contigs longer than 200 bp by Prodigal v2.6.1 (97), which penalized CDSs shorter than 250 bp in order to filter out false positives. CDSs were first annotated (–proteins Prokka parameter) using a protein bank extracted from all of the complete Enterococcus genomes in the RefSeq database (4,198 genomes retrieved in April 2020). Annotations from hits with an E value cutoff of 10⁻⁹ and 80% coverage were transferred. CDSs with no hits on this bank were annotated using the Prokka default workflow and data banks.

Pangenome analysis was performed by comparing 118 de novo-sequenced genomes of E. cecorum and 30 public genomes from the NCBI with N₅₀ values above 20 kb (January 2021). One reference genome was available in the NCBI nucleotide database, NCTC12421. Protein clustering was performed by Roary (98) v-3.12.0 with 94% identity, a “percentage of isolates a gene must be in to be core” of 100%, and the parameter “without split paralogs.” Genes classified as core were genes present in all 148 genomes. Accessory genes were all other genes present in 147 or fewer genomes. The gene accumulation curves were produced with ggplot2 (99) from Roary data. The phylogenomic analysis was performed using the E. cecorum core genome. The 1,207 core genes of E. cecorum were aligned using a codon-aware alignment produced by PRANK (v170427); an unrooted tree was then constructed using the BioNJ algorithm (100) in SeaView (v4.2) (101), using the Jukes-Cantor distance and 1,000 bootstrap replicates. Clades were determined using the Jukes-Cantor distance between aligned core genes of <0.021 and a bootstrap value of >75%. The graphical representation of the phylogenetic tree was performed with iTOL (v4) (102).

ARG research was performed using the ResFinder (v2.1) (103) tool and database (26 April 2019) for 90% gene identity and 60% coverage. BLASTN was used to detect the bcr-like gene (uppP_2) and to compare the vanA locus with pVEF3 of E. faecium 01_233 (34). Positive hits had at least 90% nucleotide identity and 60% coverage. GenoPlotR (v0.8.11) was used to visualize BLASTN results with >90% identity (34).

ICEs and IMEs were detected in the 11 complete genomes and the 2 large contigs of the genome of CIRMBP-1320 using ICEscreen (https://icescreen.migale.inrae.fr) and then inspected visually for delineation. Genomic islands corresponded to large insertions of more than 10 CDSs resulting in a synteny break between two genomes. Prophage prediction was performed using the prophage detection tools PHASTER (Phage Search Tool, Enhanced Release) (104) and VIBRANT (Virus Identification by Iterative Annotation) (105). Only predicted prophages with a prediction score of ≥100 with PHASTER or VIBRANT were retained and manually inspected to determine the attachment and integration sites in reference to the NCTC12421 genome (GenBank accession number NZ_LS483306). PHASTER was further used to identify the most similar phage genomes.

Estimated marginal means of biovolumes of biofilm and adhesion to collagen and serum growth indices were adjusted for experiment and strain factors using a linear model with the emmeans R package (1.4). Strains with similar Euclidian distances between estimated marginal means were grouped using a hierarchical clustering algorithm with average linkages. Clusters of strains were defined by cutting the dendrogram dynamically using the DynamicTreeCut R package (1.63-1) (106).

Single nucleotide polymorphism analysis was performed using GATK version 4.4.2 (107) with core genes of NCTC12421 as the reference. Default parameters were applied with the HaplotypeCaller module to call SNPs. GWASs were performed using TreeWas (1.0) (108) to identify genetic loci (SNPs obtained from GATK analysis and gene presence/absence obtained from Roary analysis) associated with clade membership (clade A, B, C, D, or E), clinical origin (clinical and nonclinical poultry isolates and clinical human isolates), biofilm robustness (binary strong values versus others, binary weak values versus others, and phenotypic groups obtained by clustering), and collagen type II adhesion and/or growth in chicken serum (similarly to biofilm robustness). Significant genetic loci corresponded to a P value of ≤0.05 according to a terminal test.

A set of criteria was applied to the TreeWas output to select clade-specific SNPs or enriched genes. The criteria applied to select SNPs were more stringent, and only clade-specific SNPs were retained (sensitivity equal to 1 or 0 and specificity equal to 1 or 0). Genes whose presence was associated with clade membership (clade A, C, or E) or clinical origin (clinical and nonclinical poultry isolates and clinical human isolates), biofilm robustness, adhesion to collagen type II, or growth in chicken serum had sensitivity and specificity scores above 0.66. In addition, genes whose absence was associated with clade membership or clinical origin, biofilm robustness, adhesion to collagen type II, or growth in chicken serum had sensitivity and specificity scores below 0.33.

All genomic data have been deposited in the EMBL ENA database under project accession number PRJEB50514. Accession numbers of raw reads and assembled genomes are available in Table S2A.