Genomic evaluation of the probiotic and pathogenic features of Enterococcus faecalis from human breast milk and comparison with the isolates from animal milk and clinical specimens

Sampling and culture-depending screening

In this study, HBM samples were collected from three healthy mothers between 7 and 10 days’ post full-term pregnancy in Saudi Arabia. These samples were called HBM_SA samples. The average age of mothers was 31  ± 2.6 years. The milk samples were collected under the sterile conditions described previously (Hunt et al., 2011). Ethical approval was obtained from the Institutional Review Board of the Dr. Soliman Fakeeh Hospital, Jeddah, under the number 77/IRB/2020. Four culture media of sheep blood agar, R2A agar, MacConkey agar, and MRS agar supplemented with 50% pasteurized milk were used to isolate bacteria from the HBM_SA samples. A serial dilution approach in peptone water was adopted, and 100 µl of the sample was spread on each culture medium. The colony-forming unit (CFU) was calculated after 48 h of incubation at 37 °C. The colonies were purified by sub-culturing and preserved in autoclaved 15% glycerol and 1% skim milk in distilled water at −80 °C.

The purified colonies were identified by MALDI-TOF based VITEK-MS (BioMérieux, France) system described previously (Yasir et al., 2022). The VITEK-MS run was validated using the reference strain Escherichia coli ATCC 25922, following the manufacturer’s standard protocol (BioMérieux, Marcy-l’Étoile, France). Antimicrobial susceptibility and minimum inhibitory concentration (MIC) of E. faecalis isolates were performed using the VITEK 2 (BioMérieux, Marcy-l’Étoile, France) system with a specific AST-GP2 card for Gram-positive bacteria. MIC results were interpreted according to the Clinical and Laboratory Standards Institute guidelines 2022.

Genome sequencing and bioinformatic analysis

For genome sequencing, one E. faecalis isolate from each HBM_SA sample was selected that was purified from the MRS agar supporting probiotic bacteria growth. Briefly, genomic DNA was extracted from the fresh overnight culture of E. faecalis isolates using UltraClean^® Microbial DNA isolation kit (Qiagen, Hilden, Germany), and genome sequencing was performed as described previously (Yasir et al., 2020). A Nextera DNA Flex Library Prep Kit (Illumina, Inc., San Diego, CA, USA) was utilized for library preparation, and sequencing was conducted with a V3 kit, employing 2 × 300 bp chemistry on the MiSeq platform (Illumina, Inc., USA). The quality assessment of raw sequence reads was carried out using FastQC 0.12.0, and sequence trimming was performed using the Trimmomatic v0.32 tool (Bolger, Lohse & Usadel, 2014). Contig assemblies were prepared using the SPAdes 3.15.3 program (Prjibelski et al., 2020).

Additionally, we retrieved 15 E. faecalis genome assemblies from the NCBI microbial genomic database (Table S1). Among them, six were isolated from HBM, seven from animal milk, and two were reference strains of Symbioflor 1 (probiotic) and V583 (vancomycin-resistant clinical strain). The E. faecalis genomes included in this study were from isolates recovered from the milk of healthy women and animal milk, excluding those recovered from animals with mastitis or other reported infections. Furthermore, analysis was performed from 43 E. faecalis isolates we previously recovered from clinical specimens in the western region of Saudi Arabia (Farman et al., 2019). Annotation of the genome was performed using BV-BRC 3.33.16 (Olson et al., 2023), and the PubMLST tool was used to determine multilocus sequence typing (MLST) (Larsen et al., 2012). Single nucleotide polymorphisms (SNPs) were identified in the core genomes and utilized to construct a maximum likelihood phylogenetic tree using CSI Phylogeny 1.4 with default parameters (Kaas et al., 2014). The phylogenetic tree was visualized using the Interactive Tree of Life (iTOL v6) tool (Letunic & Bork, 2024). The acquired antimicrobial resistance genes (ARGs) and chromosomal point mutations were determined employing ResFinder 4.6.0 (Bortolaia et al., 2020). Virulence genes were retrieved from the E. faecalis genomes using VirulenceFinder 2.0 (Malberg Tetzschner et al., 2020). Using PlasmidFinder 1.3, the origin of replications of the plasmids were identified from the assembled genomes (Carattoli & Hasman, 2020). Mobile genetic elements in connection to ARGs and virulence factors were determined using MobileElementFinder v1.0.3 (Johansson et al., 2021), and PHASTER was used to detect putative prophage elements (Wishart et al., 2023). Secondary metabolites and post-translationally modified peptides (RiPPs) associated gene clusters were identified in the genome assemblies using antiSMASH 7.0 (Blin et al., 2023). Sequence reads of isolates from this study were deposited into GenBank under BioProject ID. PRJNA1059526.

