Genomic and Phenotypic Heterogeneity of Clinical Isolates of the Human Pathogens Aspergillus fumigatus, Aspergillus lentulus, and Aspergillus fumigatiaffinis

Fungal pathogens are a global threat to human health. For example, fungi from the genus Aspergillus cause a spectrum of diseases collectively known as aspergillosis. Most of the >200,000 life-threatening aspergillosis infections per year worldwide are caused by Aspergillus fumigatus. Recently, molecular typing techniques have revealed that aspergillosis can also be caused by organisms that are phenotypically similar to A. fumigatus but genetically distinct, such as Aspergillus lentulus and Aspergillus fumigatiaffinis. Importantly, some of these so-called cryptic species are thought to exhibit different virulence and drug susceptibility profiles than A. fumigatus, however, our understanding of their biology and pathogenic potential has been stymied by the lack of genome sequences and phenotypic profiling of multiple clinical strains. To fill this gap, we phenotypically characterized the virulence and drug susceptibility of 15 clinical strains of A. fumigatus, A. lentulus, and A. fumigatiaffinis from Spain and sequenced their genomes. We found heterogeneity in drug susceptibility across species and strains. We further found heterogeneity in virulence within each species but no significant differences in the virulence profiles between the three species. Genes known to influence drug susceptibility (cyp51A and fks1) vary in paralog number and sequence among these species and strains and correlate with differences in drug susceptibility. Similarly, genes known to be important for virulence in A. fumigatus showed variability in number of paralogs across strains and across species. Characterization of the genomic similarities and differences of clinical strains of A. lentulus, A. fumigatiaffinis, and A. fumigatus that vary in disease-relevant traits will advance our understanding of the variance in pathogenicity between Aspergillus species and strains that are collectively responsible for the vast majority of aspergillosis infections in humans.


