Outbreak of invasive wound mucormycosis in a burn unit due to multiple strains of Mucor circinelloides f. circinelloides resolved by whole genome sequencing

Mucorales are ubiquitous environmental molds responsible for mucormycosis in diabetic, immunocompromised, and severely burned patients. Small outbreaks of invasive wound mucormycosis (IWM) have already been reported in burn units without extensive microbiological investigations. We faced an outbreak of IWM in our center and investigated the clinical isolates with whole genome sequencing (WGS) analysis. We analyzed M. circinelloides isolates from patients in our burn unit (BU1) together with non-outbreak isolates from burn unit 2 (BU2, Paris area) and from France over a two-year period (2013-2015). For each isolate, WGS and a de novo genome assembly was performed from read data extracted from the aligned contig sequences of the reference genome (1006PhL). A total of 21 isolates were sequenced including 14 isolates from six BU1 patients. Phylogenetic classification showed that the clinical isolates clustered in four highly divergent clades. Clade1 contained at least one of the strains from the six epidemiologically-linked BU1 patients. The clinical isolates seemed specific to each patient. Two patients were infected with more than two strains from different clades suggesting that an environmental reservoir of clonally unrelated isolates was the source of contamination. Only two patients shared one strain in BU1, suggesting direct transmission or contamination with the same environmental source. WGS coupled with precise epidemiological data and analysis of several isolates per patients revealed in our study a complex situation with both potential cross-transmission and multiple contaminations with a heterogeneous pool of strains from a cryptic environmental reservoir. Importance Invasive wound mucormycosis (IWM) is a severe infection due to the environmental molds belonging to the order Mucorales. Severely burned patients are particularly at risk for IWM. Here, we used Whole Genome Sequencing (WGS) analysis to resolve an outbreak of IWM due to Mucor circinelloides that occurred in our hospital (BU1). We sequenced 21 clinical isolates, including 14 from BU1 and 7 unrelated isolates, and compared them to the reference genome (1006PhL). This analysis revealed that the outbreak was mainly due to multiple strains that seemed patient-specific, suggesting that the patients were more likely infected from a pool of diverse strains from the environment rather than from direct transmission between the patients. This study revealed the complexity of a Mucorales outbreak in the settings of IWM in burn patients, which has been highlighted based on whole genome sequencing and careful sampling.

year period (2013)(2014)(2015). For each isolate, WGS and a de novo genome assembly was 48 performed from read data extracted from the aligned contig sequences of the reference 49 genome (1006PhL). 50 A total of 21 isolates were sequenced including 14 isolates from six BU1 patients. 51 Phylogenetic classification showed that the clinical isolates clustered in four highly 52 divergent clades. Clade1 contained at least one of the strains from the six 53 epidemiologically-linked BU1 patients. The clinical isolates seemed specific to each 54 patient. Two patients were infected with more than two strains from different clades 55 suggesting that an environmental reservoir of clonally unrelated isolates was the source 56 of contamination. Only two patients shared one strain in BU1, suggesting direct 57 transmission or contamination with the same environmental source. 58 WGS coupled with precise epidemiological data and analysis of several isolates per 59 patients revealed in our study a complex situation with both potential cross-60 transmission and multiple contaminations with a heterogeneous pool of strains from a 61 cryptic environmental reservoir. 62 63 64 Phylogenetic analyses of three loci 149 ITS, D1/D2, and RPB1 sequences were analyzed as separate (data not shown) and 150 combined datasets. The topology of the multi-locus dataset with 3 methods (NJ,ML and 151 Bayesian inference) was comparable between the individual trees of the three genes 152 analyzed. Four clades (Fig. 2), respectively denoted C1 (14 isolates including 11 from 153 BU1), C2 (four including two from BU1), C3 (two strains in addition to the reference 154 strain including one from BU1) and C4 (one isolate) were identified from the analysis of 155 the combined dataset that yielded a significant support (≥95% bootstrap for NJ and ML; 156 1.0 for Bayesian inference). Isolates recovered from BU1 were distributed in three 157 clades (C2, C3, C4). All patients from BU1 had at least one isolate included in C1. 158 159

Whole genome analysis 160
To better resolve the diversity of the strains within the four clades, and because no 161 further genotyping methods existed for this organism, whole genome sequencing was 162 performed. Because the biology and the genetics of this organism is poorly understood, 163 we first checked the reproducibility of the sequencing process and the stability of the 164 genome, to be able to define genetically-identical strains. 165

