Next Article in Journal
Endosymbiont Capture, a Repeated Process of Endosymbiont Transfer with Replacement in Trypanosomatids Angomonas spp.
Next Article in Special Issue
Transcription Factor HSF1 Suppresses the Expression of Surfactant Protein D in Cells Infected with Aspergillus fumigatus
Previous Article in Journal
Long-Term Persistence of Anti-SARS-CoV-2 Antibodies in a Pediatric Population
Previous Article in Special Issue
Aspergillus Genus and Its Various Human Superficial and Cutaneous Features
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genome-Wide Association Analysis for Triazole Resistance in Aspergillus fumigatus

Department of Biology, McMaster University, Hamilton, ON L8S 4K1, Canada
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Pathogens 2021, 10(6), 701; https://doi.org/10.3390/pathogens10060701
Submission received: 10 May 2021 / Revised: 24 May 2021 / Accepted: 2 June 2021 / Published: 4 June 2021
(This article belongs to the Special Issue Aspergillosis)

Abstract

:
Aspergillus fumigatus is a ubiquitous fungus and the main agent of aspergillosis, a common fungal infection in the immunocompromised population. Triazoles such as itraconazole and voriconazole are the common first-line drugs for treating aspergillosis. However, triazole resistance in A. fumigatus has been reported in an increasing number of countries. While most studies of triazole resistance have focused on mutations in the triazole target gene cyp51A, >70% of triazole-resistant strains in certain populations showed no mutations in cyp51A. To identify potential non-cyp51A mutations associated with triazole resistance in A. fumigatus, we analyzed the whole genome sequences and triazole susceptibilities of 195 strains from 12 countries. These strains belonged to three distinct clades. Our genome-wide association study (GWAS) identified a total of six missense mutations significantly associated with itraconazole resistance and 18 missense mutations with voriconazole resistance. In addition, to investigate itraconazole and pan-azole resistance, Fisher’s exact tests revealed 26 additional missense variants tightly linked to the top 20 SNPs obtained by GWAS, of which two were consistently associated with triazole resistance. The large number of novel mutations related to triazole resistance should help further investigations into their molecular mechanisms, their clinical importance, and the development of a comprehensive molecular diagnosis toolbox for triazole resistance in A. fumigatus.

1. Introduction

Aspergillus fumigatus is an opportunistic human fungal pathogen that is found in a broad range of substrates and is capable of surviving and growing in numerous environmental conditions. A. fumigatus is the primary cause of invasive aspergillosis, a life-threatening mold infection with high morbidity and mortality rates in immunocompromised patients. Its high sporulating capacity contributes to the environmental prevalence of A. fumigatus, leading to the high likelihood of infection in at-risk populations [1]. Globally, it is estimated that over 200,000 cases of invasive aspergillosis occur annually [2]. However, this number may represent only one-half of actual cases due to under- and mis-diagnoses [2]. Depending on factors such as population of patients, site of infection and antifungal management, mortality rates associated with invasive aspergillosis range from 60 to 90% [3].
Currently, there are four main classes of antifungals for aspergillosis treatment: azoles, polyenes, echinocandins, and allylamines. Among all antifungal agents, aspergillosis is commonly treated with triazole antifungals as the first choice because their use has been associated with better clinical response, less infusion-related toxicity, less nephrotoxicity and increased survival [4]. For aspergillosis treatment, itraconazole and voriconazole are among some of the commonly used triazole antifungals. Triazole antifungals work by inhibiting a vital enzymatic step in the synthesis of ergosterol, a major sterol and crucial part of the fungal cellular membrane [5]. Ergosterol plays a key role in membrane fluidity, membrane permeability, the activity of membrane proteins, and cell growth [5]. The triazoles work by inhibiting the demethylation of precursor sterols by binding to 14α-lanosterol demethylase (also known as Cyp51), a crucial enzyme involved in the ergosterol biosynthesis pathway. Triazoles act as competitive Cyp51 inhibitors through the binding of the N4 in their azole ring with the heme iron atom at the center of Cyp51 [5]. This binding prevents access of precursor sterols to the active site where demethylation occurs. Disruption of this enzymatic step causes significant damage to the cell membrane and results in the accumulation of toxic sterol intermediates, eventually leading to cell lysis and death [6]. However, the emergence of triazole-resistant A. fumigatus strains throughout the world has been a growing public health concern and a problem in the treatment of patients with aspergillosis.
Triazole-resistant strains have been extensively documented and characterized within multiple countries around the world. The majority of these studies have focused on the prevalence of resistant strains in a clinical setting. Furthermore, most analyses of the mechanisms of triazole resistance have focused on mutations in cyp51A, the gene coding for the triazole target enzyme [7]. The most common mutations in cyp51A among clinical-resistant strains, that develop during aspergillosis treatment, occur in amino acid sites G54, G138, M220, and G448 [7,8,9,10]. Meanwhile, the most common triazole drug resistant mutations in the global A. fumigatus population are TR34/L98H and TR46/Y121F/T289A, with many of these resistant strains originating from the environment [11,12].
The global population structure of A. fumigatus is shaped by high levels of gene flow between different populations [13,14]. Triazole-resistant A. fumigatus genotypes can rise and spread as a result of local selection due to elevated antifungal pressure within the environment. Clonal expansion of these highly fit triazole-resistant genotypes has been suggested to have led to their high abundance across the world [15,16]. Two main factors could have facilitated the spread of A. fumigatus genotypes and drug-resistant genes among geographic populations: the high dispersal ability of its asexual spores by wind and contemporary anthropogenic influences such as travel and trade [1]. Additionally, as aspergillosis is one of the leading causes of fungal deaths in avian species, bird migration may also be a factor in A. fumigatus dispersal [17,18]. The study by Ashu and colleagues further noted a large number of triazole-resistant genotypes and determined certain resistance genotypes were more commonly found in certain population genetic clusters than others [13]. Specifically, their analyses of 2026 A. fumigatus strains from 13 countries revealed that certain-resistant genotypes were mostly clustered into one genetic population and it was suggested that clonal expansion might have contributed to such a distribution.
Many research laboratories and hospitals around the world have been tracking the distribution of triazole-resistant clinical and environmental strains [19]. When examining prior epidemiological data, an increasing trend of triazole-resistant strains and infections has been reported. For example, within the Netherlands, the number of triazole-resistant infections has increased from 7.6% in 2013 to 14.7% in 2018 [20]. Another study in the Netherlands has also reported an increasing resistance rate in Radboud University Medical Center, from 0.79% between 1996 and 2001 to 7.04% between 2012 and 2016 [21]. These upward trends have also been reported in Iran (3.3% in 2013 to 6.6% in 2015), the United Kingdom (0.43% between 1998–2011 to 2.2% between 2015–2017), and in Texas, United States (7.2% between 1999–2002 to 22.6% between 2003–2015) [22,23,24]. The increasing prevalence of triazole-resistant A. fumigatus strains has become a major burden to many health institutions.
Triazole resistance in A. fumigatus is typically separated into two main categories, Cyp51A-mediated and non-Cyp51A-mediated mechanisms of resistance. In addition, studies on triazole resistance have mainly focused on three molecular mechanisms: (i) mutations in the Cyp51A protein, (ii) overexpression of the Cyp51A protein, and (iii) upregulation of drug efflux pumps. Alterations in Cyp51A are the most commonly studied mechanism for triazole resistance. Until 2008, all reported triazole resistance focused on the context of mutations in cyp51A. However, from 2008 onwards, the frequency of-resistant strains with no mutations in cyp51A was increasing [25,26]. At present, most epidemiological studies of triazole resistance in A. fumigatus only investigate mutations at the cyp51A gene. Consequently, mutations in other genes remain largely uncharacterized.
Microbial genome-wide association studies (GWAS) are a relatively new but powerful tool in understanding the relationships between genetic variations and microbial phenotypes. There have been several successful GWAS applications in identifying novel genomic markers responsible for antifungal drug resistance, with several studies focused on examining azole resistance in plant fungal pathogens [27,28]. In A. fumigatus, Zhao and colleagues recently conducted a genome-wide association study for itraconazole sensitivity in non-resistant clinical isolates from Japan [29]. In the current study, a GWAS was performed for itraconazole and voriconazole resistance in A. fumigatus based on a global population of published genomes. The aim of the study was to determine the genetic variants associated with triazole resistance using genome-wide SNP data, with the focus placed on novel non-cyp51A related mutations, as well as conduct a phylogenetic analysis with 195 strains to examine the phylogenetic distributions of itraconazole and voriconazole resistance in A. fumigatus.

2. Results

2.1. Phylogenetic Tree

Phylogenetic analysis of the whole-genome SNPs grouped the 195 samples and the reference strain Af293 into three large clades based on pairwise SNP differences between all 196 strains (Figure 1). Within each clade, whole-genome SNP differences were identified as ≤35,112 in Clade I, ≤45,160 in Clade II, and ≤48,670 in Clade III. Among the analyzed strains, 15 were in Clade I, 134 strains and the reference Af293 were in Clade II, and 46 strains were in Clade III. Geographically, the Clade II strains were from 10 countries, with 16 strains found in Canada, 4 in India, 1 in Ireland, 27 in Japan, 16 in Netherlands, 7 in Portugal, 1 in Singapore, 14 in Spain, 9 in the United Kingdom, and 37 in the United States. Furthermore, two strains were collected from the International Space Station. Clade III strains were obtained from the following seven countries: Canada (n = 1), India (n = 8), Netherlands (n = 8), Spain (n = 6), the United Kingdom (n = 16), Germany (n = 1), and the United States (n = 6). Finally, the 15 Clade I strains were from four countries: Canada (n = 1), Peru (n = 1), Portugal (n = 1), and Spain (n = 12). Within each clade, the samples were predominantly from infected patients, with 93.33% (14/15) in Clade I, 82.84% (111/134) in Clade II, and 82.61% (38/46) in Clade III of the analyzed samples were from clinical sources. The overall percentage of isolates from patients in the whole sample-set was 83.59%. According to the European Committee on Antimicrobial Susceptibility Testing (EUCAST), the MIC breakpoints for susceptible strains was set at ≤1 mg/L and the area of technical uncertainty (ATU) was set at 2 mg/L for both itraconazole and voriconazole antifungals. Among the samples with available MIC information and using an MIC ≥ 2 mg/L as the resistance cut-off value for both antifungals, 61.48% (75/122) were itraconazole resistant and 43.90% (54/123) were voriconazole resistant. At a clade-level, the percentages of itraconazole-resistant strains were as follows: 0% in Clade I, 57.75% in Clade II, and 82.93% in Clade III. For voriconazole, the percentage of resistant strains were 0% in Clade I, 37.50% in Clade II, and 65.85% in Clade III. Using MICs ≥ 4 mg/L as the resistance cut-off value for both antifungals, the itraconazole-resistant rate remained the same in our sample set, at 61.48% (75/122). However, the voriconazole-resistant rate dropped to 34.96% (43/123). The percentages of itraconazole-resistant strains in each clade remained the same while the voriconazole-resistant rates were as follows: 0% in Clade I, 27.78% in Clade II, and 56.10% in Clade III. Information on the 195 strains and clade divisions can be found in Supplementary Table S1.

2.2. Known Mutations Associated with Triazole Resistance

The MIC data for triazole drugs in this population identified 122 and 123 strains with known itraconazole and voriconazole MIC values respectively. We first examined the statistical association between mutations at 44 amino acid sites that had been previously found to be related to triazole resistance in A. fumigatus. The 44 known sites were mainly identified based on epidemiological surveys and are listed in Table 1.
Among these 44 known amino acid sites, 22 SNPs at 20 amino acid positions were found in our sample-set using the filtered vcf file, prior to multiallelic site removal (Table 1). Fisher’s Exact tests were conducted on these sites to determine their statistical associations with triazole resistance (Table 1). For these tests, using the 122 strains with known MIC values for both antifungals, we identified SNPs significantly associated with itraconazole and pan-azole resistance (Table 1).
According to EUCAST guidelines, MIC breakpoints for susceptible strains are set at ≤1 mg/L with an ATU of 2 mg/L for both itraconazole and voriconazole. To accommodate this buffer region, two resistance criteria were used and tested in this study. The first test defined resistant strains as having MIC values ≥ 2 mg/L and the second test set the resistance values at MIC ≥ 4 mg/L. A Bonferroni-corrected p-value threshold of 4.07 × 10−4 (0.05/122) was used to evaluate associations between the 22 SNPs and triazole resistance. Of the 22 known SNPs tested, only one in the Lysine-98 amino acid site, located in the gene cyp51A, was found to be significantly associated with itraconazole and pan-azole resistance in both Fisher’s Exact tests (Table 1).
We further sought to conduct Fisher’s Exact tests using subgroups consisting of solely itraconazole resistant (i.e., resistant to itraconazole but susceptible to voriconazole) or solely voriconazole resistant (i.e., resistant to voriconazole but susceptible to itraconazole) strains groups. However, the sample sizes of these subgroups were all below the requirement needed to achieve the desired Bonferroni-corrected p-value threshold. Thus, these subgroups were omitted from testing.
To unmask the effect of all these listed known mutation sites in cyp51A associated with triazole resistance, our study conducted a stepwise analysis of these sites using Fisher’s Exact tests. First, additional Fisher’s Exact tests were conducted after strains with the well-documented L98H mutation in cyp51A, which alone with its accompanying tandem repeat TR34 can confer triazole resistance, were removed. From the 122 strains with known MIC values, 21 strains contained the TR34/L98H mutation (Table S1). Using both MIC resistance thresholds and a Bonferroni-corrected threshold of 4.95 × 10−4 (0.05/101), the additional Fisher’s Exact tests identified no SNPs significantly associated with itraconazole and/or pan-azole resistance among these 22 known mutations.
To unmask the effect of other known cyp51A mutations associated with triazole resistance, additional Fisher’s Exact tests were also conducted after removal of strains containing any of these known mutations (Table S1). From the strains with known MIC values, 64 strains contained the known mutations in these cyp51A sites (Table S1). After removal of the 64 strains and using a Bonferroni-corrected threshold of 8.62 × 10−4 (0.05/58), the additional Fisher’s Exact tests identified no SNPs significantly associated with itraconazole and/or pan-azole resistance among these 22 mutation sites.
A final set of Fisher’s Exact tests were conducted focusing on a clade-level. Clade II was chosen for these additional analyses as the cluster contained the greatest number of strains and none of the Clade II strains contained the L98H mutation in cyp51A. The strains from Clade II with both itraconazole and voriconazole MIC values (n = 71) were used in this final set of the Fisher’s Exact tests. Using a Bonferroni-corrected threshold of 7.04 × 10−4 (0.05/71), no SNPs were found to be significantly associated with itraconazole and/or pan-azole resistance from these 22 mutations sites.
Together, the stepwise analyses results revealed that these well-characterized mutation sites do not account for the observed triazole resistance in our sample sets. Therefore, additional modes of action and uncharacterized novel mutations should be investigated for their possible involvement with triazole susceptibility in A. fumigatus.

2.3. Genes Overexpressed with Triazole Exposure

