Disentangling a complex nationwide Salmonella Dublin outbreak associated with raw-milk cheese consumption, France, 2015 to 2016

On 18 January 2016, the French National Reference Centre for Salmonella reported to Santé publique France an excess of Salmonella enterica serotype Dublin (S. Dublin) infections. We investigated to identify the source of infection and implement control measures. Whole genome sequencing (WGS) and multilocus variable-number tandem repeat analysis (MLVA) were performed to identify microbiological clusters and links among cases, animal and food sources. Clusters were defined as isolates with less than 15 single nucleotide polymorphisms determined by WGS and/or with identical MLVA pattern. We compared different clusters of cases with other cases (case–case study) and controls recruited from a web-based cohort (case–control study) in terms of food consumption. We interviewed 63/83 (76%) cases; 2,914 controls completed a questionnaire. Both studies’ findings indicated that successive S. Dublin outbreaks from different sources had occurred between November 2015 and March 2016. In the case–control study, cases of distinct WGS clusters were more likely to have consumed Morbier (adjusted odds ratio (aOR): 14; 95% confidence interval (CI): 4.8–42) or Vacherin Mont d’Or (aOR: 27; 95% CI: 6.8–105), two bovine raw-milk cheeses. Based on these results, the Ministry of Agriculture launched a reinforced control plan for processing plants of raw-milk cheeses in the production region, to prevent future outbreaks.

On 18 January 2016, the French National Reference Centre for Salmonella reported to Santé publique France an excess of Salmonella enterica serotype Dublin (S. Dublin) infections. We investigated to identify the source of infection and implement control measures. Whole genome sequencing (WGS) and multilocus variable-number tandem repeat analysis (MLVA) were performed to identify microbiological clusters and links among cases, animal and food sources. Clusters were defined as isolates with less than 15 single nucleotide polymorphisms determined by WGS and/or with identical MLVA pattern. We compared different clusters of cases with other cases (case-case study) and controls recruited from a web-based cohort (case-control study) in terms of food consumption. We interviewed 63/83 (76%) cases; 2,914 controls completed a questionnaire. Both studies' findings indicated that successive S. Dublin outbreaks from different sources had occurred between November 2015 and March 2016. In the case-control study, cases of distinct WGS clusters were more likely to have consumed Morbier (adjusted odds ratio (aOR): 14; 95% confidence interval (CI): 4.  or Vacherin Mont d'Or (aOR: 27; 95% CI: 6.8-105), two bovine raw-milk cheeses. Based on these results, the Ministry of Agriculture launched a reinforced control plan for processing plants of raw-milk cheeses in the production region, to prevent future outbreaks.

Background
Nontyphoidal Salmonella is a main cause of bacterial food-borne infection in Europe [1,2]. The majority of human infections is caused by a limited number of Salmonella serotypes among the 2,600 described to date [3,4]. Salmonella enterica serotype Dublin (S. Dublin) is particularly invasive in humans and more often leads to severe disease and higher mortality rates compared with other serotypes [4][5][6][7]. S. Dublin is host-adapted to bovines and is frequently isolated from cattle, with raw milk or raw-milk cheeses as a typical vehicle for food-borne outbreaks [8,9]. In 2012, a major S. Dublin outbreak occurred in France, with 103 cases linked to Saint-Nectaire (bovine raw-milk cheese) consumption [10,11]. In 2015, 34 S. Dublin cases were reported linked to the consumption of Reblochon (bovine raw-milk cheese) (data not shown; Santé publique France).

In
France, the National Reference Center for Salmonella (NRC) and the French Agency for Food, Environmental and Occupational Health and Safety (ANSES) routinely collect and serotype human and non-human Salmonella isolates, respectively [12][13][14], using the Kauffmann-White-Le Minor scheme [3]. The S. Dublin isolates collected are frequently susceptible to all antibiotics and show an indistinguishable pulsed-field gel electrophoresis (PFGE) pattern. To better distinguish S. Dublin isolates, multilocus variable-number tandem repeat analysis (MLVA) has recently been used for surveillance and outbreak investigations [11,15]. Moreover, whole genome sequencing (WGS) of Salmonella has been shown to discriminate between closely related isolates of S. Dublin [16,17].