Bacterial count and community analysis in the HBM_SA samples

In the bacterial count analysis, an average of 2.3 × 10⁴  ± 1.9 × 10⁴ cfu/ml was obtained on sheep blood agar, and 2.7 × 10⁴  ± 2.5 × 10⁴ cfu/ml was obtained on R2A agar. On probiotic supporting MRS agar, 1.6 × 10²  ± 75.1 cfu/ml was obtained, whereas no growth was observed on selective MacConkey agar to detect pathogenic E. coli. Eighty-seven isolates were purified from the HBM_SA samples and were predominantly classified into Staphylococcus epidermidis and Staphylococcus hominis. E. faecalis isolates were retrieved on MRS agar in three studied samples. Staphylococcus aureus was found in the HMD11_46M sample, and Staphylococcus haemolyticus was found in the HMD9_11M sample.

Genomic annotation of E. faecalis isolates

The three E. faecalis isolates recovered from HBM_SA samples were phenotypically resistant to erythromycin and tetracycline and sensitive to vancomycin (Table S2). Genome sequencing of three E. faecalis isolates, one from each HBM_SA sample, was performed. The sequenced E. faecalis isolates Efa_HMD9_11M, Efa_HMD11_46M, and Efa_HMD12_49M were assembled into 58, 49, and 56 contigs, respectively. The average G+C content of the isolates was 37.2  ± 0.02%, and the average total length was 3,068,846  ± 36,594 bp (Fig. 1). According to the BV-BRC 3.33.16 (Olson et al., 2023) analysis, the genome quality was good, with a CheckM completeness of 100, and no CheckM contamination was detected in two of the isolates (Table S3). However, 0.5 CheckM contamination was detected in Efa_HMD9_11M. Protein-coding sequences (CDS) were detected in the range of 3,024 to 3,121. At the genomic level, the average nucleotide identity (ANI) value among the three E. faecalis isolates was 99.98 ± 0.02%, and 99.00 ± 0.02% with the reference probiotic strains of E. faecalis Symbioflor 1 (Fig. 1). The ANI value was 98.8 ± 0.3% with the clinical reference strain of E. faecalis V583. No substantial difference was detected in the G+C content of the E. faecalis isolates’ genomes from human milk (37.5 ± 0.1%), animal milk (37.5 ± 0.2%) retrieved from the NCBI genome database, and clinical isolates (37.4 ± 0.1%). The CDS were detected in the range of 2,580 to 3,330 (Table S3).

E. faecalis genomes analysis for secondary metabolites genes

Several secondary metabolite genes were detected in the genomes of the E. faecalis isolates analyzed in this study and were mainly carried by the post-translationally modified peptides (RiPPs) (Table S4). Two genomic regions were found in four isolates from HBM, coding for RiPPs, including three isolates from the HBM_SA samples and the strain APC 3825 (Table S4). The lanthipeptide-class-II RiPPs were commonly found in these isolates, revealing 100% similarity with the cytolysin ClyLl/cytolysin ClyLs biosynthetic gene cluster. RiPP-like biosynthetic gene clusters (BGCs) were also detected in these four HBM isolates, showing 22% similarity with the carnobacteriocin XY (Table S4). No BGCs were detected in the genome assemblies of HBM isolates from China. Compared to isolates from human milk, cyclic-lactone-autoinducers were found in four E. faecalis isolates from animal milk and 20 clinical isolates. Lassopeptide (RiPPs) was also found in both animal and clinical isolates. RiPP-like (carnobacteriocin XY) was explicitly found in isolates from HBM (Table S4). Overall, a limited number of BGCs were detected in each genome, and no consistent pattern was found among the identified BGCs in relation to the source of isolation of these E. faecalis isolates.

Comparative functional analysis of E. faecalis genomes