166
Establishing the genetic threshold to determine genetically-identical strains 167 For the three strains isolated from single spore colonies (i.e. P05_600_BU1_SLS, 168 P04_603_BU1_SLS and P03_594_BU1_SLS), the pairwise evolutionary distance was 169 estimated between de novo assemblies of the parental and the single spore colony (i.e. 170 0.00035, 0.00044 and 0.00029, respectively). This information gave us the expected 171 distances between pairs of genomes arising from identical strains and independently 172 sequenced distinct isolates. The largest of the three distances (i.e. 0.00044) was 9 therefore selected as a cutoff below which two compared isolate genomes were defined 174 A phylogenetic classification of the whole genome of the 21 isolates was then performed 184 ( Fig. 3). This phylogenetic tree allows classification of the genomes in four main clades 185 corresponding exactly to the same clades (C1 to C4) described based on the analysis of 186 three loci (Fig. 2). Clades C2, C3, and C4 contained isolates that are clearly distinct from 187 those inside clade C1, e.g. estimated pairwise evolutionary distances between isolates 188 from C1 and those inside C2, C3, and C4 are 0.0187, 0.0385, and 0.0390 on average, 189 respectively, whereas C1 pairwise intra-distance is 0.0014 on average. 190 Furthermore, increased resolution of WGS allowed to robustly identify strains and 191 understand which clinical isolates belong to which strain. Indeed, the 21 investigated 192 isolates and the reference 1006PhL could be partitioned into 14 distinct strains (S1 to 193 S14; Fig. 3 In BU1, the isolate of the potential index case P01 (P01_617_BU1_SLS, S12) was different 197 from the isolates subsequently recovered in BU1. The isolate from P02, corresponding to 198 the specific strain S5 clustered in C1 together with isolates from P03, P04, P05, and P06. 199 P03 and P06 were also infected with two isolates from S1 (#594 and #592) and S6 200 (#032 and #023), respectively. P04 was infected over 10 days with two strains S9 and 201 S7 belonging to C2 (isolates #601 and #559) and C1 (isolates #602 and #603), 202 respectively. P05 was infected over 11 days with isolates belonging to three strains from 203 C1, S8 (isolates #598 and #600) and S3 (#622) and S1 (#599). Indeed, P04 and P05 had 204 mixed infections during the course of their disease, suggesting initial contamination 205 with a mixture of strains, this latter strain also recovered in P03. This suggests cross 206 contamination or common infection in P3 and P5. This have also been observed in BU2,207 where two patients (P09 and P10) shared the strain S2 (#703 and #704) that clustered 208 in C1. Two patients (P11 and P08) were also infected with strains that belong to C2 and 209 C3, respectively. The environmental (#615) and the colonization (#621) isolates from 210 our hospital clustered in C2 and C4 respectively. The patient from Eastern France 211 clustered in C1 (S4) but with a specific and different strain than the other C1 strains, as 212 expected for a geographically unrelated infection. 213

DISCUSSION 214
Because isolation of Mucorales is rare in the hospital, the observation of the same 215 species in two independent samples/patients has long been considered as a sufficient 216 criterion to suspect and assess transmission or common contamination. Here, we 217 investigated further outbreak-related and -unrelated isolates based on WGS analysis. To 218 our knowledge, this is the first time WGS analysis has been utilized to resolve an 219 outbreak of invasive Mucorales infection in the context of nosocomial acquisition. 220 Because WGS was applied for the first time in this setting, we first sought to evaluate 221 reproducibility of the sequencing process and intra-culture variation/stability during in 222 vivo passage. As the genome of three selected strains was sequenced and assembled 223 twice, we were able to compare contig sets belonging a priori to the same strain. This 224 method led to the definition of a pairwise distance cutoff. Therefore, if two genome 225 sequences belonging to different isolates have a pairwise distance below this cutoff, it 226 was inferred a posteriori that the two isolates corresponded to the same strain. This cut-227 off should vary as a function of the organism, the method of sequencing, the 228 bioinformatics pipeline and the pathophysiology of the disease, reinforcing that such 229 data should be obtained each time an investigation of an outbreak due to rare organisms 230 is undertaken. 231 Of note, using pairwise evolutionary distances could be considered as a fast but accurate 232 alternative to the well-known Average Nucleotide Identity (ANI) approach to compare 233 genomes (26, 27) because both were shown to be linearly correlated (28). 234 Contrary to the initial hypothesis of a single strain transmission in BU1, we observed 235 that all of the patients from the BU1 outbreak (P01 to P06) were infected by different 236 strains. Surprisingly, our data revealed that each strain was patient-specific in BU1, 237 except for S1, suggesting that the outbreak in BU1 was due to multiple strains present in 12 and acquired from a local environmental "reservoir" containing clonally unrelated 239 isolates. This hypothesis is reinforced by two patients (P04 and P05) with IWM co-240 infected by more than 2 genetically distinct strains. Another hypothesis is that the 241 patients could have been exposed to specific strains or a mixture of strains before 242 arriving in BU1. However, a delay between admission and the first positive culture was 243 16 days (median), making the hypothesis that exposure occurred in the environment of 244 BU1 more likely. In the settings of severe burns where invasive fungi do sporulate on the 245 wounds, transmission by air or by the hands of healthcare workers to other patients are 246 both possible. A major point to emphasize is that only two patients from BU1 (P03 and 247 P05) shared the same strain S1. This feature has also been observed in BU2 suggesting 248 that transmission between patients is possible. However, this does not rule out the 249 hypothesis of a contamination by the same strain from the environment. 250 Environmental investigation of the outbreak in BU1 failed to identify the source of 251 infection using culture of multiple samples from the environment, as well as by PCR. 252 DNA amplified with the Mucor/Rhizopus PCR (25) was detected only in the Bair Hugger 253 filters that were used during the hospitalization of the IWM patients P03, P04 and P05. 254 Despite the negative result of this investigation, it is likely that the contamination came 255 from a local source because this has already been described with linens or 12lastoplast 256 in burn units (5,18,19). 257 Our findings about the genetic structure of M. circinelloides f. circinelloides are 258 reminiscent of the WGS investigations of the Apophysomyces spp. outbreak in Joplin (16), 259 or of the S. clavata outbreak in France (29), which revealed that several genetic groups 260 can be responsible for infections over the same period of time. In our case, in a given 261 restricted area (BU1), we identified a large diversity of isolates responsible for IWM and 262 were not able to find isolates belonging to unique strains recovered in different places.

13
At the other end of this spectrum, for Apophysomyces trapeziformis several genetically-264 identical isolates were recovered in different places at a distance of several miles (16). 265 In the case of E. rostratum, all outbreak isolates have closely-related genomes suggesting 266 that a unique strain was responsible for the outbreak (30). These differences could be 267 In order to study the genetic relatedness between the BU1 clinical isolates, additional 300 isolates identified as M. circinelloides from other sources were selected (Table 1)

Polyphasic identification of isolates 329
The 21 isolates were sent to the NRCMA, where the purity was verified and 330 identification to the species-level performed using phenotypic and molecular 331 identification. In details, microscopic examination was performed on 5 to 7 day old 332 cultures growth on 2% malt agar at 30°C. Amplification and sequencing of the ITS1-5.8S-333 ITS2 region and the D1/D2 region of the LSU rDNA were performed as described 334 previously (35). The amplification of the RPB1 gene (RNA polymerase II largest subunit) 335 was made with primers RPB1Ac and RPB1Cr (36). The PCR products were then 336 sequenced and the consensus sequences were obtained as already described (35). 337 Sequences were subjected to pairwise alignments against curated fungal reference 338 databases available at the on-line MycoBank database (http://www.mycobank.org/). 339 340

Reproducibility of the sequencing process 373
In a specific experiment dedicated to determine the reproducibility of all the sequencing 374 process (from extraction to sequence analysis), additional single spore isolation was 375 performed on malt extract agar 2% plates for 3 isolates recovered from patients P03, 376 P04 and P05 in BU1 (P05_600_BU1_SLS, P04_603_BU1_SLS, and P03_594_BU1_SLS, see 377 Table 1). One colony was thus selected from each isolate for additional sequencing. 378 379

Genome stability experiments 380
Genome stability during vegetative growths and host infections was analyzed with 381 1006PhL strain. In brief, the strain was grown on PDA for one day. The colony was 382 streaked to isolate a single colony, which then transferred onto a new PDA. After four 383 days of incubation at 30°C under the light, the spores were collected and subjected to 384 the same procedures up to three passages (VEG1, VEG2 and VEG3 isolates). 385 For infection passages, spores were suspended in sterile PBS. Male 8 weeks-old BALB/c 386 mice were infected with 10 6 spores in 200 μL of sterile PBS via tail vein injection. At day 387 3 post-inoculation, the mice were sacrificed, their brain collected and placed onto PDA 388 after brain homogeneization. After one day of incubation, the fungal colony emerging 389 from the brain was streaked to isolate single colonies. One colony was transferred onto 390 PDA and incubated at 30°C under light for four days. Spores were then used as inoculum 391 Epidemiological map of thirteen patients whose isolates were selected in our study. 420 Isolates from several geographic areas in France (from dark to light gray) have been 421 studied: burn unit 1 of Hospital 1, Paris France; wards of SLS, Paris, France; Burn unit 2 422 of PER in Paris area, France and eastern France (one isolate from STR).