Molecular Epidemiology of Azole-Resistant Aspergillus fumigatus in France Shows Patient and Healthcare Links to Environmentally Occurring Genotypes

Resistance of the human pathogenic fungus Aspergillus fumigatus to antifungal agents is on the rise. However, links between patient infections, their potential acquisition from local environmental sources, and links to global diversity remain cryptic. Here, we used genotyping analyses using nine microsatellites in A. fumigatus, in order to study patterns of diversity in France. In this study, we genotyped 225 local A. fumigatus isolates, 112 azole susceptible and 113 azole resistant, collected from the Bourgogne-Franche-Comté region (Eastern France) and sampled from both clinical (n = 34) and environmental (n = 191) sources. Azole-resistant clinical isolates (n = 29) were recovered mainly from cystic fibrosis patients and environmental isolates (n = 84) from market gardens and sawmills. In common with previous studies, the TR34/L98H allele predominated and comprised 80% of resistant isolates. The genotypes obtained for these local TR34/L98H isolates were integrated into a broader analysis including all genotypes for which data are available worldwide. We found that dominant local TR34/L98H genotypes were isolated in different sample types at different dates (different patients and types of environments) with hospital air and patient’s isolates linked. Therefore, we are not able to rule out the possibility of some nosocomial transmission. We also found genotypes in these same environments to be highly diverse, emphasizing the highly mixed nature of A. fumigatus populations. Identical clonal genotypes were found to occur both in the French Eastern region and in the rest of the world (notably Australia), while others have not yet been observed and could be specific to our region. Our study demonstrates the need to integrate patient, healthcare, and environmental sampling with global databases in order to contextualize the local-scale epidemiology of antifungal resistant aspergillosis.


INTRODUCTION
Fungal pathogens pose a growing threat to the health of humans, animals, ecosystems, food security, and the global economy, making their effective control necessary (Fisher et al., 2012). Azoles fungicides, with more than 30 compounds available , are widely used in agriculture to protect cereal, vegetables, and vines from phytopathogenic fungi, as well as in the cultivation of ornamental plants and to preserve materials such as timber. These fungicides are used worldwide, particularly in Europe and Asia, where they represent one of the most commonly used classes of pesticides Chen et al., 2016). Due to their broad-spectrum action across the fungal kingdom, the azoles are also widely used as first-line drugs in medicine (human and animals) for the treatment of superficial and invasive fungal infections. Within this context of widespread use of antifungals across different compartments, the occurrence of multi-resistant pathogenic fungi represents an emerging threat to the control of plant, animal, and human diseases . Recently reviewed (Verweij et al., 2020), substantial progress has been made towards building the argument that ecological "hotspots" occur where biotic and abiotic conditions converge to allow the growth of the thermophilic Aspergillus fumigatus in contact with sub-MIC concentrations of agricultural triazoles, generating conditions that are suitable for adaptation to drug pressure. Specifically, highly related genotypes of azole-resistant A. fumigatus have now been described in clinical and environmental samples, suggesting that humans are increasingly exposed to drugresistant aerosolized A. fumigatus spores with broad public health consequences (Sewell et al., 2019;Rhodes et al., 2021).
Studying the population genetics and evolution of pathogens can help further our understanding of the evolution and spread of antimicrobial resistance, including tracing potential sources of inocula using molecular epidemiology. Specifically, for azoleresistant A fumigatus (ARAf), it is important to understand whether genotypes are randomly distributed or cluster geographically, and whether a common source of colonization may exist (van der Torre et al., 2021). In 2012, we described a first case of invasive aspergillosis caused by TR 34 /L98H ARAf in a farmer (Rocchi et al., 2014). In 2013, we described a second case, in a woodworker, who developed sinus aspergillosis due to TR 34 / L98H A. fumigatus (Jeanvoine et al., 2016). Motivated by these findings, the mycology team of Besancon (Eastern France) undertook a systematic screening of A. fumigatus to assess the burden of azole resistance across both clinical and environmental samples. In 2017, we observed an increase in the number of cystic fibrosis (CF) patients colonized with ARAf, with 16% of patients with A. fumigatus-positive cultures observed carrying the TR 34 / L98H resistance allele (Godeau et al., 2019) whereas previously, only one patient had this type of isolate in 2015. This percentage of resistance mirrors that reported by a British team (Abdolrasouli et al., 2018) and is among the highest occurrence in CF patients observed in European studies (Amorim et al., 2010;Mortensen et al., 2011;Burgel et al., 2012;Morio et al., 2012;Fischer et al., 2014;Prigitano et al., 2017;Guegan et al., 2018;Seufert et al., 2018;Lavergne et al., 2019).
Since 2017, a number of ARAf baring the TR 34 /L98H allele were detected in the air of Besancon hospital following fungal aerocontamination monitoring. Undertaken weekly using seven air sample locations in general hospital corridors, our surveillance has detected ARAf periodically between 2017 and 2019 (Rocchi et al., 2021). In parallel, two environmental studies were conducted in agricultural fields where azoles had been sprayed, in sawmills (Jeanvoine et al., 2017) and market gardens (Rocchi et al., 2018).
Across our team, we have collected 1,100 strains of the Aspergilli complex (A. fumigatus and cryptic species) isolated in our region (rural landscape in Eastern France). Here, we assess the genetic relatedness among a retrospective sample of 225 susceptible and resistant A. fumigatus from our collection for both clinical and environmental samples by comparing them to data (worldwide scale) held within pre-existing genotype databases.