From the subsystem functions, 217 were commonly present among the genome assemblies of the three E. faecalis isolated from HBM_SA samples and reference strains Symbioflor 1 and V583 (Table S5). Nine subsystem functions were unique to V583, such as resistance to vancomycin and teicoplanin, potassium-transporting ATPase, toxin-antitoxin replicon stabilization, ibrA and ibrB co-activators of prophage gene expression, and aminoglycoside-modifying enzymes (Table S5). Bile hydrolysis and type I restriction-modification systems were not found in V583 but were present in Symbioflor 1 and E. faecalis isolates from HBM_SA samples. All detected subsystem functions were common in the HBM_SA E. faecalis isolates (Table S5). Moreover, inositol catabolism, formaldehyde assimilation (ribulose monophosphate pathway), and tetracycline resistance were unique to isolates from HBM_SA samples and were not detected in both reference strains (Table S5). In comparing HBM_SA isolates with other HBM isolates from China, 224 out of 230 subsystem functions were commonly present, whereas 229 subsystem functions found in the APC 3825 isolate from Ireland were commonly found in the HBM_SA isolates (Table S5). Interestingly, the five E. faecalis isolates from HBM in China commonly shared all the 226 detected subsystem functions. Subsystem functions like hyaluronate utilization, inositol catabolism, formaldehyde assimilation, and tetracycline resistance were uniquely detected in the HBM_SA and APC 3825 isolates compared to HBM isolates from China. The hyaluronic acid capsule and hyaluronic acid-containing cell wall-associated genes were found in HBM isolates from China and not in the HBM_SA isolates. In total, 190 out of 248 subsystem functions were common in the isolates recovered from human and animal milk, whereas 214 out of 259 subsystem functions were common between HBM_SA isolates and clinical isolates recovered from patient specimens in Saudi Arabia (Table S5). Overall, 188 out of 263 subsystem functions were common in all the analyzed genome assemblies (Table S5).

In comparing gene families from the CAZy database, which catalogs microbial enzymes catalyzing carbohydrates, we identified 43 CAZy families and sub-families that were common in the HBM_SA isolates. Among these, 35 were shared with both reference strains, Symbioflor 1 and V583 (Table S6). The GH24 family was specifically found in the HBM_SA isolates. In contrast, seven CAZy families and sub-families (CE4, GH136, GH154, GH31, GH88, GT26, PL12_1) were common among the HBM_SA isolates and V583 strains but were not detected in the Symbioflor 1 genome (Table S6). As a result of the functional annotation conducted with RAST, carbohydrate metabolism was the most enriched metabolic category in the E. faecalis isolates from the HBM_SA samples and in both reference strains (Fig. 2A). These isolates exhibited the ability to metabolize di- and oligosaccharides, monosaccharides, and organic acids (Fig. 2B). The second-highest percentage of PEGs was in the protein metabolism category, followed by amino acids and derivatives, and nucleosides and nucleotides (Fig. 2A). Genes encoding for lipoic acid, pyridoxine, riboflavin, FMN, FAD, and tetrapyrroles from the subclass cofactors, vitamins, prosthetic groups, pigments, as well as genes for fatty acids, isoprenoids, and triacylglycerols, were found in the HBM_SA isolates and both reference strains (Fig. 2B).

Phylogenetic and multilocus sequence typing analysis

In a phylogenetic tree based on SNPs, HBM_SA isolates were clustered together, and several clinical isolates previously recovered from patient specimens in Saudi Arabia (Fig. 3). The isolates from HBM in China were clustered along the QH5 isolate from Yalk milk with Symbioflor 1. Overall, the isolates from HBM and animal milk clustered distinctly from the highly characterized virulent strain V583 (Fig. 3). The MLST analysis classified three E. faecalis isolates from HBM_SA samples and the APC 3825 strain isolated from HBM in Ireland as sequence type ST179 (Fig. 4A). Five E. faecalis isolates from HBM in China were classified as ST25. Variations were found in the sequence types of E. faecalis isolates from animal milk and clinical samples. None of the isolates from animal milk were classified as ST179 or ST25, whereas 11 isolates from clinical specimens were classified as ST179 (Fig. 4A). The probiotic reference strain Symbioflor 1 was classified as ST248, and the clinical reference strain for vancomycin resistance, V583, was classified as ST6 (Fig. 4A).

Comparative antimicrobial resistance genes analysis