Outbreak detection
On 18 January 2016, the French NRC reported to Santé publique France (SpFrance, the French national public health agency) an excess of S. Dublin infections across the country, with 37 S. Dublin isolates identified between mid-November 2015 and mid-January 2016, compared with 10 S. Dublin isolates during the same period in the two previous years. An outbreak investigation team with experts from SpFrance, NRC, ANSES and the French Directorate General for Food (DGAL) launched extensive epidemiological, microbiological and food investigations to confirm the outbreak, identify the vehicle of transmission and propose appropriate control measures.

Methods
We carried out both epidemiological and microbiological investigations on subsets of S. Dublin cases and isolates, respectively.

Microbiological investigations
Salmonella Dublin isolates During the years 2015 and 2016, a total of 324 S. Dublin isolates were collected, 223 from clinical NRC isolates (108 in 2015 and 115 in 2016) and 101 from non-human ANSES isolates (62 and 39, in each year respectively). Of those, a total of 235, including 147 (83 in 2015 and 64 in 2016) clinical and 88 (62 and 26, respectively) non-human isolates, were extensively studied by WGS and/or MLVA. We also analysed 54 'historical' subtyped S. Dublin isolates collected between 1929 and 2014, 31 from humans and 22 from non-human specimens, as well as the human isolate which gave the name to the serotype in 1929 (number 65k) [18]. In total, we included 289 isolates in this study ( Figure 1).