Current Systematic Screening in Besancon Hospital
For both clinical and environmental samples, A. fumigatus were first isolated on malt agar media supplemented by itraconazole (2 mg L −1 ) or voriconazole (1 mg L −1 ) (van der Linden et al., 2013;Bader et al., 2015) at 37°C for clinical ones and 48°C for environmental ones. Subsequently, resistance was measured by concentration gradient Etest strip technique (Etest, bioMeŕieux, France) and/or checked by European Committee on Antimicrobial Susceptibility Testing (EUCAST) broth microdilution method (Arendrup et al., 2013). In case of azole-resistant isolate, a 495bp section of the beta-tubulin gene (Glass and Donaldson, 1995) and the cyp51A gene and its promoter were amplified and sequenced (Mellado et al., 2001;Alanio et al., 2011).
PCRs were carried out in a 50-µl volume, containing 25 µl of GoTaq ® G2 Hot Start Green Master Mix, 2 µl of each set of primers (10 µM, Table 1), 50-100 ng of genomic DNA, and water (molecular biology quality) to 50 µl. Amplification was performed in a C1000 Touch ™ thermal cycler (BIORAD) with 2 min at 95°C, and 34 cycles of 15 s at 95°C, 15 s at 61°C (for betatubulin primer sets), or 15 s at 55°C (for cyp51A primer sets), 10 s at 72°C, and one final cycle of 5 min at 72°C. The PCR products were analyzed by electrophoresis on 1% agarose gels, visualized with SYBR ™ Safe DNA Gel Stain (ThermoFisher Scientific), and then purified (NucleoSpin ® Gel and PCR clean-up, Macherey-Nagel) to Sanger sequencing (Applied Biosystems 3130 Genetic Analyzer). and to have groups with a comparable number of isolates between the different environments. For clinical isolates, especially those from CF patients, sensitivity to azoles was assayed using Etest strip technique and confirmed by EUCAST method for only one colony growing on medium with azole from sputum until 2018. For environmental analysis, sensitivity of all strains growing on medium with azole were assayed with EUCAST method and we inoculated the samples in parallel on azole free media (DG18). Due to the large number of A. fumigatus isolated on DG18, we selected only one isolate for every 10 growing on media to test their sensitivity. Isolates growing on medium with azoles but with EUCAST MIC > 2 µg/ml were stored as susceptible strains.
The majority of environmental resistant ARAf isolates came from market gardens (n=43), sawmills (n=24), and air of hospital (n=12). For sawmills and market gardens, we chose ASAf to match the location of ARAf isolates for each professional field (preferably within the same sample when possible, or in the same field or building for sawmills or, if not possible, on the same farm). In total, we analyzed a total of 63 isolates from market gardens (43 ARAf and 20 ASAf) and 65 from sawmills (24 ARAf and 42 ASAf).
With respect to isolates from hospital air sampling, as the detection of ARAf did not always coincide with an increase in the total number of A. fumigatus in hospital corridors, we selected 45 ASAf in the following way: ASAf isolated from a contamination peak in June 2017, localized in one hospitalization ward (after a heating system malfunction), with a selection of 24 ASAf in six rooms and bathrooms of the department. ASAf were isolated during a general contamination peak in December 2018 (68 colony forming units/m 3 of air) in the general corridors of the hospital where we selected 12 ASAf. Six ASAf were isolated in June 2018, on the same day as a TR 34 /L98H ARAf and three ASAf were isolated the following week.
In total, we analyzed 112 ASAf coming from 5 clinical and 107 environmental samples and 113 ARAf coming from 29 clinical and 84 environmental samples ( Figure 1).

