A previous prevalence study conducted in French slaughterhouses led to the identification of 17 farms that provided cattle that were positive for EHEC O157:H7 at slaughter (14). The objectives of the present study were to identify if EHEC O157:H7 shedding has persisted in these farms, and to explore factors that might play a role in the maintenance of O157:H7. Farmers were asked to participate on a voluntary basis and were assured of confidentiality. Two had ceased with cattle production and two did not want to participate. The study was, thus, conducted in 13 farms, named from A to M (Table 1). A first visit (spring/summer 2013) was conducted in order to identify the EHEC O157:H7 herd status. Twenty-one fecal samples were collected from young bulls and dairy cows in each farm (except for farms F and K, where only 20 and 14 samples could be collected, respectively). In the bovine species, young bulls and cows are, respectively, defined by the European legislation [Regulation (EU) No 1308/2013] as uncastrated males aged from 12 months to less to 24 months, and as females that have calved. One environmental sample per farm was also collected. Two farms were confirmed to be positive for EHEC O157:H7 and were consequently subjected to other herd visits. During 2013/2014, farms G and M were subjected to three and four supplementary herd visits, respectively (Table 1). Farms G and M housed about 440 and 160 animals (i.e., cows and young bulls), respectively.

At each visit, fecal and environmental samples were collected. Feces were taken just after defecation. Environmental sampling was performed by attaching moistened gauze to the outside of disposable plastic boots (overshoe) and by walking 40 steps around the areas where the animals were kept, avoiding fecal pats. All samples were kept chilled and sent to the laboratory by overnight courier for analysis.

Upon arrival, each fecal sample (10 g) was 10-fold diluted (w/v) in 90 ml of modified tryptone soya broth (Oxoid, Dardilly, France) supplemented with novobiocin (Oxoid, Dardilly, France) at 16 mg L⁻¹ in Stomacher bags. Environmental samples (two pieces of gauze) were placed in Stomacher bags containing 225 ml of modified tryptone soya broth supplemented with novobiocin at 16 mg L⁻¹. Stomacher bags were mixed and incubated overnight at 37°C. Bacterial DNA was extracted from 1 ml of each enriched broth using a lysis tube (Pall GeneDisc Technologies, Bruz, France), as already described (15).

DNA extracts were subjected to a sequential real-time PCR-based approach for the detection of O157:H7 EHEC-associated genetic markers (stx1, stx2, eae-γ1, and rfbE_O157) as previously described (14). Isolations were performed for samples that tested positive by PCR for the simultaneous presence of these markers. Isolation procedure consisted of immunomagnetic separation (IMS) assays using Dynabeads anti-E.coli O157 (Invitrogen, Cergy Pontoise, France), as recommended by the manufacturer. Ten microliters of immunoconcentrated bacteria were plated onto cefixime–tellurite–sorbitol–MacConkey (CT-SMAC) agar (Oxoid, Dardilly, France) and ChromoID O157:H7 agar (bioMérieux, Marcy l'Etoile, France). For each sample, five plates of each agar were plated, and all media were incubated overnight at 37°C. Up to 10 suspect colonies were tested by slide agglutination with serogroup-specific antisera (Statens Serum Institut, Copenhagen, Danemark). The serogroup O157 was confirmed by real-time PCR and the presence of stx1, stx2, and eae-γ1 was screened as described above. The presence of the fliC_H7 alleles was also investigated as previously described (16). Isolates were also confirmed as E. coli using an API 20E test (bioMérieux, Marcy l'Etoile, France). Based on PCR results, E. coli isolates positive for stx, eae-γ1, rfbE_O157, and fliC_H7 genes were classified as Shiga toxin-producing EHEC O157:H7. E. coli isolates positive for eae-γ1, rfbE_O157, and fliC_H7, and negative for stx genes were classified as EPEC O157:H7.