INTRODUCTION
Aspergillosis is a major health problem, with rapidly evolving epidemiology and new groups of at-risk patients (Patterson et al., 2016). Aspergillosis infections are usually caused by inhalation of airborne asexual spores (conidia) of Aspergillus fumigatus and a few other Aspergillus species (Rokas et al., 2020). Aspergillosis covers a spectrum of diseases (Latgé and Chamilos, 2020). For example, non-invasive diseases caused by Aspergillus, such as aspergilloma, are currently classified as chronic pulmonary aspergillosis and are commonly associated to pulmonary tuberculosis . In atopic patients, the most severe form of aspergillosis is allergic bronchopulmonary aspergillosis (ABPA), which develops following sensibilization to A. fumigatus allergens in atopic patients with cystic fibrosis or individuals with genetic predisposition to ABPA (Agarwal et al., 2013). However, the most common invasive type of infection is invasive pulmonary aspergillosis (IPA), whose risk is significantly increased in immunocompromised individuals, in patients with acute leukemia and recipients of hematopoietic stem cells transplantation, or in solid-organ transplant recipients (Brown et al., 2012). Importantly, IPA has recently been described in new groups of traditionally low-risk patients, such as patients in intensive care units recovering from bacterial sepsis (Latgé and Chamilos, 2020).
Although A. fumigatus is the major etiologic agent of aspergillosis, a few other Aspergillus species, such as Aspergillus flavus, Aspergillus terreus, Aspergillus niger, and Aspergillus nidulans, can also cause infections (Zakaria et al., 2020). While most of these pathogens can be phenotypically easily distinguished, infections can also be caused by Aspergillus species that are morphologically very similar to A. fumigatus (Rokas et al., 2020). These close pathogenic relatives of A. fumigatus are considered sibling species or cryptic species because they are undistinguishable from each other and from A. fumigatus by classical identification methods (Alastruey-Izquierdo et al., 2014); these species vary mostly in their colony growth, robustness of the production of conidia, conidial surface markings, presence and absence of septation in phialides, and maximum growth temperatures (Taylor et al., 2000;Balajee et al., 2005;Katz et al., 2005). As a result of their near identical morphological characteristics, most of these cryptic species have only recently been described. For example, Aspergillus lentulus was first described in 2005 in a case of human aspergillosis (Balajee et al., 2005). Similarly, A. fumigatiaffinis, another pathogenic species that is closely related to A. fumigatus, was first described in 2005 (Hong et al., 2005). Even though cryptic species were only discovered relatively recently, understanding their genetic and phenotypic similarities and differences from the major pathogen A. fumigatus is important for two reasons. First, their prevalence in the clinic has been estimated to be between 11 and 19% (Balajee et al., 2009;Alastruey-Izquierdo et al., 2014;Negri et al., 2014). Second, several of these species, including A. lentulus and A. fumigatiaffinis, have been shown to differ in their drug susceptibility to amphotericin B and azoles compared to A. fumigatus (Alastruey-Izquierdo et al., 2014).
Antifungal resistance is of worldwide concern in human pathogenic Aspergillus species as well as in many other human, animal, and plant fungal pathogens (Parker et al., 2014;Sharma and Chowdhary, 2017). Several antifungalresistance mechanisms have been proposed in fungi (Sharma and Chowdhary, 2017;Perez-Cantero et al., 2020). In azoleresistant Aspergillus strains, known mechanisms are particularly well-described in genes of the cytochrome P450 sterol 14 α-demethylase family (cyp51), and include sequence variants in diverse positions of the Cyp51A protein sequence (e.g., G54, G138, M220, G448, Y121, P216, F219, A284, Y431, G432, and G434; reviewed in Wei et al., 2015;Perez-Cantero et al., 2020), as well as combinations of the aforementioned protein sequence changes with tandem repeat (TR) variants in the promoter region, such as TR34/L98H or TR46/Y121F/T289A (reviewed in Wei et al., 2015). Non-cyp51 based mechanisms of antifungal resistance, such as in multidrug efflux pumps and pathways such as ergosterol biosynthesis and stress response, have also been proposed (Perez-Cantero et al., 2020). Mechanisms of echinocandin resistance have mostly been attributed to FKS subunits of glucan synthase (Sharma and Chowdhary, 2017). While most of these studies are in Candida species (Desnos-Ollivier et al., 2008;Garcia-Effron et al., 2008), a recent study in A. fumigatus also observed mutations associated with echinocandin resistance (Jiménez-Ortigosa et al., 2017).
An emerging realization in the study of Aspergillus pathogens is the presence of phenotypic heterogeneity among strains of the same species (Keller, 2017). For example, recent studies have shown how variation in hypoxic growth phenotypes is associated with virulence among A. fumigatus strains (Kowalski et al., 2016. Similarly, A. fumigatus strains have previously been shown to exhibit great quantitative and qualitative heterogeneity in light response ; in this case, heterogeneity in light response was not associated with heterogeneity in virulence. Finally, Ries et al. (2019) found a high heterogeneity among A. fumigatus strains with regard to nitrogen acquisition and metabolism during infection and correlation between nitrogen catabolite repressionrelated protease secretion and virulence. These studies highlight the biological and clinical relevance of understanding strain heterogeneity in Aspergillus pathogens, especially with respect to virulence and antifungal drug susceptibility. However, comparisons of strain heterogeneity in virulence and drug resistance profiles among clinical strains in A. fumigatus and closely related cryptic species, such as A. lentulus and A. fumigatiaffinis, are lacking.
To address this gap in the field, we phenotypically characterized and sequenced the genomes of 15 clinical strains of A. fumigatus, A. lentulus, and A. fumigatiaffinis from Spain. At the phenotypic level, we found strain heterogeneity in both virulence and drug susceptibility profiles within each species as well as differences in drug susceptibility profiles between the three species. Interestingly, we found that the virulence profiles of the three species were similar. At the genomic level, we found that gene families known to influence drug susceptibility, such as cyp51, exhibit variation in their numbers of paralogs and sequence among these species and strains. Similarly, we found variability in the number of paralogs within and between species in many genes known to be important for virulence in A. fumigatus. Characterization of the genomic similarities and differences of clinical strains of A. lentulus, A. fumigatiaffinis, and A. fumigatus that vary in disease-relevant traits will advance our understanding of the variation in pathogenicity between Aspergillus species and strains that are collectively responsible for the vast majority of aspergillosis infections in humans.

Strains and Species Identification
To understand the degree of genomic heterogeneity among strains, we sequenced six clinical strains of A. fumigatus, five of A. lentulus, and four of A. fumigatiaffinis available in the Mycology Reference Laboratory of the National Center for Microbiology (CNM) in Instituto de Salud Carlos III in Spain (Supplementary Table S1). For initial species identification, we sequenced the Internal Transcribed Spacer region (ITS) and beta-tubulin (benA) gene amplicons (primer pairs in Supplementary Table S2). We downloaded reference sequences for the type strains of A. fumigatiaffinis IBT12703 and A. lentulus IFM54703, and of Aspergillus clavatus NRRL1 (section Clavati), which we used as the outgroup. We aligned DNA sequences with MAFFT v.7.397 (Katoh and Standley, 2013), followed by model selection and phylogenetic inference in IQ-TREE v.1.6.7 (Nguyen et al., 2015).

Characterization of Virulence and Antifungal Susceptibility Profiles
To understand the pathogenic potential of the 15 clinical strains, we carried out virulence assays using the moth Galleria mellonella model of fungal disease (Fuchs et al., 2010;Slater et al., 2011). Briefly, we obtained moth larvae by breeding adult moths that were kept for 24 h prior to infection under starvation, in the dark, and at a temperature of 37 • C. We selected only larvae that were in the sixth and final stage of larval development. We harvested fresh asexual spores (conidia) from each strain from yeast extract-agar-glucose (YAG) plates in PBS solution and filtered through a Miracloth (Calbiochem). For each strain, we counted the spores using a hemocytometer and created a 2 × 10 8 conidia/ml stock suspension. We determined the viability of the administered inoculum by plating a serial dilution of the conidia on YAG medium at 37 • C. We inoculated 5 µl (1 × 10 6 conidia/larvae) to each larva (n = 10). We used as the control a group composed of larvae inoculated with 5 µl of PBS. We performed inoculations via the last left proleg using a Hamilton syringe (7000.5KH). After infection, we maintained the larvae in petri dishes at 37 • C in the dark and scored them daily (i.e., recorded the number of dead larvae each day) during a 10-day period. We considered larvae that did not move in response to touch as dead.
We tested the virulence of each clinical strain by infecting 10 larvae, i.e., for each strain tested we have one experimental replicate with an sample size n of 10. We performed two sets of analyses. First, we statistically assessed if the survival curves of different strains in a given species are identical (null hypothesis of strain homogeneity) or different (alternative hypothesis of strain heterogeneity). Second, we used strains within each species as "biological replicates" and statistically assessed if the survival curves between species were similar or different. We performed these statistical assessments using the log-rank test implemented in the survival R package (Therneau, 2014), followed by multiple test correction of p-values (Benjamini and Hochberg). Scripts used to perform these analyses are available on the GitLab repository 1 under 'experimentalData/'.
To measure the antifungal susceptibility of the clinical strains, we applied the EUCAST (European Committee for Antimicrobial Susceptibility Testing) reference microdilution method version 9.3.1 (Arendrup et al., 2017), in which fungi are grown on plates with increasing concentrations of antifungals and the first concentration in which fungal growth is inhibited (MIC) is recorded. For all strains, we tested their susceptibility to four antifungal drug classes: (a) Polyenes: amphotericin B (Sigma-Aldrich Quimica, Madrid, Spain); (b) Azoles: itraconazole (Janssen Pharmaceutica, Madrid, Spain), voriconazole (Pfizer SA, Madrid, Spain), and posaconazole (Schering-Plough Research Institute, Kenilworth, NJ, United States); (c) Echinocandins: caspofungin (Merck & Co. Inc., Rahway, NJ, United States), micafungin (Astellas Pharma Inc., Tokyo, Japan), and anidulafungin (Pfizer SA, Madrid, Spain); and (d) Allylamines: Terbinafine (Novartis, Basel, Switzerland). The final concentrations tested ranged from 0.03 to 16 mg/L for amphotericin B, terbinafine, and caspofungin; from 0.015 to 8 mg/L for itraconazole, voriconazole and posaconazole; from 0.007 to 4 mg/L for anidulafungin; and from 0.004 to 2 mg/L for micafungin. A. flavus ATCC 204304 and A. fumigatus ATCC 204305 were used as quality control strains in all tests performed. MICs for amphotericin B, itraconazole, voriconazole, posaconazole, and terbinafine, and minimal effective concentrations (MECs) for anidulafungin, caspofungin, and micafungin were visually read after 24 and 48 h of incubation at 35 • C in a humid atmosphere. To assess the relationship between antifungal susceptibility and strain/species identification, we carried out principal component analysis (PCA) with scaled MIC/MEC values with the R package FactoMineR (Lê et al., 2008), and data visualization with the factoextra v.1.0.6 package. Scripts used to perform these analyses are available on the GitLab repository (see text footnote 1) under 'experimentalData/'.

Genome Sequencing
To understand the genomic similarities and differences within and between these pathogenic Aspergillus species and how they are associated with differences in drug susceptibility and virulence profiles, we sequenced the genomes of all 15 strains. Each strain was grown in glucose-yeast extract-peptone (GYEP) liquid medium (0.3% yeast extract and 1% peptone; Difco, Soria Melguizo) with 2% glucose (Sigma-Aldrich, Spain) for 24 h to 48 h at 30 • C. After mechanical disruption of the mycelium by vortex mixing with glass beads, genomic DNA of isolates was extracted using the phenol-chloroform method (Holden, 1994). The preparation of DNA libraries was performed using the Nextera R TM DNA Library PrepKit (Illumina Inc., San Diego, CA, United States) according to manufacturer's guidelines. DNA quantification was carried out using the QuantiFluor R dsDNA System and the QuantiFluor R ST Fluorometer (Promega, Madison, WI, United States) and its quality was checked with the Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, United States). Sequencing was performed in the Illumina platform NextSeq500, following the manufacturer's protocols (Illumina Inc., San Diego, CA, United States). We performed an initial quality analysis of the sequence reads using FastQC, v.0.11.7 2 . We inspected sequence reads for contaminants using BLAST (Altschul et al., 1990) and MEGAN5 (Huson and Weber, 2013). We trimmed low quality bases (LEADING = 3; TRAILING = 3; SLIDINGWINDOW: windowSize = 4 and requiredQuality = 15), removing both short sequences (<90 bp) and Nextera adaptors, with Trimmomatic v.0.38 (Bolger et al., 2014).

Orthogroup Identification
To identify orthologs (and closely related paralogs) across strains, we performed all-vs.-all searches with blastp 2.7.1+ (Altschul et al., 1990) using the strains' predicted proteomes. We used OrthoFinder v.2.3.3 (Emms and Kelly, 2019) to generate orthogroups using pre-computed BLAST results (-og option) and a Markov Clustering (MCL) inflation value of 1.5. We considered an orthogroup "species-specific" if it possessed one or more protein sequences from only one species. Information on performing these analyses is available on the GitLab wiki page 'orthology-calling 3 '.

Identification of Single Nucleotide Polymorphisms and Insertions/Deletions
To characterize genetic variation within and between the three pathogenic Aspergillus species, we assessed single nucleotide polymorphisms (SNPs) and insertions/deletions (indels). We used BWA-MEM v.0.7.17 (Li and Durbin, 2009) with default parameters to map reads to the reference genome sequences for A. fumigatus, A. lentulus, and A. fumigatiaffinis (CNM-CM8686, CNM-CM7927, and CNM-CM6805, respectively). We did not use type strains as reference genomes for the species under study, because they are not from Spain. Duplicate reads were identified using PICARD MarkDuplicates, v.2.9.2 4 . We indexed genomes using SAMTOOLS v.1.8 ) for subsequent variant detection analyses.
We used GenomeAnalysisTK (GATK) v.3.6 for SNP calling with the recommended hard filtering parameters (McKenna et al., 2010;Depristo et al., 2011). We used SnpEff v.4.3t (Cingolani et al., 2013) to annotate and predict the functional effect of SNPs and indels. Variants assumed to have high (disruptive) impact in the protein, probably causing protein truncation, loss of function or triggering nonsense mediated decay were classified as "high, " variants assumed that might change protein effectiveness but were non-disruptive were classified as "moderate, " and variants most likely to be harmless or unlikely to change protein behavior were classified as "low." Finally, non-coding variants or variants affecting non-coding genes, where predictions are difficult or there is no evidence of impact, were classified as "modifier." Details can be found on the SnpEff manual 5 .
We aligned protein and coding sequences for genes of interest with MAFFT v.7.397 (Katoh and Standley, 2013), using the -auto mode. We used Jalview v.2.10.3 (Waterhouse et al., 2009) to visualize SNPs, and a Python script to recover nonsynonymous mutations compared to the reference, A. fumigatus A1163. Enrichment analysis of GO terms in genes with high impact SNPs and indels for each species was carried out with GOATOOLS v.0.9.9 (Klopfenstein et al., 2018). Scripts used to perform these analyses are available on the GitLab repository (see text footnote 1) under 'genomePolymorphisms/' and 'goatools/'.