Genotyping Method
Microsatellite genotyping by short tandem repeats for A. fumigatus (STRAf) assay (de Valk et al., 2005) was performed to assess the genetic relatedness between 225 susceptible and resistant A. fumigatus.
Peaks were visualized with Microsatellite analysis application on ThermoFisher cloud dashboard (https://apps.thermofisher. com/apps/spa/#/dashboard). For each marker, the correct peak was selected and transformed into a repeat number based on the Af293 STRAf genotype. The number of repeats was inferred from the available genome for this strain by an automated script using R Studio software.

Data Analysis
Microsatellite genotype distances were calculated using Bruvo's distance (Bruvo et al., 2004). The analyses were done in three steps. Firstly, all susceptible and resistant isolates from the mycology team of Besancon were used. Genotypic richness was assessed by the total number of unique multilocus genotypes (MLGs) or expected MLGs (eMLGs), which approximates the number of genotypes that would be expected at the largest, shared sample size based on rarefaction. Diversity (Simpson's, Shannon-Wiener, or Stoddart and Taylor's Indexes) and evenness were also calculated for each type of isolate (azole susceptible and resistant) with the R package poppr (Kamvar et al., 2015).
Secondly, only local TR 34 /L98H isolates were analyzed in order to see the proximity according to their origin. Thirdly, we compared our local TR 34 /L98H isolates to those with available data worldwide. For this, we used a global dataset housed within the online AfumID application (https://afumid.shinyapps.io/ afumID/) (Sewell et al., 2019). We thus included a further 194 isolates that had previously been isolated and genotyped from 16 countries. In order to ensure harmonization, the number of alleles was recalibrated between the two databases using the reference strain Af293.
Data were represented by circular genetic tree build by a Neighbor-Joining method clustering from a matrix of bruvo distance (Bruvo et al., 2004) using ape library (Paradis and  1 | Sequences of primers used for beta-tubulin and cyp51A and its promoter.

Genotyping Results: Besancon Local Scale
The total number of MLGs was higher in the sensitive group than in the resistant group with 102 MLGs found in the 112 ASAf isolates and 53 MLGs observed among the 113 barcoded ARAf.
When comparing the two most abundant groups (ASAf and TR 34 /L98H ARAf), it was apparent that genotypic richness, represented by the total number of MLG or expected MLG (eMLG), was lower in the TR 34 /L98H ARAf group ( Table 2). Regardless of the index used (Simpson's, Shannon-Wiener, or Stoddart and Taylor's Indexes), the genotype diversity was also lower in the TR 34 /L98H ARAf group. The evenness measure ( Table 2) indicates that the MLGs observed in the ASAf population are closer to equal abundance (which would be equal to 1) than those in the TR 34 /L98H ARAf population. Indeed, when we look at the distribution of MLGs for each group, there are many more unique ASAf profiles than ARAf profiles.
When one MLG was shared by isolates of ASAf, it was only shared by two isolates and from a similar source (clinical or environmental). Only one MLG was found both in a clinical isolate (1E005, CF patient June 2017) and in an environmental isolate (2E042, hospital air June 2018).
One MLG was shared by sensitive and resistant A. fumigatus isolated from one CF patient: one susceptible (1E007) and two resistant (2E022 and 1E008) isolates. These three strains had F46Y/M172V/N248T/D255E/E427K polymorphism in cyp51A.     For ARAf, five MLGs corresponding to TR 34 /L98H allele were share by several isolates (Table 3). Two main MLGs (A and B in Figure 2), differing only by a single repeat at the STRAf 3A loci and corresponding to 17 and 12 isolates, respectively, were shared by isolates recovered from both clinical and environmental samples. In addition to being isolated both in the clinic and in the environment, these isolates were also widely found in different environments across our region and also from different categories of patients or date of sampling ( Figure 3 and Table 3).
Two highly dissimilar MLGs of eight and six isolates (MLG C and D, Figure 3) were also found in the same market garden, in six different soil samples coming from two different sites where difenoconazole was used to treat vegetable fields. These fields were found 1.5 km apart. In this market garden, we also recovered 14 different MLGs, with one MLG (MLG E, Figure 3) shared with one sawmill 40 km away (all green MLG in Figure 3, except MLG B and 1E073 MLG, which were found in another market garden). Genotyping Results for TR 34 /L98H Isolates: Global Scale A subset of the dominant MLGs in the Franche-Comtéarea were commonly shared with isolates found worldwide: MLG B (found in CF patients, hospital air, sawmills, and market gardens) was also detected in Australia and Ireland and MLG E (detected in market gardens and sawmills) was also isolated in Denmark (Figure 4). In comparison, MLG A, C, and D, while common in our dataset, were not observed in the global dataset, suggesting that they are geographically spatially constrained.
Other MLGs found in Eastern France were also closely related to those occurring more widely including those found in the Netherlands (links between the MLG on the Figure 4). One MLG was shared by four isolates in Besancon (a sawmill, twice in hospital air and once in a CF patient) and one isolate in the Netherlands. Another MLG is common with Ireland and Iran isolates and one is common with the Netherlands and another country (not known in the AfumID application). Nine MLG corresponding to nine isolates (five from market gardens, two from sawmills, one from hospital air, and one from CF patient) appear to be completely different (no link in Figure 4 with fixed cutoff) and were only found in the Besancon collection of isolates

DISCUSSION
In our study, we describe linked patient sampling, healthcare bioaerosol monitoring, and environmental sampling of A. fumigatus in order to explore potential links using molecular epidemiology. STRAf genotyping that we undertook demonstrated spatial heterogeneity, with some clones occurring far more widely than others. These data highlight the importance of incorporating spatial scale towards understanding of the local population genetic structure of A. fumigatus and the epidemiology ARAf. FIGURE 3 | Minimum spanning network for TR34/L98H isolates (nine microsatellite loci STRAf using Bruvo's distance, with distance cutoff < 0.2). Colors represent the origin of sampling and numbers are ID isolates. Link thickness is proportional to genotype similarity. The most common multilocus genotypes (MLG) are denoted by the letters A, B, C, D, and E.
Currently, the modes of dissemination, the occurrence of resistant genotypes in the environment, and the likelihood of exposure to at-risk patients remain poorly known. When applied locally, genotyping methods can shed light into the source of infection with various studies showing possible nosocomial acquisition of A. fumigatus and potentially ARAf (van der Torre et al., 2021). Three studies have suggested that dissemination from patient to patient occurs, as well as from patient to the healthcare environment (Pegues et al., 2002;Lemaire et al., 2018;Engel et al., 2019). For animals, the dynamics of resistant animal populations in livestock in relation to the environment (except for poultry) are poorly described and are unknown in wildlife. In both cases, it is likely that exposure to ARAf genotypes is higher for individuals in environment that are frequently exposed to fungicides and/or also to organic matter (leaf litter, compost, and cultivated soils rich in organic matter). Clearly, studies to assess the extent that resistant genotypes are stratified by region and environment will lend nuanced understanding to the ecological processes that are thought to generate ARAf.
In this study, as in others, we observe a reduced genetic diversity of TR 34 /L98H isolates (the most common resistance mechanism) compared to susceptible strains, suggesting a more recent emergence of this genotype . This allele is suggested to have emerged around 1997, perhaps in the Netherlands and dispersed across Europe, although other regions of origin are possible (Snelders et al., 2012;Abdolrasouli et al., 2015;Bader et al., 2015;Garcia-Rubio et al., 2018). Currently, TR 34 /L98H is increasingly described in other global settings and the occurrence of TR 34 /L98H genotypes does not seem to be limited to Europe. Aerial dispersal is most likely a major contributor to the long-distance migration of airborne A. fumigatus conidia owing to the widespread detection of aerosolized A. fumigatus spores (Shelton et al., 2020), which can also be moved in traded products (Dunne et al., 2017). FIGURE 4 | Minimum spanning network for TR 34 /L98H isolates (nine microsatellite loci STRAf using Bruvo's distance, with distance cutoff < 0.2) from Besancon mycology team and from the world (AfumID application). Colors represent the origin of sampling and numbers are ID isolates. Link thickness is proportional to genotype similarity. The most common multilocus genotypes (MLGs) found in Besancon laboratory are denoted by the letters A, B, C, D, and E.
Indeed, identical genotypes were described in several centers and countries, sometimes separated by thousands of kilometres (Chowdhary et al., 2012;Sewell et al., 2019). In 2012, a study of environmental isolates from Iran reported a TR 34 /L98H MLG closely related to Dutch and Indian isolates (Badali et al., 2013). A recent Japanese publication described TR 34 /L98H isolates clustering with strains from the Netherlands and France (Tsuchido et al., 2019).
In our study, one of the two dominant MLGs (MLG B) isolated in different environments in the Besancon region and isolated from different patients over several years was identical to an MLG in Australia and Ireland. Owing to the highly diverse nature of the STRAf barcoding loci, the likelihood by these MLGs being identical by chance rather than by descent is vanishingly unlikely (Sewell et al., 2019). Conversely, we have also found MLGs that have so far not been found in other countries. This is particularly the case for MLG A, which has high frequency in varied environments (sawmills and hospital air) as well as in different patients over several years. In Taiwan, a team has observed a potential global spread of TR 34 /L98H isolates but also a local genetic clustering of some clinical and environmental TR 34 /L98H (Wang et al., 2018). Other works have also reported local clusters, as in China with clinical and environmental TR 34 /L98H/S297T/F495I isolates (Deng et al., 2017;Chen et al., 2018). The two isolates carrying this mutation found in the French sawmill do not have the same STRAf profiles as those described in these works. In Germany, local factors have been evoked as the reason for the increased prevalence of specific genotypes (Bader et al., 2015). Resistance could therefore also result from the appearance of different clones located in different countries or regions when there is sufficient selection pressure. This could explain the greater diversity of genotypes found in vegetable crops in the market gardens, for example.
In our study, we have shown that ARAf occurs in the air of hospital and, in this case, could be acquired by patients during episodes of immunosuppression. In our hospital, resistant isolates were found in the general corridors, not in the hematology departments, thus reducing the risk of acquiring invasive aspergillosis caused by a resistant genotype. However, we described MLGs found in CF patients and in environments such as hospital air, and also with hospital visits coinciding with the isolation of certain resistant genotypes; this is the case, for example, for MLGs A and B. However, the wide distribution of these MLGs, both locally for MLG A (different environments) and globally for MLG B (different environments and in other countries), prevents us from conclusively asserting nosocomial acquisition of these isolates in CF patients. This is because, on a probabilistic basis, there is a non-negligible likelihood that these dominant genotypes may have been acquired through exposures outside of the hospital. However, a single MLG found only in Besancon was isolated from hospital air (1E023) and from sputum of a CF patient (1E017) 1.5 months apart. Even if we do not have the analytical power to definitively prove this, the potential for nosocomial contamination of the healthcare environment cannot currently be excluded.
Microsatellite genotyping is widely used in eukaryotes. In A. fumigatus, it has been used to identify the potential origins of resistant genotypes, and how genotypes are shared between patients, geographical regions, countries, and continents (Korfanty et al., 2019). Genotyping analyses with the ninemicrosatellite-based method in A. fumigatus (de Valk et al., 2005) provide high discriminatory power and reproducibility (Klaassen and Osherov, 2007). For our local isolate comparisons, we can be sure that MLG A and B of the Besancon region isolates are different, as they were analyzed with the same methodology (PCR, fragment analysis, and copy number transformation). However, when we compare the profiles of isolates analyzed in different laboratories, how can we be sure that two isolates with only one microsatellite difference, for example, are indeed different or that the difference is due to interlaboratory variations? MLG A (only local isolates) and MLG B (local and found in Australia and Ireland) only differ by one repeat on a microsatellite, STRAf 3A, which is the loci showing the greatest diversity, as has already been shown for agricultural strains (Ashu et al., 2018). It is therefore important to secure or find ways to harmonize methods of analysis on a large scale so that results can be reliably compared, and ensure in our case that it is indeed the MLG A that is local and that the MLG B has a global distribution. Interlaboratory controls should also be considered with, in particular, the implementation of standardization of sizing data by using allelic ladders (de Valk et al., 2009). Moreover, microsatellite genotyping does not provide sufficiently fine-grained information on the degree of proximity of two strains with two MLG, for example, whereas whole genome sequencing (WGS) can (Hagiwara et al., 2014;Rhodes et al., 2021). High-throughput WGS techniques are now poised to surpass other typing methods because of the high resolution of the data generated and the decreasing costs as performance increases. Better genetic knowledge of resistant strains through genomic approaches will help define the population structure of A. fumigatus and inform the future evolutionary trajectory of resistant phenotypes.
While WGS techniques have a higher analysis resolution than other genotyping techniques, the STRAf method, coupled with mutation screening and even mating type characterization, appears to be a good compromise to rapidly and simply screen relatively large numbers of isolates that have resulted from surveillance programs (van der Torre et al., 2021). Then, as in our case, where a "super clone" emerges from the analysis, it would be judicious to use WGS analysis to refine analyses in order to further explore patient-environmental links. Going forward, there is a need to continue to study the environmental conditions (both in host and in the environment) that facilitate the development, selection, and spread of resistance genotypes with large-scale multidisciplinary collaborations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
Conceptualization: SR and MF. Molecular analysis: SR, CG, and AL. Data analysis, SR, TS, and BV. Writing-original draft preparation: SR. Writing, review, and editing: SR, MF, TS, CG, BV, and LM. Funding acquisition: SR and LM. All authors contributed to the article and approved the submitted version.

FUNDING
The routine analyses were carried out within the framework of hospital activity, without external funding. The in-depth analysis of resistant strains was financially supported by LTSER "Zone Atelier Arc Jurassien", "Observatoire des Sciences de l'Univers Terre Homme Environnement Temps Astronomie", Bourgogne Franche-Comte-University, and RECOTOX network supported by AllEnvi and co-carried by INRAE and CNRS-INEE. MF and TS were supported by the UK Medical Research Council and the UK Natural Environmental Research Council.