Enumeration of EHEC O157:H7 was performed in a liquid medium by using a most probable number method. One gram of feces, stored at 4°C, was 10-fold serial diluted in tryptone salt broth (Oxoid, Dardilly, France). One milliliter of each dilution was added, in duplicate, to 9 ml of modified tryptone soya broth supplemented with novobiocin at 16 mg L⁻¹. Broths were incubated overnight at 37°C. In order to isolate EHEC O157:H7 from positive (turbid) broths, IMS-based isolations were performed. Twenty-five microliters of immunoconcentrated bacteria were plated onto CT-SMAC agar and ChromoID O157:H7 agar. For each positive broth, five plates of each agar were plated, and all media were incubated overnight at 37°C. Suspect colonies were tested by PCR for the presence of stx, eae-γ1, and rfbE_O157 genes. Enumeration was finally calculated from the number of EHEC O157:H7-positive duplicates or each dilution according to McCrady's tables, and expressed as the most probable number per gram of feces (MPN/g) (17). In order to perform generic E. coli counts, 10 g of feces was homogenized in 90 ml of saline solution, and 10-fold serial dilutions were prepared. Decimal dilutions were plated onto Petrifilm^TM Select E. coli (Grosseron, Saint Herblain, France). Enumerations were performed after incubation at 42°C during 24 h.

Subtyping of stx₁ and stx₂ genes was performed as previously described (18). The presence of additional EHEC-virulence markers (enterohemolysin ehxA gene and OI-122 associated genes, namely, pagC, sen, efa1, and efa2 genes) was screened by PCR as previously described (19, 20). The presence of typical EPEC markers, i.e., bfpA and EPEC adherence factor (EAF) genes, was also tested by PCR (21, 22). E. coli isolates positive for the eae gene, and negative for the bfpA gene and EAF plasmid were classified as atypical enteropathogenic E. coli (aEPEC) (1). Finally, EHEC and EPEC strains were typed using the Standard PulseNet PFGE protocol for E. coli O157, as described previously (14).

Samples were processed for DNA extraction and purification using the QIAamp DNA Stool Mini kit (Qiagen), following the manufacturer recommendations. The V4 region of the 16S rDNA gene was amplified and sequenced at the Génome Québec Innovation Center of McGill University (Montréal, Canada). The FastStart High Fidelity PCR System kit (Roche Diagnostics), together with primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) were used to amplify the V4 region of the 16S rDNA gene (23). Sequencing was done using an Illumina MiSeq system (Illumina, USA). Resulting FASTQ files can be found in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA599584.

The 16S Illumina data have been analyzed by following the ≪Bioconductor Workflow for microbiome data analyses≫ (24) based on the dada2 algorithm developed to improve the detection of Amplicon Sequence Variants (ASVs) (25). The dada2 package has been used on R 3.5.1 (26) on a 32-core with 32-Gb RAM desktop computer for the pipeline's step of filtering, trimming, dereplication, to infer the sample composition, and to remove chimera. To improve the detection of rare variants from 16S amplicons and maintain a linear memory requirement, the dada2 procedure using the pseudo-pool algorithm of samples was used. The multiple-sequence alignment was performed using the DECIPHER R package (27). The Silva nr v.132 database (28) and a naïve Bayesian classifier (29) was used to assign taxonomy from kingdom to species.

The non-bacteria kingdom and singleton sequences (ASV, which only appears once) were removed before analyses. All the diversity analyses were performed with the Phyloseq R package (30). Raw ASV abundances were used for the alpha-diversity analyses. Beta-diversity analyses were performed using transformed abundance table with the DESeq2's variance stabilizing transformation (31, 32). Non-metric multidimensional scaling (NMDS) was used for the multivariate analysis of the community using the “vegan” package in R. Permutational multivariate analysis of variance (33) from the “vegan” R package were performed on ASV distance matrix to fit a regression model between the abundance distance matrix and variables. The microbiome differential abundance testing and log2foldChange estimate (34) were performed using the default multiple-inference correction of DESeq2 (Benjamini–Hochberg). Indicator species were identified using the “indicspecies” package in R (35, 36). We chose to use this method which focuses on identifying species that are both restricted to one group and with high fidelity (most samples in that group have the species), and thus, it offers a good complementary alternative method to the differential representation analysis using DESeq2. For this, a genus-level identity ASV table (similar to the one used for the DESeq2 analysis) was used as input. Each ASV ecological niche preference [shedders (S) or non-shedders (NS)]was identified using the Pearson's phi coefficient of association (corrected for unequal sample sizes) using the “indicspecies” package and 10,000 permutations. Only data from young bulls were used for all comparisons between S and NS to avoid confounding effects. All samples were considered as independent.