We further examined the potential overlap between the genome-wide population level SNPs identified here with previously identified genes not listed in Table 1 but were related with triazole resistance in A. fumigatus. Specifically, we extracted information about specific genes that were overexpressed in A. fumigatus during exposure to itraconazole and/or voriconazole. Table 2 summarizes the genes that were overexpressed upon exposure to each antifungal. The overexpression of these genes under triazole stress were determined using RT-qPCR and RNA-seq information [25,38,39]. Supplementary Table S2 describes the details on the experimental conditions and setup associated with each gene listed in Table 2. Specifically, previous work demonstrated that ten ATP-binding cassette (ABC) transporters (abcA-1, abcA-2, abcB, abcC, abcD, abcE, atrF, mdr1, mdr4, and AFUA_5G02260), four major facilitator superfamily (MFS) transporters (AFUA_2G11580, mfs56, mfsA and mfsC), the 14-alpha sterol demethylase cyp51A, and 16 transcription factors (ace1, AFUA_1G02870, AFUA_1G04140, AFUA_1G16460, AFUA_2G01190, AFUA_3G09130, AFUA_4G06170, AFUA_4G13600, AFUA_5G02655, AFUA_5G06350, AFUA_5G07510, AFUA_6G01960, AFUA_6G03430, AFUA_7G03910, AFUA_8G07360, and fumR) were overexpressed following itraconazole exposure [21,35]. Similarly, five ABC transporters (mdr1, abcB, abcC, abcD and abcE), three MFS multidrug transporters (mfsA, mfsB and mfsC), a F-box domain protein (fbpA), an AAA-family ATPase (aaaA), a C6 zinc finger domain protein (finA), a BZIP transcription factor (cpcA), and a putative C2H2 zinc-finger transcription factor (zfpA) were overexpressed with voriconazole exposure [38].
A summary of the overexpressed genes and specific fold changes, revealed in previous studies by RT-qPCR and RNA-seq during triazole exposure, are detailed in Table 2. Using these studies, a total of 37 overexpressed genes with triazole exposure were identified for further investigation. However, subsequent analysis excluded cyp51A-related mutations as they have already been extensively searched and discussed in Section 2.2.
We identified SNPs in these 37 overexpressed genes and their neighbouring intergenic regions using the soft-filtered vcf file. A total of 3230 SNP sites were identified in these overexpressed genes from our dataset. Using the same procedure as that used for the 22 known mutation sites, we identified SNPs significantly associated with itraconazole-resistance and pan-azole resistance in these 36 overexpressed genes. Multiple Fisher’s Exact tests, with a Bonferroni-corrected threshold of 1.55 × 10−5 (0.05/3230), were conducted on these sites.
Using a MIC threshold of 2 mg/L and all 122 strains, we found 57 SNPs and 11 SNPs to be significantly associated with itraconazole and pan-azole resistance, respectively (Table S3). For itraconazole, these SNPs were located in or beside 14 genes: abcC (n = 9), abcD (n = 1), abcE (n = 8), fbpA (n = 1), fumR (n = 1), mfsA (n = 5), mfsB (n = 1), mfsC (n = 1), AFUA_1G16460 (n = 6), AFUA_2G01190 (n = 1), AFUA_4G13600 (n = 3), AFUA_5G02655 (n = 5), AFUA_6G01960 (n = 9), and AFUA_6G03430 (n = 6). Among these 57 SNPs, 46 were found in intergenic or intronic regions, two were non-coding transcript variants, eight were synonymous variants, and one was a missense variant (Table S3). For pan-azole resistance, the 11 associated SNPs were located in or beside six genes: mfsA (n = 1), mfsB (n = 1), AFUA_2G01190 (n = 1), AFUA_4G06170 (n = 1), AFUA_4G13600 (n = 3), and AFUA_6G03430 (n = 4). The 11 SNPs comprised of 10 intergenic variants and one missense variant. Next, using the MIC threshold of 4 mg/L as the resistance cut-off, we found 57 SNPs and 10 SNPs to be significantly associated with itraconazole and pan-azole resistance, respectively (Table S3). When compared to the previous results obtained using a MIC threshold of 2 mg/L, two variants were no longer significantly associated with pan-azole resistance: a missense variant in mfsA and an intergenic variant in mfsB. Furthermore, a newly found synonymous variant in AFUA_1G04140 was significantly associated with pan-azole resistance using the MIC threshold of 4 mg/L (Table S3).
Fisher’s Exact tests were also conducted after removal of the 21 strains containing the L98H mutation in cyp51A. Using both MIC resistance thresholds of 2 mg/L and 4 mg/L, three SNPs were found to be significantly associated with itraconazole resistance (Table S3). The three SNPs consisted of the previously found intergenic variant in mfsB, and two novel intergenic variants—one in AFUA_6G01960 and the second in AFUA_6G01960 (Table S3). No SNPs were found to be significantly associated with pan-azole resistance.
A third set of Fisher’s Exact tests were conducted after removal of the 64 strains containing the mutations in cyp51A and using both MIC thresholds. Using the MIC resistance thresholds of 4 mg/L, one SNP was found to be significantly associated with pan-azole resistance. This SNP was found in the intergenic region of abcA (Table S3).
A final set of Fisher’s Exact tests was completed and focused solely on strains from Clade II (n = 71). Using both MIC resistance thresholds, 2 mg/L and 4 mg/L, no SNPs were found to be significantly associated with itraconazole and/or pan-azole resistance in this sample set.

2.4. Genome-Wide Association Study

In addition to examining known triazole resistance mutations and SNPs in genes overexpressed during triazole exposure, a genome-wide association study (GWAS) was performed on the 122 and 123 strains with known itraconazole and voriconazole MIC values to investigate potential novel mutations associated with triazole sensitivity. The results of our analyses are summarized in Figure 2. Specifically, the itraconazole GWAS Manhattan plot can be found in Figure 2A and for voriconazole, in Figure 2B. The generated quantile–quantile plots for the GWAS results displayed no systematic inflation in our samples (Figure S1A,B).
We further examined the top 20 significant SNPs identified by the GWAS analysis. Among the 20 SNPs obtained from the itraconazole GWAS, 13 (65%) were located in intergenic regions and 7 (35%) within protein-coding regions (Table 3). These seven SNPs consisted of five missense variants, one synonymous variant, and one non-coding transcript variant (Table 3). In terms of the top 20 SNPs found from the voriconazole GWAS, 10 (50%) were found in intergenic regions and the remaining 10 in coding regions (Table 4). These 10 coding-region SNPs consist of four missense variants, five synonymous variants, and one non-coding transcript variant (Table 4). Among the top 20 SNPs associated with each of the two drugs, only one was shared. This variant was a missense A to C mutation at the position 2,538,614 on chromosome 1, in the gene AFUA_1G09780. The remaining 38 SNPs were unique to each of the two triazole drugs.
Additional GWAS analyses were conducted in the same stepwise manner seen in the previous Fisher’s Exact tests. Firstly, to alleviate any potential masking effect caused by the L98H mutation in cyp51A, the 21 strains with the L98H mutation were removed. A second GWAS, using the same previous pipeline, was then conducted. The results of the second GWAS are summarized in Figure 3A,B as Manhattan plots for itraconazole and voriconazole, respectively. The generated quantile–quantile plots for both GWAS results displayed no genomic inflation (Figure S2A,B).
The top 20 significant SNPs identified by the second GWAS analyses were examined. Among the 20 SNPs obtained from the itraconazole GWAS, 13 (65%) were located in intergenic regions and 7 (35%) within protein-coding regions (Table 5). These seven SNPs comprised of four missense variants, one synonymous variant, and two non-coding transcript variants (Table 5). In terms of the top 20 SNPs obtained from the second voriconazole GWAS, 10 (50%) were found in intergenic regions and the remaining 10 in coding regions (Table 6). These 10 coding-region SNPs consist of 5 missense variants, 3 synonymous variants, and 2 non-coding transcript variants (Table 6). Among the top 20 SNPs associated with each of the two drugs, none were shared between the two triazole drugs.
A third set of GWAS analyses was also done to alleviate any potential masking effect caused by the known mutations in cyp51A, previously listed in Table 1. The 64 strains with cyp51A mutations were removed and the third GWAS, using the same previous pipeline, was then conducted. The results of the third GWAS are summarized in Figure 4A,B as Manhattan plots for itraconazole and voriconazole, respectively. The generated quantile–quantile plots for both GWAS results displayed no genomic inflation (Figure S3A,B).
The top 20 significant SNPs identified by the third GWAS analyses were examined. Among the top 20 SNPs obtained from the itraconazole GWAS, 11 (55%) were located in intergenic regions and 9 (45%) within protein-coding regions (Table 7). These nine SNPs comprised of three missense variants, two synonymous variants, three non-coding transcript variants and one intragenic variant (Table 7). In terms of the top 20 SNPs obtained from the third voriconazole GWAS, 10 (50%) were found in intergenic regions and the remaining 10 in coding regions (Table 8). These 10 coding-region SNPs consist of six missense variants and four synonymous variants (Table 8). Among the top 20 SNPs associated with each of the two drugs, two SNPS were shared between the two triazole drugs. The first variant was a synonymous C to A mutation at the position 2,539,714 on chromosome 4, in the gene AFUA_4G09770. The second mutation was a synonymous T to C mutation at position 2,131,740 of chromosome 5, in the gene AFUA_5G08390.
A final set of GWAS was completed to focus our analysis on a clade-level, using strains from Clade II. The strains from Clade II with itraconazole (n = 71) and voriconazole (n = 72) MIC values were used for the fourth GWAS, using the same previous pipelines. The results of this GWAS are summarized in Figure 5A,B as Manhattan plots for itraconazole and voriconazole, respectively. The generated quantile–quantile plots for both GWAS results displayed no genomic inflation (Figure S4A,B).
The top 20 significant SNPs identified by the GWAS analyses on strains from Clade II were examined. Among the top 20 SNPs obtained from the itraconazole GWAS, 15 (75%) were located in intergenic regions and 5 (25%) within protein-coding regions (Table 9). These five SNPs comprised of two missense variants, two synonymous variants, and one non-coding transcript variant (Table 9). In terms of the top 20 SNPs obtained from the voriconazole GWAS, 6 (30%) were found in intergenic regions and the remaining 14 in coding regions (Table 10). These 14 coding-region SNPs consist of seven missense variants, five synonymous variants, one non-coding transcript variant and one intragenic variant (Table 10). Among the top 20 SNPs associated with each of the two drugs, no mutation sites were shared between the two triazole drugs.

2.5. Linkage Disequilibrium Analysis

Linkage disequilibrium analyses were conducted using the top 20 SNPs obtained by the four GWAS analyses and all 314,999 SNPs in the soft-filtered vcf file to search for SNPs highly linked (R2 > 0.85) to these significantly associated SNPs. Specifically, we focused on highly linked non-synonymous mutations. The results of this association analysis are presented in Table 11 for itraconazole and in Table 12 for voriconazole. In total, for itraconazole resistance, we identified 15 additional highly linked missense variants located in 13 (putative) protein-coding genes (Table 11). For voriconazole resistance, this analysis revealed 11 additional missense SNPs located in 11 different (putative) protein coding genes (Table 12). None of these additional missense SNPs were shared between the two drugs.
Fisher’s Exact tests, with a Bonferroni-corrected p-value threshold of 4.07 × 10−4 (0.05/122), were conducted to examine associations among these highly linked mutations to itraconazole and pan-azole resistance (Table 13). MIC resistance thresholds of 2 mg/L and 4 mg/L were both tested for these 26 sites and using all 122 strains. Both MIC thresholds identified four SNPs to be significantly associated with itraconazole resistance as well as two of these SNPs also being associated with pan-azole resistance (Table 13).
Additional Fisher’s Exact tests were also conducted after the removal of the 21 strains with the L98H mutation in cyp51A and using a Bonferroni-corrected p-value threshold of 4.95 × 10−4 (0.05/101) (Table 14). For both MIC resistance thresholds, the results showed that the three previously noted SNPs, in AFUA_1G17380, AFUA_3G09040 and AFUA_3G09070, were again significantly associated with itraconazole resistance. After removal of the 21 strains and using both MIC thresholds, two of these SNPs, AFUA_3G09040 and AFUA_3G09070, were now also significantly associated with pan-azole resistance. Another shared SNP with the previous analysis is the missense variant in AFUA_1G03370, which was found to be significantly associated with pan-azole resistance at both MIC resistance thresholds. Furthermore, using the MIC threshold of 2 mg/L, a novel missense variant in AFUA_7G06290 was found to be associated with pan-azole resistance (Table 14).
A third set of Fisher’s exact tests were conducted after removal of the 64 strains containing known cyp51A mutations and using a Bonferroni-corrected p-value threshold of 8.62 × 10−4 (0.05/58) (Table 15). For both MIC resistance thresholds, the tests determined three previously identified SNPs to be significantly associated with both itraconazole and pan-azole resistance. These three SNPs were a missense variant in AFUA_1G17380, AFUA_3G09040, and AFUA_3G09070. Using both MIC thresholds, the tests also identified the previous AFUA_7G06290 missense variant to be significantly associated with pan-azole resistance.
Lastly, another set of Fisher’s Exact test was conducted to focus solely on strains from Clade II and used a Bonferroni-corrected threshold of 7.04 × 10−4 (0.05/71) (Table 16). For both MIC resistance thresholds, the tests determined three previously identified SNPs to be significantly associated with both itraconazole and pan-azole resistance. These SNPs were a missense variant in AFUA_1G17380, AFUA_3G09040 and AFUA_3G09070. Furthermore, using both MIC thresholds, the tests also identified the previously noted missense variant in AFUA_7G06290 to be significantly associated with pan−azole resistance (Table 16).

3. Discussion

