During the development phase, 12 S. aurantiacum strains were selected as tester strains from the Molecular Mycology Research Laboratory Culture Collection at Sydney Medical School - Westmead Hospital (as indicated by the strain numbers in bold and italics in Supplementary Table S2 ). These strains were selected to represent diverse geographic regions, had known clinical associations, and known genetic characteristics (identical and diverse genotypes) as established by PCR fingerprinting and AFLP analysis (Delhaes et al., 2008). In the application phase, the developed scheme was applied to a total of 188 strains that comprised 106 clinical, 1 veterinary and 81 environmental strains. Details of these strains are provided in Supplementary Table S2 . Most of strains originated from Australia (n = 84), followed by France (n = 48), Austria (n = 15), Germany (n = 13), Spain (n = 5), New Zealand (n = 4), the UK (n = 3), Nepal (n = 2), Thailand (n = 2), Ireland (n = 1), and the USA (n = 1). All isolates were grown on Sabouraud dextrose agar (Oxoid, Hampshire, UK) and incubated at 30°C for 5 to 7 days before DNA extraction.

Extraction of genomic DNA was performed according to a previously published protocol (Ferrer et al., 2001), with minor modification. The mycelia from 5-day-old cultures were harvested and placed in 1.5 ml tubes. After washing in deionized water, mycelia were frozen in liquid nitrogen. Using a sterile miniature pestle, frozen mycelia were finely ground to disrupt the fungal cell walls; 500 µl of SDS lysis buffer and 5 µl of 2-mercaptoethanol were then added and the mixture was mixed vigorously by vortexing. After incubation at 65°C for 1 hour, with 2-3 times vortexing in between, 500 µl of phenol:chloroform:isoamyl alcohol (25:24:1) (Sigma, St. Louis, USA) were added. The tubes were flipped for 2 minutes to ensure thorough mixing, followed by centrifugation at 14,000 rpm for 15 minutes. DNA from the aqueous phase was transferred to a fresh 1.5 ml tube. An equal amount of isopropanol (Merck, Kilsyth, Australia) was then added to precipitate the DNA. The tubes were then incubated at - 20°C for a minimum of 1 hour but usually overnight to increase DNA yield. The precipitated DNA was pelleted by centrifugation at 14,000 rpm for 15 minutes. After washing with 500 µl 70% ethanol (Merck) and centrifugation, the DNA pellet was dried at room temperature and reconstituted in 100 µl of sterile distilled water. DNA concentration was determined spectrophotometrically.

In the preliminary development stage, 24 gene loci, namely AAT, ACT, ANXC4, ATP6, BGT, BT2, EF1a, CAL, CAT, CHS, D1D2, FKS, LIP, GLN, IGS, MDH1, mtSSU, MP1, RPB1, RPB2, SOD2, TUB, VPS13 and ZRF, were amplified from the 12 tester strains (see above) to identify the most polymorphic loci. These loci had been previously utilized in either phylogenetic and/or genotyping studies of other Scedosporium species and/or other fungi, such as Candida, Aspergillus and Penicillium species (O’Donnell et al., 1998; Liu et al., 1999; Bougnoux et al., 2002; Cruse et al., 2002; Dodgson et al., 2003; Gilgado et al., 2005; Bain et al., 2007; Hoffman et al., 2007). For each genetic locus, the 12 sequences obtained from the tester strains were aligned using BioEdit™ Sequence Alignment Editor (Tom Hall, Carlsbad, USA) to determine the sequence variation. The genetic loci, which demonstrated a high polymorphism and yielded the largest number of sequence types in combination, were selected to form the new S. aurantiacum MLST scheme.