The presence of EHEC O157:H7 was screened in 13 farms (A to M), which provided cattle carrying EHEC O157:H7 at slaughterhouse 2 years before (14). The first visit allowed to identify two positive farms, namely, G and M, having three to seven young bulls shedding EHEC O157:H7 (Table 1). Two to three additional visits were performed in these two farms over a 1-year period, and each of them led to the detection of at least one shedder among young bulls. The percentage of young bulls shedding EHEC O157:H7 ranged from 14.3 to 31.8% and from 3.2 to 54.2% of the animals tested in farms G and M, respectively, depending on the time of visit. In each farm, an environmental sample led to the isolation of EHEC O157:H7 (during the second visit for farm G and the third visit for farm M). Finally, attempts to isolate EHEC O157:H7 also led to identification of young bulls harboring EPEC O157:H7 in farms G and M. None of the dairy cows tested from the farms visited in this study (including farms G and M) shed EHEC O157:H7 (Table 1).

For farm G, the 11 EHEC O157:H7 strains isolated during the first and second visits carried the stx_1a subtype, whereas the four strains isolated during the third visit carried the stx_2c subtype (Table 2). For farm M, the 25 EHEC O157:H7 isolates possessed simultaneously the stx_1a and stx_2c subtypes (Table 3). All EHEC O157:H7 strains were positive for eae-γ1, enterohemolysin ehxA, and OI-122 associated genes (pagC, sen, efa1, and efa2 genes). These genetic characteristics were shared by EPEC O157:H7 isolates (two from farm G and three from farm M), except that these EPEC lacked the stx genes (Tables 2, 3). All EPEC were negative for bfp and EPEC adherence factor (EAF) plasmid, thereby justifying their classification as atypical EPEC (aEPEC).

EHEC O157:H7 counts were low (Tables 2, 3). These counts were below 10 MPN/g of feces in 60% of the young bulls. Only one young bull from farm M (second visit) showed an EHEC count above 10⁴ MPN/g of feces and could be considered as a supershedder (Table 3). By contrast, generic E. coli counts ranged from to 1 × 10⁵ to 4.9 × 10⁷ CFU/g of feces indicating that EHEC O157:H7 represented minor subpopulations of E. coli.

The 40 EHEC and 5 aEPEC O157:H7 isolated from farms G and M were subtyped by pulse-field gel electrophoresis (PFGE) (Tables 2, 3 and Supplementary Figure 1). Two additional EHEC O157:H7 strains isolated 2 years before in the same two farms were added to the dendrogram (see Supplementary Figure 1). With regard to farm M, 23 EHEC O157:H7 and 2 EPEC O157:H7 showed high genetic relatedness as most of them clustered together and displayed PFGE types (PT) (referred to as PT 1 to PT 12) with more than 82% similarity (Supplementary Figure 1). Four other strains, including the EHEC O157:H7 strain isolated in 2010, did not belong to this cluster. Among the cluster, several strains shared the same PT, including strains isolated from different animals indicating transmission between these animals or exposure to a common source. Some strains with an identical PT were also isolated from different herd visits illustrating their persistence over time in farm M (Table 3). Strain 1bisP5-1 isolated during the second visit from a supershedder young bull (i.e., 1.3 × 10⁴ MPN/g) displayed the same PT (PT 8) as strain 1M9-1 shed during the first visit at a lower level (i.e., 60 MPN/g). Strain 2MP4-1 isolated from an environmental sample displayed the same PT (PT 1) as strains isolated from seven young bulls. With regard to farm G, a higher diversity of PFGE profiles was observed (Supplementary Figure 1). Of the 18 O157:H7 strains isolated from this farm, only two strains shared the same PT (PT 13), and these were isolated during the same visit (Table 2).

Fecal samples from 52 animals corresponding to 9 dairy cows (C) and 43 young bulls (YB), including 33 beef bulls and 10 dairy bulls, from farm M were analyzed by 16S rDNA amplicon MiSeq sequencing. These included the samples of 24 young bulls shedding EHEC O157:H7 (S) (Table 3) and the samples of 28 non-shedders (NS) consisting of 16 young beef bulls, 9 dairy cows, and 3 young bulls shedding EPEC O157:H7. The 16S rDNA amplicon analysis of the fecal samples generated 4,757,143 high-quality sequences, with an average of 179,514 sequences per sample after quality filtering. The overall number of ASVs detected by analysis was 7,264, based on a 97% nucleotide sequence identity between reads and after removing Archaea (37) and non-attributed (5) sequences. Rarefaction curves showed that the sampling was sufficient for all the 52 samples (not shown).

To compare the microbiota of S vs. NS, only the data from young bulls (43 samples) were used, in order to avoid confounding effects. Alpha diversity indices are shown in Supplementary Figure 2. By repeated measure of analysis of variance (ANOVA), it was apparent that there were no significant differences in alpha diversity between S and NS (43 samples) or between YB and C (52 samples) (p > 0.1). We analyzed the relative abundance of bacteria present in the feces of the animals. Animals were grouped by S vs. NS (Figure 1), and by YB vs. C (Figure 2). At the phylum, family, and genus level, there were no obvious differences between groups. The major phyla found in all groups were Bacteroidetes, Firmicutes, and Proteobacteria (Figures 1, 2). While the relative abundance of Bacteroidetes varied from 35 to 57% depending on the animals, with a mean at 45.9%, that of Firmicutes was more variable, from 14.5 to 59.2%. Proteobacteria relative abundance was also highly variable, from 0.1 to 31.8% depending on the animals. Five animals with a very high proportion of Proteobacteria (1-M17, 1b-P21, 1b-P1, 1-M7, and 1b-P14) were also those with the lowest proportion of Firmicutes and were all young bulls (Figures 1A, 2A). Three of these animals were EHEC O157:H7 shedders (1-M17, 1b-P1, and 1-M7). In samples from animals 1b-P1, 1b-P14, and 1b-P21, the relative abundance of Tenericutes was also higher than for the other animals (from 6 to 10%). Analysis at the genus level (Figures 1B, 2B) showed that in animals with a high proportion of Proteobacteria, this phylum was mostly composed of Pseudomonas, Acinetobacter, Thauera, or non-attributed genera. In these samples, the Enterobacteriaceae family was not represented (not shown).

Overall differences in community structure were assessed using NMDS (Supplementary Figure 3). No clear separation of samples was observed for bacteria associated with S and NS or C and YB. We further examined the differences in bacterial community structure between S and NS samples (43 YB samples) using Permanova (on Bray Curtis distance), which showed that samples did not cluster by shedding status (p > 0.05). However, when adding the sampling time (farm visit) as variable, the community structure appeared different between S and NS feces (p = 0.0496). In addition, the sampling time had clearly a significant effect on the microbiota beta diversity (p = 0.0001) but not on alpha diversity (not shown).

The analysis of ASVs differential abundance between YB and C fecal samples (using DESeq2) identified 24 significantly differential ASVs (Supplementary Figure 4), but the DESeq2 analysis did not find differentially abundant ASVs between S and NS. However, using the indicator analysis “indicspecies” package, we identified six ASVs as significant indicators of the S group, and four ASVs as indicators of the NS group (Table 4). Six of these indicator species belong to the Firmicutes phylum. The bacterial indicators of S belong to family XIII UCG-001, Slackia, and Campylobacter genera, and Ruminococcaceae NK4A21A, Lachnospiraceae-UGC-010, and Lachnospiraceae-GCA-900066575 groups. The NS group indicator ASVs were affiliated to Pirellulaceae-1088-a5 gut group, Anaerovibrio, Victivallis, and Sellimonas genera.