In this study, we analyzed the genomic polymorphisms among 195 A. fumigatus isolates collected from 12 countries as well as the International Space Station to investigate the potential associations between genomic SNPs and triazole resistance. Phylogenetic analyses of the whole-genome SNPs identified three main clades in this sample, with Clade I being very divergent from the other two clades. Most strains in this clade were from Spain and they likely represent a cryptic species within A. fumigatus sensu stricto. Among these 195 strains, the minimum inhibitory concentrations of two triazoles, itraconazole and voriconazole, were reported for 122 and 123 strains, respectively. Over the past two decades, an increasing number of studies have been conducted to investigate the genetic diversity and population structure of A. fumigatus using different molecular markers [14,40,41,42,43]. A previous study exploring global population genetic variation by Ashu et al. identified 8 genetic clusters by examining nine short tandem repeats in 2026 A. fumigatus isolates from 13 countries [13]. However, a more recent study analyzing the same short tandem repeats of 4049 A. fumigatus isolates identified two broad genetic clusters [14]. The whole-genome SNP analyses here revealed three divergent clades and within both Clades II and III, several sub-clades with significant bootstrap supports were also found. Therefore, the true number and composition of the genetic clusters in the global A. fumigatus population remain uncertain and depend on how clades and genetic clusters are defined. However, based on previous studies, most genetic clusters and clades contain geographically and ecologically diverse strains, consistent with frequent gene flow and great adaptability of A. fumigatus genotypes [14,26].
Among geographic and ecological populations, different frequencies of triazole resistance have been reported, likely reflecting their variations in strain source, clinical antifungal usage, agricultural fungicide usage, and surveillance techniques [44,45,46,47,48,49,50,51]. In 2017, Garcia-Rubio et al. reviewed previously published literature and reported that the global triazole-resistant rate ranged from 0.55% to 30% [30]. In the samples analyzed here and using an MIC threshold of 2 mg/L, 61.48% of all isolates with available MIC data were itraconazole resistant and 43.90% were voriconazole resistant. Furthermore, 63.46% and 43.81% of the clinical isolates were itraconazole and voriconazole resistant, respectively. Similarly, there was a high frequency of environmental isolates resistant to itraconazole and voriconazole, at 50.00% and 44.44%, respectively. Using the MIC threshold of 4 mg/L, resistance frequencies for itraconazole remained the same, however, these values changed for voriconazole. The resistance rate for voriconazole in clinical strains decreased to 35.24% and for environmental strains, it changed to 33.33%. The high rates of resistance among strains analyzed here could be attributed to the biases among research groups in preferentially submitting drug-resistant strains for whole-genome sequencing. However, the broad range of triazole MIC values among the large number of sequenced strains allowed us to infer potential novel genetic variants not identified in previous studies.
Indeed, in this study, GWAS for both itraconazole and voriconazole resistance identified novel genes and mutations linked to triazole resistance. We conducted four GWAS analyses. The first analysis used all strains with known MIC values for itraconazole (n = 122) and voriconazole (n = 123). Meanwhile, the following two analyses investigated novel mutations associated with itraconazole and voriconazole resistance at other SNP loci, unrelated to cyp51A. Specifically, the second GWAS removed the 21 strains containing the L98H mutation in cyp51A and the third GWAS removed the 64 strains containing known mutations in cyp51A related to triazole resistance. The last GWAS was done on a clade-level, focusing the analysis on strains from Clade II with known triazole MIC values (n = 71). For each GWAS, we focused our investigation on the top 20 SNPs obtained via the GWAS analyses for each of the two drugs. We identified a total of six missense variants to be putatively associated with itraconazole resistance. These six missense variants were located in six genes: AFUA_1G09780, AFUA_2G08060, AFUA_3G09090, AFUA_5G08150, AFUA_6G01860, and AFUA_8G02350. The first mutation was in AFUA_1G09780, which encodes for a stomatin family protein. The stomatin proteins belong to a highly conserved family of integral membrane proteins. In humans, stomatin interacts with various ion channels and modulates their activity. The proteins are also thought to perform specific scaffolding functions in membranes. However, the functional information of stomatin proteins in fungi is scarce. The ortholog of AFUA_1G09780 in the closely related Aspergillus nidulans species is AN1287 (stoB). The AN1287 protein is located in the inner mitochondrial membrane [52]. In addition to targeting ergosterol biosynthesis, triazole exposure also causes production of deleterious mitochondrial reactive oxygen species (ROS). Therefore, since triazoles promote ROS accumulation, the mitochondrial membrane complexes represent another group of targets when studying resistance [53]. Of note, this mutation in AFUA_1G09780, was found to be significantly associated with itraconazole resistance in all four GWAS analyses. The second mutation was in AFUA_5G08150, which encodes for a putative ABC bile acid transporter. Although this specific gene has not been previously linked to triazole resistance, other ABC transporter members are known to modulate triazole extrusion. Previous studies have found multiple ABC transporter members to be overexpressed with triazole exposure and that expression levels of various ABC transporter members were higher among triazole-resistant isolates [23,37,54]. In most cases, overexpression in this family of transporters seems to prevent intracellular drug concentrations from reaching levels needed to be effective at impacting ergosterol biosynthesis. The next variant was found in AFUA_8G02350, encoding a putative polyketide synthase. The polyketide synthases are involved in the biosynthesis of polyketides, which are a large and structurally diverse group of secondary metabolites [55]. These compounds have a wide range of biological activities that are important for ecological and evolutionary adaptation in fungi [55]. Furthermore, the gene AFUA_8G02350 is within a terpene hybrid cluster [56]. The fourth gene, AFUA_3G09090, encodes for an uncharacterized protein. The fifth missense mutation related to itraconazole resistance was found in AFUA_2G08060, which encodes an involucrin repeat protein. This protein has been found to be involved in tethering Woronin bodies to septal pores [57]. In multicellular fungi such as A. fumigatus, cells are connected to each other via intercellular bridges called septal pores and Woronin bodies plug these pores upon injury to avoid excessive loss of cell content. Deletion mutants with impaired Woronin bodies have shown to have impaired stress resistance and delayed hyphal wounding response [57]. The missense mutation identified here may enhance the strains’ ability to respond quickly to triazole drugs. The last missense mutation was in AFUA_6G01860, which encodes a putative MFS lactose permease. This class of protein plays a role in transmembrane transport and is an MFS family member, a superfamily of transport proteins that are mainly responsible for antifungal resistance development through drug efflux activity.
Interestingly, our itraconazole GWAS analysis results differ from those in a GWAS conducted by Zhao and colleagues who examined SNPs associated with itraconazole sensitivity [29]. They completed a GWAS using 76 clinical A. fumigatus isolates collected from Japan. Our comparisons revealed no overlap between our top 20 SNPs and the SNPs they found to be highly associated with itraconazole sensitivity. Two factors might have contributed to the different observations. In the first, the study by Zhao et al. focused on itraconazole sensitivity in non-resistant clinical isolates of A. fumigatus, with itraconazole MIC ranging from 0.125 to 1 mg/L among their 76 isolates [29]. In contrast, the itraconazole MICs for our 122 strains with itraconazole MICs ranged from 0.13 to 32 mg/L. Secondly, all the strains analyzed by Zhao et al. were from one country, Japan [29]. In contrast, our itraconazole GWAS included 122 strains from eight countries, Canada (n = 12), India (n = 12), Japan (n = 8), Netherlands (n = 21), Spain (n = 19), Germany (n = 1), the United Kingdom (n = 24), and the United States (n = 25). Together, these results suggest that additional novel SNPs associated with itraconazole sensitivity or resistance will likely be present in other geographic and/or ecological populations of A. fumigatus. Another recent preprint by Rhodes and colleagues also conducted a GWAS on itraconazole resistance using treeWAS, a phylogenetic tree-based GWAS approach [58]. A comparison between the top 20 SNPs from our GWAS analyses and their significantly associated SNPs to itraconazole resistance found no overlap in SNP sites. This difference in results is most likely related to sample selection, as their sample set focused on strains obtained from the United Kingdom and Republic of Ireland. Another potential reason for the discrepancy is that our study used a quantitative phenotype, based on MIC values, for our GWAS while Rhodes and colleagues used a binary phenotype, separating strains into susceptible and resistant strains by defining resistance as MIC ≥ 2 mg/L. However, their results are consistent with our general conclusion that a large number of additional novel mutations, un-related to cyp51A, are significantly associated with triazole resistance in A. fumigatus.
The voriconazole GWAS analyses identified a total of 17 missense variants to be putatively associated with resistance. These variants were found in 17 different genes: AFUA_1G03370, AFUA_1G09780, AFUA_1G12540, AFUA_1G17410, AFUA_2G01700, AFUA_2G13030, AFUA_4G09580, AFUA_5G01000, AFUA_5G02210, AFUA_5G14610, AFUA_6G09870, AFUA_7G00740, AFUA_7G05960, AFUA_8G01250, AFUA_8G01340, AFUA_8G01940, and AFUA_8G02280. Two of these missense variants were found in AFUA_1G03370 and AFUA_5G02210, which encodes for putative proteins of unknown functions. A third gene, AFUA_1G12540, encodes a putative TMEM1 family protein with uncharacterized function. The next seven variants were found in members of enzyme families with roles across a large number of biological processes, which comprised of the genes AFUA_1G17410 that encodes a putative beta-glucosidase, AFUA_2G01700 that encodes a putative serine/threonine protein kinase, AFUA_2G13030 that encodes a phenylalanyl-tRNA synthetase, AFUA_5G01000 that encodes a putative oxidoreductase of the 2-oxoglutarate (2OG)-Fe(II) oxygenase superfamily, AFUA_5G14610 that encodes a putative carboxypeptidase Y, AFUA_7G00740 that encodes a putative protein kinase, and AFUA_8G01250 that encodes a putative GNAT family acetyltransferase. The next variant was found in AFUA_7G05960, which encodes a putative C2H2 finger domain protein. Three significantly associated missense variants also encoded for putative C6 zinc cluster transcription factors, which were found in AFUA_6G09870, AFUA_8G01940 and AFUA_8G02280. Other members of this transcription factor family have been linked to triazole resistance. For example, a previous transcriptome study had found finA, a C6 zinc finger domain protein, displayed increased mRNA levels during adaptation to voriconazole exposure [38]. In addition, another C6 zinc-cluster transcription factor, AtrR, had been found to be associated with triazole resistance by regulating expression of genes related to ergosterol biosynthesis [59]. However, the Zn cluster family is the largest family of transcription factors known in eukaryotes and thus additional testing is required [60]. Interestingly, a missense mutation in AFUA_1G09780 was also found significantly associated to voriconazole and was the same SNP found to be associated with itraconazole resistance. The remaining two genes were AFUA_8G01340 that encodes a putative MFS sugar transporter and AFUA_4G09580, which encodes the major allergen Aspf2. Although our study focused on examining missense variants, significantly associated SNPs obtained by GWAS also included synonymous, intergenic and intronic variants. These variants can have biological consequences and contribute to functional changes in the protein. For example, synonymous mutations can affect critical cis-regulatory sequences, alter mRNA structure, and impact translational speed [61]. Furthermore, non-coding variants can be found within potential regulatory sequences such as enhancers, promoters, and 5′ and 3′ UTRs [62]. Through these regulatory roles, non-coding variants can influence processes such as transcription, translation and splicing.
The results of the itraconazole and voriconazole GWAS showed few overlaps between significant SNPs; the first GWAS having one SNP overlap in the gene AFUA_1G09780 and the third GWAS having two shared SNPs, a synonymous mutation in AFUA_4G09770 and a synonymous mutation in AFUA_5G08390. The remaining two GWAS found no overlapping SNPs between the two antifungal drugs. Several reasons could have contributed to the low number of shared SNPs. Although all azoles operate using the same common mode of action, decreasing ergosterol synthesis by inhibiting the fungal enzyme 14α-sterol demethylase, there are differences between itraconazole and voriconazole in terms of their mechanisms of action. Voriconazole also inhibits 24-methylene dihydrolanasterol demethylation in Aspergillus and its antifungal activity is likely a result of a combination of effects in addition to inhibition of ergosterol synthesis [63]. Secondly, our sample set consisted of a large number of strains that had different susceptibilities to the two drugs and were only resistant to one of the two antifungals. Finally, our sample set was not a natural randomly mating population but were selected strains sequenced by different laboratories based on their own specific objective and purpose. As a result of the diversity of strains and their originating populations, some of the shared azole-resistance related mutations may have been less frequent in our studied sample set and were, thus, likely filtered out during quality control.
Linkage disequilibrium was also evaluated using the top 20 SNPs obtained by each of the four GWAS analyses to identify additional highly linked missense SNPs. Four sets of Fisher’s Exact tests were conducted: using all 122 strains, after removal of the 21 strains with the L98H mutation, after removal of the 64 strains with the cyp51A mutations, and using only the 71 Clade II strains. Interestingly, the last three tests identified two SNPs, missense variants in AFUA_3G09040 and AFUA_3G09070, to be significantly associated with itraconazole as well as pan-azole resistance using both MIC thresholds (Table 14, Table 15 and Table 16). In the first test, conducted using all 122 strains and at both MIC thresholds, these two missense variants were also found to be significantly associated with itraconazole resistance (Table 13). Another SNP in AFUA_1G17380 was also found in three tests to be significantly associated with itraconazole and pan-azole resistance using both MIC thresholds (Table 13, Table 15, and Table 16). Furthermore, in the remaining test, the SNP site was significantly associated with itraconazole resistance (Table 14). In terms of the function of these three genes, AFUA_1G17380 encodes a putative 3-oxoacyl-(acyl-carrier-protein) reductase, AFUA_3G0904 encodes for an uncharacterized protein and AFUA_3G09070 encodes a putative carboxylesterase.
In addition, our study examined 20 previously known amino acid sites associated with triazole resistance. Fifteen of these known amino acid sites were in the cyp51A gene. Of these 15 sites, several have been functionally validated in previous studies and found to directly contribute to triazole resistance. These validated mutation sites consisted of G54, L98, Y121, G138, M220, T289 and G448 [8,64]. However, hot-spot SNP sites that confer triazole resistance by itself only include five of the seven sites, at G54, Y121, G138, M220 and G448. Meanwhile, for L98 and T289, a combination with a tandem repeat is required for triazole resistance, specifically TR34/L98H and TR46/Y121F/T289A [65]. Mutations in the G138 site has also been validated to cause multi-azole resistance [66]. Lastly, SNP sites P216 and F219 both confer resistance to itraconazole when mutated [67]. Specifically, there have been two amino acid substitutions in F219 in triazole-resistant strains, F219I and F219L [67]. However, interestingly, we have instead identified a different substitution F219S in our sample set; although it was not significantly associated with resistance (Table 1). We have also included the amino acid sites F46, M172, N248, D255, and E427; although only 3 of these sites could be found in our filtered genotype file, specifically F46, D255, and E427. Strains with a combination of these five mutations, as F46Y/M172V/D255E or F46Y/M172V/D255E/N248T/E427K, have shown higher triazole MICs than the wild-type strains [30]. The three remaining cyp51A mutation sites, H147, S297 and F495, have been identified to be associated with resistant isolates but have not been functionally validated. Mutations in H147 were found to coincide with isolates with G448 mutations and were not found to be associated with resistance by itself and is thought to only increase protein stability [9,68]. Similarly, S297 and F495 mutations are found in some TR34/L98H mutant A. fumigatus strains but have not been proven to be sufficient for triazole resistance by themselves [12,69]. However, these two amino acid positions, S297 and F495, are located near the triazole binding pocket of Cyp51A [70]. The remaining five known amino acid sites associated with triazole resistance are non-cyp51A mutations. Three of these SNP sites (I412, P309, and S305) are in hmg1, which encodes a 3-hydroxy-3-methyl-glutaryl-coenzyme A (HMG-CoA) reductase. Mutation at these three sites have been identified in triazole-resistant isolates and their association with triazole resistance has been functionally validated by inserting each SNP into the hmg1 gene of the laboratory strain akuBKU80 [33]. These three sites, I412, P309, and S305, are predicted to be within the conserved sterol-sensing domain of hmg1 [33]. Furthermore, akuBKU80 mutants with the substitution I412S and S305P possessed significantly different cellular sterol profiles compared to the unaltered akuBKU80 strain; independent of cyp51A and cyp51B expression [33]. Another known mutation site is L167 in the uncharacterized AFUA_7G01960 gene. The nonsense mutation L167* was first identified in the multi-azole-resistant clinical isolate V157-62, and subsequently functionally validated by inserting the specific SNP into an azole-susceptible clinical isolate, V130-15 [36]. At this amino acid position, a nonsense mutation was generated, which increased resistance to itraconazole [36]. Furthermore, this SNP was also associated with decreased ergosterol in the fungal membrane [36]. Interestingly, overexpression of the AFUA_7G01960 gene itself has also been correlated with increased voriconazole resistance [38]. Taken together with bioinformatic analysis, AFUA_7G01960 is predicted to be a putative transcription factor involved in ergosterol biosynthesis and mutation at L167 likely prevents its activity, thus leading to increased resistance [36]. The last known SNP site we investigated was E180 in AFUA_2G10600, a gene encoding the mitochondrial 29.9 KD NADH oxidoreductase subunit of respiratory complex I. The amino acid substitution E180D is present in itraconazole-resistant clinical isolates of A. fumigatus [37]. Furthermore, restriction enzyme-mediated insertion of this mutation in the AFUA_2G10600 gene led to increased itraconazole resistance [71]. This insertion was at a Xhol site, 534 bp from the start codon of the gene. This increase in resistance indicates that intact AFUA_2G10600 may confer azole susceptibility through mitochondrial NADH metabolism or NAD/NADH redox stress [71]. Complete deletion of the coding region for this 29.9KD subunit was also found to result in itraconazole resistance in the laboratory strain A1163 KU80, increasing from an MIC of 0.25 mg/L to >8 mg/L, which further supports its contribution to itraconazole resistance [37].
Using a Fisher’s Exact test on these 22 mutations, with all 122 strains and using both MIC resistance thresholds (2 mg/L and 4 mg/L), we found only one mutation, L98H in cyp51A, to be highly associated with itraconazole and pan-azole resistance (Table 1). Examining all samples with known triazole MIC data, this L98H mutation was found in 21 strains. We also determined using coverage data across the promoter region of cyp51A that all 21 strains with the L98H mutation were accompanied with the common 34-bp tandem repeat (Figure S5). A subsequent Fisher’s Exact test was, thus, done after removing strains with the L98H mutation in cyp51A (n = 21). Additional Fisher’s exact tests were also conducted after removal of all strains containing the cyp51A mutations (n = 64) and conducted again with only Clade II strains (n = 71). However, these additional tests identified no SNPs significantly associated with itraconazole and/or pan-azole resistance. A potential reason why the previously functionally validated sites were not highly associated with triazole resistance in our study is that strain counts for mutation genotypes at these sites were low in our 122-strain sample set and only ranged between 1 and 6 strains, making them unable to meet our criteria (>5% frequency in the population) for inclusion (Table S4).
We further examined the distribution of mutation phenotypes in these functionally validated sites in our three clades, using all 195 strains. Interestingly, for cyp51A, the mutation L98H was only present in strains of Clade III (n = 22). For the T289 site, three strains had the mutation T289A and this mutation always accompanied with the substitution Y121F. Furthermore, these three strains were all from Clade III. Mutations in the site G54, specifically G54V (n = 5), G54E (n = 1), G54W (n = 3) and G54R (n = 2), were found mostly in Clade II strains with a rate of 90.91% (10/11). Mutations M220I (n = 2), M220V (n = 1), P216L (n = 4), and F219S (n = 1) were also only found in Clade II strains. For the gene hmg1, mutations in this gene were only seen in Clade II strains (n = 7).
Fisher’s Exact tests were also conducted on 37 genes previously found to be overexpressed with triazole exposure using 3230 SNP sites (Table S3). The tests were conducted using both MIC resistance thresholds, 2 mg/L and 4 mg/L. The test conducted with all 122 strains and the MIC threshold set at 2 mg/L found 57 SNPs in or beside 14 genes to be associated with itraconazole resistance. For pan-azole resistance, 11 significantly associated SNPs were found. These SNPs were located in or beside six genes. The test conducted with the MIC threshold of 4 mg/L found the same 57 SNPs in or beside 14 genes to be significantly associated with itraconazole resistance. However, there were slight changes in the SNPs significantly associated with pan-azole resistance when using the 4 mg/L MIC threshold. The 4 mg/L MIC threshold determined 10 SNPs in or beside five genes to be associated with pan-azole resistance. Furthermore, to focus on novel mutations associated with triazole resistance not linked to cyp51A, two additional Fisher’s Exact tests were also completed after removing the 21 strains containing the L98H mutation and again after removing the 64 strains with well-known mutations in cyp51A (Table S3). Using both MIC thresholds, the results after removal of the 21 strains found three SNPs to be significantly associated with itraconazole resistance. One SNP had already been noted as associated with itraconazole resistance by the first set of Fisher’s Exact test, done prior to the 21-strain removal, which was an intergenic variant in mfsB. However, two novel intergenic variants, one in AFUA_6G01960 and the second in AFUA_1G16460, were found to be significantly associated with itraconazole resistance as well. Furthermore, after removal of the 21 strains, no SNPs were found to be significantly associated with pan-azole resistance (Table S3). In terms of the results after removal of the 64 strains, one novel SNP site was found to be significantly associated with triazole differences. Using the MIC threshold of 2 mg/L, the test identified an intergenic variant in abcA to be associated with itraconazole resistance (Table S3). A final set of Fisher’s Exact tests were conducted on a clade-level, using only strains in Clade II (n = 71). However, no SNPs sites were significantly associated with triazole resistance using this sample set. The result differences between tests could be due to genetic hitchhiking alongside the resistance polymorphism L98H in cyp51A as well as to the reduced sample size, thus decreasing the sample count of certain SNPs and making them unable to meet the Bonferroni-corrected critical p-value threshold criteria in these tests.
In about 20% to 70% of triazole-resistant clinical A. fumigatus strains, no mutations related to cyp51A were observed [29,72]. Molecular assays for the detection of A. fumigatus and its cyp51A alterations have been produced to provide rapid detection of cyp51A-mediated triazole resistance in clinical samples of A. fumigatus [72]. However, as shown in our analyses, in many of the triazole-resistant strains, relying on assays targeting only the cyp51A mutations would lead to misidentification of these strains as triazole susceptible and cause inappropriate treatment strategies. Indeed, over the last decade, novel triazole resistance mechanisms have been increasingly reported including mutations in hapE, hmg1, yap1, and cox10 genes [10,34,35,73]. The wide and growing range of resistance mechanisms seen in A. fumigatus demonstrates the high potential this fungus has for stress adaptation, including adaptation to antifungal drug resistance. Delays in the initiation of appropriate antifungal therapy are associated with overall increased mortality. Thus, tools enabling direct detection of resistance using rapid molecular methods can greatly facilitate optimal therapy for individual patients. Together, the putative variants found in this study represent promising candidates for future studies to investigate emerging mechanisms of triazole resistance in A. fumigatus. Moreover, these candidate SNPs hold great potential for developing additional diagnostic markers for accurate and rapid identification of triazole resistance in a clinical setting.

