Ruminants are the primary reservoir of Shiga toxin-producing Escherichia coli (STEC) O157 (1), a significant foodborne pathogen that has the potential to cause disease in humans (2). STEC infections can have severe clinical outcomes, including death in young children (3). An estimated 2,801,000 global STEC O157 acute human illnesses occur annually (4), and in the US, 68% of STEC O157 infections are considered to be foodborne (5). Ground beef or fertilizer derived from cattle effluent is implicated as the vehicle in many foodborne outbreaks (6, 7) and sporadic illnesses (8), prompting large scale control efforts to minimize STEC contamination along the production chain in the USA (9) and Australia (10, 11). Environmental and animal contact are also considered significant exposure pathways for sporadic STEC O157 human infections (12). Shedding of STEC O157 by beef cattle has been described both in feedlot and pasture based settings (13), with variation in prevalence estimates between herds and between regions (10, 14, 15). Super shedding has been implicated as a significant driver of pathogen transmission (16, 17); however, its role in the epidemiology and transmission dynamics of E. coli O157 is still being explored (18).

The molecular epidemiology of STEC O157 in cattle is complex, as multiple subtypes are often circulating simultaneously on individual farms (19). Further, identical subtypes can be found on farms separated by hundreds of kilometers without any evidence of an epidemiological link between farms (9). Traditionally, molecular typing of STEC O157 has primarily been based on pulse field gel electrophoresis, multilocus variable number of tandem repeat analysis (MLVA), or phage typing (20–22). However, these comparative typing methods have shortcomings that may be exacerbated when supported by insufficient epidemiological study designs (23).

With the increased application of whole-genome sequencing technology in outbreak investigations (24), progress is being made in better capturing the evolutionary and transmission dynamics of STEC O157. For example, different lineages have been associated with human disease and asymptomatic carriage in the bovine host and genomic signatures have been identified that are associated with altered pathogenic potential (25). The identification of phylogenetically informative single nucleotide polymorphisms (SNPs) has further supported epidemiological investigations of STEC O157 transmission, offering discriminatory power unmatched by macrogenomic techniques (26, 27).

A better understanding of the temporal variation of STEC O157 shedding can help target control efforts. Strain-specific characteristics, such as pathogenicity and adaptation mechanisms, can influence shedding (28) and ultimately control strategies (16). Thus, the application of molecular tools in a spatiotemporal context can greatly enhance the understanding of infectious disease dynamics (29–31). For example, persistence of a single strain type versus shedding of multiple types over time may be an indicator for different epidemiological patterns, exposures, and transmission routes.

To address the challenges described, the shedding dynamics and molecular diversity of STEC O157 in an Australian beef herd was recently investigated (32, 33). Variation in shedding was observed both within and between cattle over time, with a total of 33 MLVA profiles identified during the course of the 9-month longitudinal study. Clonal replacement was found to occur over time as well as shedding of identical strains between cows. The objective of the work presented here was to integrate epidemiological data with molecular typing results from said studies to further analyze the spatial and temporal trends in genotypes shed and to better understand the previously identified patterns of shedding synchronization.

Previously reported molecular typing data from Jones et al. (32), including MLVA types and WGS results, were used for this study. To analyze the spatial and temporal trends of the MLVA genotypes, we used a new approach to define relationships of MLVA types and GCs. We first computed the relationship of MLVA types that had WGS data and determined a VNTR difference of two loci as the cutoff for GCs. Using this cutoff, we used bootstrapped MSTs to estimate the relationship between 33 different MLVA types. The MST with the highest bootstrap percent (43 ± 5.8) is composed of 33 nodes, each representing a unique MLVA type (Figure 1). A total of three GCs were identified based on a maximum distance of two loci to be considered the same cluster. A consensus network MST was also computed (not shown), which supported the GC definitions obtained from the bootstrapped MST. Nineteen, seven, and three MLVA types were found in GC1, GC2, and GC3, respectively. Four MLVA types (1, 9, 13, and 26) did not belong to a GC, as they differed by at least three locus variants from any other type.

