Sixty-five Japanese A. fumigatus clinical strains were provided through the National Bio-Resource Project (NBRP), Japan (http://nbrp.jp/) (Supplementary Table 1). These samples originated from different patients, several different sources and infections, and 15 cities. Whole genome paired-end Illumina sequence data for an additional 11 A. fumigatus isolates that were previously sequenced and have ITCZ MIC data (Takahashi-Nakaguchi et al., 2015) (Supplementary Table 1) were downloaded from NCBI Sequence Read Archive (SRA) (Leinonen et al., 2011) using the SRA toolkit (https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=show&f=software&m=software&s=software).

Minimal inhibitory concentration (MIC) of ITCZ for each isolate was determined following the Clinical and Laboratory Standards Institute (CLSI) M38-A2 broth microdilution method (John, 2008). Before MIC calculations, each strain was cultured using a potato dextrose agar plate (Becton Dickinson, Sparks, MD, US) for 5 days at 30°C degrees to produce the fungal conidia. Harvested conidia were suspended in standard RPMI 1640 broth (pH = 7) (Sigma Aldrich, St. Louis, US-MO). For each isolate, 2.5 × 10⁴ conidia per l mL were incubated in RPMI 1640 broth (pH=7) with a range of ITCZ concentrations (8, 4, 2, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625 μg/ml) at 35°C for 48 h. MIC values represent the lowest ITCZ concentrations that completely inhibited growth.

Genomic DNA (gDNA) isolation was performed as previously described (Zhao et al., 2019). gDNA was directly isolated from conidia stocks using the MasterPure™ Yeast DNA Purification Kit (Lucigen/Epicenter) following the manufacturer's instructions, with several minor modifications. Conidia stocks were centrifuged at 14,000 RPM for 5 min to obtain a pellet. Next, 300 ml of yeast cell lysis solution was added to the pellet along with 0.4 ml of sterile 1.0 mm diameter silica beads. Lysis was carried out on a Biospec Mini-BeadBeater-8 at medium intensity for 8 min. One μl of RNase was added to the cell lysis solution and incubated at 65°C for 30 min. DNA isolation and purification were conducted according to the manufacturer's instructions for the remainder of the protocol. PCR-free 150-bp paired-end libraries were constructed and sequenced by Novogene (https://en.novogene.com/) on an Illumina NovaSeq 6000.

Raw reads were first deduplicated using tally (Davis et al., 2013) with the parameters “–with-quality” and “–pair-by-offset” to remove potential PCR duplication during library construction. Next, we used trim_galore v0.4.2 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) to trim residual adapter sequences from reads, and trim reads where quality score was below 30, with the parameters “–stringency 5” and “-q 30,” respectively. Trimmed reads shorter than 50 bp were then discarded using the option “–length 50.” Next, the deduplicated and trimmed read set was mapped to the A fumigatus Af293 reference genome (Nierman et al., 2005) using BWA-MEM v0.7.15 aligner (Li and Durbin, 2009). The resulting SAM files were converted into sorted BAM files using the “view” and “sort” functions in samtools 1.4.1 (Li et al., 2009).

Because A. fumigatus is haploid, we followed the best practice pipeline for “Germline short variant discovery” (Van Der Auwera et al., 2013) in Genome Analysis ToolKit (GATK) v4.0.6.0 (Mckenna et al., 2010). The function “HaplotypeCaller” was used to call short variants (SNPs and INDELs) with the sorted BAM file for each sample. The resulting g.vcf files of all 76 samples were then combined to generate a joint-called variant file using the function “GenotypeGVCFs.” Next only SNPs were extracted from the joint-called variant file using the function “SelectVariants.” To limit false positive variant calling, the function “VariantFiltration” was used to carry out “hard filtering” with the following parameters: “QD < 25.0 || FS > 5.0 || MQ < 55.0 || MQRankSum < −0.5 || ReadPosRankSum < −2.0 || SOR > 2.5”. 206,055 polymorphic loci were predicted after hard filtering.

To investigate the population structure of the A. fumigatus isolates we used a subset of population genetic informative SNPs. We used VCFtools v0.1.14 (Danecek et al., 2011) (http://vcftools.sourceforge.net/) with options “–maf 0.05 –max-missing 1 –thin 3500,” to filter the full set of SNPs and require a minor allele frequency ≥ 5%, no missing data across all samples, and at least 3.5 Kb distance between SNPs. 6,324 SNPs remained after filtering, and subsequent population structure analysis was conducted with this marker set. In addition, to test the consistency of population assignments with different number of SNPs, population structure analysis was conducted with a dense SNP set where thinning was not applied (59,433 SNP sites) and an additional thinned SNP set where markers were spaced apart by at least 35 Kb (756 SNPs).

To conduct population structure analysis, we first used the model-based program ADMIXTURE v1.3 (Alexander et al., 2009) for K = 1–10, where K indicates the number of populations. The 5-fold cross-validation (CV) procedure was calculated to find the most likely K with option “–cv = 5.” For each K the CV error was calculated and the K with lowest CV error indicated the most likely population number. Additionally, we used the non-model based population structure software DAPC (Jombart et al., 2010) in the “adegenet” package v2.1.2 (Jombart, 2008) in R v3.5.3 (Team, 2013) to the predict the number and assignment of individuals into populations. DAPC applies a Bayesian clustering method to identify populations without evolutionary models. The most likely number of populations was inferred by calculating the Bayesian Information Criterion (BIC) for each K.

Lastly, we also constructed a phylogenetic network with the alignment of 6,324 SNPs. The phylogenetic network was built using SplitsTree v4.14.4 (Huson and Bryant, 2006) with the neighbor joining method and 1,000 replicates for bootstrap analysis.

Genome Wide Association (GWA) analysis was conducted to identify genetic variants that were significantly correlated with ITCZ MIC. For GWA analysis, we filtered our complete set of SNPs with VCFtools to include SNPs with a minor allele frequency ≥5%, SNP sites with ≤ 10% missing data, and SNPs that were biallelic. This filtering procedure resulted in 68,853 SNPs that were used for GWA.

Two models were used to perform GWA between each of the 68,853 SNPs and ITCZ MIC. When ITCZ MIC data was treated as a quantitative trait (Supplementary Table 1), we used a linear mixed model with a genetic distance matrix for population structure correction in Tassel (Bradbury et al., 2007). GWA was also performed when ITCZ MIC was treated as a binary trait (MIC ≤ 0.5 = more sensitive, and MIC > 0.5 = less sensitive). In this GWA analysis, we used a mixed effect logistic model with an empirical covariance matrix as a population structure correction in RoadTrips (Thornton and Mcpeek, 2010). Quantile–quantile(Q-Q) plots were generated using the R package “qqman” (Turner, 2014) in order to evaluate potential p-value inflation. The potential functional effects of candidate SNPs were predicted using SnpEff v4.3t (Cingolani et al., 2012) with the A. fumigatus Af293 reference genome annotation.

To investigate the expression patterns of our candidate genes Afu2g02220 and Afu2g02140, we obtained FPKM values as well as fold-change and p-values for pairwise comparisons from FungiDB (https://fungidb.org/fungidb/) (Stajich et al., 2012) for oxidative stress, iron depletion, growth in blood and minimal media, and ITCZ exposure (Irmer et al., 2015; Kurucz et al., 2018).

A. fumigatus strain CEA10 was used as the genetic background for the deletion of Afu2g02220. The deletion was carried out using a clustered regularly interspaced short palindromic repeats (CRISPR)/Cas9-mediated protocol for gene editing, as previously described (Al Abdallah et al., 2017). Briefly, two Protospacer Adjacent Motif (PAM) sites, at both upstream and downstream of Afu2g02220, were selected using the EuPaGDT tool (Peng and Tarleton, 2015) and custom crRNAs were designed using the 20 base pairs of sequence immediately upstream of the PAM site. The crRNAs used are as follows: 5′ crRNA of Afu2g02220 = CTGTTATTTTCTTCGGGTCT and 3′ crRNA of Afu2g02220 = TGGACCAGGAAGAAACTGAG. Both crRNAs were purchased from IDT (Integrated DNA Technologies, Inc.). Complete guideRNAs (gRNAs) were then assembled in vitro using the custom designed crRNA coupled with a commercially acquired tracrRNA. The assembled gRNAs were then combined with commercially purchased Cas9 to form ribonucleoproteins for transformation, as previously described (Al Abdallah et al., 2017). Repair templates carrying a hygromycin resistance (HygR) cassette were PCR amplified to contain 40-basepair regions of microhomology on either side for homologous integration at the double strand DNA break induced by the Cas9 nuclease. Protoplast-mediated transformations were then carried out using the hygromycin repair templates and Cas-ribonucleoproteins for gene targeting. Homologous integrations were confirmed by PCR. The primers used are as follows:

Colony diameter was used as an estimate of growth rate to compare KO and WT strains in the presence and absence of ITCZ. For each strain, 10⁴ conidia were inoculated onto glucose minimal media (GMM) agar plates without ITCZ or with 0.15 μg/ml ITCZ. We used 0.15 μg/ml ITCZ because the parent strain (CEA10) is sensitive to ITCZ and did not grow in ITCZ concentrations ≥ 0.30 μg/ml. GMM was prepared as previously described (Shimizu and Keller, 2001). Colony diameter was measured with a digital caliper after 72 h at 37°C. Experiments were performed with ten replicates. To compare the growth rate between WT and KO strains, an ANOVA was performed between WT, ΔAfu2g02220-1, and ΔAfu2g02220-2 followed by a post-hoc Dunnett's test using the WT as the control group. Statistical analysis was conducted using JMP®, PRO 14 (SAS Institute Inc., Cary, NC, 1989–2019).

We conducted whole genome sequencing (WGS) for 65 isolates of A. fumigatus from Japan and analyzed them in combination with an additional 11 previously sequenced isolates (Takahashi-Nakaguchi et al., 2015). Deduplicated, quality trimmed, and adapter trimmed WGS data of the 76 isolates were used for joint SNP calling with GATK (Mckenna et al., 2010) and yielded 206,055 SNPs. To reduce the linkage between adjacent SNPs for population structure analysis, we subsampled SNPs so that they were separated by at least 3.5 kb, which yielded 6,324 SNPs. This subsampled dataset was used for population structure and phylogenetic analysis.

Population structure is a main confounding factor in GWA studies that can lead to false positive associations (Sul et al., 2018). Therefore, we investigated the population structure of the 76 A. fumigatus isolates using the model-based approach implemented in ADMIXTURE (Alexander et al., 2009), as well as a non-model approach where population structure is inferred using discriminant analysis of principal components (DAPC) (Jombart et al., 2010). In ADMIXTURE, cross-validation (CV) error was estimated for each K from K = 1–10. The CV error is calculated by systematically withholding data points, and the lowest value represents the best estimate of the number of ancestral populations (Alexander and Lange, 2011). Using this approach K = 4 was the most likely population number (Figure 1A). DAPC uses the Bayesian Information Criterion (BIC) to evaluate the optimal number of clusters (K). K = 4 was also the most likely scenario as evaluated by BIC in DAPC (Figure 1B). Population assignment was highly consistent when the entire SNP set was used, or when subsampled datasets consisting of 6,324 or 756 markers were used to limit linkage between markers (Supplementary Figure 1). At K = 4, DAPC assigned the 76 isolates into four distinct populations with no admixture, while ADMIXTURE assigned 30 of 76 individuals to more than one population. For population assignment, we placed isolates into their respective population based on their largest membership coefficient. Using this approach, only two isolates, IFM51978 and IFM61610 (Figure 1C, indicated by black arrows), were assigned into different populations between the two methods. Phylogenetic network analysis further supports the presence of four main populations and individual population assignment into these populations (Supplementary Figure 2).

The ITCZ MIC of all isolates ranged from 0.125 to 1 μg/ml (MIC_0.125 = 3, MIC_0.25 = 17, MIC_0.50 = 35, and MIC₁ = 21). For reference, ITCZ resistance is typically defined by MIC ≥ 4 (Tashiro et al., 2012). GWA was independently conducted when MIC data was treated as a quantitative trait, and when MIC was treated as a binary trait (“more-sensitive” = MIC < 0.5 or “less-sensitive” = MIC ≥ 0.5). Populations 1, 2, 3, and 4 had 1, 0, 2, and 0 individuals with MIC = 0.125, 5, 5, 4, and 3 individuals with MIC = 0.25, 10, 14, 10, and 1 individuals with MIC = 0.5, and 11, 7, 3, and 0 individuals with MIC = 1, respectively (Figure 1C).

We hypothesized that GWA would allow us to identify genes and/or genetic variants with minor contributions to ITCZ sensitivity. To test this hypothesis, we performed GWA with a set of 68,853 SNPs that have a minor allele frequency >5% and < 10% missing data, and the matched ITCZ MICs. Because these isolates have clear population structure (Figure 1) we used a mixed effect model GWA, which can reduce the inflated false-positive effect stemming from population structure (Yu et al., 2006; Price et al., 2010; Power et al., 2017) and has previously been applied in microbial GWA (Alam et al., 2014; Earle et al., 2016). We performed this mixed-model GWA with a covariance matrix as population correction for ITCZ MIC when treated as a quantitative trait (Figure 2A) and as a binary trait (Figure 2B) using Tassel 5 (Bradbury et al., 2007) and RoadTrips (Thornton and Mcpeek, 2010), respectively. We generated quantile-quantile (Q-Q) plots of expected vs. observed p-values to inspect p-value inflation, which could be the product of inadequate population structure correction. The Q-Q plots indicate that the distribution of p-values for both analyses are not inflated (Supplementary Figure 3).

We considered the 20 SNPs with the lowest p-values (lower 0.03 percentile) in each analysis as significant (Table 1). Of the 20 SNPs significantly associated with ITCZ MIC when MIC was treated as a quantitative trait (Figure 2A), five SNPs were located in genes (four in exons and one in an intron), 7 SNPs were located in 3′ UTR regions, two SNPs were located in 5′ UTR regions, and six SNPs were located in intergenic regions (Table 1). Of the four SNPs located in exons, one was synonymous (in Afu2g02220) while the remaining three SNPs were non-synonymous (in Afu2g02140, Afu4g00350, and Afu6g11980) (Table 1). Significant SNPs mapped to chromosomes 2 (N = 5), 3 (N = 8), 4 (N = 2), 5 (N = 2), 6 (N = 1), and 8 (N = 1) (Figure 2A).

Of the 20 SNPs significantly associated with ITCZ MIC when MIC was treated as a binary trait (Figure 2B), 12 SNPs were located in genes (11 in exons and one in an intron), two SNPs were located in 3′ UTR regions, 1 SNP was located in a 5′ UTR regions, and 4 SNPs were located in intergenic regions (Table 1). Of the 11 SNPs located in exons, six were synonymous (in Afu2g02220, Afu2g02140, Afu2g02290, Afu2g02170, and Afu2g01910) while the remaining five were non-synonymous (in Afu2g01930, Afu2g02140, and Afu2g01910) (Table 1). Interestingly, in this analysis, 19 of the 20 SNPs with lowest p-values were located to a 165 KB region on chromosome 2 (position 413,387 – 579,284) (Figure 2B).

Two significant SNPs overlapped between the quantitative trait and binary trait GWA analyses (Figure 2C). The SNP located in Afu2g02220 encodes a synonymous variant and had the ninth lowest and lowest p-values in the quantitative trait and binary trait analyses, respectively (Figures 2A,B). Afu2g02220 is annotated as a sterol 3-β-glucosyltransferase (Table 1). The SNP located in Afu2g02140 encodes a non-synonymous variant (A233G) and had the 10th lowest and seventh lowest p-values in the quantitative trait and binary trait analyses, respectively (Figures 2A,B). Afu2g02140 contains a CUE domain (as predicted by PFAM) (El-Gebali et al., 2019), which has been shown to bind to ubiquitin (Donaldson et al., 2003; Shih et al., 2003). For both Afu2g02220 and Afu2g02140, the major allele was associated with higher MIC values and the minor allele was absent in all isolates with ITCZ MIC = 1, and nearly absent in isolates with ITCZ MIC = 0.5 (Figure 2D, Supplementary Figure 4).

To investigate whether gene expression of Afu2g02220 and Afu2g02140 could be modulated by environmental stress, we analyzed A. fumigatus RNA-seq data publicly available on FungiDB (Stajich et al., 2012), during oxidative stress, iron depletion, ITCZ exposure, and growth in blood and minimal media (Irmer et al., 2015; Kurucz et al., 2018). Afu2g02220 was up-regulated during iron starvation (FPKM_control = 20.33, FPKM_FeStarvation = 32.70, and p-value = 5.7e⁻⁴), oxidative stress induced by H₂O₂ (FPKM_control = 20.33, FPKM_H2O2 = 30.61, and p-value = 2.7e⁻³), iron starvation + H₂O₂ (FPKM_control = 20.33, FPKM_{FeStarvation+H2O2} = 39.93, and p-value = 6.7e⁻²³), and during exposure to ITCZ in strain A1160 (FPKM_−ITCZ = 48.20, FPKM_+ITCZ = 67.37, and p-value = 1.7e⁻⁴) (Supplementary Figures 5A,B). Afu2g02140 was not significantly up-regulated during any condition, and expressed at lower levels across all conditions compared to Afu2g02220 (Supplementary Figure 5).

We chose to functionally examine the role of Afu2g02220 because (i) the SNP located in this gene had highly significant p-values in both GWA analyses (ii) Afu2g02220 has a predicted functional role in sterol metabolism, and ITCZ targets the ergosterol biosynthesis pathway and (iii) Afu2g02220 was up-regulated during ITCZ exposure (Supplementary Figure 5). Thus, we used an established CRISPR/Cas9 method (Al Abdallah et al., 2017) to knockout (KO) Afu2g02220 by replacing it with the indicator gene hygromycin B phosphotransferase (hygR) in the A. fumigatus CEA10 genetic background (Figure 3A). We generated two independent KOs of Afu2g02220 which we validated by via PCR (Figure 3B).

To test the effect of Afu2g02220 on ITCZ sensitivity, we grew the wild type (WT) and ΔAfu2g02220 strains in the presence of 0.15 μg/ml of ITCZ and measured colony diameter after 72 h of incubation at 37°C. We observed a qualitative reduction in conidia production in KO strains (Supplementary Figure 6). In minimal media without ITCZ ΔAfu2g02220-1 and ΔAfu2g02220-2 growth rates did not significantly differ from the WT (ΔAfu2g02220-1 = 45.016 ± 0.027 mm, ΔAfu2g02220-2 = 45.018 ± 0.030 mm, WT = 44.994 ± 0.024 mm) (Figure 3C). This result suggests that the background growth rate of ΔAfu2g02220 is not impacted by the gene deletion. However, at ITCZ concentrations of 0.15 μg/ml we observed a minor but consistent reduction in growth in KO strains compared to WT (ΔAfu2g02220-1 = 18.594 ± 0.105 mm, ΔAfu2g02220-1 = 18.615 ± 0.022 mm, WT = 19.239 ± 0.021 mm) (p-value = 2e⁻¹⁶ for both KOs) (Figure 3D). These results suggest that Afu2g02220 plays a minor role in ITCZ sensitivity.