4. Materials and Methods

4.1. Whole Genome Sequences and Strains

Whole-genome sequences for 195 A. fumigatus isolates were used in this study. This sample set comprised of 184 whole-genome sequences obtained from the National Center for Biotechnology Information (NCBI) Sequence Read Archive and an additional 12 isolates that were sequenced from our previous study [63]. This strain collection spans 12 countries, across four continents consisting of 61 strains from North America, 1 from South America, 91 from Europe, 40 from Asia, as well as two strains from the International Space Station. In total, 163 of the 195 strains were isolated from a clinical environment, 29 from the natural environment and 3 of unknown sources. Among them, 122 and 123 samples had antifungal susceptibility profiles to itraconazole and voriconazole, respectively. These profiles are recorded as the minimum inhibitory concentrations (MICs) and are presented in Table S1.

4.2. Variant Calling

For genome sequence analysis, a modified pipeline from our previous study was used [74]. In brief, FastQC v0.11.5 was used to check for read quality and low-quality sequences were trimmed using Trimmomatic v0.36 [75,76]. The reads were then mapped to the reference A. fumigatus strain Af293 (GenBank accession GCA_000002655.1) using the BWA-MEM algorithm v0.7.17 [77]. Duplicate reads were removed using MarkDuplicates in the Picard tool and variants were called using FreeBayes v0.9.21-19 [78,79]. The initial variant filtering was done via vcftools to remove indels, variants with a quality score below 15, and variants with a call rate less than 0.90 [80]. A second filtering step removing multiallelic sites was also conducted using vcftools and this resulting vcf file was named the “soft-filtered” file, which contained 314,999 SNP sites.

4.3. Phylogenetic Analysis

To infer evolutionary relationships among the 195 samples, nucleotides of SNP sites were concatenated for each sample and the invariant sites of sequence alignment were removed using RAxML ascertainment bias correction [81]. The maximum likelihood phylogenetic tree was constructed based on 314,999 SNP sites, using the ASC_GTRCAT nucleotide substitution model and 500 bootstrap replicates in RAxML v8.0.25 [81]. The phylogeny was then visualized using iTOL [82]. Strains were assigned into clades based on pairwise SNP comparisons, with a threshold set at 50,000 SNPs.

4.4. Genome-Wide Association Study and Linkage Disequilibrium

Variants were annotated with SnpEff v5.0 using the Af293 reference genome annotation to determine functional effects of genetic variants [83]. Highly linked (VIF > 2) SNP markers were removed using PLINK 1.90 beta to ensure uniform sampling of the genome [84]. Association analysis via a mixed linear model was done in TASSEL 5 using two parameters: a population structure defined by 5 principal component vectors, determined based on the scree plot, and a kinship matrix calculated using the Identity by State method (Centered IBS) [85]. To avoid biases due to imbalanced allele frequencies, the minimum allele frequency was also set to 0.05 using TASSEL 5. A total of 21,432 SNP sites remained for conducting the itraconazole GWAS and a total of 21,226 SNP sites for the voriconazole GWAS. A second GWAS was conducted after the removal of strains that contained the L98H mutation in cyp51A (n = 21). For the second association analysis, a total of 22,411 and 21,214 SNP sites remained for conducting the itraconazole and voriconazole GWAS, respectively. A third GWAS was also conducted after the removal of strains that contained well-known mutations associated with in cyp51A (n = 64). For the third association analysis, a total of 20,176 and 20,278 SNP sites remained for conducting the itraconazole and voriconazole GWAS, respectively. Lastly, we conducted a GWAS examining only the strains from Clade II for itraconazole (n = 71) and voriconazole (n = 72). For the analysis focusing solely on strains from Clade II, a total of 16,702 SNP sites remained for conducting the itraconazole GWAS and a total of 16,782 SNP sites remained for the voriconazole GWAS. Using the results of the GWAS, further association mapping between the top 20 SNPs and all SNPs in the soft-filtered vcf file was conducted using TASSEL 5 to determine additional highly linked SNPs of interest.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/pathogens10060701/s1, Figure S1: Quantile–quantile (Q-Q) plots from the GWAS for (A) itraconazole and (B) voriconazole, Figure S2: Quantile–quantile (Q-Q) plots from the second GWAS analysis, after removal of the 21 strains containing the L98H mutation in cyp51A, for (A) itraconazole and (B) voriconazole, Figure S3: Quantile–quantile (Q-Q) plots from the third GWAS analysis, after removal of the 64 strains containing the well-known mutations in cyp51A, for (A) itraconazole and (B) voriconazole, Figure S4: Quantile–quantile (Q-Q) plots from the fourth GWAS analysis, focusing on strains from Clade II, for (A) itraconazole and (B) voriconazole, Figure S5: Coverage at promoter region of cyp51A from position 1,782,000 to 1,782,200 bp for the 21 strains with L98H mutation and known triazole MIC values, Table S1: Additional information on the total 195 Aspergillus fumigatus genomes, Table S2: Overexpressed genes with itraconazole and voriconazole exposure determined through previous RT-qPCR and RNA-seq analyses, Table S3: Significant single-nucleotide polymorphisms located in or near genes overexpressed with triazole exposure.

Author Contributions

Conceptualization, J.X.; methodology, Y.F. and Y.W.; software, Y.F. and Y.W.; formal analysis, Y.F. and Y.W.; investigation, Y.F. and Y.W.; resources, J.X.; writing—original draft preparation, Y.F., Y.W., G.A.K. and M.A.; writing—review and editing, J.X., Y.F., Y.W., G.A.K. and M.A.; visualization, Y.F., Y.W., G.A.K. and M.A.; supervision, J.X.; project administration, J.X.; funding acquisition, J.X. Y.F. and Y.W. share co-first authorship for this manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by grants from the Natural Sciences and Engineering Research Council of Canada (Grant No. CRDPJ 474638-14), and by the Institute of Infectious Diseases Research (IIDR) Antibiotic Resistance Initiative and the Faculty of Science’s Global Science Initiative of McMaster University. The APC was funded by IIDR. Y.F. and M.A. have been supported by Ontario Graduate Scholarships. Y.W. is supported by a MacData Fellowship. G.A.K. is supported by MITACS and NSERC PGS-D Scholarships.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Accession numbers for all 195 isolates used in this study are available and listed in Supplementary Table S1.

Acknowledgments