MT14, MT2, and MT3 were each the first shed MLVA type of GC1, GC2, and GC3, respectively, with MT2 and MT3 each found at a central node in the MST in relation to the remaining members of the cluster. MT14 was the first GC1 MT to be isolated (on 5th November 2012) and was predominantly found in the beginning of the study period, whereas MT20 was the most commonly isolated MLVA type and the only GC1 type found after movement to the final paddock on February 28, 21 weeks after the start of the study.

All cows shed isolates from at least one of the GCs during the 9-month study period, with 9% (n = 2) shedding isolates from a single GC, 61% (n = 14) shedding isolates from two different GCs, and 30% (n = 7) shedding isolates from three different GCs (Figure 2). Five cows shed isolates from a previously isolated GC after shedding a different GC (e.g., shed Cluster 1, then Cluster 2, then Cluster 1 again). Further, in six instances, an isolate from a different GC was recovered in consecutive sampling points. Figure 2 illustrates the number of GCs isolated from each cow in each paddock. Cows spent the most amount of time in Paddock 6, which is also the only paddock where all three GCs were found. A significant difference was found in the distribution of GCs between paddocks (P < 0.001), but not within cows (P = 0.142), based on Fisher’s exact tests.

The number of isolates identified from each GC at each sampling point is shown in Figure 3. A trend can be seen, with a higher number of GC2 isolates found in the early part of the study period, followed by an apparently bimodal distribution of GC1 shedding in the middle, and increasing numbers of cows shedding GC3 isolates at the end of the study period. GC1 was isolated from cows in all paddocks, with the exception of Paddock 1. GC2 was the most frequently isolated GC in Paddock 2 and 3 but was not detected after cattle were moved to Paddock 4. GC3 was only isolated from cattle in Paddock 6.

The analysis presented here expands upon previous results from a 9-month longitudinal study of 23 co-grazed beef cattle in Australia (32, 33), supporting a comprehensive investigation of the shedding dynamics of STEC O157 genotypes. Three STEC O157 GCs were detected, as well as four unique MLVA types that did not belong to any of the GCs identified. Most cows shed MLVA types from more than one cluster. The investigation revealed a stronger association between GC and paddock number than between GC and individual cow ID, which supports the hypothesis that shedding of E. coli strain types is more strongly influenced by environmental exposure, than individual cow attributes (39).

Interestingly, the number of cows shedding GC1 isolates appeared to follow a bimodal epidemic curve, with the first peak occurring in mid-January and the second in early March. Toward the end of the study period, GC3 showed a similar trend with an increasing number of cows shedding that cluster until the study terminated in June. This might indicate that cows developed an immune response to the different GCs, resulting in a peak shedding period, followed by a period of decreased shedding. This is supported by the fact that despite not showing clinical signs of illness, cattle do mount an immune response to STEC O157 (40) and decreased shedding has previously been observed in cows vaccinated against STEC O157 (41, 42), with evidence of strain dependence (43).

Co-colonization of individual cows with multiple genotypes has been documented previously (44) and is likely to have occurred here. Shedding of a previously isolated GC after shedding of different GCs suggests possible co-colonization, which was observed in five cows over the course of the study period. However, coinfection with multiple genotypes cannot be confirmed, as only one isolate was saved from each cow per sampling point.

Interestingly, a change from Paddock 3 to 4 and a change from Paddock 5 to 6 marked noticeable shifts in the STEC O157 population structure identified in this analysis. In the first instance, GC2 seemed to disappear from cattle upon changing to a new paddock. Conversely, GC3 was not detected until cattle moved to Paddock 6. This again suggests that environmental forces may play a substantial role in the colonization of cattle, such as movement to an environment where new STEC O157 types are present or the introduction of a new type (e.g., by wildlife vectors). For example, starlings have been implicated in the transmission of STEC O157 between dairy farms (45), and free-ranging deer sharing rangeland with cattle have been found to shed STEC O157 (46). Isolation and genetic analysis of STEC O157 from the environment and wildlife populations in this study would help elucidate the sources and transmission dynamics between different hosts and environments.