In isolates from HBM, the fewest number of ARGs were identified. The lsaA gene, intrinsic to E. faecalis and conferring resistance to lincomycin, clindamycin, dalfopristin, pristinamycin IIA, and virginiamycin M, was commonly found in the studied isolates (Fig. 4A). The tetM and ermB genes, responsible for tetracycline and macrolide antibiotic resistance, were detected in the HBM_SA isolates. E. faecalis HBM isolates from China carried only the lsaA gene, similar to the Symbioflor 1 strain (Fig. 4A). The APC 3825 isolate carried tetM gene. However, the ARG pattern of E. faecalis isolates from animal milk of sheep and yalk resembled that of HBM isolates (Fig. 4A). Two isolates from cattle milk carried multiple resistance genes. In the isolates from clinical samples, 30 isolates were carrying ≥6 ARGs, conferring resistance to clinically important antibiotics from the classes of aminoglycoside, phenicol, diaminopyrimidine, and macrolide. Furthermore, no mutation was detected in the gyrA and parC genes, causing quinolone resistance in the HBM_SA isolates, similar to the Symbioflor 1 strain (Fig. 4B). However, the known parC p.S517N mutation was detected in HBM isolates from China (Fig. 4B). Most isolates from animal milk did not carry mutations in the gyrA and parC genes (Fig. 4B). The LJM05 strain from cattle carried an unknown mutation parC p.V13I (GTA − > ATA, V − > I). In clinical isolates, parC p.S80I mutation was detected in 18 isolates, and gyrA p.S83I mutation was detected in 11 isolates (Fig. 4B). Additionally, ten unknown mutations were detected in several clinical isolates (Fig. 4B).

Virulence factors associated genes analysis

In comparing virulence genes in E. faecalis isolates from HBM, 13 genes were commonly found in isolates from HBM_SA samples, APC 3825, and HBM isolates from China (Fig. 5). However, the cylM, cylL, cylA, and acm genes were detected in the HBM_SA isolates and APC 3825 but not in the HBM isolates from China (Fig. 5, Table S1). Compared to the reference probiotic strain, 12 virulence genes detected in HBM isolates were common to Symbioflor 1. The virulence gene profiles in the isolates from animal milk matched those of HBM isolates, with ten virulence genes commonly detected among the isolates from HBM and animal milk (Fig. 5). In clinical isolates, more variation was noticed in the virulence genes detected in the 13 to 22 genes range (Fig. 5). The virulence genes fsrB and hylB were found in clinical reference strains V583 but were not detected in the HBM_SA isolates. An operon of three genes, ebpA, ebpB, and ebpC, encoding the pilus subunits of PilB, was found in the genome of E. faecalis clinical isolates, but in the HBM isolates, ebpB was not detected. Interestingly, genes encoding exoenzymes and toxins like gelatinase (gelE) were commonly found in most isolates, but the hyaluronidase gene hylB was not found in the HBM isolates and was detected in 16 clinical isolates and V583 (Fig. 5). The gene espfm, involved in the virulence and biofilm-forming capacity, and fsrB, involved in the quorum-sensing system of E. faecalis, were not detected in the HBM isolates (Fig. 5).

Prophages, and mobile genetic elements analysis

In the analysis of putative prophage elements using PHASTER, intact phiFL1A was commonly found in HBM_SA E. faecalis isolates but was not present in intact form in either of the reference strains. The intact phage EFC_1 and phiFL4A were found in Efa_HMD11_46M. The phiFL3A and phiFL4A were detected in Efa_HMD12_49M. The V538 strain was carrying phiFL4A and phBC6A52. The probiotic Symbioflor 1 genome carried four intact phages of PBL1c, phiFL2A, phiEf11, and phBC6A52. The mobile genetic elements of insertion sequence ISLgar5, ISSsu5, Tn6009, and composite transposon cn_11752_ISLgar5 were commonly found in the three HBM_SA isolates, but those were not found in the Symbioflor 1 strain, and ISLgar5 was detected in the V583 strains along with five other MGEs (Table S7).

Plasmids and associated antibiotic resistance genes analysis

A diversity of plasmids was found in the analyzed 61 genomes of E. faecalis isolates, classified into 15 different replicons (Fig. 6). The HBM_SA isolates commonly carried the replicons repUS43 and rep9b. The rep9a replicon was not detected in the Efa_HMD9_11M isolate. No plasmids were detected in the HBM isolates from China and Symbioflor 1 (Fig. 6). The repUS43 replicon was the most common replicon type identified, also prevalent in animal milk (four out of eight, 50%) and clinical isolates (36 out of 43, 83.7%), followed by rep9b detected in 17 out of 43 clinical isolates (Fig. 6). The tetM gene from HBM_SA isolates were found on the repUS43 replicon carrying plasmid, which carried transposon Tn6009. The ermB gene was found on a contig carrying the IS1380 insertion sequence in the Efa_HMD11_46M and Efa_HMD12_49M isolates.