We thank Compute Canada and Brian Golding for access to their servers for our analyses.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Kwon-Chung, K.J.; Sugui, J.A. Aspergillus Fumigatus—What Makes the Species a Ubiquitous Human Fungal Pathogen? PLoS Pathog. 2013, 9, e1003743. [Google Scholar] [CrossRef]
  2. Brown, G.D.; Denning, D.W.; Gow, N.A.R.; Levitz, S.M.; Netea, M.G.; White, T.C. Hidden Killers: Human Fungal Infections. Sci. Transl. Med. 2012, 4, 165rv13. [Google Scholar] [CrossRef] [Green Version]
  3. Tekaia, F.; Latgé, J.-P. Aspergillus Fumigatus: Saprophyte or Pathogen? Curr. Opin. Microbiol. 2005, 8, 385–392. [Google Scholar] [CrossRef] [PubMed]
  4. Latgé, J.-P.; Chamilos, G. Aspergillus Fumigatus and Aspergillosis in 2019. Clin. Microbiol. Rev. 2019, 33, e00140-18. [Google Scholar] [CrossRef] [PubMed]
  5. Cowen, L.E.; Sanglard, D.; Howard, S.J.; Rogers, P.D.; Perlin, D.S. Mechanisms of Antifungal Drug Resistance. Cold Spring Harb. Perspect. Med. 2015, 5. [Google Scholar] [CrossRef]
  6. Alcazar-Fuoli, L.; Mellado, E. Ergosterol Biosynthesis in Aspergillus Fumigatus: Its Relevance as an Antifungal Target and Role in Antifungal Drug Resistance. Front. Microbiol. 2012, 3, 439. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Verweij, P.E.; Zhang, J.; Debets, A.J.M.; Meis, J.F.; van de Veerdonk, F.L.; Schoustra, S.E.; Zwaan, B.J.; Melchers, W.J.G. In-Host Adaptation and Acquired Triazole Resistance in Aspergillus Fumigatus: A Dilemma for Clinical Management. Lancet Infect. Dis. 2016, 16, e251–e260. [Google Scholar] [CrossRef]
  8. Lestrade, P.P.A.; Meis, J.F.; Melchers, W.J.G.; Verweij, P.E. Triazole Resistance in Aspergillus Fumigatus: Recent Insights and Challenges for Patient Management. Clin. Microbiol. Infect. 2019, 25, 799–806. [Google Scholar] [CrossRef]
  9. Howard, S.J.; Cerar, D.; Anderson, M.J.; Albarrag, A.; Fisher, M.C.; Pasqualotto, A.C.; Laverdiere, M.; Arendrup, M.C.; Perlin, D.S.; Denning, D.W. Frequency and Evolution of Azole Resistance in Aspergillus Fumigatus Associated with Treatment Failure1. Emerg. Infect. Dis. 2009, 15, 1068–1076. [Google Scholar] [CrossRef]
  10. Camps, S.M.T.; Dutilh, B.E.; Arendrup, M.C.; Rijs, A.J.M.M.; Snelders, E.; Huynen, M.A.; Verweij, P.E.; Melchers, W.J.G. Discovery of a HapE Mutation That Causes Azole Resistance in Aspergillus Fumigatus through Whole Genome Sequencing and Sexual Crossing. PLoS ONE 2012, 7, e50034. [Google Scholar] [CrossRef] [Green Version]
  11. Chowdhary, A.; Kathuria, S.; Xu, J.; Meis, J.F. Emergence of Azole-Resistant Aspergillus Fumigatus Strains Due to Agricultural Azole Use Creates an Increasing Threat to Human Health. PLoS Pathog. 2013, 9, e1003633. [Google Scholar] [CrossRef]
  12. Snelders, E.; van der Lee, H.A.L.; Kuijpers, J.; Rijs, A.J.M.M.; Varga, J.; Samson, R.A.; Mellado, E.; Donders, A.R.T.; Melchers, W.J.G.; Verweij, P.E. Emergence of Azole Resistance in Aspergillus Fumigatus and Spread of a Single Resistance Mechanism. PLoS Med. 2008, 5, e219. [Google Scholar] [CrossRef] [PubMed]
  13. Ashu, E.E.; Hagen, F.; Chowdhary, A.; Meis, J.F.; Xu, J. Global Population Genetic Analysis of Aspergillus Fumigatus. mSphere 2017, 2, e00019-17. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Sewell, T.R.; Zhu, J.; Rhodes, J.; Hagen, F.; Meis, J.F.; Fisher, M.C.; Jombart, T. Nonrandom Distribution of Azole Resistance across the Global Population of Aspergillus Fumigatus. mBio 2019, 10. [Google Scholar] [CrossRef] [Green Version]
  15. Chowdhary, A.; Kathuria, S.; Xu, J.; Sharma, C.; Sundar, G.; Singh, P.K.; Gaur, S.N.; Hagen, F.; Klaassen, C.H.; Meis, J.F. Clonal Expansion and Emergence of Environmental Multiple-Triazole-Resistant Aspergillus Fumigatus Strains Carrying the TR34/L98H Mutations in the Cyp51A Gene in India. PLoS ONE 2012, 7, e52871. [Google Scholar] [CrossRef]
  16. Ahangarkani, F.; Badali, H.; Abbasi, K.; Nabili, M.; Khodavaisy, S.; de Groot, T.; Meis, J.F. Clonal Expansion of Environmental Triazole Resistant Aspergillus Fumigatus in Iran. J. Fungi 2020, 6, 199. [Google Scholar] [CrossRef]
  17. Melo, A.M.; Stevens, D.A.; Tell, L.A.; Veríssimo, C.; Sabino, R.; Xavier, M.O. Aspergillosis, Avian Species and the One Health Perspective: The Possible Importance of Birds in Azole Resistance. Microorganisms 2020, 8, 37. [Google Scholar] [CrossRef]
  18. Cacciuttolo, E.; Rossi, G.; Nardoni, S.; Legrottaglie, R.; Mani, P. Anatomopathological Aspects of Avian Aspergillosis. Vet. Res. Commun. 2009, 33, 521–527. [Google Scholar] [CrossRef]
  19. Resendiz Sharpe, A.; Lagrou, K.; Meis, J.F.; Chowdhary, A.; Lockhart, S.R.; Verweij, P.E.; ISHAM/ECMM Aspergillus Resistance Surveillance working group. Triazole Resistance Surveillance in Aspergillus Fumigatus. Med. Mycol. 2018, 56, 83–92. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Lestrade, P.P.A.; Buil, J.B.; van der Beek, M.T.; Kuijper, E.J.; van Dijk, K.; Kampinga, G.A.; Rijnders, B.J.A.; Vonk, A.G.; de Greeff, S.C.; Schoffelen, A.F.; et al. Paradoxal Trends in Azole-Resistant Aspergillus Fumigatus in a National Multicenter Surveillance Program, the Netherlands, 2013–2018. Emerg. Infect. Dis. 2020, 26, 1447–1455. [Google Scholar] [CrossRef]
  21. Buil, J.B.; Snelders, E.; Denardi, L.B.; Melchers, W.J.G.; Verweij, P.E. Trends in Azole Resistance in Aspergillus Fumigatus, the Netherlands, 1994–2016. Emerg. Infect. Dis. 2019, 25, 176–178. [Google Scholar] [CrossRef] [Green Version]
  22. Nabili, M.; Shokohi, T.; Moazeni, M.; Khodavaisy, S.; Aliyali, M.; Badiee, P.; Zarrinfar, H.; Hagen, F.; Badali, H. High Prevalence of Clinical and Environmental Triazole-Resistant Aspergillus Fumigatus in Iran: Is It a Challenging Issue? J. Med. Microbiol. 2016, 65, 468–475. [Google Scholar] [CrossRef]
  23. Heo, S.T.; Tatara, A.M.; Jiménez-Ortigosa, C.; Jiang, Y.; Lewis, R.E.; Tarrand, J.; Tverdek, F.; Albert, N.D.; Verweij, P.E.; Meis, J.F.; et al. Changes in In Vitro Susceptibility Patterns of Aspergillus to Triazoles and Correlation With Aspergillosis Outcome in a Tertiary Care Cancer Center, 1999–2015. Clin. Infect. Dis. 2017, 65, 216–225. [Google Scholar] [CrossRef] [Green Version]
  24. Abdolrasouli, A.; Petrou, M.A.; Park, H.; Rhodes, J.L.; Rawson, T.M.; Moore, L.S.P.; Donaldson, H.; Holmes, A.H.; Fisher, M.C.; Armstrong-James, D. Surveillance for Azole-Resistant Aspergillus Fumigatus in a Centralized Diagnostic Mycology Service, London, United Kingdom, 1998–2017. Front. Microbiol. 2018, 9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Fraczek, M.G.; Bromley, M.; Buied, A.; Moore, C.B.; Rajendran, R.; Rautemaa, R.; Ramage, G.; Denning, D.W.; Bowyer, P. The Cdr1B Efflux Transporter Is Associated with Non-Cyp51a-Mediated Itraconazole Resistance in Aspergillus Fumigatus. J. Antimicrob. Chemother. 2013, 68, 1486–1496. [Google Scholar] [CrossRef] [PubMed]
  26. Zhou, D.; Korfanty, G.A.; Mo, M.; Wang, R.; Li, X.; Li, H.; Li, S.; Wu, J.-Y.; Zhang, K.-Q.; Zhang, Y.; et al. Extensive Genetic Diversity and Widespread Azole Resistance in Greenhouse Populations of Aspergillus Fumigatus in Yunnan, China. mSphere 2021, 6. [Google Scholar] [CrossRef] [PubMed]
  27. Mohd-Assaad, N.; McDonald, B.A.; Croll, D. Multilocus Resistance Evolution to Azole Fungicides in Fungal Plant Pathogen Populations. Mol. Ecol. 2016, 25, 6124–6142. [Google Scholar] [CrossRef]
  28. Talas, F.; Kalih, R.; Miedaner, T.; McDonald, B.A. Genome-Wide Association Study Identifies Novel Candidate Genes for Aggressiveness, Deoxynivalenol Production, and Azole Sensitivity in Natural Field Populations of Fusarium Graminearum. Mol. Plant. Microbe Interact. 2016, 29, 417–430. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Zhao, S.; Ge, W.; Watanabe, A.; Fortwendel, J.R.; Gibbons, J.G. Genome-Wide Association for Itraconazole Sensitivity in Non-Resistant Clinical Isolates of Aspergillus Fumigatus. Front. Fungal. Biol. 2021, 1. [Google Scholar] [CrossRef]
  30. Garcia-Rubio, R.; Cuenca-Estrella, M.; Mellado, E. Triazole Resistance in Aspergillus Species: An Emerging Problem. Drugs 2017, 77, 599–613. [Google Scholar] [CrossRef] [PubMed]
  31. Wei, X.; Zhang, Y.; Lu, L. The Molecular Mechanism of Azole Resistance in Aspergillus Fumigatus: From Bedside to Bench and Back. J. Microbiol. 2015, 53, 91–99. [Google Scholar] [CrossRef] [PubMed]
  32. Gonzalez-Jimenez, I.; Lucio, J.; Amich, J.; Cuesta, I.; Sanchez Arroyo, R.; Alcazar-Fuoli, L.; Mellado, E. A Cyp51B Mutation Contributes to Azole Resistance in Aspergillus Fumigatus. J. Fungi 2020, 6, 315. [Google Scholar] [CrossRef] [PubMed]
  33. Rybak, J.M.; Ge, W.; Wiederhold, N.P.; Parker, J.E.; Kelly, S.L.; Rogers, P.D.; Fortwendel, J.R. Mutations in Hmg1, Challenging the Paradigm of Clinical Triazole Resistance in Aspergillus Fumigatus. mBio 2019, 10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Hagiwara, D.; Arai, T.; Takahashi, H.; Kusuya, Y.; Watanabe, A.; Kamei, K. Non-Cyp51A Azole-Resistant Aspergillus Fumigatus Isolates with Mutation in HMG-CoA Reductase. Emerg. Infect. Dis. 2018, 24, 1889–1897. [Google Scholar] [CrossRef] [Green Version]
  35. Wei, X.; Chen, P.; Gao, R.; Li, Y.; Zhang, A.; Liu, F.; Lu, L. Screening and Characterization of a Non-Cyp51A Mutation in an Aspergillus Fumigatus Cox10 Strain Conferring Azole Resistance. Antimicrob. Agents Chemother. 2017, 61. [Google Scholar] [CrossRef] [Green Version]
  36. Ballard, E.; Weber, J.; Melchers, W.J.G.; Tammireddy, S.; Whitfield, P.D.; Brakhage, A.A.; Brown, A.J.P.; Verweij, P.E.; Warris, A. Recreation of In-Host Acquired Single Nucleotide Polymorphisms by CRISPR-Cas9 Reveals an Uncharacterised Gene Playing a Role in Aspergillus Fumigatus Azole Resistance via a Non-Cyp51A Mediated Resistance Mechanism. Fungal Genet. Biol. 2019, 130, 98–106. [Google Scholar] [CrossRef] [PubMed]
  37. Bromley, M.; Johns, A.; Davies, E.; Fraczek, M.; Mabey Gilsenan, J.; Kurbatova, N.; Keays, M.; Kapushesky, M.; Gut, M.; Gut, I.; et al. Mitochondrial Complex I Is a Global Regulator of Secondary Metabolism, Virulence and Azole Sensitivity in Fungi. PLoS ONE 2016, 11, e0158724. [Google Scholar] [CrossRef] [Green Version]
  38. Da Silva Ferreira, M.E.; Malavazi, I.; Savoldi, M.; Brakhage, A.A.; Goldman, M.H.S.; Kim, H.S.; Nierman, W.C.; Goldman, G.H. Transcriptome Analysis of Aspergillus Fumigatus Exposed to Voriconazole. Curr. Genet. 2006, 50, 32–44. [Google Scholar] [CrossRef] [PubMed]
  39. Du, W.; Zhai, P.; Wang, T.; Bromley, M.J.; Zhang, Y.; Lu, L. The C 2 H 2 Transcription Factor SltA Contributes to Azole Resistance by Coregulating the Expression of the Drug Target Erg11A and the Drug Efflux Pump Mdr1 in Aspergillus Fumigatus. Antimicrob. Agents Chemother. 2021, 65, e01839-20. [Google Scholar] [CrossRef]
  40. Meis, J.F.; Chowdhary, A.; Rhodes, J.L.; Fisher, M.C.; Verweij, P.E. Clinical Implications of Globally Emerging Azole Resistance in Aspergillus Fumigatus. Philos. Trans. R Soc. Lond B Biol. Sci. 2016, 371. [Google Scholar] [CrossRef] [Green Version]
  41. Klaassen, C.H.W.; Gibbons, J.G.; Fedorova, N.D.; Meis, J.F.; Rokas, A. Evidence for Genetic Differentiation and Variable Recombination Rates among Dutch Populations of the Opportunistic Human Pathogen Aspergillus Fumigatus. Mol. Ecol. 2012, 21, 57–70. [Google Scholar] [CrossRef] [Green Version]
  42. Molecular Epidemiology of Aspergillus Fumigatus Isolates Recovered from Water, Air, and Patients Shows Two Clusters of Genetically Distinct Strains. Available online: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC193792/ (accessed on 31 March 2021).
  43. Bain, J.M.; Tavanti, A.; Davidson, A.D.; Jacobsen, M.D.; Shaw, D.; Gow, N.A.R.; Odds, F.C. Multilocus Sequence Typing of the Pathogenic Fungus Aspergillus Fumigatus. J. Clin. Microbiol. 2007, 45, 1469–1477. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Seufert, R.; Sedlacek, L.; Kahl, B.; Hogardt, M.; Hamprecht, A.; Haase, G.; Gunzer, F.; Haas, A.; Grauling-Halama, S.; MacKenzie, C.R.; et al. Prevalence and Characterization of Azole-Resistant Aspergillus Fumigatus in Patients with Cystic Fibrosis: A Prospective Multicentre Study in Germany. J. Antimicrob. Chemother. 2018, 73, 2047–2053. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Romero, M.; Messina, F.; Marin, E.; Arechavala, A.; Depardo, R.; Walker, L.; Negroni, R.; Santiso, G. Antifungal Resistance in Clinical Isolates of Aspergillus Spp.: When Local Epidemiology Breaks the Norm. J. Fungi 2019, 5, 41. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Verweij, P.E.; van de Sande-Bruisma, N.; Kema, G.H.J.; Melchers, W.J.G. Azole resistance in Aspergillus fumigatus in the Netherlands--increase due to environmental fungicides? Ned. Tijdschr. Geneeskd 2012, 156, A4458. [Google Scholar]
  47. Özmerdiven, G.E.; Ak, S.; Ener, B.; Ağca, H.; Cilo, B.D.; Tunca, B.; Akalın, H. First Determination of Azole Resistance in Aspergillus Fumigatus Strains Carrying the TR34/L98H Mutations in Turkey. J. Infect. Chemother. 2015, 21, 581–586. [Google Scholar] [CrossRef]
  48. Vermeulen, E.; Maertens, J.; De Bel, A.; Nulens, E.; Boelens, J.; Surmont, I.; Mertens, A.; Boel, A.; Lagrou, K. Nationwide Surveillance of Azole Resistance in Aspergillus Diseases. Antimicrob. Agents Chemother. 2015, 59, 4569–4576. [Google Scholar] [CrossRef] [Green Version]
  49. Sewell, T.R.; Zhang, Y.; Brackin, A.P.; Shelton, J.M.G.; Rhodes, J.; Fisher, M.C. Elevated Prevalence of Azole-Resistant Aspergillus Fumigatus in Urban versus Rural Environments in the United Kingdom. Antimicrob. Agents Chemother. 2019, 63. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Van der Linden, J.W.M.; Arendrup, M.C.; Warris, A.; Lagrou, K.; Pelloux, H.; Hauser, P.M.; Chryssanthou, E.; Mellado, E.; Kidd, S.E.; Tortorano, A.M.; et al. Prospective Multicenter International Surveillance of Azole Resistance in Aspergillus Fumigatus. Emerg. Infect. Dis. 2015, 21, 1041–1044. [Google Scholar] [CrossRef] [PubMed]
  51. Chowdhary, A.; Sharma, C.; Meis, J.F. Azole-Resistant Aspergillosis: Epidemiology, Molecular Mechanisms, and Treatment. J. Infect. Dis. 2017, 216, S436–S444. [Google Scholar] [CrossRef] [Green Version]
  52. Takeshita, N.; Diallinas, G.; Fischer, R. The Role of Flotillin FloA and Stomatin StoA in the Maintenance of Apical Sterol-Rich Membrane Domains and Polarity in the Filamentous Fungus Aspergillus Nidulans. Mol. Microbiol. 2012, 83, 1136–1152. [Google Scholar] [CrossRef]
  53. Shekhova, E.; Kniemeyer, O.; Brakhage, A.A. Induction of Mitochondrial Reactive Oxygen Species Production by Itraconazole, Terbinafine, and Amphotericin B as a Mode of Action against Aspergillus Fumigatus. Antimicrob. Agents Chemother. 2017, 61. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Esquivel, B.D.; Rybak, J.M.; Barker, K.S.; Fortwendel, J.R.; Rogers, P.D.; White, T.C. Characterization of the Efflux Capability and Substrate Specificity of Aspergillus Fumigatus PDR5-like ABC Transporters Expressed in Saccharomyces Cerevisiae. mBio 2020, 11, e00338-20. [Google Scholar] [CrossRef] [Green Version]
  55. Fatema, U.; Broberg, A.; Jensen, D.F.; Karlsson, M.; Dubey, M. Functional Analysis of Polyketide Synthase Genes in the Biocontrol Fungus Clonostachys Rosea. Sci. Rep. 2018, 8, 15009. [Google Scholar] [CrossRef]
  56. Bignell, E.; Cairns, T.C.; Throckmorton, K.; Nierman, W.C.; Keller, N.P. Secondary Metabolite Arsenal of an Opportunistic Pathogenic Fungus. Philos. Trans. R Soc. Lond. B Biol. Sci. 2016, 371. [Google Scholar] [CrossRef] [PubMed]
  57. Beck, J.; Echtenacher, B.; Ebel, F. Woronin Bodies, Their Impact on Stress Resistance and Virulence of the Pathogenic Mould Aspergillus Fumigatus and Their Anchoring at the Septal Pore of Filamentous Ascomycota. Mol. Microbiol. 2013, 89, 857–871. [Google Scholar] [CrossRef]
  58. Rhodes, J.; Abdolrasouli, A.; Dunne, K.; Sewell, T.R.; Zhang, Y.; Ballard, E.; Brackin, A.P.; van Rhijn, N.; Tsitsopoulou, A.; Posso, R.B.; et al. Tracing Patterns of Evolution and Acquisition of Drug Resistant Aspergillus Fumigatus Infection from the Environment Using Population Genomics. bioRxiv 2021. [Google Scholar] [CrossRef]
  59. Paul, S.; Stamnes, M.; Thomas, G.H.; Liu, H.; Hagiwara, D.; Gomi, K.; Filler, S.G.; Moye-Rowley, W.S. AtrR Is an Essential Determinant of Azole Resistance in Aspergillus Fumigatus. mBio 2019, 10. [Google Scholar] [CrossRef] [Green Version]
  60. Shelest, E. Transcription Factors in Fungi: TFome Dynamics, Three Major Families, and Dual-Specificity TFs. Front. Genet. 2017, 8. [Google Scholar] [CrossRef] [Green Version]
  61. Hunt, R.C.; Simhadri, V.L.; Iandoli, M.; Sauna, Z.E.; Kimchi-Sarfaty, C. Exposing Synonymous Mutations. Trends Genet. 2014, 30, 308–321. [Google Scholar] [CrossRef] [PubMed]
  62. Barrett, L.W.; Fletcher, S.; Wilton, S.D. Regulation of Eukaryotic Gene Expression by the Untranslated Gene Regions and Other Non-Coding Elements. Cell Mol. Life Sci. 2012, 69, 3613–3634. [Google Scholar] [CrossRef] [Green Version]
  63. Pemán, J.; Salavert, M.; Cantón, E.; Jarque, I.; Romá, E.; Zaragoza, R.; Viudes, Á.; Gobernado, M. Voriconazole in the Management of Nosocomial Invasive Fungal Infections. Ther. Clin. Risk Manag. 2006, 2, 129–158. [Google Scholar] [CrossRef] [Green Version]
  64. Chen, P.; Liu, M.; Zeng, Q.; Zhang, Z.; Liu, W.; Sang, H.; Lu, L. Uncovering New Mutations Conferring Azole Resistance in the Aspergillus Fumigatus Cyp51A Gene. Front. Microbiol. 2020, 10. [Google Scholar] [CrossRef] [Green Version]
  65. Bernal-Martínez, L.; Gil, H.; Rivero-Menéndez, O.; Gago, S.; Cuenca-Estrella, M.; Mellado, E.; Alastruey-Izquierdo, A. Development and Validation of a High-Resolution Melting Assay To Detect Azole Resistance in Aspergillus Fumigatus. Antimicrob. Agents Chemother. 2017, 61. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Snelders, E.; Karawajczyk, A.; Schaftenaar, G.; Verweij, P.E.; Melchers, W.J.G. Azole Resistance Profile of Amino Acid Changes in Aspergillus Fumigatus CYP51A Based on Protein Homology Modeling. Antimicrob. Agents Chemother. 2010, 54, 2425–2430. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Camps, S.M.T.; van der Linden, J.W.M.; Li, Y.; Kuijper, E.J.; van Dissel, J.T.; Verweij, P.E.; Melchers, W.J.G. Rapid Induction of Multiple Resistance Mechanisms in Aspergillus Fumigatus during Azole Therapy: A Case Study and Review of the Literature. Antimicrob. Agents Chemother. 2012, 56, 10–16. [Google Scholar] [CrossRef] [Green Version]
  68. Chowdhary, A.; Sharma, C.; Hagen, F.; Meis, J.F. Exploring Azole Antifungal Drug Resistance in Aspergillus Fumigatus with Special Reference to Resistance Mechanisms. Future Microbiol. 2014, 9, 697–711. [Google Scholar] [CrossRef]
  69. Liu, M.; Zeng, R.; Zhang, L.; Li, D.; Lv, G.; Shen, Y.; Zheng, H.; Zhang, Q.; Zhao, J.; Zheng, N.; et al. Multiple Cyp51A-Based Mechanisms Identified in Azole-Resistant Isolates of Aspergillus Fumigatus from China. Antimicrob. Agents Chemother. 2015, 59, 4321–4325. [Google Scholar] [CrossRef] [Green Version]
  70. Liu, M.; Zheng, N.; Li, D.; Zheng, H.; Zhang, L.; Ge, H.; Liu, W. Cyp51A -Based Mechanism of Azole Resistance in Aspergillus Fumigatus: Illustration by a New 3D Structural Model of Aspergillus Fumigatus CYP51A Protein. Med. Mycol. 2016, 54, 400–408. [Google Scholar] [CrossRef] [Green Version]
  71. Bowyer, P.; Mosquera, J.; Anderson, M.; Birch, M.; Bromley, M.; Denning, D.W. Identification of Novel Genes Conferring Altered Azole Susceptibility in Aspergillus Fumigatus. FEMS Microbiol. Lett. 2012, 332, 10–19. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Dudakova, A.; Spiess, B.; Tangwattanachuleeporn, M.; Sasse, C.; Buchheidt, D.; Weig, M.; Groß, U.; Bader, O. Molecular Tools for the Detection and Deduction of Azole Antifungal Drug Resistance Phenotypes in Aspergillus Species. Clin. Microbiol. Rev. 2017, 30, 1065–1091. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Qiao, J.; Liu, W.; Li, R. Truncated Afyap1 Attenuates Antifungal Susceptibility of Aspergillus Fumigatus to Voriconazole and Confers Adaptation of the Fungus to Oxidative Stress. Mycopathologia 2010, 170, 155–160. [Google Scholar] [CrossRef]
  74. Fan, Y.; Wang, Y.; Xu, J. Comparative Genome Sequence Analyses of Geographic Samples of Aspergillus Fumigatus-Relevance for Amphotericin B Resistance. Microorganisms 2020, 8, 1673. [Google Scholar] [CrossRef] [PubMed]
  75. Babraham Bioinformatics—FastQC A Quality Control Tool for High Throughput Sequence Data. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed on 31 March 2021).
  76. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  77. Li, H. Aligning Sequence Reads, Clone Sequences and Assembly Contigs with BWA-MEM. arXiv 2013, arXiv:1303.3997. [Google Scholar]
  78. Picard Tools—By Broad Institute. Available online: https://broadinstitute.github.io/picard/ (accessed on 31 March 2021).
  79. Garrison, E.; Marth, G. Haplotype-Based Variant Detection from Short-Read Sequencing. arXiv 2012, arXiv:1207.3907. [Google Scholar]
  80. Danecek, P.; Auton, A.; Abecasis, G.; Albers, C.A.; Banks, E.; DePristo, M.A.; Handsaker, R.E.; Lunter, G.; Marth, G.T.; Sherry, S.T.; et al. The Variant Call Format and VCFtools. Bioinformatics 2011, 27, 2156–2158. [Google Scholar] [CrossRef] [PubMed]
  81. Stamatakis, A. RAxML Version 8: A Tool for Phylogenetic Analysis and Post-Analysis of Large Phylogenies. Bioinformatics 2014, 30, 1312–1313. [Google Scholar] [CrossRef]
  82. Letunic, I.; Bork, P. Interactive Tree of Life (ITOL) v4: Recent Updates and New Developments. Nucleic Acids Res. 2019, 47, W256–W259. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. Cingolani, P.; Platts, A.; Wang, L.L.; Coon, M.; Nguyen, T.; Wang, L.; Land, S.J.; Lu, X.; Ruden, D.M. A Program for Annotating and Predicting the Effects of Single Nucleotide Polymorphisms, SnpEff: SNPs in the Genome of Drosophila Melanogaster Strain W1118; Iso-2; Iso-3. Fly 2012, 6, 80–92. [Google Scholar] [CrossRef] [Green Version]
  84. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.R.; Bender, D.; Maller, J.; Sklar, P.; de Bakker, P.I.W.; Daly, M.J.; et al. PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am. J. Hum. Genet. 2007, 81, 559–575. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  85. TASSEL: Software for Association Mapping of Complex Traits in Diverse Samples | Bioinformatics | Oxford Academic. Available online: https://academic.oup.com/bioinformatics/article/23/19/2633/185151 (accessed on 31 March 2021).