Twenty-four loci were selected for initial screening of genetic polymorphisms (see above). For the PCR amplification of these loci, primers were used as previously published except for SOD2 for which primers were specifically designed in the current study (O’Donnell et al., 1998; Liu et al., 1999; Bougnoux et al., 2002; Cruse et al., 2002; Dodgson et al., 2003; Gilgado et al., 2005; Bain et al., 2007; Hoffman et al., 2007) ( Supplementary Table S1 ). Partial or full SOD2 gene sequences were obtained from the GenBank database (http://www.ncbi.nlm.nih.gov). Sequences from as many fungal genera as possible for SOD2 were downloaded and then aligned using the program BioEdit™ Sequence Alignment Editor. Initial primers for SOD2 were designed based on the aligned sequences. Following the initial amplification and selection of the six most polymorphic loci, S. aurantiacum specific primers were subsequently designed for those loci ( Table 1 ) based on the obtained sequences and on the nucleotide sequence of these loci identified via BLAST searches in the draft genome of S. aurantiacum strain WM 09.24 (available in the DDBJ/EMBL/GenBank under the accession number JUDQ00000000) (Pérez-Bercoff et al., 2015).

PCR amplifications were performed in a total volume of 50 µl. Generally, each reaction mixture contained: 1x PCR buffer (20 mM Tris-HCl, 50 mM KCl), 1.5-2.5 mM MgCl₂, 100-200 µM of each deoxyribonucleotide triphosphate (dATP, dCTP, dGTP, dTTP) (Invitrogen, Carlsbad, USA), 0.2-0.4 mM each of forward and reverse primer, and 1.25 U of DNA polymerase (Bioline™, London, UK). 20-60 ng of template DNA was added to the reaction mixture. Sterile distilled water in place of DNA was used as a negative control. The primer sequences used in the development phase are listed in Supplementary Table S1 and for the final consensus MLST scheme in Table 1 . Initial PCR amplifications were performed in a thermal cycler (Perkin Elmer Cetus, Norwalk, USA) under the following conditions: an initial denaturation at 94°C for 5-10 minutes, followed by 35 cycles of 94°C for 45 seconds at a temperature ranging from 50-60°C for annealing depending on the amplified genetic loci (see Supplementary Table S1 ), followed by an extension step 1 min at 72°C and a denaturation step of 45 seconds at 94°C, with a final extension step at 72°C for 10 min. Optimized annealing temperatures for the primer sets used in the final MLST scheme are as listed in Table 1 . PCR products were separated on 1.4% agarose gels in Tris-borate-EDTA (TBE) buffer, stained with ethidium bromide (Sigma) and visualized by UV transillumination. The products were purified using PureLink™ PCR Purification Kit (Invitrogen) following the manufacturers protocol and the concentration of purified DNA was measured spectrophotometrically. DNA sequencing was performed by Macrogen Inc., Seoul, Korea (http://www.macrogen.co.kr/eng/sequencing), and the Australian Genome Research Facility (AGRF) Pty. Ltd., St. Lucia, Queensland, Australia (http://www.agrf.org.au). The quality of nucleotide sequences was verified by aligning both forward and reverse strands using the software Sequencher™ 5.4 (Gene Codes Corp., Ann Arbor, USA).

During the development phase, the consensus sequences of the 12 selected strains were aligned for each genetic locus. The genetic loci that were most polymorphic were selected for inclusion in the final S. aurantiacum MLST scheme. In the application phase, all obtained sequences from the 188 strains were aligned and analyzed. Sequence alignments were performed using BioEdit™ Sequence Alignment Editor. For each locus, numbers were assigned to designate unique allelic variants, with a single bp difference resulting in a new allele type. These numbers were subsequently combined to yield unique sequence types (see Supplementary Table S2 ). A S. aurantiacum MLST database using the BioloMICS software version 21.07.9.324 (BioAware, Hannut, Belgium) was constructed for the six loci at the Molecular Mycology Research Laboratory and can be accessed at http://mlst.mycologylab.org. GenBank accession numbers for all MLST sequences generated in this study are listed in Supplementary Table S3 . All allele and sequence types can be accessed via the specific MLST S. aurantiacum website at http://mlst.mycologylab.org. Phylogenetic trees were constructed using the software MEGA version 11 (The Biodesign Institute, Tempe, USA) (Tamura et al., 2021). To further investigate the geographical relationship between genotypes, a goeBURST minimum spanning tree was generated from the aligned concatenated sequences of the strains studied using the PHYLOViZ 2.0 analysis software (http://www.phyloviz.net/).

Assessment of the likelihood of selective pressure at each locus was estimated by the ratio of non-synonymous to synonymous nucleotide substitutions (d_N/d_S ) (Nei and Gojobori, 1986). To test for purifying selection the codon-based Z-test using evolutionary pathway by Nei and Gojobori was performed (Nei and Gojobori, 1986) using MEGA version 11 (Tamura et al., 2021). To evaluate the variability of the selected loci, the haplotype diversity (Hd), nucleotide diversity (π) and the average number of nucleotide differences (k) were determined using the software DNA Sequence Polymorphism DnaSP version 6.12.03 (Rozas et al., 2017). Testing for neutrality utilizing the Tajima’s D test (equal to zero at neutral equilibrium) (Tajima, 1989) was performed in MEGA version 11 (Tamura et al., 2021).

Intragenic linkage disequilibrium (LD), and intragenic recombination rates were calculated by using DNA Sequence Polymorphism DnaSP version 6.12.03 (Rozas et al., 2017). Evidence of recombination was shown using the 4-gamete test to infer the minimum number of recombination events (Rm).

Patients’ demographic and clinical data were recorded. Associations between genotypes and clinical variables were explored. The variables examined included age, sex, and geographical origin, source of isolates, infection status and predisposing factors ( Supplementary Table S2 ). We investigated potential associations between sequence type and the variables using Pearson Chi Square test. Statistical analysis was performed using PASW Statistics 18 (SPSS Inc., Chicago, IL, USA) and STATA 15 (StataCorp LP, College Station, USA). P-values of <0.05 were considered statistically significant.

Virulence studies based upon in vivo survival in a murine model (Harun et al., 2010b) were performed. Eighteen strains listed in Supplementary Table S2 (indicated by underlined strain numbers) and in Figure 4 were selected to represent a wide global spectrum. Five seven-week-old female Balb/C mice were used, in which 0.2 ml of a conidial suspension (10⁶ conidia/ml) was inoculated intravenously via the lateral tail vein. The mice were monitored daily till 30 days post-inoculation for signs of infection, including ruffling of fur, inactivity, loss of weight, difficulty in breathing, and neurological signs such as ataxia. In accordance with the protocol approved by Western Sydney Local Health District Animal Ethics Committee (WSLHD AEC), mice that were deteriorating prior to that end point were sacrificed. Mean survival time (MST), estimated by Kaplan-Meier method, was used as a parameter to compare the relative pathogenicity among the selected strains. Comparison between each group was performed by the log-rank test as part of the software package PASW Statistics 18. Graphs were plotted using GraphPad Prism version 5.0b (GraphPad Software Inc., San Diego, USA). P-values of < 0.05 was considered as statistically significant.

DNA sequences of 24 genetic loci (AAT, ACT, ANXC4, ATP6, BGT, BT2, EF1a, CAL, CAT, CHS, D1D2, FKS, LIP, GLN, IGS, MDH1, mtSSU, MP1, RPB1, RPB2, SOD2, TUB, VPS13, and ZRF,) ( Supplementary Table S1 ) were screened to assess their polymorphisms and hence suitability as candidate genetic loci for the S. aurantiacum MLST scheme. Following analysis of the DNA sequences obtained from 12 tester strains, the following six loci, actin (ACT), calmodulin (CAL), elongation factor 1-alpha (EF1a), RNA polymerase II subunit (RPB2), manganese superoxide dismutase (SOD2) and beta tubulin (TUB), were found to be the most variable, and were therefore selected to form the S. aurantiacum MLST scheme (http://mlst.mycologylab.org) ( Table 1 ), and to identify the allele and sequence types of S. aurantiacum strains.

The sizes of the six MLST gene fragments obtained including all gaps were: 830 bp for the ACT locus, 689 bp for the CAL locus, 715 bp for the EF1α locus, 952 bp for the RPB2 locus, 437 bp for the SOD2 locus and 401 bp for the TUB locus ( Table 2 ). Seventy-seven (1.89%) polymorphic sites were identified across all six genes combined, which represents a total of 4,024 bp. The number of variable nucleotide sites per locus ranged between 5 (0.73%, CAL) and 18 (4.12%, SOD2). The variability among loci is shown in Table 2 .

The ratio of non-synonymous to synonymous nucleotide substitutions (d_N/d_S ) was < 1 for five of the six MLST loci studied ( Table 2 ). In all loci, the probability (p value was > 0.05 and therefore the null hypothesis of strict neutrality (dN=dS) was not rejected. All but one locus showed an dN-dS of negative value (dN<dS, dN/dS ratio of <1), thus indicating that with exception to SOD2 none of the studied loci were under positive selective pressure ( Table 2 ). The nucleotide diversity (π) ranged from 0.00196 for CAL to 0.01335 for SOD2, the haplotype diversity (Hd) ranged from 0.582 for TUB to 0.864 for RPB2, and the average number of nucleotide differences (k) ranged from 1.35305 for CAL to 5.78177 for SOD2, indicating an equal range for all loci, except for SOD2 ( Table 2 ). In the neutrality test, none of the Tajima’s D values obtained for the studied loci deviated significantly from zero, except for SOD2, suggesting the occurrence of balancing selection (Tajima, 1989) ( Table 2 ).

The intragenic recombination test identified 1-4 recombination events (Rm) at the ACT, EF1a, RPB2, SOD2, and TUB loci, but no recombination at the CAL locus ( Table 3 ). Based on the concatenated multilocus sequence data, the interlocus LD was assessed over all segregating sites using pairwise comparisons. The LD (|D′| Y = 0.8257–0.0970X) was detected with a negative slope, indicating a decrease in linkage with increased nucleotide distance. Of the 2556 pairwise comparisons, 798 were significant by the Fisher exact test, and 229 were significant after Bonferroni correction ( Table 3 ).

The sequence alignments showed polymorphisms in all six loci denoting the presence of different alleles ( Table 2 ). Each specific sequence was considered a unique allele type, which was then assigned a unique allele type number. For example, 22 alleles were defined for the ACT locus, which were assigned as allele types AT1 to AT22, accordingly. The CAL locus showed 8 alleles, the EF1α locus 23 alleles, the RPB2 locus 15 alleles, the SOD2 locus 18 alleles and the TUB locus 12 alleles ( Table 2 ). A total of 159 sequence types were obtained by combining the allele types of the six MLST loci studied ( Supplementary Table S2 ). Most of the strains exhibited unique sequence types ( Supplementary Table S2 and Figures 2 , 3 ). Only a few strains shared the same sequence type, but most of them originated from the same patient or from closely linked soil samples. For example, strains IHEM 23081 and IHEM 23092, which both exhibited the sequence type ST140, were recovered from respiratory secretions of the same patient at a one-year interval; likewise, six strains shared the sequence type ST108, i.e. 110349103-01/1, 110349103-01/2 and 110349103-01/3, as well as 110349211-01/1, 110349103-01/2 and 110349103/3, but they were recovered from two soil samples collected at the same location on the banks of the Loire river, in France. On the contrary, strains IHEM 23080 and IHEM 23081, which were recovered from the same clinical sample, exhibited distinct sequence types (ST 139 and ST140, respectively).

Figure 1 shows the distribution of sequence types according to the geographical origin. 69 of the 159 (43%) sequence types were only found in Australia, but none of them were predominant. Most of the Australian strains were obtained from New South Wales (NSW), suggesting a possible study bias. Among the 61 studied strains from NSW, 49 different sequence types were identified, most being specifically identified from this state since only two of these sequence types were also found outside NSW (ST12 and ST15, which were also found in Western Australia (WA)). Twenty strains collected in WA were studied, which revealed 19 sequence types, 17 of them being specifically identified from WA. ST26 and ST27 were found only in South Australia (SA), and ST28 was found only in Queensland (QLD) ( Figure 1 and Supplementary Table S2 ). The sequence types ST39, ST40, ST41 and ST42 were found exclusively in New Zealand.

Ninety-five strains collected in Europe were analyzed in this study. MLST analysis of these strains revealed a total number of 81 sequence types. Almost all of them were country specific, with the highest number of sequence types (39 STs) being found in France ( Figure 1 and Supplementary Table S2 ). A unique sequence type was identified from two distinct countries, i.e. the sequence type ST82, which was isolated from both Germany and Austria, with no obvious connection of the patients to each other.

Other sequence types obtained in this study include ST55, which was unique to the United States, ST99 and ST100, which were specific to Thailand, and ST105 and ST106, which were found in Nepal ( Figure 1 and Supplementary Table S2 ).

Parsimonious trees were constructed for each MLST locus ( Supplementary Figures S1 – S6 ). These analyses resulted in different tree topologies, indicating variable rates of gene evolution for each genetic locus ( Supplementary Figures S1 –S6 ).

Maximum parsimony analysis of the combined gene sequences obtained from the six loci studied revealed a high genetic diversity amongst the 188 S. aurantiacum strains investigated ( Figure 2 ), with most strains forming unique sequence types, which resulted in the widespread topology of the combined tree ( Figure 2 and Supplementary Table S2 ). Sequence types of the strains from different Australian states demonstrated no tendency to group with each other. When comparing the two major clusters ( Figure 2 ), there was a significant difference according to the geographic origin of the strains (p < 0.0005). Most Australian strains belonged to cluster 2, while almost all European strains, except one German strain RKI94-0197, were grouped together in cluster 1 ( Figure 2 ). The four isolates from New Zealand grouped either basal to cluster 1, which may be called ‘global cluster’ (WM 07.96, WM 07.97) or to cluster 2, the ‘Australian cluster’ (WM 07.101, WM 07.108) ( Figure 2 ). In between the two major clusters the central basal group contains strains from Australia, Austria, France, Germany, Ireland, Netherlands, Spain, and the United Kingdom.

The goeBURST analysis of the obtained MLST data confirmed the high genetic diversity seen in the Australian S. aurantiacum population, as well as in the non-Australian population ( Figure 3 ). Most of the Australian strains grouped separate to the European strains, forming a closely connected gene network. However, a subset of Australian strains was intermixed with and are closely related to the global strains, indicating certain genetic links between Australian and non-Australian strains ( Figure 3 ). The results of the goeBURST analysis ( Figure 3 ) confirmed the genetic relationships identified in the phylogenetic analysis of the combined genes as reflected in the concatenated gene tree ( Figures 2 , 3 ).

Both clinical and environmental strains were widely distributed throughout the tree showing no sequence types to be indicative of a human clinical, veterinary, or environmental origin. Similar findings were found for the distribution of strains recovered from patients with either “invasive” disease or “colonization” (hereafter referred to as “invasive” or “colonizing” strains) ( Figure 2 ).

Overall, the statistical analyses (PASW Statistics and STATA II) performed showed no significant associations between the different sequence types, which occurred largely as singletons, and clinical variables. In addition, there were no significant associations between a particular cluster and patient age (p = 0.078), sex (p = 0.076) or infection site (p > 0.05). Most analyzed strains were recovered from respiratory secretions from patients with chronic lung disease, such as cystic fibrosis; these strains that colonized the airways, were distributed throughout both clusters. There was only a single sequence type which was shared between colonizing and invasive strains, ST15, with WM 06.476 and WM 07.555 being colonizing strains, and WM 07.452 being an invasive strain.

When the two major clusters were compared with the patients predisposing factors (e.g., chronic lung disease, malignancy, diabetes, corticosteroid administration, chemotherapy, trauma, and drowning), there was no association between predisposing factors and any particular cluster (p > 0.05). An exception was a correlation of certain clusters and chronic lung disease, which was noted to show a significant difference (p = 0.005). In the “global cluster” there was a group of colonizing strains including WM 08.2114 and WM 08.215 (both ST53), WM 06.571 (ST37), 10-03-12 (ST149) and 10.03.10.92 (ST154), whereas the “Australian cluster” comprised another group of colonizing strains composed of WM 06.481, WM 08.209, and WM 06.480 (all ST17), WM 06.546, WM 08.210, and WM 08.211 (all ST52), WM 06.492 (ST21) and WM 06.479 (ST16) ( Figure 2 ).

The results of the survival analysis in mice, expressed as percent survival, for 18 S. aurantiacum strains are shown in Figure 4 . All infected mice showed evidence of active infection (ruffling of fur, severe weight loss and neurological abnormalities, such as ataxia), as early as day 3 post-inoculation ( Figure 4 ). An overall comparison between survival curves showed a significant difference (p = 0.007), with WM 06.482 being the most virulent strains followed by WM 09.24 and WM 08.52, and strains WM 08.269 and WM 08.202 being the least virulent strains. Of note, pairwise comparisons among all tested strains showed variable results. Pairwise comparison between invasive and colonizing strains, and between clinical and environmental strains revealed no significant difference.