The use of the distance method to generate an MST was found to be the most appropriate for the data available, as single-repeat changes were previously observed to occur most frequently (47); however, in that same study, 25% of the mutational events included multiple-repeat changes. Such large changes would likely not have affected the assignment of GCs in this study as, by definition, GC was inclusive of single-locus variants. Here, the definition of GC was supported by WGS data of 18 isolates representing eight MLVA types with genome lineages 1, 2, and 3 corresponding to the GC1, GC2, and GC3, respectively. Consequently, genotypes were clustered into groups differing by no more than two loci, which is the same threshold previously suggested by Mellor et al. (48). Previous analysis of the same MLVA data by Jones et al. (32) used a threshold of a single VNTR difference to define a clonal complex; however, this conflicts with the WGS data that showed that isolates with as many as five VNTR differences were closely related. Thus, WGS provides necessary resolution to infer relatedness of strains in this context and should be incorporated into the interpretation of other typing methods when such data are available.

While the results generated by this study have provided additional insight into STEC O157 transmission dynamics on the study site, there are some shortcomings. In this study, paddock changes were convenience based and driven by herd management needs, such as pasture quality and nutritional requirements, and thus cows did not spend equal amounts of time in each location. As such, additional analysis of STEC O157 shedding patterns in other herds and environments are needed to substantiate our findings. Furthermore, as with any epidemiological study, unaccounted for selection bias could have occurred. The strong sampling design did, however, offer insight into the spatiotemporal nature of STEC O157 shedding of cattle in this pasture based environment.

This study used WGS data to assist defining MLVA GCs. In the absence of WGS information, the true evolutionary relationship of isolates within and between MLVA types cannot be confirmed. Care should be taken when making inferences about the population dynamics based on MLVA data, and results should be considered to be of preliminary nature. This is highlighted by the fact that isolates with up to five VNTR locus differences can still be nearly identical based on WGS. This is not entirely surprising, as Vogler et al. (47) determined that the mutation rate of VNTR loci in STEC O157 is high and variable between loci, confounding the analytical power of MLVA. Similarly, MLVA typing of Mycobacterium avium subspecies paratuberculosis was found to be unreliable to detect related and unrelated isolates, as it both over- and underestimated relatedness in comparison with WGS (49).

Importantly, results from the current analysis could inform potential future typing studies that were out of scope for the investigation presented here and could provide a rationale for the prioritization of additional isolates for WGS using existing samples. For example, building onto these initial studies, STEC O157 isolates (n = 295) from an intensive sampling study performed by the same research group (50) have been collected and could be used to investigate the substantial diurnal variation in shedding that was identified in 24 cows sampled one or two times per day for 2 weeks. Such high-resolution genetic analyses would provide further insight into STEC O157 daily shedding dynamics in Australian beef cattle if based on molecular typing methods compatible with the spatiotemporal range under investigation. We suggest that MLVA typing in isolation would not be appropriate in this context and that WGS would be needed to resolve daily trends in STEC O157 shedding.

In conclusion, this study has provided further insight into STEC O157 shedding dynamics in an Australian beef herd over time and identified environmental forces as a likely driver of pathogen transmission. The knowledge gained can be used as a stepping stone to design higher resolution studies to answer additional questions related to cow- and herd-level shedding synchronization of specific strains over time.

PM, CA, and JH drafted the manuscript and made a lead contribution to the conception, design, analysis, and interpretation of the work. GL, MJ, SO, and RL critically reviewed the manuscript and made contributions to field, laboratory, analysis, and interpretation of the work.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.