Figure 1. Maximum likelihood phylogenetic tree detailing the strain characteristics. Branches with a red dot represent those with over 75% bootstrap support, based on 500 bootstrap iterations. The inner-most circle denotes the clade affiliation of strains with strain names corresponding to those in Supplementary Table S1. The second inner-most circle represents country of origin for individual strains with different colors representing different countries as shown in the left “Country” panel. The third circle from the inside denotes strain ecological niche, with hollow squares representing strains from the natural environment, solid red squares representing strains from the clinical environment, and the source for the remaining strains (unmarked) were unknown. The itraconazole and voriconazole minimum inhibitory concentrations (MIC) were represented in the two outer circles with different colors representing different MIC values as shown in the left “TriazoleMIC” panel. The white boxes in the two outer circles represent strains with no MIC data. The branch length separating Clade I from the two other clades was manually truncated to make relationships in the other two clades more visible.
Figure 1. Maximum likelihood phylogenetic tree detailing the strain characteristics. Branches with a red dot represent those with over 75% bootstrap support, based on 500 bootstrap iterations. The inner-most circle denotes the clade affiliation of strains with strain names corresponding to those in Supplementary Table S1. The second inner-most circle represents country of origin for individual strains with different colors representing different countries as shown in the left “Country” panel. The third circle from the inside denotes strain ecological niche, with hollow squares representing strains from the natural environment, solid red squares representing strains from the clinical environment, and the source for the remaining strains (unmarked) were unknown. The itraconazole and voriconazole minimum inhibitory concentrations (MIC) were represented in the two outer circles with different colors representing different MIC values as shown in the left “TriazoleMIC” panel. The white boxes in the two outer circles represent strains with no MIC data. The branch length separating Clade I from the two other clades was manually truncated to make relationships in the other two clades more visible.
Pathogens 10 00701 g001
Figure 2. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 122) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 123). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Figure 2. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 122) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 123). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Pathogens 10 00701 g002aPathogens 10 00701 g002b
Figure 3. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus after removal of strains containing the L98H mutation in cyp51A. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 101) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 102). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Figure 3. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus after removal of strains containing the L98H mutation in cyp51A. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 101) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 102). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Pathogens 10 00701 g003aPathogens 10 00701 g003b
Figure 4. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus after removal of strains containing the mutations in cyp51A. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 58) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 59). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Figure 4. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus after removal of strains containing the mutations in cyp51A. (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 58) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 59). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Pathogens 10 00701 g004aPathogens 10 00701 g004b
Figure 5. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus in Clade II (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 71) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 72). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Figure 5. The Manhattan plot showing genome-wide SNPs associated with triazole resistance in A. fumigatus in Clade II (A) SNPs associated with itraconazole resistance in A. fumigatus isolates (n = 71) and (B) SNPs associated with voriconazole resistance in A. fumigatus isolates (n = 72). The top 20 SNPs in each analysis are separated out by the red dashed line. The plot is depicted with chromosome position on the X-axis and the −log10(p-value) on the Y-axis.
Pathogens 10 00701 g005
Table 1. The 44 known mutation sites previously reported to be associated with triazole resistance and results of the Fisher’s Exact tests using 122 A. fumigatus strains with known itraconazole and voriconazole MICs.
Table 1. The 44 known mutation sites previously reported to be associated with triazole resistance and results of the Fisher’s Exact tests using 122 A. fumigatus strains with known itraconazole and voriconazole MICs.
GeneCodonAmino Acid ChangeChromosome—Position (bp)Fisher’s Exact Test (p-Values), MIC ≥ 2 mg/LFisher’s Exact Test (p-Values), MIC ≥ 4 mg/LReferences
ItraconazolePan-AzoleItraconazolePan-Azole
cyp51A
(AFUA_4G06890)
N22NA 1[30]
* F46YCHR 4—1,781,6864.50 × 10−33.54 × 10−24.50 × 10−32.96 × 10−2[30]
S52NA 1[31]
G54V, ECHR 4—1,781,6621.55 × 10−11.001.55 × 10−11.00[30]
W, RCHR 4—1,781,6638.16 × 10−22.72 × 10−28.16 × 10−28.90 × 10−3
Q88NA 1[31]
L98HCHR 4—1,781,4591.19 × 10−53.33 × 10−61.19 × 10−51.47 × 10−5[30]
V101NA 1[31]
Y121FCHR 4—1,781,3902.81 × 10−12.42 × 10−12.81 × 10−19.71 × 10−2[30]
N125NA 1[31]
G138CCHR 4—1,781,3408.09 × 10−21.008.09 × 10−21.00[31]
Q141NA 1[31]
H147YCHR 4—1,781,3131.001.001.001.00[31]
F165NA 1[30]
* M172NA 1[30]
P216LCHR 4—1,781,1055.23 × 10−14.95 × 10−15.23 × 10−12.19 × 10−1[30]
F219SCHR 4—1,781,0962.91 × 10−12.44 × 10−12.91 × 10−11.04 × 10−1[30]
M220ICHR 4—1,781,0921.001.001.001.00[30]
VCHR 4—1,781,0942.87 × 10−11.002.87 × 10−11.00
M236NA 1[31]
* N248NA 1[30]
* D255ECHR 4—1,780,9876.30 × 10−11.006.30 × 10−11.00[30]
D262NA 1[30]
A284NA 1[30]
T289ACHR 4—1,780,8872.91 × 10−12.44 × 10−12.91 × 10−11.04 × 10−1[30]
S297TCHR 4—1,780,8635.26 × 10−11.005.26 × 10−11.00[31]
P394NA 1[31]
* E427KCHR 4 -
1,780,473
5.00 × 10−33.69 × 10−25.00 × 10−33.11 × 10−2[30]
Y431NA 1[30]
G432NA 1[30]
G434NA 1[30]
T440NA 1[30]
G448SCHR 4—1,780,4101.001.001.004.71 × 10−1[30]
N479NA 1[30]
Y491NA 1[30]
F495ICHR 4—1,780,2695.22 × 10−11.005.22 × 10−11.00[31]
cyp51B
(AFUA_7G03740)
G457NA 1[32]
hapE
(AFUA_6G05300)
P88NA 1[30]
hmg1
(AFUA_2G03700)
F262NA 1[33]
S305PCHR 2—985,9595.22 × 10−14.95 × 10−15.22 × 10−12.14 × 10−1[33]
P309LCHR 2—985,9721.001.001.001.00[33]
I412T, SCHR 2—986,2811.56 × 10−11.17 × 10−11.56 × 10−14.34 × 10−2[33]
erg6
(AFUA_4G03630)
A350NA 1[34]
cox10
(AFUA_4G08340)
R243NA 1[35]
AFUA_7G01960L167Stop GainedCHR 7—531,5821.001.001.004.66 × 10−1[36]
AFUA_2G10600E180DCHR 2—2,714,1886.39 × 10−22.33 × 10−26.39 × 10−28.79 × 10−3[37]
* The reference strain Af293 contains the cyp51A mutations F46Y, M172V, N248T, D255E, and E427K. 1 The mutation sites were not found in the soft filtered genotype file, prior to multiallelic site removal.
Table 2. Overexpressed genes associated with triazole exposure in A. fumigatus from previous RT-qPCR and RNA-seq studies.
Table 2. Overexpressed genes associated with triazole exposure in A. fumigatus from previous RT-qPCR and RNA-seq studies.
Overexpressed Gene NameEncoded ProteinFold Change When Exposed to ItraconazoleFold Change When Exposed to VoriconazoleReferences
abcA-1
(AFUA_1G17440)
ABC multidrug transporter7.1NA[25]
abcA-2
(AFUA_2G15130)
~6.50NA[25]
abcB
(AFUA_1G10390)
~4.50~5.00–13.00[25,38]
abcC
(AFUA_1G14330)
~5.50~5.00–>20.00[25,38]
abcD
(AFUA_6G03470)
~4.50~2.00–>20.00[25,38]
abcE
(AFUA_7G00480)
~1.00~2.00–>20.00[25,38]
atrF
(AFUA_6G04360)
31.7NA[25]
mdr1
(AFUA_5G06070)
~5.00~2.00–5.00[25,38]
mdr4
(AFUA_1G12690)
~4.70NA[25]
AFUA_5G02260ABC multidrug transporter, putative~4.90NA[25]
AFUA_2G11580MFS multidrug transporter, putative14.2NA[25]
mfs56
(AFUA_1G05010)
~4.50–700.00NA[25]
mfsA
(AFUA_8G05710)
MFS multidrug transporter
~4.70~1.50–11.00[25,38]
mfsB
(AFUA_1G15490)
NA~4.00–18.00[38]
mfsC
(AFUA_1G03200)
~7.90~2.50–30.00[25,38]
cyp51A
(AFUA_4G06890)
14-alpha sterol demethylase21.00–550.90NA[25]
fbpA
(AFUA_1G14050)
F-box domain proteinNA~ >50.00–600.00[38]
aaaA
(AFUA_7G06680)
AAA-family ATPase, putativeNA~2.00–90.00[38]
finA
(AFUA_8G05800)
C6 zinc finger domain proteinNA~4.00–40.00[38]
AFUA_1G02870Transcription factor involved in oxidative stress response, putative2.48–2.61NA[39]
AFUA_1G04140C6 finger domain protein, putative2.04–2.94NA[39]
AFUA_6G019602.01–3.02NA[39]
AFUA_6G034302.78–2.93NA[39]
fumR
(AFUA_8G00420)
C6 zinc finger transcription factor4.00–4.70NA[39]
AFUA_5G075102.39–3.50NA[39]
AFUA_3G09130C6 transcription factor, putative1.73–2.22NA[39]
AFUA_8G073601.90–1.92NA[39]
cpcA
(AFUA_4G12470)
BZIP transcription factorNA>1.50–~5.50[38]
AFUA_1G16460BZIP transcription factor (LziP), putative1.75–2.12NA[39]
AFUA_7G03910C2H2 zinc finger protein2.50–2.86NA[39]
ace1
(AFUA_3G08010)
C2H2 zinc-finger transcription factor, putative1.66–2.32NA[39]
AFUA_4G136002.30–2.71NA[39]
zfpA
(AFUA_8G05010)
NA~1.50–60.00[38]
AFUA_2G01190Cu-dependent DNA-binding protein, putative1.30–2.10NA[39]
AFUA_4G06170Predicted DNA-binding transcription factor3.79–3.89NA[39]
AFUA_5G026552.75–3.84NA[39]
ada
(AFUA_5G06350)
DNA repair and transcription factor, putative1.23–2.05NA[39]
Table 3. Top 20 significant SNPs obtained from the GWAS that were associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 3. Top 20 significant SNPs obtained from the GWAS that were associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-value)Gene IDAnnotationPredicted Effect
12,538,614A to C6.15AFUA_1G09780Stomatin family proteinMissense Variant
(Asp418Ala)
21,845,323C to T5.14AFUA_2G06330—AFUA_2G07340Ubiquitin C-terminal hydrolase, putative—COP9 subunit 3, putativeIntergenic Region
21,899,353C to T4.58AFUA_2G07430—AFUA_2G07440DDHD domain protein—Thioesterase family proteinIntergenic Region
32,408,041T to C4.29AFUA_3G09400—AFUA_3G09450MFS transporter (Hol1), putative—Alpha/beta fold family hydrolase, putativeIntergenic Region
8623,331G to T4.09AFUA_8G02330Endoglucanase, putativeNon-coding Transcript Variant
43,737,973C to T3.66AFUA_4G14300—AFUA_4G14310Dynamin family GTPase, putative—APH domain-containing proteinIntergenic Region
52,063,521C to A3.42AFUA_5G08150ABC bile acid transporter, putativeMissense Variant
(His105Gln)
52,069,483G to A3.23AFUA_5G08160—AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
31,953,910G to A3.14AFUA_3G07730—AFUA_3G07740Uncharacterized protein—Uncharacterized proteinIntergenic Region
63,054,001C to G3.08AFUA_6G12145—AFUA_6G12150Uncharacterized protein—BZIP transcription factor (Atf7), putativeIntergenic Region
31,266,358A to G3.08AFUA_3G04310—AFUA_3G05320SnoRNA binding protein, putative—C2H2 finger domain protein, putativeIntergenic Region
52,069,698A to G3.03AFUA_5G08160—AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
21,781,938G to A2.99AFUA_2G06205—AFUA_2G06220Yippee family protein—Zinc knuckle domain proteinIntergenic Region
61,353,971T to C2.98AFUA_6G06350—AFUA_6G06360Proteasome subunit alpha type 3, putative—Mating alpha-pheromone (PpgA)Intergenic Region
41,363,615T to C2.96AFUA_4G04820—AFUA_4G05830C-4 methyl sterol oxidase (Erg25), putative—Methylthioribose-1-phosphate isomerase (Mri1)Intergenic Region
8635,137A to G2.94AFUA_8G02350Polyketide synthase, putativeMissense Variant
(Thr1206Ala)
32,316,978A to G2.91AFUA_3G09090RING finger domain proteinMissense Variant
(Glu298Gly)
62,508,121A to G2.91AFUA_6G10140—AFUA_6G10150C6 transcription factor, putative—Uncharacterized proteinIntergenic Region
22,074,852A to C2.89AFUA_2G08060Involucrin repeat proteinMissense Variant
(Lys779Thr)
22,080,579T to C2.89AFUA_2G08060Involucrin repeat proteinSynonymous Variant
(His2640His)
Table 4. Top 20 significant SNPs obtained from the GWAS that were associated with voriconazole resistance.
Table 4. Top 20 significant SNPs obtained from the GWAS that were associated with voriconazole resistance.
ChromosomePosition (bp)Change−log10(p-value)Gene IDAnnotationPredicted Effect
21,870,902G to A4.69AFUA_2G06330—AFUA_2G07340Ubiquitin carboxyl-terminal hydrolase—COP9 subunit 3, putativeIntergenic Region
1975,914G to A4.02AFUA_1G03370Uncharacterized proteinMissense Variant
(Ser174Asn)
8613,458G to A3.99AFUA_8G02290—AFUA_8G02300Uncharacterized protein—FMN-dependent dehydrogenase family proteinIntergenic Region
34,040,199T to C3.64AFUA_3G15350—AFUA_3G15380Short chain dehydrogenase family protein, putative—MFS multidrug transporter, putativeIntergenic Region
24,689,008C to T3.33AFUA_2G17600Conidial pigment polyketide synthase (Alb1)Synonymous Variant
(Val357Val)
5564,519A to C3.29AFUA_5G02210Uncharacterized proteinMissense Variant
(Met287Arg)
12,538,614A to C3.29AFUA_1G09780Stomatin family proteinMissense Variant
(Asp418Ala)
21,851,010G to A2.91AFUA_2G06330—AFUA_2G07340Ubiquitin carboxyl-terminal hydrolase—COP9 subunit 3, putativeIntergenic Region
8611,467C to A2.85AFUA_8G02280C6 transcription factor, putativeMissense Variant
(Glu79Asp)
22,087,757C to A2.79AFUA_2G08060Involucrin repeat proteinNon-coding Transcript Variant
63,648,516T to C2.78AFUA_6G143305-oxo-L-prolinase, putativeSynonymous Variant
(Glu131Glu)
11,337,273A to G2.77AFUA_1G04700—AFUA_1G04710Ras guanyl-nucleotide exchange factor (RasGEF), putative—Cytoplasmic tRNA 2-thiolation protein 1Intergenic Region
24,805,099C to T2.77AFUA_2G18070—AFUA_2G18100Neutral protease 2—Telomere-associated RecQ helicase, putativeIntergenic Region
2426,803C to T2.75AFUA_2G01740Sulfate transporter, putativeSynonymous Variant
(Ala141Ala)
3269,388G to T2.74AFUA_3G01150—AFUA_3G01160GPI anchored cell wall protein, putative—Choline monooxygenase, chloroplasticIntergenic Region
62,383,015C to T2.70AFUA_6G09745—AFUA_6G09760Uncharacterized protein—Cytochrome P450 monooxygenase, putativeIntergenic Region
2420,712T to C2.68AFUA_2G01710GPI anchored protein, putativeSynonymous Variant
(Ile294Ile)
71,182,007A to C2.68AFUA_7G05020—AFUA_7G05030Uncharacterized protein—Pectin lyase BIntergenic Region
5184,363G to A2.68AFUA_5G00650—AFUA_5G00660Uncharacterized protein—Uncharacterized proteinIntergenic Region
2441,695C to T2.67AFUA_2G01780Small nucleolar ribonucleoprotein complex subunit (Utp15), putativeSynonymous Variant
(Val184Val)
Table 5. Top 20 significant SNPs obtained from the second GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 5. Top 20 significant SNPs obtained from the second GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
12,538,614A to C5.37AFUA_1G09780Stomatin family proteinMissense Variant
(Asp418Ala)
21,845,323C to T4.14AFUA_2G06330-AFUA_2G07340Ubiquitin C-terminal hydrolase, putative—COP9 subunit 3, putativeIntergenic Region
8623,331G to T3.96AFUA_8G02330Endoglucanase, putativeNon-coding Transcript Variant
21,899,353C to T3.91AFUA_2G07430-AFUA_2G07440DDHD domain protein—Thioesterase family proteinIntergenic Region
43,737,973C to T3.69AFUA_4G14300-AFUA_4G14310Dynamin family GTPase, putative—APH domain-containing proteinIntergenic Region
52,063,521C to A3.52AFUA_5G08150ABC bile acid transporter, putativeMissense Variant
(His105Gln)
* 3267,884T to G3.36AFUA_3G01140-AFUA_3G01150Uncharacterized protein—GPI anchored cell wall protein, putativeIntergenic Region
52,069,483G to A3.27AFUA_5G08160-AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
52,069,698A to G3.13AFUA_5G08160-AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
* 32,389,222G to A3.12AFUA_3G09400-AFUA_3G09450MFS transporter (Hol1), putative—Alpha/beta fold family hydrolase, putativeIntergenic Region
63,054,001C to G2.99AFUA_6G12145-AFUA_6G12150Uncharacterized protein—BZIP transcription factor (Atf7), putativeIntergenic Region
* 32,414,011A to G2.87AFUA_3G0948015-hydroxyprostaglandin dehydrogenase (NAD(+))Synonymous Variant
(Ser60Ser)
31,953,910G to A2.84AFUA_3G07730-AFUA_3G07740Uncharacterized protein—Uncharacterized proteinIntergenic Region
* 6145,947T to C2.84AFUA_6G00570-AFUA_6G00580Uncharacterized protein—Ankyrin repeat proteinIntergenic Region
* 6262,795G to A2.79AFUA_6G01860MFS lactose permease, putativeMissense Variant
(Val106Met)
* 31,883,390C to A2.78AFUA_3G07510-AFUA_3G07520Uncharacterized protein—
Exo-beta-1,3-glucanase, putative
Intergenic Region
32,316,978A to G2.77AFUA_3G09090RING finger domain proteinMissense Variant
(Glu298Gly)
31,266,358A to G2.76AFUA_3G04310-AFUA_3G05320SnoRNA binding protein, putative—C2H2 finger domain protein, putativeIntergenic Region
* 13,885,980G to A2.74AFUA_1G14540Oxidoreductase, short-chain dehydrogenase/reductase familyNon-coding Transcript Variant
* 6734,136G to T2.72AFUA_6G03400-AFUA_6G03430Uncharacterized protein—C6 finger transcription factor (FsqA)Intergenic Region
Unique SNP sites are denoted by asterisks “*” (n = 8).
Table 6. Top 20 significant SNPs obtained from the second GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 6. Top 20 significant SNPs obtained from the second GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
1975,914G to A3.97AFUA_1G03370Uncharacterized proteinMissense Variant
(Ser174Asn)
21,870,902G to A3.83AFUA_2G06330-AFUA_2G07340Ubiquitin carboxyl-terminal hydrolase—COP9 subunit 3, putativeIntergenic Region
8613,458G to A3.71AFUA_8G02290-AFUA_8G02300Uncharacterized protein—FMN-dependent dehydrogenase family proteinIntergenic Region
* 7195,144A to G3.67AFUA_7G00740Protein kinase, putativeMissense Variant
(Ile188Val)
* 23,345,583A to G3.53AFUA_2G13030Phenylalanyl-tRNA synthetaseMissense Variant
(Asp343Gly)
* 2416,242C to T3.49AFUA_2G01700Carbon catabolite derepressing protein kinase (Snf1), putativeNon-coding Transcript Variant
63,648,516T to C3.41AFUA_6G143305-oxo-L-prolinase, putativeSynonymous Variant
(Glu131Glu)
62,383,015C to T3.29AFUA_6G09745-AFUA_6G09760Uncharacterized protein—Cytochrome P450 monooxygenase, putativeIntergenic Region
* 8237,297T to G3.28AFUA_8G01030-AFUA_8G01050Uncharacterized protein—Lipase/esterase, putativeIntergenic Region
* 8379,123C to T3.28AFUA_8G01480-AFUA_8G01490Potassium channel, putative—Endoglucanase, putativeIntergenic Region
* 4776,628A to G3.14AFUA_4G02800-AFUA_4G02805Haemolysin-III family protein—Asp hemolysin-like proteinIntergenic Region
24,689,008C to T3.12AFUA_2G17600Conidial pigment polyketide synthase (Alb1)Synonymous Variant
(Val357Val)
22,087,757C to A3.12AFUA_2G08060Involucrin repeat proteinNon-coding Transcript Variant
* 8292,607C to T2.99AFUA_8G01250GNAT family acetyltransferase, putativeMissense Variant
(Arg134Cys)
5564,519A to C2.94AFUA_5G02210Uncharacterized proteinMissense Variant
(Met287Arg)
* 71,019,801A to G2.94AFUA_7G04470-AFUA_7G04480Uncharacterized protein—DNA mismatch repair protein (Msh3)Intergenic Region
* 3341,035T to A2.92AFUA_3G01370-AFUA_3G01400MFS transporter, putative—ABC multidrug transporter, putativeIntergenic Region
* 3386,560T to C2.90AFUA_3G01520MFS multidrug transporter, putativeSynonymous Variant
(Val170Val)
* 81,631,284G to A2.87AFUA_8G06690-AFUA_8G06700Cytochrome P450 alkane hydroxylase—AnnexinIntergenic Region
* 62,940,890T to C2.86AFUA_6G11780-AFUA_6G11790Uncharacterized protein—Uncharacterized proteinIntergenic Region
Unique SNP sites are denoted by asterisks “*” (n = 12).
Table 7. Top 20 significant SNPs obtained from the third GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 7. Top 20 significant SNPs obtained from the third GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
8635,137A to G3.72AFUA_8G02350Polyketide synthase (PKS), putativeMissense Variant
(Thr1206Ala)
8623,331G to T3.50AFUA_8G02330Endoglucanase, putativeNon-coding Transcript Variant
52,069,698A to G3.26AFUA_5G08160-AFUA_5G08170Cyclin, putative—Autophagy-related protein 3 (Atg3)Intergenic Region
* 5419,750A to G3.11AFUA_5G01640-AFUA_5G01650Ankyrin repeat protein—bZIP transcription factor (JlbA), putativeIntergenic Region
12,538,614A to C3.06AFUA_1G09780Stomatin family proteinMissense Variant
(Asp418Ala)
* 8629,524G to T3.01AFUA_8G02340-AFUA_8G02350Uncharacterized protein—Polyketide synthase, putativeIntergenic Region
* 8576,158T to C2.96AFUA_8G02210-AFUA_8G02220Alpha-ketoglutarate-dependent taurine dioxygenase—Uncharacterized proteinIntergenic Region
* 52,131,740T to C2.94AFUA_5G08390Response regulator, putative (Ssk1)Synonymous Variant
(Lys532Lys)
43,737,973C to T2.83AFUA_4G14300-AFUA_4G14310Dynamin family GTPase, putative—APH domain-containing proteinIntergenic Region
* 42,539,714C to A2.77AFUA_4G09770Velvet domain-containing proteinSynonymous Variant
(Leu193Leu)
13,885,980G to A2.74AFUA_1G14540Oxidoreductase, short-chain dehydrogenase/reductase familyNon-coding Transcript Variant
* 31,256,445T to A2.74AFUA_3G04310-AFUA_3G05320SnoRNA binding protein, putative—
C2H2 finger domain protein, putative
Intergenic Region
* 62,141,290T to C2.72AFUA_6G09000-AFUA_6G09010PHD finger domain protein, putative—U1 snRNP splicing complex subunit (Luc7), putativeIntergenic Region
32,389,222G to A2.70AFUA_3G09400-AFUA_3G09450MFS transporter (Hol1), putative—Alpha/beta fold family hydrolase, putativeIntergenic Region
* 5810,835T to C2.67AFUA_5G03020-AFUA_5G030306 0S ribosomal protein L4, putative—
C6 transcription factor, putative
Intergenic Region
* 62,891,637A to G2.65AFUA_6G11620-AFUA_6G11630Formyltetrahydrofolate deformylase, putative—FAD-dependent isoamyl alcohol oxidase, putativeIntergenic Region
52,063,521C to A2.64AFUA_5G08150ABC bile acid transporter, putativeMissense Variant
(His105Gln)
* 81,069,676A to G2.63AFUA_8G04680Oxidoreductase, short-chain dehydrogenase/reductase family, putativeNon-coding Transcript Variant
* 11,585,001C to T2.62AFUA_1G00410C6 transcription factor, putativeIntragenic Variant
* 43,891,318A to C2.59AFUA_4G14751-AFUA_4G14770Uncharacterized protein—
Protostadienol synthase (HelA)
Intergenic Region
Unique SNP sites compared to the previous two GWAS analyses are denoted by asterisks “*” (n = 12).
Table 8. Top 20 significant SNPs obtained from the third GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 8. Top 20 significant SNPs obtained from the third GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
8613,458G to A2.99AFUA_8G02290-AFUA_8G02300Uncharacterized protein—FMN-dependent dehydrogenase family proteinIntergenic Region
* 53,732,385G to A2.59AFUA_5G14315Uncharacterized proteinSynonymous Variant
(Phe212Phe)
* 3246,050C to A2.59AFUA_3G01060-AFUA_3G01070Uncharacterized protein—Tyrosinase, putativeIntergenic Region
* 8388,274G to A2.58AFUA_8G01510-AFUA_8G01520Uncharacterized protein—PectinesteraseIntergenic Region
* 5794,519G to T2.57AFUA_5G02970LCCL domain proteinSynonymous Variant
(Thr24Thr)
* 42,494,977G to C2.55AFUA_4G09580Major allergen (Aspf2)Missense Variant
(Gly276Ala)
* 8293,836G to A2.52AFUA_8G01260Uncharacterized proteinSynonymous Variant
(Pro383Pro)
* 62,379,483T to C2.45AFUA_6G09745-AFUA_6G09760Uncharacterized protein—Cytochrome P450 monooxygenase, putativeIntergenic Region
* 62,424,223C to A2.44AFUA_6G09870C6 transcription factor, putativeMissense Variant
(Val360Phe)
* 8791,268A to G2.44AFUA_8G02870-AFUA_8G03870Uncharacterized protein—
Uncharacterized protein
Intergenic Region
* 62,480,554C to T2.43AFUA_6G10050-AFUA_6G10060Small oligopeptide transporter, OPT family—F-actin-capping protein subunit alphaIntergenic Region
* 21,785,216G to A2.42AFUA_2G06205-AFUA_2G06220Yippee family protein—Zinc knuckle domain proteinIntergenic Region
* 71,458,738C to G2.37AFUA_7G05960C2H2 finger domain protein, putativeMissense Variant
(Arg759Pro)
* 81,548,514C to T2.37AFUA_8G06410MFS multidrug transporter, putativeSynonymous Variant
(Arg17Arg)
* 52,131,740T to C2.37AFUA_5G08390Response regulator, putative (Ssk1)Synonymous Variant
(Lys532Lys)
* 21,774,354T to C2.35AFUA_2G06205-AFUA_2G06220Yippee family protein—Zinc knuckle domain proteinIntergenic Region
* 42,539,714C to A2.35AFUA_4G09770Velvet domain-containing proteinSynonymous Variant
(Leu193Leu)
* 21,787,001C to T2.35AFUA_2G06205-AFUA_2G06220Yippee family protein—Zinc knuckle domain proteinIntergenic Region
21,870,902G to A2.34AFUA_2G06330-AFUA_2G07340Ubiquitin carboxyl-terminal hydrolase—COP9 subunit 3, putativeIntergenic Region
* 14,762,609A to G2.32AFUA_1G17410Beta-glucosidase, putativeMissense Variant
(Val287Ala)
Unique SNP sites compared to the previous two GWAS analyses are denoted by asterisks “*” (n = 18).
Table 9. Top 20 significant SNPs obtained from the fourth GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 9. Top 20 significant SNPs obtained from the fourth GWAS associated with itraconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
12,538,614A to C4.86AFUA_1G09780Stomatin family proteinMissense Variant
(Asp418Ala)
21,845,323C to T3.21AFUA_2G06330—AFUA_2G07340Ubiquitin C-terminal hydrolase, putative—COP9 subunit 3, putativeIntergenic Region
8623,331G to T3.04AFUA_8G02330Endoglucanase, putativeNon-coding Transcript Variant
21,899,353C to T3.03AFUA_2G07430—AFUA_2G07440DDHD domain protein—Thioesterase family proteinIntergenic Region
52,063,521C to A2.91AFUA_5G08150ABC bile acid transporter, putativeMissense Variant
(His105Gln)
32,389,222G to A2.88AFUA_3G09400-AFUA_3G09450MFS transporter (Hol1), putative—Alpha/beta fold family hydrolase, putativeIntergenic Region
52,069,483G to A2.76AFUA_5G08160—AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
63,054,001C to G2.76AFUA_6G12145—AFUA_6G12150Uncharacterized protein—BZIP transcription factor (Atf7), putativeIntergenic Region
52,069,698A to G2.71AFUA_5G08160—AFUA_5G08170Cyclin, putative—
Autophagy-related protein 3 (Atg3)
Intergenic Region
3267,884T to G2.70AFUA_3G01140-AFUA_3G01150Uncharacterized protein—GPI anchored cell wall protein, putativeIntergenic Region
* 62,895,225T to C2.56AFUA_6G11620-AFUA_6G11630Formyltetrahydrofolate deformylase, putative—FAD-dependent isoamyl alcohol oxidase, putativeIntergenic Region
31,953,910G to A2.54AFUA_3G07730—AFUA_3G07740Uncharacterized protein—Uncharacterized proteinIntergenic Region
43,737,973C to T2.54AFUA_4G14300—AFUA_4G14310Dynamin family GTPase, putative—APH domain-containing proteinIntergenic Region
* 3228,628C to T2.54AFUA_3G00970-AFUA_3G00980Uncharacterized protein—MFS transporter Liz1/Seo1, putativeIntergenic Region
* 52,042,856G to A2.53AFUA_5G08050-AFUA_5G08060Aminopeptidase P, putative—Importin 13, putativeIntergenic Region
*3 2,456,111A to G2.53AFUA_3G09630-AFUA_3G09640Asparaginyl-tRNA synthetase Slm5, putative—Camp independent regulatory proteinIntergenic Region
* 3247,848G to A2.53AFUA_3G01060-AFUA_3G01070Uncharacterized protein—Tyrosinase, putativeIntergenic Region
* 3220,452G to A2.52AFUA_3G00930C6 transcription factor, putativeSynonymous Variant
(Ile259Ile)
32,408,041T to C2.50AFUA_3G09400—AFUA_3G09450MFS transporter (Hol1), putative—Alpha/beta fold family hydrolase, putativeIntergenic Region
32,414,011A to G2.48AFUA_3G0948015-hydroxyprostaglandin dehydrogenase (NAD(+))Synonymous Variant
(Ser60Ser)
Unique SNP sites compared to the previous three GWAS analyses are denoted by asterisks “*” (n = 6).
Table 10. Top 20 significant SNPs obtained from the fourth GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
Table 10. Top 20 significant SNPs obtained from the fourth GWAS associated with voriconazole resistance, arranged based on their −log10(p-values) from the highest to lowest.
ChromosomePosition (bp)Change−log10(p-Value)Gene IDAnnotationPredicted Effect
8613,458G to A3.82AFUA_8G02290-AFUA_8G02300Uncharacterized protein—FMN-dependent dehydrogenase family proteinIntergenic Region
8611,467C to A3.22AFUA_8G02280C6 transcription factor, putativeMissense Variant
(Glu79Asp)
3246,050C to A2.93AFUA_3G01060-AFUA_3G01070Uncharacterized protein—Tyrosinase, putativeIntergenic Region
* 412,352G to A2.80Chr Start—AFUA_4G00100Rhamnogalacturonase, putativeIntergenic Region
* 8641,537T to C2.79AFUA_8G02380-AFUA_8G02390FAD-dependent monooxygenase, putative—Uncharacterized proteinIntergenic Region
42,539,714C to A2.75AFUA_4G09770Velvet domain-containing proteinSynonymous Variant
(Leu193Leu)
* 8331,435C to A2.62AFUA_8G01340MFS sugar transporter, putativeMissense Variant
(Leu239Met)
* 5295,677C to T2.62AFUA_5G01180RAN small monomeric GTPase (Ran), putativeSynonymous Variant(Ser74Ser)
* 5256,650T to C2.61AFUA_5G01000Oxidoreductase, 2OG-Fe(II) oxygenase family, putativeMissense Variant
(Ser110Pro)
* 2417,623C to T2.60AFUA_2G01700Carbon catabolite derepressing protein kinase (Snf1), putativeMissense Variant
(Arg188Gln)
2426,803C to T2.59AFUA_2G01740Sulfate transporter, putativeSynonymous Variant
(Ala141Ala)
* 8503,790C to G2.59AFUA_8G01940C6 finger domain protein, putativeMissense Variant
(Pro261Arg)
2420,712T to C2.55AFUA_2G01710GPI anchored protein, putativeSynonymous Variant
(Ile294Ile)
* 11,138,713A to G2.53AFUA_1G00410C6 transcription factor, putativeIntragenic Variant
2441,695C to T2.53AFUA_2G01780Small nucleolar ribonucleoprotein complex subunit (Utp15), putativeSynonymous Variant
(Val184Val)
* 53,788,892C to T2.53AFUA_5G14610Carboxypeptidase Y, putativeMissense Variant
(Val254Met)
* 51,815,994A to G2.50AFUA_5G07300- AFUA_5G07310Electron transfer flavoprotein, beta subunit—DUF500 domain proteinIntergenic Region
* 13,306,670A to C2.50AFUA_1G12540TMEM1 family protein, putativeMissense Variant
(Phe879Cys)
* 41,285,247G to A2.50AFUA_4G04570Uncharacterized proteinNon-coding Transcript Variant
8388,274G to A2.50AFUA_8G01510-AFUA_8G01520Uncharacterized protein—PectinesteraseIntergenic Region
Unique SNP sites compared to the previous three GWAS analyses are denoted by asterisks “*” (n = 12).
Table 11. Additional non-synonymous SNPs found to be highly linked to the 46 SNP sites obtained by GWAS analyses for itraconazole.
Table 11. Additional non-synonymous SNPs found to be highly linked to the 46 SNP sites obtained by GWAS analyses for itraconazole.
ChromosomePositionGene IDPredicted Effect
(Amino Acid Substitution)
Description
22,079,605AFUA_2G08060Missense Variant (Ala2316Ser)Involucrin repeat protein
22,083,296AFUA_2G08060Missense Variant (Asn3546Ser)Involucrin repeat protein
22,086,695AFUA_2G08060Missense Variant (Val4679Ala)Involucrin repeat protein
3587,378AFUA_3G02360Missense Variant (Leu413Gln)Carboxylic ester hydrolase
31,604,491AFUA_3G06490Missense Variant (Gln531Arg)Uncharacterized protein
31,629,278AFUA_3G06570Missense Variant (Gln77Pro)Uncharacterized protein
31,693,467AFUA_3G06800Missense Variant (Arg615Thr)Uncharacterized protein
31,700,605AFUA_3G06820Missense Variant (Lys540Arg)Oxidoreductase, FAD-binding
32,132,951AFUA_3G08280Missense Variant (Glu28Lys)Cell cycle regulatory protein (Srw1), putative
32,155,356AFUA_3G08400Missense Variant (Glu393Lys)SNF2 family helicase/ATPase, putative
32,304,691AFUA_3G09040Missense Variant (Ser13Leu)Uncharacterized protein
32,311,362AFUA_3G09070Missense Variant (Ile406Thr)Carboxylesterase, putative
32,409,306AFUA_3G09450Missense Variant (Pro220Leu)Alpha/beta fold family hydrolase, putative
43,875,753AFUA_4G14712Missense Variant (Pro208Ser)C6 transcription factor, putative
62,583,985AFUA_6G10420Missense Variant (Gln309Glu)Uncharacterized protein
Table 12. Additional non-synonymous SNPs found to be highly linked to the 62 SNP sites obtained by GWAS analyses for voriconazole.
Table 12. Additional non-synonymous SNPs found to be highly linked to the 62 SNP sites obtained by GWAS analyses for voriconazole.
ChromosomePositionGene IDPredicted Effect
(Amino Acid Substitution)
Description
1976,070AFUA_1G03370Missense Variant (Ser226Leu)Uncharacterized protein
14,754,138AFUA_1G17380Missense Variant (Leu226Pro)3-oxoacyl-(Acyl-carrier-protein) reductase, putative
2437,241AFUA_2G01760Missense Variant (Thr1812Ala)NACHT domain protein
2541,777AFUA_2G02170Missense Variant (Ser67Pro)Nuclear condensin complex subunit (Smc4), putative
5205,924AFUA_5G00730Missense Variant (Val814Phe)H /K ATPase alpha subunit, putative
53,290,025AFUA_5G12670Missense Variant (Phe390Ser)Nucleoporin (Nup192), putative
63,252,789AFUA_6G12890Missense Variant (Arg878Gly)Vacuole-associated enzyme activator complex component (Vac14), putative
63,330,314AFUA_6G13180Missense Variant (Ala529Thr)CECR1 family adenosine deaminase, putative
71,457,904AFUA_7G05960Missense Variant (Arg1037Gln)C2H2 finger domain protein, putative
71,541,519AFUA_7G06290Missense Variant (Gln666Leu)NACHT domain protein, putative
8332,292AFUA_8G01340Missense Variant (Met524Ile)MFS sugar transporter, putative
Table 13. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests (n = 122).
Table 13. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests (n = 122).
ChromosomePosition (bp)Gene IDPredicted Effect
(Amino Acid Substitution)
Fisher’s Exact Test (p-Values), MIC ≥ 2 mg/LFisher’s Exact Test (p-Values), MIC ≥ 4 mg/L
ItraconazolePan-AzoleItraconazolePan-Azole
1976,070AFUA_1G03370Missense Variant (Ser226Leu)2.37 × 10−5 *8.10 × 10−6 *2.37 × 10−5 *5.48 × 10−6 *
14,754,138AFUA_1G17380Missense Variant (Leu226Pro)8.40 × 10−6 *5.57 × 10−5 *8.40 × 10−6 *9.21 × 10−5 *
32,304,691AFUA_3G09040Missense Variant (Ser13Leu)3.92 × 10−4 *3.82 × 10−33.92 × 10−4 *2.41 × 10−3
32,311,362AFUA_3G09070Missense Variant (Ile406Thr)3.74 × 10−4 *1.99 × 10−33.74 × 10−4 *2.41 × 10−3
32,409,306AFUA_3G09450Missense Variant (Pro220Leu)2.79 × 10−4 *7.05 × 10−5 *2.79 × 10−4 *1.64 × 10−4 *
* Statistically significant association between SNP and antifungal resistance.
Table 14. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests after removing the 21 strains with the L98H mutation in cyp51A (n = 101).
Table 14. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests after removing the 21 strains with the L98H mutation in cyp51A (n = 101).
ChromosomePosition (bp)Gene IDPredicted Effect
(Amino Acid Substitution)
Fisher’s Exact Test (p-Values), MIC ≥ 2 mg/LFisher’s Exact Test (p-Values), MIC ≥ 4 mg/L
ItraconazolePan-AzoleItraconazolePan-Azole
1976,070AFUA_1G03370Missense Variant (Ser226Leu)3.20 × 10−34.79 × 10−4 *3.20 × 10−31.84 × 10−4 *
14,754,138AFUA_1G17380Missense Variant (Leu226Pro)3.15 × 10−4 *2.12 × 10−33.15 × 10−4 *1.41 × 10−3
32,304,691AFUA_3G09040Missense Variant (Ser13Leu)1.06 × 10−5 *7.39 × 10−5 *1.06 × 10−5 *6.19 × 10−5 *
32,311,362AFUA_3G09070Missense Variant (Ile406Thr)5.07 × 10−6 *3.47 × 10−5 *5.07 × 10−6 *6.19 × 10−5 *
71,541,519AFUA_7G06290Missense Variant (Gln666Leu)1.08 × 10−34.35 × 10−4 *1.08 × 10−31.29 × 10−3
* Statistically significant association between SNP and antifungal resistance.
Table 15. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests after removing the 64 strains with the mutations in cyp51A (n = 58).
Table 15. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests after removing the 64 strains with the mutations in cyp51A (n = 58).
ChromosomePosition (bp)Gene IDPredicted Effect
(Amino Acid Substitution)
Fisher’s Exact Test (p-Values), MIC ≥ 2 mg/LFisher’s Exact Test (p-Values), MIC ≥ 4 mg/L
ItraconazolePan-AzoleItraconazolePan-Azole
14,754,138AFUA_1G17380Missense Variant (Leu226Pro)1.87 × 10−5 *1.25 × 10−4 *1.87 × 10−5 *2.91 × 10−4 *
32,304,691AFUA_3G09040Missense Variant (Ser13Leu)8.33 × 10−5 *2.33 × 10−5 *8.33 × 10−5 *3.10 × 10−5 *
32,311,362AFUA_3G09070Missense Variant (Ile406Thr)8.33 × 10−5 *2.33 × 10−5 *8.33 × 10−5 *3.10 × 10−5 *
71,541,519AFUA_7G06290Missense Variant (Gln666Leu)1.18 × 10−33.64 × 10−4 *1.18 × 10−35.03 × 10−4 *
* Statistically significant association between SNP and antifungal resistance.
Table 16. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests and strains in Clade II (n = 71).
Table 16. Highly linked significant SNP sites associated with triazole resistance determined using Fisher’s Exact tests and strains in Clade II (n = 71).
ChromosomePosition (bp)Gene IDPredicted Effect
(Amino Acid Substitution)
Fisher’s Exact Test (p-Values), MIC ≥ 2 mg/LFisher’s Exact Test (p-Values), MIC ≥ 4 mg/L
ItraconazolePan-AzoleItraconazolePan-Azole
14,754,138AFUA_1G17380Missense Variant (Leu226Pro)1.68 × 10−6 *4.81 × 10−5 *1.68 × 10−6 *4.83 × 10−5 *
32,304,691AFUA_3G09040Missense Variant (Ser13Leu)2.59 × 10−5 *8.51 × 10−5 *2.59 × 10−5 *1.73 × 10−5 *
32,311,362AFUA_3G09070Missense Variant (Ile406Thr)1.17 × 10−5 *3.48 × 10−5 *1.17 × 10−5 *1.73 × 10−5 *
71,541,519AFUA_7G06290Missense Variant (Gln666Leu)3.21 × 10−35.88 × 10−4*3.21 × 10−32.85 × 10−4*
* Statistically significant association between SNP and antifungal resistance.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fan, Y.; Wang, Y.; Korfanty, G.A.; Archer, M.; Xu, J. Genome-Wide Association Analysis for Triazole Resistance in Aspergillus fumigatus. Pathogens 2021, 10, 701. https://doi.org/10.3390/pathogens10060701

AMA Style

Fan Y, Wang Y, Korfanty GA, Archer M, Xu J. Genome-Wide Association Analysis for Triazole Resistance in Aspergillus fumigatus. Pathogens. 2021; 10(6):701. https://doi.org/10.3390/pathogens10060701

Chicago/Turabian Style

Fan, Yuying, Yue Wang, Gregory A. Korfanty, Meagan Archer, and Jianping Xu. 2021. "Genome-Wide Association Analysis for Triazole Resistance in Aspergillus fumigatus" Pathogens 10, no. 6: 701. https://doi.org/10.3390/pathogens10060701

APA Style

Fan, Y., Wang, Y., Korfanty, G. A., Archer, M., & Xu, J. (2021). Genome-Wide Association Analysis for Triazole Resistance in Aspergillus fumigatus. Pathogens, 10(6), 701. https://doi.org/10.3390/pathogens10060701

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop