A. fumigatus IP_23 (originally known as CBS 386.75) and IP_24 (originally known as CBS 110.46) were obtained from the Westerdijk Fungal Biodiversity Institute (KNAW) of the Netherlands (Sarfati et al., 2002). IFM47072 and IFM59985 are clinical isolates originally isolated from Japan (Zhao et al., 2021). Isolates were grown on potato dextrose agar (PDA) plates at 37°C for 72 hours. DNA was isolated directly from spores using the MasterPure yeast DNA purification kit following the manufacturer’s instructions, with several minor modifications, as previously described (Zhao et al., 2019), and as follows. Spores were collected with wash-store solution (49.9% glycerol, 50% potato dextrose broth and 0.1% Tween-20), centrifuged at 14,000 RPM for 5 minutes, and the supernatant was discarded. Next, 300 ml of yeast cell lysis solution was added to the spores along with 0.4 ml of 1.0 μm diameter silica beads. Cell lysis was carried out on a Biospec Mini-BeadBeater-8 at medium intensity for 8 minutes. One μl of RNase was added to the cell lysis solution and incubated at 65°C for 30 minutes. DNA purification was conducted according to the manufacturer’s instructions for the remainder of the protocol.

IFM47072 and IFM59985 were previously sequenced and whole-genome Illumina data is available through the NCBI Sequence Read Archive under run accession numbers SRR11977809 and SRR11977822, respectively (Zhao et al., 2021). 150-bp paired-end libraries were constructed and sequenced for IP_23 and IP_24 by Novogene on an Illumina NovaSeq 6000 sequencer. Trim_Galore v0.4.2 was used to remove residual adaptor sequences and trim reads at low-quality sites using the parameters “–paired”, “–stringency 1”, “–quality 30” and “–length 50” (Martin, 2011). Read sets were then error corrected and assembled using SPAdes v3.13.1 with the “–careful” parameter and K-mer sizes of 55, 77, and 99 (Bankevich et al., 2012). Genome assemblies for IP_23 and IP_24 are available through NCBI accession numbers JALLAH000000000 and JALLAI000000000, respectively. Assembly quality, using the genome assembly as input, was assessed with BUSCO v5 using the “eurotiomycetes_odb10” ortholog set through the gVolante v2.0.0 server (Simao et al., 2015; Nishimura et al., 2017). Gene prediction was performed with Augustus v3.3.2 (Stanke and Waack, 2003) using the following parameters “–strand=both”, “–genemodel=complete”, “–gff3=on”, “–codingseq=1”, “–protein=1”, “–uniqueGeneId=true”, and “–species=aspergillus_fumigatus”.