Genetic Determinants Important for Virulence
To examine whether SNPs, indels, and number of paralogs in a given orthogroup were associated with virulence, we recovered 215 genes in A. fumigatus Af293 considered genetic determinants of virulence based on their presence in PHIbase (Winnenburg, 2006) and in previously published studies (Abad et al., 2010;Kjaerbølling et al., 2018). We obtained functional annotation of these virulence-related genes from FungiDB (Basenko et al., 2018).
To identify single-copy orthologous genes among the 24 genomes, we implemented the BUSCO, v.2.0.1 pipeline (Waterhouse et al., 2013;Simão et al., 2015). Specifically, we used the BUSCO pipeline to identify single-copy orthologous genes from genomes using the Eurotiomycetes database of 4,046 orthologs from OrthoDB, v9 (Waterhouse et al., 2013). Among the 4,096 orthologs, we identified 3,954 orthologs with at least 18 taxa represented and aligned the protein sequence each ortholog individually using Mafft, v7.294b (Katoh and Standley, 2013), with the same parameters as described elsewhere . We then forced nucleotide sequences onto the protein alignment with a custom Python, v3.5.2 script (indicated on the Gitlab repository README.md file) using BioPython, v1.7 (Cock et al., 2009). The resulting nucleotide alignments were trimmed using trimAl, v1.4 (Capella-Gutierrez et al., 2009), with the 'gappyout' parameter. The trimmed alignments were then concatenated into a single matrix with 7,147,728 sites. We then used the concatenated data matrix as input into IQ-TREE, v1.6.11 (Nguyen et al., 2015), with the 'nbest' parameter set to 10. The best-fitting model of substitutions was automatically determined using the Bayesian information criterion. The best-fitting model was a general time general time-reversible model with empirical base frequencies, a discrete Gamma model with 4 rate categories, and a proportion of invariable sites (GTR+I+F+G4) (Tavaré, 1986;Yang, 1994Yang, , 1996Vinet and Zhedanov, 2011). Lastly, we evaluated bipartition support using 5,000 ultrafast bootstrap approximations (Hoang et al., 2018).
In order to build the phylogeny with Cyp51 paralogs, we recovered protein sequences from two orthogroups that included Cyp51A and Cyp51B from A. fumigatus Af293 (Afu4g06890 and Afu7g03740, respectively). We generated a maximum-likelihood phylogeny in IQ-Tree v. 1.6.12 (Nguyen et al., 2015), using 1000 Ultrafast Bootstrap Approximation (UFBoot) replicates. The LG+G4 model was chosen as the best according to Bayesian Information Criterion. The protein sequences and tree files are available on the GitLab repository (see text footnote 1) under ' AntifungalGenes/'.