EHEC O157:H7 were isolated from young bulls, but never from cows. This is in agreement with other studies showing that younger animals are associated with an increased risk of shedding EHEC O157:H7 (38, 39). Moreover, the fact that EHEC O157:H7 strains were isolated from only a fraction of young bulls tested at each visit suggests that some positive animals might have been missed due to certain limitations of our studies related to the sporadic pattern and short duration of EHEC O157:H7 fecal shedding by animals. Indeed, a growing number of longitudinal studies have led to the hypothesis that E. coli O157 prevalence may be heterogeneous (40–42). Finally, two EHEC O157:H7 strains were isolated from the farm environment, including one strain in farm M that displayed the same PFGE type as several isolates retrieved from fecal samples. As already proposed by others, our results suggest that the analysis of environmental samples might be useful to determine the herd status with respect to EHEC shedding by animals (43). All EHEC O157:H7 strains isolated were positive for stx_1a and/or stx_2c. Additionally, they were all positive for the eae-γ1 gene, the enterohemolysin ehxA gene, and OI-122-associated genes. These EHEC virulence markers have already been shown to be associated with O157:H7 strains that cause severe outbreaks (19, 44).

As EHEC O157:H7 supershedders have been identified to play a key role in the transmission of EHEC O157:H7 at the farm level (10), we performed counts of EHEC O157:H7 to identify such animals. For 60% of the shedders, these counts were very low (≤10 MPN/g of feces), and only one animal was a supershedder. This is in agreement with longitudinal studies that identified low numbers of supershedders in the herds examined (11, 45, 46). Interestingly, the strain from the supershedder shared the same PT as one strain shed in an earlier visit at a low concentration, suggesting that strain-specific characteristics are not sufficient for causing supershedding by cattle, as observed previously (47). Finally, the comparison of the EHEC O157:H7 counts to the generic E. coli counts from the same animals showed that EHEC O157:H7 represented a minority of the E. coli population. This observation confirmed the interest of researches focused on mechanisms of competition among E. coli strains in order to identify factors controlling the proliferation of EHEC in the digestive tract, such as colicins, for example (48).

Further characterization of EHEC O157:H7 strains revealed different patterns of shedding dynamics in the two farms investigated. In farm M, the 25 EHEC O157:H7 strains isolated during four herd visits showed identical virulence profiles. Notably, they were all positive for the stx_1a and stx_2c subtypes. PFGE analysis revealed that the majority of these strains were highly related and clustered together, including strains with either identical PT or showing 82% of similarity. Although up to five band losses have been reported following subculturing in the laboratory (49), only a few studies have considered the stability of PFGE profiles following host passage or exposure to environmental stimuli. Yoshii N. et al. reported a one to six band loss following PFGE analysis of E. coli O157:H7 after calf inoculation (50). In the light of these elements, our results showed that the E. coli O157:H7 isolates were highly related, suggesting they might originate from the same strain. Nevertheless, two EHEC O157:H7 strains from farm M belonged to PT quite different from that of the clustered isolates suggesting the circulation, although at a lower extent, of additional EHEC O157:H7 strains within this farm. The situation in farm G differed from that in farm M on several aspects. First, strains carrying different stx subtypes were isolated, including stx_1a-positive strains during the two first visits followed by stx_2c-positive strains during the last visit. Second, these strains showed a high diversity by PFGE analysis. Altogether, these results indicated that a higher number of EHEC O157:H7 strains that are unrelated circulated within this farm. The persistence of an O157:H7-positive herd status was, therefore, related to distinct events as observed by others (40, 42, 51, 52). On the one hand, it may result from continuous exposure of cattle to various strains as observed in farm G and, to a lower extent, in farm M. This was also the case for both farms G and M when sampled 2 years apart as the two strains isolated previously at slaughter were not the same as those isolated here. On the other hand, persistence of a positive herd status may be caused by the circulation of a predominant strain (and its derivatives resulting from genetic variations associated with minor PFGE variations) as observed in farm M. One strain (2M9-1) belonging to the cluster of related strains from farm M, was subjected to further studies in order to explore factors involved in the persistence (53). Our results showed that this strain (renamed MC2) had a large number of colonization factors known to bind to intestinal epithelium and to biotic or abiotic surfaces. We also showed that MC2 had the capacity to produce biofilms and to activate stress fitness genes in bovine feces. Taken together, these results could explain the persistence of this strain and its derivatives in the farm environment, acting as a source of contamination for cattle.