We previously sequenced the genomes of IFM47072 and IFM59985 and inferred their evolutionary relationship with 74 additional isolates (Zhao et al., 2021). This analysis revealed four major A. fumigatus populations. We selected a minimum of four individuals from each population and included IFM47072, IFM59985, IP_23 (386.75) and IP_24 (110.46) to infer the phylogenetic relationship between samples. Briefly, GATK v4.0.6.0 was used to call SNVs across all samples relative to the Af293 reference genome using the best practice pipeline for “Germline short variant discovery” (Van der Auwera et al., 2013). Next, hard filtering was carried out to reduce false positives using the “VariantFiltration” function with the following parameters: “QD < 25.0 || FS > 5.0 || MQ < 55.0 || MQRankSum < −0.5 || ReadPosRankSum < −2.0 || SOR > 2.5”. Hard filtering parameters were determined following the “Hard-filtering germline short variants” protocol (https://gatk.broadinstitute.org/hc/en-us/articles/360035890471-Hard-filtering-germline-short-variants). To reduce linkage between SNVs, we used VCFtools v0.1.14 to space SNV markers by a minimum of 3.5 Kb using the following option “–thin 3500”. Phylogenetic analysis of an alignment of 7,386 SNVs was conducted in MEGAX v10.0.5 (Kumar et al., 2018). We used model test in MEGAX to predict the best fit nucleotide substitution model, which was determined as the model with the lowest Bayesian Information Criterion (BIC) value. The general time reversible (GTR) gamma substitution model, with 100 bootstrap replicates, was used to construct the maximum likelihood phylogenetic tree.

Quality and adapter trimmed read sets from A. fumigatus IP_23, IP_24, IFM47072 and IFM59985 were mapped against the A. fumigatus Af293 reference genome (Genbank Assembly Accession number: GCA_000002655.1) (Nierman et al., 2005) using BWA-MEM v0.7.15 (Li and Durbin, 2009). Joint variant calling (SNVs and indels) was conducted using freebayes v1.3.1 with the default settings with the exception of setting ploidy to haploid (–ploidy = 1) (Garrison and Marth, 2012). While the vast majority of indels consisted of 1 or 2 bps (68%), indels as large as 49 bp were identified.

We implemented several filtering steps to prioritize the identification of variants putatively associated with the white-spore phenotype (depicted in Figure 2). We identified 51,551 SNVs and 15,183 indels that differed in at least one sample relative to the Af293 reference genome. Next, we narrowed our candidate list by identifying variants that were (i) identical in IP_23 and IP_24, (ii) different from the Af293 genotype, and (iii) different from the IFM47072 and IFM59985 genotypes (4,279 SNVs and 1,785 indels) (Data Sheet S1). Finally, we further filtered our candidate variant list by prioritizing variants annotated as “high impact” using SnpEff v4.3t with the A. fumigatus Af293 reference genome annotation (Cingolani et al., 2012) (19 SNVs and 45 indels). Saccharomyces cerevisae and Schizosaccharomyces pombe orthologs of the A. fumigatus candidate genes were identified using FungiDB (Stajich et al., 2012; Basenko et al., 2018).

Because DHN melanin is the molecule responsible for the greenish pigment in wild-type A. fumigatus spores, we focused on variants in the DHN melanin encoding gene cluster (Heinekamp et al., 2012). This gene cluster contains 6 genes (Afu2g17530, Afu2g17540, Afu2g17550, Afu2g17560, Afu2g17580 and Afu2g17600). We used the “intersect” function in bedtools v2.29.2 (Quinlan and Hall, 2010) to identify SNVs and indels for which genotypes were identical in the white-spore isolates, different from Af293, and different from the green-spore isolates. We focused on pksP (Afu2g17600) because it contained 3 candidate SNVs (all synonymous variants in the pksP gene on chromosome 2: 4687866, 4688908 and 4691782) and 2 candidate indels (one intron variant chromosome 2: 4692683 and one frameshift variant chromosome 2: 4692995) (Data Sheet 1).

The pksP gene was extracted from the IP_23, IP_24, IFM47072, IFM59985 and Af293 genomes using the samtools “faidx” function (Li et al., 2009). pksP sequences were aligned with MAFFT using the default settings (Katoh and Standley, 2013). Alignments were visualized in BioEdit (Hall et al., 2011), variants were visually inspected, and the coding region was extracted and translated. PksP protein domains were predicted in white-spore and green-spore isolates using the PFAM webserver (http://pfam.xfam.org/search/sequence) (El-Gebali et al., 2019).

The Alphafold2 model (Jumper et al., 2021) for A. fumigatus PksP was downloaded from EMBL-EBI (https://alphafold.ebi.ac.uk/entry/Q4WZA8) along with the corresponding predicted aligned error (PAE) file. The model was pruned of low confidence residues (pLDDT< 50). The PAE of the deleted domains were then evaluated for inter-residue confidence. Within the truncated region, each of the three terminal domains (two phosphopantetheine attachment site (PP-binding) and one thioesterase) are well supported in their internal residue placement (pLDDT > 70) and structural arrangement (PAE < 5 Å).

Crystal structure of the thioesterase domain of Aspergillus parasiticus PksA was downloaded from PDB (3ILS) (Korman et al., 2010). Structural superposition between the crystal structure and the AlphaFold2 predicted thioesterase domain of A. fumigatus PksP was performed using Chimera Matchmaker (Meng et al., 2006). Hydrogen bond identification was done in UCSC Chimera (version 1.15, FindHBond), using the rotamer orientation predicted by AlphaFold2. Hydrogen bonds were predicted only between residues having both high absolute confidence (pLDDT > 70) and precise relative placement (PAE < 5 Å).

We compared growth rate on minimal media (MM) at 37°C and 48°C, to examine whether our natural albino variants displayed thermotolerance phenotypes, as previously described for albino strains isolated from the Brazilian rainforest (Couger et al., 2018). Aspergillus MM was prepared as previously described (Cove, 1966). For each isolate, ~10⁵ spores were inoculated onto the center of MM agar plates and incubated at 37°C for 2 days and 48°C for 3 days in darkness. Experiments were conducted in triplicate and at the conclusion of the experiment colony diameter was measured with digital calipers. Because colony shape is often slightly irregular, we took two random measurements of diameter per plate and averaged these values. Additionally, to assess whether the putative non-functionality of IP_23 and IP_24 PksP contributes to growth patterns at optimal or thermal stress conditions, we performed the same experiment described above, with A. fumigatus Af293 and the Af293 ΔpksP null mutant (which produces white conidia) generated as previously described (Al Abdallah et al., 2017).

Genome assemblies for IP_23 (386.75), IP_24 (110.46), IFM47072 and IFM59985 were generated with SPAdes (Bankevich et al., 2012). Genome size, as indicated by the cumulative length of the genome assemblies, was highly similar across the four isolates and ranged from 28.36 Mb – 28.80 Mb (Table S1). We used Augustus to predict gene models in each assembly, and similarly, found highly comparable gene numbers (8,855 – 8,908) (Table S1). We used BUSCO to assess genome completeness for each assembly and found that >97% of genes were completely recovered using the eurotiomycetes_odb10 dataset (Simao et al., 2015). These results reflect high quality genome assemblies.

We previously performed population genomic analysis of 76 clinical A. fumigatus isolates from Japan and identified 4 major populations (Zhao et al., 2021). To examine the relationship of IP_23 and IP_24 to these isolates, we conducted phylogenetic analysis with at least 4 individuals from each population. Our analysis again confirmed the presence of 4 populations (Figure 1). IP_23 and IP_24 clustered within population 1 and were closely related to IFM47072 and IFM59985 (Figure 1). Thus, we focused our comparative genomic analysis on the white-spore isolates IP_23 and IP_24, and the closely related green-spore isolates IFM47072 and IFM59985 (Figure 1). We identified 10,950 SNV sites in which IFM47072 and IFM59985 differed in genotype and 2,129 SNV sites in which IP_23 and IP_24 differed in genotype, suggesting IFM47072 and IFM59985 and IP_23 and IP_24 are closely related but not clones.

We used FreeBayes to perform joint genotyping of IP_23, IP_24, IFM47072 and IFM59985 relative to the Af293 genome. We identified 51,551 SNVs and 15,183 indels that differed in at least one sample relative to the Af293 reference genome. Next, we narrowed our candidate list by identifying variants that were (i) identical in IP_23 and IP_24, (ii) different from the Af293 genotype, and (iii) different from the IFM47072 and IFM59985 genotypes. This approach yielded 4,279 SNVs and 1,785 indels (Data Sheet S1). Next, we predicted the putative functional effects of these candidate variants using SnpEff (Cingolani et al., 2012). We hypothesized that the white-spore phenotype was the result of a loss-of-function mutation and further prioritized candidate variants that were annotated as “High Impact” by SnpEff (i.e., loss of stop codon, gain of stop codon, loss of start codon, splice donor variant, and splice acceptor variant). We identified 19 SNVs and 45 indels annotated as putative high impact mutations (Figure 2, Data Sheet S2). One of the 64 candidate high impact variants was present in the biosynthetic gene in the DHN-melanin gene cluster (pksP, Afu2g17600), and was annotated as a frameshift variant (Data Sheet S2). This variant, which is only present in IP_23 and IP_24, is a single nucleotide deletion in the fifth exon of pksP and results in a premature stop codon yielding a predicted protein length of 1,686 amino acids in the white-spore isolates compared to 2,148 amino acids in the green-spore isolates (Figure 3). To confirm this variant, which was identified through a mapping-based approach, we extracted the pksP gene from each genome assembly and aligned the sequences. Again, we observed the single nucleotide deletion at position 5,035.

To gain insight into the functional consequences of the PksP protein truncation, we predicted PFAM protein domains in the green-spore and white-spore PksP proteins. The green-spore PksP contains 8 protein domains (Figure 3D). However, the last two-and-a-half PksP protein domains are absent in the white-spore protein. These domains include two consecutive PP-binding domains (PF00550) and a thioesterase domain (PF00975) (Figure 3D).

Structural modeling of the PksP protein using AlphaFold2 (Cramer, 2021) confirms that the truncation removes several highly structured, C-terminal motifs (Figure 4A). The largest of the structural motifs predicted by AlphaFold2 superimposes over the crystal structure of the thioesterase domain of Aspergillus parasiticus PksA (Korman et al., 2010) to sub-angstrom RMSD (RMSD=0.987Å; Figure 4B). The truncation also removes six residues (Lysine 1711, Isoleucine 1710, Threonine 1709, Proline 1708, Tyrosine 1691, Arginine 1690) predicted to participate in hydrogen bonding with residues within the remaining portion of the protein (Figure 4C). Thus, the AlphaFold2 model not only corroborates the thioesterase predictions, it also suggests that the truncation may affect conformational stability within the remaining protein through the removal of hydrogen bond partners.

Couger et al. (2018) isolated an albino strain of A. fumigatus from Brazilian rainforest soil that was thermotolerant (Couger et al., 2018). To assess whether IP_23 and IP_24 were also thermotolerant, we compared growth rate on minimal media (MM) at 37°C and 48°C. A one-way ANOVA was performed to compare the effect of temperature on growth rate across green-spore and white-spore isolates, which revealed a significant difference at both temperatures (37°C: F-ratio = 86.8, d.f. = 3, p-value = 6.7e-6 and 48°C: F-ratio = 899.9, d.f. = 3, p-value = 1.9e-10). A post-hoc Tukey’s HSD test revealed significant growth rate differences between all green-spore and white-spore comparisons (all p-values ≤ 0.0001) but not between green-spore isolates or white-spore isolates (all p-values ≥ 0.09). The white-spore isolates grew significantly slower at 37°C and 48°C (Figures 5A, B). White-spore isolates did not grow at 48°C.

To assess whether PksP functionality directly contributes to the thermosensitivity phenotype in IP_23 and IP_24, we assessed growth rate patterns of A. fumigatus Af293 and Af293 ΔpksP (Al Abdallah et al., 2017), at 37°C and 48°C (Figures 5C, D). We did not observe significant growth rate differences between Af293 and ΔpksP for either temperature (Figures 5C, D), suggesting that pksP is not directly involved in the thermosensitivity phenotype observed in IP_23 and IP_24.

To provide additional insight into variants and/or genes that may contribute to the thermosensitivity phenotype, we identified S. cerevisiae orthologs for the 57 candidate genes containing the 64 high impact mutations described earlier (Figure 2). We observed a 1 bp insertion resulting in a frameshift in Afu4g04740 in IP_23 and IP_24 (Chr4: 1,336,071) (Datasheet 2). The S. cerevisiae ortholog of Afu4g04740 is SRM1, and mutants of SRM1 (also known as prp20) are temperature sensitive (Vijayraghavan et al., 1989; Fleischmann et al., 1991; Clement et al., 2006). In addition, null mutants of pim1, the Schizosaccharomyces pombe ortholog of Afu4g04740, are also thermosensitive (Matsuo et al., 2011).