Clinical Strains Show Varying Antifungal Drug Susceptibility
To study susceptibility to antifungals across all strains of the three Aspergillus pathogens, we employed the EUCAST reference microdilution method with the four different known classes of antifungal drugs ( Table 1). By performing PCA on the antifungal drug susceptibility values of all 15 strains, we found that the strains exhibited high heterogeneity in their drug resistance profiles ( Figure 1A). In many cases, we found that strains from different species were more similar to each other (e.g., strain CNM-CM8686 from A. fumigatus with strain CNM-CM6069 from A. lentulus) than to other strains from the same species (e.g., strain CNM-CM8686 with strain CNM-CM8057 from A. fumigatus), highlighting the magnitude of heterogeneity in drug susceptibility of these species and strains. Principal component 1 (PC1) explained 37.2% of the variation and separated almost all A. fumigatus strains from those of the other two species. Principal component 2 (PC2) explained 21% of the variation, but did not separate species. The individual contributions of each antifungal drug to each PC are shown in Supplementary Figure S1. Finally, we found that the susceptibility of amphotericin B (in the polyenes class) was negatively correlated with micafungin (echinocandins) and terbinafine (allylamines), whereas anidulafungin (echinocandins) and voriconazole (azoles) were positively correlated (Supplementary Figure S2). Interestingly, the drugs exhibiting these negative or positive correlations are from different classes (e.g., polyenes versus allylamines or echinocandins versus azoles).
We also looked at the differences in susceptibility between strains for each antifungal drug ( Figure 1B). Our data show that clinical strains of A. fumigatus exhibit lower MICs to amphotericin B compared to A. lentulus and A. fumigatiaffinis, albeit different levels are observed among different strains (one-way ANOVA; α < 0.05; Tukey multiple comparisons of means for amphotericin B) ( Table 1). With the exception of susceptibility of A. fumigatus and A. lentulus to amphotericin B, for which a significant difference is observed between these two species, we observed high heterogeneity among strains of different species for the other drugs (Table 1). Among azoles, itraconazole and voriconazole displayed higher levels of variability across strains. With respect to terbinafine, the four A. fumigatiaffinis strains exhibited low MICs, whereas four A. fumigatus strains displayed higher MICs (MIC values >1 mg/L) and the other two A. fumigatus strains even higher; finally, one A. lentulus strain (CNM-CM8694) displayed the highest MICs across all strains (albeit other strains showed in general lower MICs). Among echinocandin drugs, caspofungin showed high MECs for the three species. In particular, one strain of A. fumigatiaffinis and three of A. lentulus were notable in exhibiting very high MECs (MECs ≥ 1 mg/L). MECs for micafungin and anidulafungin were low (≤0.125 mg/L) for all strains.

Clinical Strains Within Each Species Show Varying Levels of Virulence
Given functional similarities of the greater wax moth Galleria mellonella innate immune system with that of mammals, and prior work showing that moth larvae and mice exhibit similar survival rates when infected with A. fumigatus (Slater et al., 2011;Mead et al., 2019), we infected G. mellonella larvae with all 15 strains to assess their virulence profiles (Figure 2). Survival curves revealed high heterogeneity in virulence across clinical strains within each of the three species (Figure 2). We observed highly virulent strains for which all ten larvae were dead at day 10, such as A. fumigatus Af293 (one of our reference strains), A. fumigatiaffinis CNM-CM5878 and A. lentulus CNM-CM8927. In contrast, other strains were less virulent and >25% larvae survived to the last day of data collection, such as A. lentulus CNM-CM6069 and CNM-CM8060. Moreover, we found significant heterogeneity in the survival curves between strains within each species (Benjamini and Hochberg adjusted p-values: 0.00285 in A. fumigatus, 0.00054 in A. fumigatiaffinis, and 0.014 in A. lentulus; log-rank test) (Figures 2A-C). We also tested differences between species (considering each strain as a biological replicate), and observed no significant difference between the kill curves of the various species (p-value = 0.17; log-rank test) -that is, we found that both A. lentulus and A. fumigatiaffinis were as virulent as A. fumigatus (Figure 2D).

Genomic Variation Within and Between Spanish Strains of A. fumigatus, A. lentulus, and A. fumigatiaffinis
To begin exploring the potential genetic underpinnings of species and strain variation in drug susceptibility and virulence, we conducted comparative genomic analyses. The genomes of all 15 strains were of high quality and contained 97-98% of expected complete and single-copy BUSCOs (Supplementary Table S3). A. lentulus and A. fumigatiaffinis genomes had larger gene repertoires (9,717-9,842 and 10,329-10,677, respectively) than A. fumigatus (8,837-8,938), consistent with previous genome studies of A. lentulus and A. fumigatus (Nierman et al., 2005;Fedorova et al., 2008;Kusuya et al., 2016). A genome-scale phylogenetic analysis using the nucleotide sequences of BUSCOs with previously sequenced strains ( Figure 3A) supports the close relationship between A. lentulus and A. fumigatiaffinis.

Genome Diversity Among and Within Species Across Clinical Strains
Examination of orthogroups across the 15 strains and three species revealed that most genes (7,938) are shared by all three species (Figure 3B). A. fumigatiaffinis has a larger set of species-specific genes (1,062) than A. lentulus (656) or A. fumigatus (645), consistent with its larger genome size and gene number. The numbers of shared genes between A. lentulus and A. fumigatiaffinis are also higher than intersections between each of them with A. fumigatus, consistent with their closer evolutionary relationship ( Figure 3A). Within each species, most orthogroups are found in all strains (9,008, 8,321, and 9,423 in A. lentulus, A. fumigatus, and A. fumigatiaffinis, respectively); approximately 5.4-6.13% of genes in each species appear to vary in their presence between strains (Supplementary Figure S3). Among these, we noted that orthogroups that are present all but one strain are usually the most frequent ( Figure 3C).
We identified a total of 114,378, 160,194, and 313,029 SNPs in A. fumigatus, A. fumigatiaffinis, and A. lentulus, respectively. We identified 406, 493, and 747 SNPs in A. fumigatus, A. fumigatiaffinis, and A. lentulus, respectively, as high-impact polymorphisms; these polymorphisms are those whose mutation is presumed to be highly deleterious to protein function. Similarly, out of a total of 11,698 (A. fumigatus), 20,135 (A. fumigatiaffinis) and 34,506 (A. lentulus) indels segregating within each species, we identified 615, 1,739, and 1,830 high-impact indels in A. fumigatus, A. fumigatiaffinis, and A. lentulus, respectively.

Polymorphisms in Major Antifungal Target Genes Correlate With Antifungal Susceptibility
Given the observed variation within and between species in antifungal drug susceptibility, we examined DNA sequence polymorphisms in genes known to be involved in antifungal susceptibility to azoles and echinocandins. In particular, we examined patterns of sequence variation in the 14α-sterol demethylase gene cyp51A (Afu4g06890) and in the 1,3-betaglucan synthase catalytic subunit gene fks1 (Afu6g12400). Using A. fumigatus A1163 as reference, we identified important speciesand strain-specific polymorphisms in both cyp51A and fks1 ( Figure 4A and Table 2 shows a detailed breakdown of all SNP and indel polymorphisms per strain).
An alignment of Cyp51A protein sequences in the three species shows possible insertions in different sites (Figure 4A red arrow). We observed substitutions in at least one of the clinical strains in the three species in 42 positions that might be correlated with the strains' varying drug susceptibility levels. For instance, Cyp51A in A. fumigatus CNM-CM8714 revealed a welldocumented substitution related to azole resistance at position 98 FIGURE 2 | High heterogeneity of virulence levels among Spanish strains of three closely related Aspergillus pathogens. We found significant heterogeneity in the survival curves between strains within A. fumigatiaffinis (A), A. lentulus (B), and A. fumigatus (C) (Benjamini and Hochberg adjusted p-values using the log-rank test are shown). We also tested differences between species (considering each strain as a biological replicate), and found that the virulence profiles of Spanish strains of all three species are not significantly different (D) (p-value using the log-rank test is shown).
(L98H) (Figure 4A -blue arrows), which might be correlated to its lower susceptibility to itraconazole or voriconazole compared to other A. fumigatus strains (Figure 1B).
We also looked at the promoter region of the cyp51A gene ( Figure 4B) and identified the TR insertions TR34 and TR46 (region highlighted between blue arrows), previously reported in antifungal resistant strains (Dudakova et al., 2017). These changes were specific to certain clinical strains of A. fumigatus and were previously reported in combination with specific point mutations leading to amino acid substitutions. For example, A. fumigatus CNM-CM8714 carries the TR34 promoter insertions combined with L98H ( Figure 4A -blue arrow), whereas A. fumigatus CNM-CM8057 has a TR46 insertion combined with Y121F/T289A (Figure 4A -blue arrow). There are other variants (short indels) that were exclusive to either A. lentulus or A. fumigatiaffinis, or both. The newly sequenced A. fumigatiaffinis strains form a separated group that is closely related to A. novofumigatus. All A. lentulus strains in this work group together and share an ancestor with A. lentulus IFM54703, the only sequenced strain in this species to date. The A. fumigatus strains sequenced in this work form different internal groups in the clade with other strains in the species (e.g., strains CNM-CM8714 and CNM-CM8812 group together and strains CNM-CM8686 and CNM-CM8689 form another group). (B) A. fumigatiaffinis and A. lentulus shares the highest number of common orthogroups and A. fumigatiaffinis displays the highest number of species-specific orthogroups. We considered species-specific orthologs those that were present in at least one strain of a given species, with no representative from another species. (C) Orthogroups shared by all and "all but one" strains are the most frequent in three closely related Aspergillus pathogens. A. lentulus, A. fumigatus, and A. fumigatiaffinis have 9,008, 8,321, and 9,423 orthologous genes present in all strains, respectively. The five largest combinations of orthogroups are shown. As expected, the most frequent combination of orthogroups are those in all strains but one.
Examination of the Fks1 protein sequence alignment from strains of the three species also revealed substitutions in 39 sites ( Figure 4A). We also observed an insertion at position 1,626 of A. lentulus CNM-CM8927 (red arrow). Fks1 also showed substitutions at positions comprising an important hot-spot 2 (HS2) (blue arrows): all A. lentulus strains have a substitution at position 1,349 (I1349V) and all A. fumigatiaffinis have a substitution at position 1,360 (T1360I).
Examination of orthogroups revealed that the orthogroup that includes the cyp51A gene (Afu4g06890) contained additional paralogs of the cyp51 family in A. fumigatiaffinis. Thus, we carried out a phylogenetic analysis with the amino acid sequences with the orthogroups containing cyp51A and cyp51B genes in A. fumigatus Af293 (Figure 4C) that comprises the three species in this work. We observed three well-defined clades. The A. fumigatiaffinis paralog related to cyp51A is likely to represent cyp51C, which has been previously reported in other Aspergillus species, such as A. flavus and A. oryzae (Hagiwara et al., 2016;Perez-Cantero et al., 2020). Sequence identity between the putative Cyp51C protein in A. fumigatiaffinis CNM-CM6805 and  Cyp51C (XM_002383890.1) and Cyp51A (XM_002375082.1) of A. flavus (Liu et al., 2012) is 471/512 (92%) and 391/508 (77%), respectively.

Genetic Determinants Involved in Virulence: Single-Nucleotide Polymorphisms, Insertions/Deletions Across Strains and Within Species Conservation
To explore the genetic underpinnings of the observed strain heterogeneity in virulence we next examined the SNPs and indels in 215 genes that have previously been characterized as genetic determinants of virulence in A. fumigatus (Supplementary Table S5).
Most virulence genetic determinants (146 genes) were found in single-copy in all strains (Supplementary Table S6), whereas 57 genes varied in their number of paralogs across clinical strains (Figure 5). We also identified four virulence determinants that had no orthologs in either A. lentulus or A. fumigatiaffinis, such as Afu6g07120 (nudC), which is an essential protein involved in nuclear movement (Morris et al., 1998), and considered an essential gene in A. fumigatus (Hu et al., 2007). Interestingly, we noted 17 virulence determinants that are present in A. fumigatus and A. fumigatiaffinis but absent in A. lentulus (Figure 5 -top panel), such as Afu8g00200 (ftmD), one of the genes in the fumitremorgin biosynthetic gene cluster (Abad et al., 2010).
Several virulence determinants exhibited larger numbers of paralogs in one or more species. For example, the conidial pigment polyketide synthase alb1 (Afu2g17600), which is involved in conidial morphology and virulence (Tsai et al., 1998), is one of the determinants with highest number of paralogs in A. lentulus and A. fumigatiaffinis (n = 7) when compared to A. fumigatus strains (n = 4). For determinants that contained a gene in at least one strain, we tested correlations between number of paralogs and virulence (lethal time 50: day at which 50% of the larvae were dead, or "ND-end": the number of dead larvae at the end of the experiment) and we observed no significant correlation suggesting paralog number does not associate with virulence.