Although different factors have been associated with EHEC shedding (37), several studies have suggested that microbiota composition plays a critical role in the establishment and/or ecology of E. coli O157:H7 within the intestinal tract of cattle. Some recent results supported this hypothesis (12, 13, 54–57), but the results differ between studies and the genera associated with the fecal microbiota of EHEC shedders and non-shedders varied from one study to another (58). In the present study, when analyzing the young bulls' fecal microbiota, we did not find any significant difference in alpha diversity between EHEC O157:H7 shedders vs. non-shedders. Furthermore, there was no significant difference in the overall community structure of these two groups, but the sampling time had a clear effect on the community structure. The indicator species analysis identified 10 indicators of the S or NS groups. Six ASVs were indicators of the S group, including four sequences affiliated to the Firmicutes phylum (two Lachnospiracea, one Ruminococcacea, and one Family XIII). Four ASVs were indicators of the NS group, including two Firmicutes sequences (Lachnospiracea and Veillonellacea), one Planctomycetes, and one Lentisphaerae ASVs. Identifying genera associated with non-shedding status could help in the design of potential new direct-fed microbial candidates to limit cattle STEC shedding (58).

In a previous study on 240 calves, the Shannon diversity index was lower in fecal samples from animals colonized by EHEC (54). In another study on 40 dairy cattle, the average microbiota richness of cows shedding EHEC was also lower than that of non-shedders, but statistical analysis did not find differentially abundant taxa between them (55). In a similar study on feedlot steers, Xu and colleagues found that supershedders exhibited higher bacterial richness and diversity than non-shedders, and 72 OTUs were found as differentially abundant between the two groups (13). In another work on feedlot steers, the overall bacterial community structure did not differ by E. coli O157:H7 shedding status, but 4 OTUs (Ruminococcus, Selenomonas, Campylobacter, and Streptococcus genera) were more abundant in supershedders (57). Campylobacter was also found to have a stronger preference for the gastrointestinal tract (GIT) of EHEC shedders in our work. A recent study on several farms in the USA showed that the fecal microbiota of farms with a high-STEC prevalence had greater richness compared with those of farms with a low-STEC prevalence, but no alpha diversity difference was observed when comparing EHEC shedders (stx and eae positive) and controls (59). These authors found 30 genera to be differentially abundant between STEC carriers and non-carriers, which were different from the ones identified in the present work. Wang and colleagues compared the mucosal microbiota of the rectoanal junction of supershedders and non-shedders, and they identified OTUs unique to SS or NS (56), but none of them was similar to the ones identified in the present study. As illustrated by these examples, correlating EHEC-shedding status with fecal or intestinal microbiota composition, richness, or specific taxa abundance is not straightforward. Indeed, the inconsistency in the results may reflect the variability in the inter-animal microbiota composition, which may be due to several factors including the animal age (54, 60) or diet (12, 59), but also to the farm environment or the methodology used in the different analyses. It is also known that the season has a strong impact on the shedding status of cattle, and a seasonal effect has also been shown on the microbiota composition of the GIT (feces) (9, 61, 62). Further studies should therefore compare the microbiota of a high number of animals of the same age and diet, and similar rearing mode. Finally, our study also compared the fecal microbiota of young bulls vs. lactating cows, and found 24 genera differentially abundant between these two groups. Although several studies monitored the evolution of the fecal microbiota of calves from birth to weaning (60, 63–65), only a few of them compared young bulls and lactating cows. Dill-Mcfarland et al. showed that the bacterial fecal community of 1- and 2-year-old animals, clustered clearly differently from that of younger ones (2–8 months) (60). As in our study, the genus Faecalibacterium was found much more abundant in feces of young animals.