Figure 1
Flowchart of human and non-human Salmonella Dublin isolates a subtyped and included in the study, France, 2015-2016 (n = 289) System (Roche, Basel, Switzerland), DNA was further processed for sequencing with Illumina systems (libraries using the Nextera XT DNA Library Prep kit and the sequencing with the NextSeq 500 system) generating 100 to 146 bp paired-end reads. Reads were trimmed and assembled as previously described [19]. Genomic data as multilocus sequence typing (MLST) type and resistance genes were detected from assembled sequences using web-tools (http://www.genomicepidemiology.org/). For each isolate, the paired-end reads were aligned against the S. Dublin str. 3246 reference genome (GenBank accession number: CM001151.1) using Bowtie2 with default parameters [20]. A coregenome multi-alignment of assembled genomes was also done using Harvest v1.0.1 f ParSNP function [21]. For each approach, the resulting single nucleotide polymorphisms (SNPs) were concatenated to generate a filtered multiple alignment that was used as input for the construction of a phylogenetic tree using Molecular Evolutionary Genetics Analysis (MEGA)6 [22] with a maximum-likelihood (ML) approach. The final trees were visualised in the interactive Tree Of Life [23]. All reads generated in this study have been deposited in project PRJEB28817.
Multilocus variable-number tandem repeat analysis ANSES used MLVA as described elsewhere [11], to analyse 241 isolates (including 148 human and 93 non-human), of which 110 corresponded to the outbreak period. The measured lengths for each fragment were obtained using an ABI3500 capillary electrophoresis system (Applied Biosystems, France). Data were imported into GeneMapper software (Applied Biosystems, France) where each fragment was identified according to colour and size. A normalisation of the results was done with the free access MLVA_ Normalizer software [24].

Epidemiological investigation
Case

Cluster definition
We defined a cluster as isolate sequences with < 15 SNP divergence obtained by core-genome comparison and/ or with identical MLVA pattern. Among these clusters, we also defined subclusters as isolate sequences having < 5 intra SNPs.

Study design
We compared cases belonging to a specific cluster/subcluster with other cases belonging to all other clusters/ subclusters in a case-case study and with controls in a case-control study. For the case-control study, controls were recruited from a cohort of individuals with children and adults registered on GrippeNet.fr (https:// www.grippenet.fr), an online population-based surveillance system for influenza-like illness [25,26]. During winter 2015/16, 6,515 participants reported online the presence or absence of basic symptoms on a weekly basis using a list of 19 predefined symptoms commonly or rarely related to influenza. We excluded controls reporting travel abroad during the Salmonella outbreak period and those with digestive symptoms.

Data collection
Epidemiologists from the regional offices of SpFrance interviewed the cases by telephone using a trawling questionnaire on clinical symptoms, medical history, detailed food consumption, contact with other persons experiencing diarrhoea, travel history and contact with animals. Loyalty card information was collected to trace-back supermarket purchases.

Statistical analyses
We calculated proportions, using the number of nonmissing values as denominators. For the case-case study, we calculated crude odds ratios (OR). For the case-control study, we calculated adjusted odds ratios (aOR) for age and sex using multivariable logistic regression. The initial regression models included age, sex and food items consumed by at least 50% of the cases. We performed this analysis for WGS clusters/ subclusters and MLVA clusters with at least 10 cases. We used STATA version 12.0 (Stata Corporation, Texas, United States) for this analysis.

Ethical considerations
The study was approved by the French Commission for Data Protection (Commission Nationale de l'Informatique et des Libertés). Interviewees or next of kin provided verbal consent. Only anonymised data were analysed and used for the purpose of the study.

Food production chain and animal trace-back investigations
DGAL conducted a trace-back investigation on potential contaminated products identified by the epidemiological investigations. Points of purchase, like supermarkets or cheese retailers, were reported by cases. Where possible, customer loyalty card numbers were used to identify the exact point of purchase and product batch numbers and specific production facilities. The traceback investigation was conducted for a retrospective period of up to a month before symptom onset of given cases, if date of onset was available, or if not, up to two months before the date of isolation of S. Dublin in the patient specimen. Suspected food products were tested for Salmonella if food samples were available. The number of cases peaked during week 53 (28 December-3 January) ( Figure 2). Week of onset

Description of controls
Among the 6,200 GrippeNet.fr participants (i.e. controls) who could be reached, 2,914 (47%) completed the questionnaire. Of those, 2,690 (92%) did this within 2 days; 1,916 (66%) were female; median age was 56 years (range: 2-90). The whole process, i.e. from the beginning of the recruitment of the controls to the end of data collection, took 12 days.
There were significant differences between cases and controls in terms of proportion of female (53% in cases vs 66% in controls; p value = 0.01) and age (58% of cases vs 29% of controls were > 65 years-old; p value = 0.00).

Microbiological findings
Among 241 MLVA performed for S. Dublin isolates, 49 different MLVA patterns were obtained and among these, seven gathered 10 or more isolates ( Table 2).
The WGS analysis of all the 289 S. Dublin isolates, including historical ones, indicated a unique MLST, ST10 that had previously been related to serotype Dublin in both the literature and public databases (https://enterobase.warwick.ac.uk/) [27]. As the mapping against the S. Dublin str. 3246 reference genome (GenBank accession number: CM001151.1) showed high divergence with 19,213 SNPs, the core-genome multialignment of assembled genomes using ParSNP function was preferred. A generated matrix file revealed that maximum divergence for all isolates sequenced during the outbreak including the reference strain was 791 SNPs suggesting a relative homogeneous population of S. Dublin isolates that were circulating in France. We identified 28 different clusters with < 15 SNPs, including five clusters with > 10 isolates (clusters A, B, C, F and K). Three of those, A, B and C, accounted for the majority of human (70%, 125/179) and non-human isolates (46%, 51/110) (Figure 3).
A total of 18 isolates (17 human and 1 non-human) with date of isolation between October 2015 and March 2016 were identified as belonging to the WGS cluster A and were mainly associated to the MLVA cluster 20-8-10-7-5-4 (Table 2, Figure 3). In terms of food items, only one WGS cluster A isolate from October 2015 was found in raw milk.
The WGS cluster B comprised a total of 35 isolates (20 human and 15 non-human). Two main subclusters were identified; one B3 (mainly associated to the MLVA cluster 17-8-10-7-5-4) was found in relation to cattle, milk and raw-milk cheeses (Reblochon and Morbier) with date of isolation mainly in January 2016 and the other B4 (mainly associated to MLVA cluster 18-8-10-7-5-4) was found in cases with date of isolation mainly in December 2015 and January 2016.
The WGS cluster C was the most prevalent with 123 isolates (88 human and 35 non-human). It was subdivided into nine subclusters (< 5 SNPs), C3 and C other (grouping eight smaller subclusters including C8). All the 35 isolates belonging to the C3 subcluster presented a very limited intra SNP difference (< 2) and a sufficient divergence (> 5 SNPs) to the other eight C subclusters (C other ), indicating high level of genetic relationship due to a putative common source of contamination. The WGS subcluster C3 was mainly associated with MLVA cluster 18-8-10-7-5-4 and was found in clinical isolates from patients between January and April 2016. This molecular signature was also found in one isolate from Morbier cheese in February 2016 during this investigation. It was found as well in several isolates in Morbier cheeses tested in 2015 through company internal microbiological monitoring system. The WGS subcluster C8 was mainly associated with the MLVA cluster 19-8-10-7-5-3 and harboured isolates from patients at the beginning of 2015. In that period, another S. Dublin outbreak had been investigated between February and April 2015 in seven French regions, and was found to be possibly associated with consumption of Reblochon cheese (data not shown; Santé publique France).
Among the non-human isolates, we revealed several other WGS cluster groups (in particular F and K) but no or few linked human isolates (data not shown). The review of the historical human and non-human      frequently coincided with WGS cluster C3) appeared as consuming more often Vacherin Mont d'Or cheese (OR: 5.1; 95% CI: 0.9-35). No associations were found with other raw-milk cheeses, nor with any other food items.

Food trace-back investigations
We collected 39 (62%) loyalty card numbers from 63 cases. Based on the available information, trace-back investigations were conducted among 10 supermarket brands. Twelve cheese producers were identified as potential origin of the cheeses consumed by the cases. The trace-back investigations linked one Morbier producer and three different Vacherin Mont d'Or producers, to 11, five, four and three cases, respectively. All those producers were located in the same region, i.e. Bourgogne-Franche-Comté (Eastern part of France).

Discussion
We reported one of the largest S. Dublin outbreaks in France in the past few years. Two different bovine rawmilk cheeses, Morbier and Vacherin Mont d'Or, were the most likely vehicles of transmission for this food-borne outbreak. For the present outbreak investigation we used two different typing methods on a large panel of strains (both historical and obtained during the 2015-16 outbreak period, as well as from human and nonhuman origins). The first method was MLVA, which had already been used in previous investigations of other outbreaks in France [11]. The second, WGS, was used for the first time in the current S. Dublin investigation and demonstrated increased capacity to discriminate clusters. We also used two different epidemiological methods: a case-case study that allowed a rapid analysis and identification of suspected sources and a case-control study that was more statistically powerful to confirm the suspected associations.
In this investigation, MLVA was deemed sufficient to identify a link between human cases, food and animal sources. However, MLVA could not distinguish some of the clusters identified by WGS. Our results suggested that at least two outbreaks of S. Dublin occurred during the same period, and potentially originated from two different sources. WGS cluster A and subcluster C3 occurred in different periods indicating that they might belong to distinct outbreaks. The retrospective use of WGS also confirmed the occurrence of different S. Dublin outbreaks in 2012 [10,11]. High resolution molecular tools like WGS may facilitate linkage of human cases to sources, especially in serotypes with limited intrinsic genetic variation, and may also provide a more detailed picture of the extent and context of the outbreak.
Recruiting controls from an online health cohort survey for an ongoing outbreak investigation was novel in France and served as a pilot to evaluate the suitability of this method in future food-borne outbreak investigations. This method allowed conducting the case-control study in a timely manner with minimum resources, achieving a high response rate [28].
Our investigation pointed towards several cheese producers from the same region as sources of the outbreaks. In this region, an increase in salmonellosis incidence was observed in cattle at the end of summer 2015 (data not shown; Santé publique France). This could explain the increase of contaminated cheese batches in autumn and winter 2015. Veterinary and food investigations were challenging due to (i) high number of possibly implicated processing plants of raw-milk cheeses, and (ii) the high frequency of cheese consumption by cases and the variety of cheeses and places of purchase. It was difficult to identify the exact batches that cases consumed because some cheeses were sold at the deli counter, sliced on demand. Furthermore, the probable low levels of contamination of the implicated cheeses may have led to false negative test results, possibly allowing some contaminated batches to enter the market.
Previous studies indicated that S. Dublin is frequently isolated in live cattle, and that cheeses made with unpasteurised milk may be contaminated with S. Dublin [8,9]. S. Dublin infection is also responsible for substantial losses in the dairy industry [29]. A modelling exercise in Denmark [29] estimated the gross margin losses due to S.Dublin infection in a 200-cow stall-herd to be up to EUR 188 per stall annually averaged over the 10-year period following introduction of infection. In that study, relative simple and cheap control measures such as improving calving and colostrum management could lead to significant decreases in prevalence of S. Dublin in some herds. In other herds, it was reported that these measures might have to be supplemented by changes in hygiene and feeding practices. It might be worthwhile conducting studies in France to evaluate the impact of recent or future measures on S. Dublin prevalence at herd level. In addition, infected cattle might carry chronic and possibly asymptomatic infections while still contributing to onwards transmission by excreting pathogens in faeces [30,31]. SNP-typing based on WGS is a promising tool to monitor the routes and the spread of S. Dublin between herds in traditional regions of cheese production, as already reported in previous studies [16,17]. The combination of epidemiological studies in human and nonhuman sectors and the use of WGS may improve the cost effectiveness of control measures for S. Dublin in France, by targeting contaminated herds.

Limitations
Our investigations suffered from several limitations. First, cases were interviewed by phone, while controls completed a shorter online questionnaire, which could have led to obtaining exposure data with different degrees of accuracy. To minimise this bias, we used the same questions for cases and controls. Second, age and sex distribution of controls differed from that of the cases. We thus included age and sex in the multivariable analysis to adjust for those characteristics. Third, it was difficult to identify the exact sources of contamination due to probably low levels of contamination by S. Dublin of the cheese batches. As the cattle contamination was diffuse, it was difficult to incriminate specific cheese producers in the Bourgogne-Franche-Comté region as sources of contamination. Furthermore, suspected batches of aOR: adjusted odds ratio; CI: confidence interval; MLVA: Multilocus variable-number tandem repeat analysis; NA: not applicable; NS: did not remain significant in the final logistic regression model; OR: odds ratio; WGS: whole genome sequencing. a Only cheeses consumed by at least 50% of the cases are included in the analysis; cases and controls could report more than one cheese consumed. b Cases who belong to the other known subtypes. c Compares cluster cases with other cases with other known subtypes. d Compares cluster cases with controls. e Adjusted for age and sex. cheese identified through trace-back investigations were no longer available for testing. Then, even if epidemiological investigations were carried out within a very constrained period of time to allow the ad hoc microbiological analyses to be launched to support the generated hypotheses, we had to deal with the impossibility to get through the whole process due to the lack of products.

Conclusions and recommendations
Microbiological, epidemiological and environmental evidence pointed towards two raw-milk cheeses, Morbier and Vacherin Mont d'Or, as vehicles of the S. Dublin infections. The use of MLVA and WGS subtyping methods allowed the identification of different clusters and of the potential vehicles of infection, highlighting the importance of adequate subtyping methods during Salmonella outbreaks and the relevance of company internal microbiological monitoring system. As a result, WGS has now been routinely implemented at the French NRC and findings of this multi-disciplinary investigation led to a reinforced control plan for processing plants of raw-milk cheeses to prevent future outbreaks.