Shu Zhao
Wenbo Ge
Akira Watanabe4
Jarrod R. Fortwendel
John G. Gibbons- 1Molecular and Cellular Biology Graduate Program, University of Massachusetts, Amherst, MA, United States
- 2Department of Food Science, University of Massachusetts, Amherst, MA, United States
- 3Department of Clinical Pharmacy and Translational Science, University of Tennessee Health Science Center, Memphis, TN, United States
- 4Division of Clinical Research, Medical Mycology Research Center, Chiba University, Chiba, Japan
- 5Organismic and Evolutionary Biology Graduate Program, University of Massachusetts, Amherst, MA, United States
Aspergillus fumigatus is a potentially lethal opportunistic pathogen that infects over ~200,000 people and causes ~100,000 deaths per year globally. Treating A. fumigatus infections is particularly challenging because of the recent emergence of azole-resistance. The majority of studies focusing on the molecular mechanisms underlying azole resistance have examined azole-resistant isolates. However, isolates that are susceptible to azoles also display variation in their sensitivity, presenting a unique opportunity to identify genes contributing to azole sensitivity. Here, we used genome-wide association (GWA) analysis to identify loci involved in azole sensitivity by analyzing the association between 68,853 SNPs and itraconazole (ITCZ) minimum inhibitory concentration (MIC) in 76 clinical isolates of A. fumigatus from Japan. Population structure analysis suggests the presence of four distinct populations, with ITCZ MICs distributed relatively evenly across populations. We independently conducted GWA when treating ITCZ MIC as a quantitative trait and a binary trait, and identified two SNPs with strong associations in both analyses. These SNPs fell within the coding regions of Afu2g02220 and Afu2g02140. We functionally validated Afu2g02220 by knocking it out using a CRISPR/Cas9 approach, because orthologs of this gene are involved in sterol modification and ITCZ targets the ergosterol biosynthesis pathway. Knockout strains displayed no difference in growth compared to the parent strain in minimal media, yet a minor but consistent inhibition of growth in the presence of 0.15 μg/ml ITCZ. Our results suggest that GWA paired with efficient gene deletion is a powerful and unbiased strategy for identifying the genetic basis of complex traits in A. fumigatus.
Introduction
Fungal infections result in more global deaths per year than deaths from tuberculosis or malaria (Brown et al., 2012). Aspergillus fumigatus is one of the most deadly fungal pathogens and results in more than 100,000 deaths per year (Brown et al., 2012). Invasive aspergillosis (IA) is the most severe infection caused by A. fumigatus and occurs when fungal growth, most commonly originating in the lung, disseminates to other parts of the body via the bloodstream (Latge, 1999). A. fumigatus is an opportunistic pathogen primarily affecting immunocompromised individuals, and unfortunately, infections have become more common due to the increased usage of immunosuppressive drugs to treat autoimmune disorders and to increase the success of organ transplantation surgery (Robinett et al., 2013; Neofytos et al., 2018; Latge and Chamilos, 2019). Even when aggressively treated with first and second-line antifungal medication, mortality rates can exceed 50% in IA patients (Latge, 1999; Lin et al., 2001). The relatively rapid emergence of A. fumigatus antifungal resistance has made treatment of infections particularly challenging.
Antifungal drugs target components that distinguish fungal cells from mammalian cells, including the fungal cell wall as well as unique components of the fungal cell membrane. For example, the echinocandins target β 1,3 glucan, the most abundant polysaccharide in the fungal cell wall, while amphotericin B (a polyene class of antifungal drug) and triazoles (an azole class of antifungal drugs) target ergosterol (Latge et al., 2017). Ergosterol plays an essential functional role in regulating cell membrane permeability and fluidity. Triazoles, such as itraconazole (ITCZ) and voriconazole, are the most common first-line treatment for A. fumigatus infections, and target the lanosterol demethylase enzymes (Cyp51A and Cyp51B in A. fumigatus) which are directly involved in the biosynthesis of ergosterol (Alcazar-Fuoli and Mellado, 2013; Revie et al., 2018). Blocking Cyp51A and Cyp51B results in the accumulation of a toxic sterol intermediate that causes severe membrane stress, impairment of growth, and cell death (Revie et al., 2018).
Strains of A. fumigatus have gained resistance to triazoles through mutations in both the coding and regulatory regions of cyp51A, and through cyp51A independent mechanisms (Garcia-Rubio et al., 2017). The three amino acid positions that are commonly found with non-synonymous mutations in cyp51A in azole resistance strains are 54, 220, and 448 (Garcia-Rubio et al., 2017). Protein structure modeling suggests that these mutations disrupt the binding efficiency of azoles to Cyp51A (Fraczek et al., 2011; Warrilow et al., 2015). Increased expression of cyp51A through a combination of a promoter region repeat and the L98H point mutation can also confer azole resistance (Mellado et al., 2007). Additionally, several transcription factors (e.g., srbA (Hagiwara et al., 2016), hapE (Camps et al., 2012), atrR (Paul et al., 2019), transporters [e.g., cdr1B (Fraczek et al., 2013), atrF (Meneau et al., 2016), various ABC transporters (Moye-Rowley, 2015) etc.], and other functional groups of genes (e.g., genes involved in calcium signaling, iron balance, signaling pathways, and the Hsp90-calcineurin pathway) have been implicated in azole resistance or susceptibility (Chen et al., 2020).
The numerous genes identified in azole resistance other than cyp51A (Garcia-Rubio et al., 2017) suggests that additional genes with additive minor effects likely play a role in fine-scale differences in azole sensitivity and resistance. Historically, most genes involved in azole resistance in A. fumigatus were discovered through a candidate gene approach (Garcia-Rubio et al., 2017), or through gene expression differences during exposure to azoles (da Silva Ferreira et al., 2006; Hokken et al., 2019). However, candidate gene methods are biased toward genes and pathways of biological interest. Alternatively, genome-wide association (GWA) studies offer a powerful and versatile approach to identify genetic variants that contribute to complex traits, such as A. fumigatus ITCZ sensitivity. In GWA, thousands to millions of high-density genetic variants are tested for a statistical association between each variant and a phenotype of interest (Gibson, 2018). Microbial GWAS methods have recently been developed (Read and Massey, 2014; Chen and Shapiro, 2015; Power et al., 2017), and has been used in other fungal species. For instance, GWA has been used to identify genes and variants associated with virulence in Heterobasidion annosum (Dalman et al., 2013), Saccharomyces cerevisiae (Muller et al., 2011), and Parastagonospora nodorum (Gao et al., 2016), fungal communication in Neurospora crassa (Palma-Guerrero et al., 2013), aggressiveness in Fusarium graminearum (Talas et al., 2016), and Zymoseptoria tritici (Hartmann et al., 2017). Here, we hypothesized that GWA could be applied in A. fumigatus to identify genes with minor effects on ITCZ sensitivity. We performed GWA in 76 non-resistant clinical isolates of A. fumigatus and identified a gene that contributes to fine-scale ITCZ sensitivity. More broadly, we demonstrate that GWA in combination with gene disruption is a useful tool for investigating medically relevant traits in A. fumigatus.
Materials and Methods
Japanese A. fumigatus Clinical Isolates
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).
Minimum Inhibitory Concentration Testing
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 × 104 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.
DNA Extraction and Illumina Whole-Genome Sequencing
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.
Quality Control and Sequence Read Mapping
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).
SNP Genotyping
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.
Population Structure of A. fumigatus Isolates
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 Analysis for Itraconazole Sensitivity
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.
RNA-Seq Based Expression Data for Afu2g02220 and Afu2g02140
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).
Gene Deletion of Afu2g02220 in A. fumigatus CEA10
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:
Afu2g02020 KO Forward Screening Primer (P1): GGATGCGTTGTTCCTGTGCG
Afu2g02220 KO Reverse Screening Primer (P2): AACGAGGGCTGGAGTGCC
Common HygR Reverse Screening Primer (P3): ACACCCAATACGCCGGCC
Comparison of ITCZ Sensitivity Between WT and KO Strains
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, 104 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).
Results
Population Structure of Clinical A. fumigatus Isolates From Japan
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).
Figure 1. Population structure and ITCZ sensitivity of the 76 Japanese clinical Aspergillus fumigatus Isolates. (A) The optimal number of genetic clusters (K, X-axis) inferred by ADMIXTURE using cross validation procedure (CV error, Y-axis). (B) The optimal number of genetic clusters (K, X-axis) inferred by DAPC using Bayesian Information Criterion (BIC, Y-axis). (C) Membership coefficients (Y-axis) for each of the 76 isolates (X-axis) from ADMIXTURE and DAPC for K = 4. The two black arrows indicate the two isolates that are assigned into different clusters by ADMIXTURE and DAPC. Population 1, 2, 3, and 4 are colored as blue, red, yellow, and gray, respectively. Binary ITCZ MIC assignment and quantitative ITCZ MIC values are depicted in the upper and lower panels below the membership coefficient plots, respectively. For binary ITCZ MIC, individuals are coded as either more-sensitive (MIC < 0.5, gray) or less-sensitive (MIC ≥ 0.5, yellow).
Itraconazole Minimum Inhibitory Concentration
The ITCZ MIC of all isolates ranged from 0.125 to 1 μg/ml (MIC0.125 = 3, MIC0.25 = 17, MIC0.50 = 35, and MIC1 = 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).
Genome-Wide Association of Itraconazole Sensitivity in A. fumigatus
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).
Figure 2. Genome-wide association (GWA) for itraconazole (ITCZ) sensitivity. GWA for ITCZ sensitivity when MIC data is treated as a quantitative trait (A) or as a binary trait (B). For binary characterization of ITCZ sensitivity, MIC < 0.5 = more-sensitive and MIC ≥ 0.5 = less-sensitive. The genomic location of the 68,853 SNPs used for GWA is depicted on the X-axis, while the –log (P-values) are depicted on Y-axis. The dotted gray horizontal line represents the cutoff line at the 20th lowest p-value. Afu2g02140 and Afu2g02220 were within the 20 SNPs with the strongest associations in both analyses and are labeled on each plot. (C) Venn diagram of the 20 SNPs most strongly associated with ITCZ MIC that overlap genes when data is treated as a quantitative trait (blue circle) and a binary trait (red circle). (D) Allele frequency of the SNP at Chr2:561,450 that falls within Afu2g02220 (Y-axis) by ITCZ MICs (X-axis).
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).
Table 1. Characterization of SNPs associated with ITCZ sensitivity.
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).
Expression of Afu2g02220 and Afu2g02140 From Existing RNA-Seq Experiments
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 (FPKMcontrol = 20.33, FPKMFeStarvation = 32.70, and p-value = 5.7e−4), oxidative stress induced by H2O2 (FPKMcontrol = 20.33, FPKMH2O2 = 30.61, and p-value = 2.7e−3), iron starvation + H2O2 (FPKMcontrol = 20.33, FPKMFeStarvation+H2O2 = 39.93, and p-value = 6.7e−23), and during exposure to ITCZ in strain A1160 (FPKM−ITCZ = 48.20, FPKM+ITCZ = 67.37, and p-value = 1.7e−4) (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).
Validation of a GWA Candidate Gene via CRISPR/Cas9 Gene Deletion
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).
Figure 3. Deletion of Afu2g02220 impairs growth in the presence of itraconazole (ITCZ). (A) Schematic of Afu2g02220 gene deletion via CRISPR/Cas9. The blue box with arrow in the upper panel represents Afu2g02220 in the parent CEA10 genome (wild type, WT), while the gray box in the lower panel represents the indicator gene HygR that replaced Afu2g02220 in ΔAfu2g02220 strains. The two black arrows on the flanking region of the locus indicate the forward primer (P1) and reverse primer (P2) used for PCR validation. The WT amplicon size is ~4.4 Kb, while the HygR gene replacement amplicon is ~3.8 Kb. (B) Validation of Afu2g02220 gene replacement via PCR. Lanes “M,” “WT,” “-1,” and “-2” indicate ladder, PCR product from WT and PCR product from the two independent knockout strains, respectively. Boxplots for colony diameter of WT and ΔAfu2g02220 strains grown at 37°C for 72 h on minimal media (C) and minimal media with 0.15 μg/ml ITCZ (D). Measurements were collected for 10 biological replicates for each experiment. Dunnett's test p-values indicate a significant reduction in growth in the KOs compared to the WT.
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−16 for both KOs) (Figure 3D). These results suggest that Afu2g02220 plays a minor role in ITCZ sensitivity.
Discussion
Here, we analyzed the association between SNP allele frequency and ITCZ MIC data from 76 Japanese clinical isolates of A. fumigatus to identify loci involved in ITCZ sensitivity. MIC values fell within a relatively tight range of 0.125–1 μg/ml [for reference, ITCZ resistant strains are defined by MICs ≥ 4 μg/ml (Tashiro et al., 2012)]. We reasoned that GWA could be a feasible tool to identify loci that contribute to the small differences in ITCZ MIC we observed across these clinical isolates. We identified several candidate SNPs and loci associated with ITCZ sensitivity, and validated the function of the top candidate by knocking it out using a CRISPR/Cas9 based approach.
We identified a synonymous variant in Afu2g02220 that showed highly significant associations with ITCZ sensitivity across GWA analyses with different underlying statistical models (Figure 2). Synonymous mutations can be functional through their (i) effect on cis-regulatory regions (e.g., splice sites or miRNA and exonic transcription factor binding sites), (ii) alteration of mRNA structure, or (iii) influence on translation speed (e.g., codon usage) (Hunt et al., 2014). Determining the mechanism by which this variant alters phenotype would require extensive in silico and in vitro experimentation. Afu2g02220 encodes a predicted sterol glycosyltransferase. This enzyme biosynthesizes sterol glucosides, which make up the common eukaryotic membrane bound lipids. Orthologs of Afu2g02220 from the ascomycete yeasts Saccharomyces cerevisiae (Atg26), Candida albicans, Pichia pastoris, as well as the amoeba Dictyostelium discoideum can use various sterols, including ergosterol, as sugar acceptors (Warnecke et al., 1999). In S. cerevisiae, Atg26 can directly bind to and glycosylate ergosterol, which yields ergosterol-glucoside (Gallego et al., 2010). In S. cerevisiae ΔAtg26 did not impair growth when cultured in complex or minimal media, low or elevated temperatures, varying osmotic stress conditions, or in the presence of nystatin, an antifungal drug that binds to ergosterol (Warnecke et al., 1999). Similarly, we did not observe a difference in growth rate between ΔAfu2g02220 and the WT when grown in minimal media (Figure 3C).
In addition to its role in sterol modification, Afu2g02220 may also have additional functions related to autophagy (Kikuma et al., 2017). Orthologs of Afu2g02220 in Pichia pastoris (PpAtg26) (Oku et al., 2003), Colletotrichum orbiculare (CoAtg26) (Asakura et al., 2009) and Aspergillus oryzae (AoAtg26) (Kikuma et al., 2017) are required for autophagy. In A. oryzae, ΔAoAtg26 shows deficiency in degradation of peroxisomes, mitochondria, and nuclei and localizes to vacuoles (Kikuma et al., 2017). ΔAoAtg26 also shows reductions in conidiation and impairment of aerial hyphae formation (Kikuma et al., 2017). Similarly, we observed a reduction in condition in ΔAfu2g02220 compared to the WT (Supplementary Figure 6).
The fungal cell wall is rigid but also dynamic in order to respond to environmental stress. Because Afu2g02220 may directly interact with ergosterol, we hypothesized that environmental stress could alter the expression of Afu2g02220. We analyzed A. fumigatus RNA-seq data during growth under iron depletion, oxidative stress, ITCZ exposure and growth in blood and minimal media (Stajich et al., 2012). We found that Afu2g02220 expression was significantly up-regulated during oxidative stress, iron depletion and ITCZ exposure (Supplementary Figure 5). However, other studies examining gene expression (da Silva Ferreira et al., 2006; Hokken et al., 2019) or protein abundance (Amarsaikhan et al., 2017) during exposure to ITCZ and voriconazole (da Silva Ferreira et al., 2006) (another triazole with the same mechanism of action as ITCZ) did not observe differential abundance of the Afu2g02220 transcript or protein. Additional experiments are necessary to determine the precise role of Afu2g02220 in stress response and ITCZ sensitivity.
Previously, Palma-Guerrero et al. (2013) used a similar approach to identify NCU04379 as a gene that contributes to fungal communication in N. crassa. This study used RNA-seq data to identify genetic variants, Fisher's exact tests to perform GWA in a closely related group of 112 isolates, and existing deletion mutants generated by the Neurospora Genome Project (Colot et al., 2006; Dunlap et al., 2007) to validate the involvement of NCU04379 in cellular communication during germling fusion. A study in S. cerevisiae used a mixed linear model to identify correlations between genotype and tolerance to hydrolysate toxins, and used homologous recombination to knockout candidate genes in two independent genetic backgrounds (Sardi et al., 2018). Interestingly, eight of 14 gene knockouts had a significant effect on phenotype in one, but not both genetic backgrounds, suggesting that the network of genes contributing to hydrolysate toxins tolerance likely differs between genetic backgrounds. The results of these studies, and of our own, broadly suggest that GWA in combination with an efficient gene disruption technique is a powerful and unbiased approach for identifying the genetic basis of polygenic phenotypes in fungal systems.
Data Availability Statement
Raw whole-genome Illumina data for the 65 isolates are available through NCBI BioProject PRJNA638646 and the 11 previously sequenced isolates by Takahashi-Nakaguchi et al. (2015) through NCBI BioProject PRJDB1541.
Author Contributions
SZ and JG designed the study and analyzed the data. AW determined itraconazole MIC and provided A. fumigatus isolates. WG and JF conducted CRISPR and growth rate experiments. All authors contributed to writing the manuscript.
Funding
This research was supported by grant R21AI137485 from the National Institutes of Health and National Institutes of Allergy and Infectious Diseases (NIAID) to JG which supports JG and SZ. JF and WG are supported by NIAID R01AI143197 to JF. AW was supported by the Japan Agency for Medical Research and Development (AMED) under Grant Number 20jm0110015.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
This manuscript has been released as a pre-print at BioRxiv, (Zhao et al., 2020). Computational analysis was conducted on the Massachusetts Green High Performance Computing Center (MGHPCC).
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffunb.2020.617338/full#supplementary-material
References
Al Abdallah, Q., Ge, W., and Fortwendel, J. R. (2017). A simple and universal system for gene manipulation in Aspergillus fumigatus: in vitro-assembled Cas9-guide RNA ribonucleoproteins coupled with microhomology repair templates. Msphere 2:e00446–17. doi: 10.1128/mSphere.00446-17
Alam, M. T., Petit Iii, R. A., Crispell, E. K., Thornton, T. A., Conneely, K. N., Jiang, Y., et al. (2014). Dissecting vancomycin-intermediate resistance in Staphylococcus aureus using genome-wide association. Genome Biol. Evol. 6, 1174–1185. doi: 10.1093/gbe/evu092
Alcazar-Fuoli, L., and Mellado, E. (2013). Ergosterol biosynthesis in Aspergillus fumigatus: its relevance as an antifungal target and role in antifungal drug resistance. Front. Microbiol. 3:439. doi: 10.3389/fmicb.2012.00439
Alexander, D. H., and Lange, K. (2011). Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12:246. doi: 10.1186/1471-2105-12-246
Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. doi: 10.1101/gr.094052.109
Amarsaikhan, N., Albrecht-Eckardt, D., Sasse, C., Braus, G. H., Ogel, Z. B., and Kniemeyer, O. (2017). Proteomic profiling of the antifungal drug response of Aspergillus fumigatus to voriconazole. Int. J. Med. Microbiol. 307, 398–408. doi: 10.1016/j.ijmm.2017.07.011
Asakura, M., Ninomiya, S., Sugimoto, M., Oku, M., Yamashita, S.-I., Okuno, T., et al. (2009). Atg26-mediated pexophagy is required for host invasion by the plant pathogenic fungus Colletotrichum orbiculare. Plant Cell 21, 1291–1304. doi: 10.1105/tpc.108.060996
Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308
Brown, G. D., Denning, D. W., Gow, N. A., Levitz, S. M., Netea, M. G., and White, T. C. (2012). Hidden killers: human fungal infections. Sci. Transl. Med. 4:165rv113. doi: 10.1126/scitranslmed.3004404
Camps, S. M., Dutilh, B. E., Arendrup, M. C., Rijs, A. J., Snelders, E., Huynen, M. A., et al. (2012). Discovery of a HapE mutation that causes azole resistance in Aspergillus fumigatus through whole genome sequencing and sexual crossing. PLoS ONE 7:e50034. doi: 10.1371/journal.pone.0050034
Chen, P., Liu, J., Zeng, M., and Sang, H. (2020). Exploring the molecular mechanism of azole resistance in Aspergillus fumigatus. J. Mycol. Med. 30:100915. doi: 10.1016/j.mycmed.2019.100915
Chen, P. E., and Shapiro, B. J. (2015). The advent of genome-wide association studies for bacteria. Curr. Opin. Microbiol. 25, 17–24. doi: 10.1016/j.mib.2015.03.002
Cingolani, P., Platts, A., Wang Le, L., Coon, M., Nguyen, T., Wang, L., et al. (2012). 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 6, 80–92. doi: 10.4161/fly.19695
Colot, H. V., Park, G., Turner, G. E., Ringelberg, C., Crew, C. M., Litvinkova, L., et al. (2006). A high-throughput gene knockout procedure for Neurospora reveals functions for multiple transcription factors. Proc. Natl. Acad. Sci. U. S. A. 103, 10352–10357. doi: 10.1073/pnas.0601456103
da Silva Ferreira, M. E., Malavazi, I., Savoldi, M., Brakhage, A. A., Goldman, M. H., Kim, H. S., et al. (2006). Transcriptome analysis of Aspergillus fumigatus exposed to voriconazole. Curr. Genet. 50, 32–44. doi: 10.1007/s00294-006-0073-2
Dalman, K., Himmelstrand, K., Olson, A., Lind, M., Brandstrom-Durling, M., and Stenlid, J. (2013). A genome-wide association study identifies genomic regions for virulence in the non-model organism Heterobasidion annosum s.s. PLoS ONE 8:e53525. doi: 10.1371/journal.pone.0053525
Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., Depristo, M. A., et al. (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. doi: 10.1093/bioinformatics/btr330
Davis, M. P., van Dongen, S., Abreu-Goodger, C., Bartonicek, N., and Enright, A. J. (2013). Kraken: a set of tools for quality control and analysis of high-throughput sequence data. Methods 63, 41–49. doi: 10.1016/j.ymeth.2013.06.027
Donaldson, K. M., Yin, H., Gekakis, N., Supek, F., and Joazeiro, C. A. (2003). Ubiquitin signals protein trafficking via interaction with a novel ubiquitin binding domain in the membrane fusion regulator, Vps9p. Curr. Biol. 13, 258–262. doi: 10.1016/S0960-9822(03)00043-5
Dunlap, J. C., Borkovich, K. A., Henn, M. R., Turner, G. E., Sachs, M. S., Glass, N. L., et al. (2007). Enabling a community to dissect an organism: overview of the Neurospora functional genomics project. Adv. Genet. 57, 49–96. doi: 10.1016/S0065-2660(06)57002-6
Earle, S. G., Wu, C.-H., Charlesworth, J., Stoesser, N., Gordon, N. C., Walker, T. M., et al. (2016). Identifying lineage effects when controlling for population structure improves power in bacterial association studies. Nat. Microbiol. 1:16041. doi: 10.1038/nmicrobiol.2016.41
El-Gebali, S., Mistry, J., Bateman, A., Eddy, S. R., Luciani, A., Potter, S. C., et al. (2019). The Pfam protein families database in 2019. Nucleic Acids Res. 47, D427–D432. doi: 10.1093/nar/gky995
Fraczek, M. G., Bromley, M., and Bowyer, P. (2011). An improved model of the Aspergillus fumigatus CYP51A protein. Antimicrob. Agents Chemother. 55, 2483–2486. doi: 10.1128/AAC.01651-10
Fraczek, M. G., Bromley, M., Buied, A., Moore, C. B., Rajendran, R., Rautemaa, R., et al. (2013). The cdr1B efflux transporter is associated with non-cyp51a-mediated itraconazole resistance in Aspergillus fumigatus. J. Antimicrob. Chemother. 68, 1486–1496. doi: 10.1093/jac/dkt075
Gallego, O., Betts, M. J., Gvozdenovic-Jeremic, J., Maeda, K., Matetzki, C., Aguilar-Gurrieri, C., et al. (2010). A systematic screen for protein–lipid interactions in Saccharomyces cerevisiae. Mol. Syst. Biol. 6:430. doi: 10.1038/msb.2010.87
Gao, Y., Liu, Z., Faris, J. D., Richards, J., Brueggeman, R. S., Li, X., et al. (2016). Validation of genome-wide association studies as a tool to identify virulence factors in Parastagonospora nodorum. Phytopathology 106, 1177–1185. doi: 10.1094/PHYTO-02-16-0113-FI
Garcia-Rubio, R., Cuenca-Estrella, M., and Mellado, E. (2017). Triazole resistance in aspergillus species: an emerging problem. Drugs 77, 599–613. doi: 10.1007/s40265-017-0714-4
Gibson, G. (2018). Population genetics and GWAS: a primer. PLoS Biol. 16:e2005485. doi: 10.1371/journal.pbio.2005485
Hagiwara, D., Watanabe, A., and Kamei, K. (2016). Sensitisation of an azole-resistant Aspergillus fumigatus strain containing the Cyp51A-related mutation by deleting the SrbA gene. Sci. Rep. 6:38833. doi: 10.1038/srep38833
Hartmann, F. E., Sanchez-Vallet, A., Mcdonald, B. A., and Croll, D. (2017). A fungal wheat pathogen evolved host specialization by extensive chromosomal rearrangements. ISME J. 11, 1189–1204. doi: 10.1038/ismej.2016.196
Hokken, M. W. J., Zoll, J., Coolen, J. P. M., Zwaan, B. J., Verweij, P. E., and Melchers, W. J. G. (2019). Phenotypic plasticity and the evolution of azole resistance in Aspergillus fumigatus; an expression profile of clinical isolates upon exposure to itraconazole. BMC Genomics 20:28. doi: 10.1186/s12864-018-5255-z
Hunt, R. C., Simhadri, V. L., Iandoli, M., Sauna, Z. E., and Kimchi-Sarfaty, C. (2014). Exposing synonymous mutations. Trends Genet. 30, 308–321. doi: 10.1016/j.tig.2014.04.006
Huson, D. H., and Bryant, D. (2006). Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267. doi: 10.1093/molbev/msj030
Irmer, H., Tarazona, S., Sasse, C., Olbermann, P., Loeffler, J., Krappmann, S., et al. (2015). RNAseq analysis of Aspergillus fumigatus in blood reveals a just wait and see resting stage behavior. BMC Genomics 16:640. doi: 10.1186/s12864-015-1853-1
John, H. (2008). Reference method for broth dilution antifungal susceptibility testing of filamentous fungi, approved standard. M38-A2. Clin. Lab Stand Inst. 28, 1–35.
Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405. doi: 10.1093/bioinformatics/btn129
Jombart, T., Devillard, S., and Balloux, F. (2010). Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11:94. doi: 10.1186/1471-2156-11-94
Kikuma, T., Tadokoro, T., Maruyama, J.-I., and Kitamoto, K. (2017). AoAtg26, a putative sterol glucosyltransferase, is required for autophagic degradation of peroxisomes, mitochondria, and nuclei in the filamentous fungus Aspergillus oryzae. Biosci. Biotechnol. Biochem. 81, 384–395. doi: 10.1080/09168451.2016.1240603
Kurucz, V., Kruger, T., Antal, K., Dietl, A. M., Haas, H., Pocsi, I., et al. (2018). Additional oxidative stress reroutes the global response of Aspergillus fumigatus to iron depletion. BMC Genomics 19:357. doi: 10.1186/s12864-018-4730-x
Latge, J. P. (1999). Aspergillus fumigatus and aspergillosis. Clin. Microbiol. Rev. 12, 310–350. doi: 10.1128/CMR.12.2.310
Latge, J. P., Beauvais, A., and Chamilos, G. (2017). The cell wall of the human fungal pathogen aspergillus fumigatus: biosynthesis, organization, immune response, and virulence. Annu. Rev. Microbiol. 71, 99–116. doi: 10.1146/annurev-micro-030117-020406
Latge, J. P., and Chamilos, G. (2019). Aspergillus fumigatus and Aspergillosis in 2019. Clin. Microbiol. Rev. 33:e00140-18. doi: 10.1128/CMR.00140-18
Leinonen, R., Sugawara, H., Shumway, M., and International Nucleotide Sequence Database, C. (2011). The sequence read archive. Nucleic Acids Res. 39, D19–21. doi: 10.1093/nar/gkq1019
Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/btp324
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Lin, S. J., Schranz, J., and Teutsch, S. M. (2001). Aspergillosis case-fatality rate: systematic review of the literature. Clin. Infect. Dis. 32, 358–366. doi: 10.1086/318483
Mckenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., et al. (2010). The genome analysis toolkit: a mapreduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. doi: 10.1101/gr.107524.110
Mellado, E., Garcia-Effron, G., Alcazar-Fuoli, L., Melchers, W. J., Verweij, P. E., Cuenca-Estrella, M., et al. (2007). A new Aspergillus fumigatus resistance mechanism conferring in vitro cross-resistance to azole antifungals involves a combination of cyp51A alterations. Antimicrob. Agents Chemother. 51, 1897–1904. doi: 10.1128/AAC.01092-06
Meneau, I., Coste, A. T., and Sanglard, D. (2016). Identification of Aspergillus fumigatus multidrug transporter genes and their potential involvement in antifungal resistance. Med. Mycol. 54, 616–627. doi: 10.1093/mmy/myw005
Moye-Rowley, W. (2015). Multiple mechanisms contribute to the development of clinically significant azole resistance in Aspergillus fumigatus. Front. Microbiol. 6:70. doi: 10.3389/fmicb.2015.00070
Muller, L. A., Lucas, J. E., Georgianna, D. R., and Mccusker, J. H. (2011). Genome-wide association analysis of clinical vs. nonclinical origin provides insights into Saccharomyces cerevisiae pathogenesis. Mol. Ecol. 20, 4085–4097. doi: 10.1111/j.1365-294X.2011.05225.x
Neofytos, D., Chatzis, O., Nasioudis, D., Boely Janke, E., Doco Lecompte, T., Garzoni, C., et al. (2018). Epidemiology, risk factors and outcomes of invasive aspergillosis in solid organ transplant recipients in the Swiss transplant cohort study. Transpl. Infect. Dis. 20:e12898. doi: 10.1111/tid.12898
Nierman, W. C., Pain, A., Anderson, M. J., Wortman, J. R., Kim, H. S., Arroyo, J., et al. (2005). Genomic sequence of the pathogenic and allergenic filamentous fungus Aspergillus fumigatus. Nature 438, 1151–1156. doi: 10.1038/nature04332
Oku, M., Warnecke, D., Noda, T., Muller, F., Heinz, E., Mukaiyama, H., et al. (2003). Peroxisome degradation requires catalytically active sterol glucosyltransferase with a GRAM domain. EMBO J. 22, 3231–3241. doi: 10.1093/emboj/cdg331
Palma-Guerrero, J., Hall, C. R., Kowbel, D., Welch, J., Taylor, J. W., Brem, R. B., et al. (2013). Genome wide association identifies novel loci involved in fungal communication. PLoS Genet. 9:e1003669. doi: 10.1371/journal.pgen.1003669
Paul, S., Stamnes, M., Thomas, G. H., Liu, H., Hagiwara, D., Gomi, K., et al. (2019). AtrR is an essential determinant of azole resistance in Aspergillus fumigatus. MBio 10:e02563-18. doi: 10.1128/mBio.02563-18
Peng, D., and Tarleton, R. (2015). EuPaGDT: a web tool tailored to design CRISPR guide RNAs for eukaryotic pathogens. Microb. Genom. 1:e000033. doi: 10.1099/mgen.0.000033
Power, R. A., Parkhill, J., and De Oliveira, T. (2017). Microbial genome-wide association studies: lessons from human GWAS. Nat. Rev. Genet. 18, 41–50. doi: 10.1038/nrg.2016.132
Price, A. L., Zaitlen, N. A., Reich, D., and Patterson, N. (2010). New approaches to population stratification in genome-wide association studies. Nat. Rev. Genet. 11, 459–463. doi: 10.1038/nrg2813
Read, T. D., and Massey, R. C. (2014). Characterizing the genetic basis of bacterial phenotypes using genome-wide association studies: a new direction for bacteriology. Genome Med. 6:109. doi: 10.1186/s13073-014-0109-z
Revie, N. M., Iyer, K. R., Robbins, N., and Cowen, L. E. (2018). Antifungal drug resistance: evolution, mechanisms and impact. Curr. Opin. Microbiol. 45, 70–76. doi: 10.1016/j.mib.2018.02.005
Robinett, K. S., Weiler, B., and Verceles, A. C. (2013). Invasive aspergillosis masquerading as catastrophic antiphospholipid syndrome. Am. J. Critic. Care 22, 448–451. doi: 10.4037/ajcc2013659
Sardi, M., Paithane, V., Place, M., Robinson, E., Hose, J., Wohlbach, D. J., et al. (2018). Genome-wide association across Saccharomyces cerevisiae strains reveals substantial variation in underlying gene requirements for toxin tolerance. PLoS Genet. 14:e1007217. doi: 10.1371/journal.pgen.1007217
Shih, S. C., Prag, G., Francis, S. A., Sutanto, M. A., Hurley, J. H., and Hicke, L. (2003). A ubiquitin-binding motif required for intramolecular monoubiquitylation, the CUE domain. EMBO J. 22, 1273–1281. doi: 10.1093/emboj/cdg140
Shimizu, K., and Keller, N. P. (2001). Genetic involvement of a cAMP-dependent protein kinase in a G protein signaling pathway regulating morphological and chemical transitions in Aspergillus nidulans. Genetics. 157, 591–600.
Stajich, J. E., Harris, T., Brunk, B. P., Brestelli, J., Fischer, S., Harb, O. S., et al. (2012). FungiDB: an integrated functional genomics database for fungi. Nucleic Acids Res. 40, D675–681. doi: 10.1093/nar/gkr918
Sul, J. H., Martin, L. S., and Eskin, E. (2018). Population structure in genetic studies: confounding factors and mixed models. PLoS Genet. 14:e1007309. doi: 10.1371/journal.pgen.1007309
Takahashi-Nakaguchi, A., Muraosa, Y., Hagiwara, D., Sakai, K., Toyotome, T., Watanabe, A., et al. (2015). Genome sequence comparison of Aspergillus fumigatus strains isolated from patients with pulmonary aspergilloma and chronic necrotizing pulmonary aspergillosis. Med. Mycol. 53, 353–360. doi: 10.1093/mmy/myv003
Talas, F., Kalih, R., Miedaner, T., and Mcdonald, B. A. (2016). 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. 29, 417–430. doi: 10.1094/MPMI-09-15-0218-R
Tashiro, M., Izumikawa, K., Minematsu, A., Hirano, K., Iwanaga, N., Ide, S., et al. (2012). Antifungal susceptibilities of Aspergillus fumigatus clinical isolates obtained in Nagasaki, Japan. Antimicrob. Agents Chemother. 56, 584–587. doi: 10.1128/AAC.05394-11
Thornton, T., and Mcpeek, M. S. (2010). ROADTRIPS: case-control association testing with partially or completely unknown population and pedigree structure. Am. J. Hum. Genet. 86, 172–184. doi: 10.1016/j.ajhg.2010.01.001
Turner, S. D. (2014). qqman: an R package for visualizing GWAS results using QQ and manhattan plots. biorxiv 005165. doi: 10.1101/005165
Van Der Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., Del Angel, G., Levy-Moonshine, A., et al. (2013). From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr. Protoc. Bioinformatics 43:bi1110s43. doi: 10.1002/0471250953.bi1110s43
Warnecke, D., Erdmann, R., Fahl, A., Hube, B., Müller, F., Zank, T., et al. (1999). Cloning and Functional expression of UGT genes encoding sterol glucosyltransferases from saccharomyces cerevisiae, Candida albicans, Pichia pastoris, and Dictyostelium discoideum. J. Biol. Chem. 274, 13048–13059. doi: 10.1074/jbc.274.19.13048
Warrilow, A. G., Parker, J. E., Price, C. L., Nes, W. D., Kelly, S. L., and Kelly, D. E. (2015). In vitro biochemical study of CYP51-mediated azole resistance in Aspergillus fumigatus. Antimicrob. Agents Chemother. 59, 7771–7778. doi: 10.1128/AAC.01806-15
Yu, J., Pressoir, G., Briggs, W. H., Bi, I. V., Yamasaki, M., Doebley, J. F., et al. (2006). A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat. Genet. 38, 203–208. doi: 10.1038/ng1702
Zhao, S., Ge, W., Watanabe, A., Fortwendel, J. R., and Gibbons, J. G. (2020). Genome-wide association for itraconazole sensitivity in non-resistant clinical isolates of Aspergillus fumigatus. bioRxiv. doi: 10.1101/2020.08.31.275297
Keywords: mycology, fungal pathogen, genome-wide association, population genomics, azole, antifungal, Aspergillus fumigatus, itraconazole
Citation: Zhao S, Ge W, Watanabe A, Fortwendel JR and Gibbons JG (2021) Genome-Wide Association for Itraconazole Sensitivity in Non-resistant Clinical Isolates of Aspergillus fumigatus. Front. Fungal Biol. 1:617338. doi: 10.3389/ffunb.2020.617338
Received: 14 October 2020; Accepted: 15 December 2020;
Published: 14 January 2021.
Edited by:
Gianni Liti, Centre National de la Recherche Scientifique (CNRS), FranceReviewed by:
Tibor Mihaly Nemeth, University of Szeged, HungaryHung-Ji Tsai, University of Birmingham, United Kingdom
Copyright © 2021 Zhao, Ge, Watanabe, Fortwendel and Gibbons. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: John G. Gibbons, jggibbons@umass.edu