DISCUSSION
A. fumigatus and the closely related species A. lentulus and A. fumigatiaffinis are important causal agents of aspergillosis (Zbinden et al., 2012;Lamoth, 2016). Importantly, the emergence of antifungal resistance is of increasing worldwide concern (Fisher et al., 2018) and antifungal resistant strains of A. lentulus and A. fumigatiaffinis (Alastruey-Izquierdo et al., 2014) have been identified. Heterogeneity in virulence across different strains of A. fumigatus has also been known for some time (Mondon et al., 1996). Analyses of strain phenotypic and genetic heterogeneity allow us to identify correlations between phenotype and genotype in strains of Aspergillus pathogens.
We found that high heterogeneity exists in drug susceptibility and virulence across different strains of A. fumigatus, A. lentulus, and A. fumigatiaffinis (Figures 1, 2). For one specific antifungal drug, amphotericin B, our results confirmed previous findings that A. fumigatus is more susceptible to amphotericin B than strains of cryptic species (Balajee et al., 2004). Studies on the intrinsic resistance to amphotericin B reported for A. terreus highlight the importance of stress response pathways, in particular heat shock proteins (such as Hsp90 and Hsp70), as well as enzymes detoxifying reactive oxygen species (Posch et al., 2018). Future work involving genomics on the cryptic species will be able to exploit changes in genetic determinants FIGURE 5 | Orthogroups for virulence determinants reveals variable number of paralogs among the three closely related Aspergillus pathogens. We searched for 215 known genetic determinants of virulence in A. fumigatus Af293 in the species of interest and found they were grouped into 203 orthogroups. 146/203 were found in single copy across all strains and are not shown here. The cladogram above the species reflects similarities between strain presence/absence patterns. A. fumigatus Af293 shows a different pattern compared to other strains of A. fumigatus, grouping with one of the A. lentulus strains (CNM-CM8927). This may reflect the phylogenetic divergence of A. fumigatus strain Af293 from other species members. Conidial pigment polyketide synthase alb1 (Afu2g17600) is one of the genetic determinants of virulence with highest number of copies in cryptic species (n = 7) when compared to A. fumigatus strains (n = 4). Gene identifiers in A. fumigatus Af293 are highlighted in bold. Color scale indicates the number of genes found within the orthogroup. involved in amphotericin B susceptibility, although this drug is not commonly used in clinical settings. Interestingly, our PCA identified pairs of positively and negatively correlated antifungal drugs from different classes (Supplementary Figure S2), suggestive of potential synergistic effects (resistance to one drug leads to resistance to the other) and trade-offs (resistance to one drug leads to susceptibility to the other), which could be important for clinical applications.
Comparison of the three species showed that the virulence profiles of A. fumigatiaffinis and A. lentulus strains were not significantly different from the virulence profiles of A. fumigatus strains (Figure 2). This finding is in contrast to a previous comparison of survival curves of the type strains of A. lentulus and A. fumigatus, which found that A. fumigatus is significantly more virulent than A. lentulus (Sugui et al., 2014). The likely explanation for this is our finding that there is significant strain heterogeneity within each species (Figure 2), suggesting that comparisons of individual strains between species are not going to be representative of the variation in virulence within species. While additional testing using diverse models of fungal disease will be required to test the validity of these observations, our findings reinforce the emerging view (Kowalski et al., 2016Ries et al., 2019;Bastos et al., 2020) that examining within-species variation in Aspergillus pathogens is an important, yet poorly studied and understood, dimension of fungal virulence.
The advent of whole genome sequencing boosted our understanding of the biology of the genus Aspergillus (de Vries et al., 2017). Several studies have previously analyzed genomic data of A. fumigatus strains (Abdolrasouli et al., 2015;Takahashi-Nakaguchi et al., 2015), uncovering cyp51A mutations in A. fumigatus populations (Abdolrasouli et al., 2015). Some studies have also used population genomic data for strains of A. fumigatus to gain insights on antifungal drug susceptibility (Garcia-Rubio et al., 2018) or virulence potential (Puértolas-Balint et al., 2019). Correlations between phenotypic traits, such as antifungal susceptibility or virulence and genetic traits have been also studied in other well-studied pathogens, such as the opportunistic yeast Candida albicans (Hirakawa et al., 2015). However, to our knowledge, this is the first study that examines the phenotypic and genetic heterogeneity among strains of species closely related to A. fumigatus.
The three main classes of antifungal drugs comprise polyenes, azoles, and echinocandins, involved in ergosterol composition of fungal membrane, ergosterol biosynthesis, and the cell wall biopolymer (1,3)-β-D-glucan, respectively (Robbins et al., 2017). Due to toxicity to host cells, polyenes are only used in exceptional cases, and first-line prophylaxis and treatment of aspergillosis is usually carried out with azoles (Garcia-Rubio et al., 2017;Garcia-Vidal et al., 2019). Mechanisms of azole resistance involving mutations in the cyp51 genes have been identified in diverse fungi, including in multiple animal and plant pathogens (Parker et al., 2014). In A. fumigatus, research has focused on azole susceptibility testing and correlation with point mutations in the cyp51A gene and TR insertions in its promoter region (Chen et al., 2020;Zakaria et al., 2020). Major changes in protein Cyp51A that correlated with azole resistance include point mutations, such as in positions G54, G138, M220, and G448, or combination of point mutations with TRs in the promoter region, such as the TR34/L98H and the TR46/Y121F/T289A (Wei et al., 2015;Beardsley et al., 2018). In previous studies, alterations such as the insertion of TR34 and TR46 (Dudakova et al., 2017) were only found in the cyp51A promoter of A. fumigatus strains, but these have also been found in other pathogens, such as the wheat pathogen Zymoseptoria tritici (Cools, Hans et al., 2012). Our work explored the promoter region of the cyp51A gene in A. lentulus and A. fumigatiaffinis, two closely related pathogenic species, and identified these promoter region changes only in two strains of A. fumigatus and not in either of the two cryptic species. Interestingly, two of the A. fumigatus strains in this work presented the combined TR34/L98H and the TR46/Y121F/T289A. We also identified other changes in proteins encoded by cyp51A and fks1 that can be used in the future to generate mutants and test the effect of mutations in well-studied wild-type strains of A. fumigatus (Chen et al., 2020).
The evolution of the gene families that contain genes involved in drug resistance might also give us clues on how drug resistance evolves in fungal populations; previous studies (Hawkins et al., 2014;Zheng et al., 2019) report two paralogs of cyp51 in diverse species, including A. fumigatus, A. nidulans, Penicillium digitatum, and Magnaporthe oryzae, while Fusarium graminearum, A. flavus, and A. oryzae have three cyp51 genes (Dudakova et al., 2017). Recently, a study proposed the existence of cyp51C gene arising from a duplication in cyp51B in A. terreus and A. carbonarius (Perez-Cantero et al., 2020). Interestingly, our study found a paralog of the cyp51A gene in A. fumigatiaffinis that likely corresponds to cyp51C. Whole-genome sequence analysis in A. flavus reported substitutions in the three paralogous genes (cyp51A, cyp51B, and cyp51C) in the context of antifungal resistance (Sharma et al., 2018). Novel substitutions identified in cyp51C and modeling of protein changes suggested possible effects on drug binding. Next steps in studies of azole susceptibility in A. fumigatiaffinis strains could include the analysis of this putative cyp51C gene and its role in the organism's observed drug susceptibility profile.
Studies on echinocandins have focused on the (1,3)-β-Dglucan synthase enzyme, encoded by the fks1 gene (Robbins et al., 2017). Particularly, two hot-spots have been studied (Gonçalves et al., 2016). Although most studies report mutations in Candida (Desnos-Ollivier et al., 2008;Garcia-Effron et al., 2008), previous work reported point mutations in fks1 hot spot 1 associated with echinocandin resistance in A. fumigatus (Jiménez-Ortigosa et al., 2017). Our work did not find changes among A. fumigatus clinical strains, and didn't find any mutation in the hot spot 1 of sequences of the cryptic species, which is in agreement with previous study that analyzed fks1 sequences in A. lentulus (Staab et al., 2010). However, we did observe changes in hot spot 2 that were specific to the cryptic species; further examination of these changes with respect to echinocandin susceptibility is an interesting future avenue of research.
Although this study focused on polymorphisms in genes cyp51A and fks1, there is also increasing research on non-cyp51 (Zakaria et al., 2020) and non-fks1 (Szalewski et al., 2018) genetic changes. Future exploitation of genomic data on strains of A. fumigatus and closely related species could also exploit these additional genes. These future studies could also exploit new antifungal drugs, such as olorofim, which has also been tested on cryptic species of Aspergillus (Rivero-Menendez et al., 2019). Finally, future phenotypic and genomic analyses can help us to better understand more complex topics involving antifungal drugs, such as resistance, persistence, and tolerance (as well as the role of tolerance in resistance) (Berman and Krysan, 2020). Given possible emergence of antifungal resistance in agriculture (Hawkins et al., 2019), future work could also exploit correlations in antifungals and the origin of these isolates.

DATA AVAILABILITY STATEMENT
All genomes sequenced as part of this work can be accessed through BioProject PRJNA592352; the raw sequence reads are also available through the NCBI Sequence Read Archive. BioSample and Assembly identifiers are presented in Supplementary Table S4. The data and scripts used in this project are available on the Gitlab repository under https://gitlab.com/SantosRAC/afum_afma_ alen2020.