The pan-genome of Aspergillus fumigatus provides a high-resolution view of its population structure revealing high levels of lineage-specific diversity driven by recombination

Abstract

Aspergillus fumigatus is a deadly agent of human fungal disease where virulence heterogeneity is thought to be at least partially structured by genetic variation between strains. While population genomic analyses based on reference genome alignments offer valuable insights into how gene variants are distributed across populations, these approaches fail to capture intraspecific variation in genes absent from the reference genome. Pan-genomic analyses based on de novo assemblies offer a promising alternative to reference-based genomics with the potential to address the full genetic repertoire of a species. Here, we evaluate 260 genome sequences of A. fumigatus including 62 newly sequenced strains, using a combination of population genomics, phylogenomics, and pan-genomics. Our results offer a high-resolution assessment of population structure and recombination frequency, phylogenetically structured gene presence-absence variation, evidence for metabolic specificity, and the distribution of putative antifungal resistance genes. Although A. fumigatus disperses primarily via asexual conidia, we identified extraordinarily high levels of recombination with the lowest linkage disequilibrium decay value reported for any fungal species to date. We provide evidence for 3 primary populations of A. fumigatus, with recombination occurring only rarely between populations and often within them. These 3 populations are structured by both gene variation and distinct patterns of gene presence-absence with unique suites of accessory genes present exclusively in each clade. Accessory genes displayed functional enrichment for nitrogen and carbohydrate metabolism suggesting that populations may be stratified by environmental niche specialization. Similarly, the distribution of antifungal resistance genes and resistance alleles were often structured by phylogeny. Altogether, the pan-genome of A. fumigatus represents one of the largest fungal pan-genomes reported to date including many genes unrepresented in the Af293 reference genome. These results highlight the inadequacy of relying on a single-reference genome-based approach for evaluating intraspecific variation and the power of combined genomic approaches to elucidate population structure, genetic diversity, and putative ecological drivers of clinically relevant fungi.

Fig 1

(A) DAPCA plot depicting 3 primary populations with clear separation between clusters. The x-axis represents the first discriminant function, y axis represents the second discriminant function, ellipses represent 95% confidence areas. (B) STRUCTURE plots of K = 2 to K = 5 populations support low levels of admixture between clades at K = 3 (highlighted with triangle). (C) Isolates with <0.85 probability of assignment to a single clade, indicative of admixture or incomplete lineage sorting between clades. (D, E) LD decay (r²) for all isolates in each respective clade or all isolates across all clades (in gray). (D) Linear–linear plot of LD decay. (E) Zoomed-in log-linear plot of LD decay, with arrows indicating LD50 (half decay) values. (F, G) LD decay (r²) for sample size normalized isolates (n = 12 isolates, averaged over 20 draws) in each respective clade or for all isolates across all clades (in gray). (F) Linear–linear plot of LD decay. (G) Zoomed-in log-linear plot of LD decay, with arrows indicating LD50 values. The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265. BP, base pair; LD, linkage disequilibrium.

Fig 2. Pan-genome gene family distribution.

(A) The pan-genome of 260 A. fumigatus strains included 15,508 gene families in total, including 8,595 (55.42%) core genes (present in >95% of strains), 3,660 (23.60%) accessory genes (present in >2 and <248 strains), and 3,253 (20.98%) singletons (present in only 1 isolate). (B) The distribution of the number of genomes represented in each gene family. (C) Gene family accumulation curves, including (green) and excluding (yellow) singletons. (D) The distribution of unique accessory gene families by strain. (E) The distribution of unique singleton gene families by strain. The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265.

Fig 3. Distribution of clade-specific gene family gains and absences.

(A) The occurrence of clade-specific gene families in that clade (x-axis), defined as genes that exist in that clade in more than 1 genome (singletons were excluded), by prevalence in that clade (y-axes), and highlights that while many clade-specific gene families that exist in a low number of isolates, others are present in nearly all the isolates of that clade and in no other isolates. The distribution of occurrence frequency varies between the 3 clades, due to differences in clade size, with 1,256 accessory gene families exclusive to Clade 1, 95 exclusive to Clade 2, and 115 exclusive to Clade 3. (B) The occurrence of clade-specific gene family absences for clade-defining absences, where each panel depicts the occurrence of each gene in the 2 clades where the gene is not absent. The distribution of each gene is connected by a gray line. (B1) Gene families lost in all isolates of Clade 1 but present in greater than 90% of the isolates in either Clade 2 or Clade 3 (n = 5). (B2) Gene families lost in all isolates of Clade 2 but present in greater than 90% of the isolates in either Clade 1 or Clade 3 (n = 25). (B3) Gene families lost all isolates of Clade 3 but present in greater than 90% of the isolates in either Clade 1 or Clade 2 (n = 125). The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265.

Fig 4. Phylogeny of A. fumigatus mapped against pan-genome distribution and MAT type.

Overall abundance of accessory and singleton gene families was not structured by clade. The distribution of MAT type was random across the phylogeny, with both mating types present in all 3 clades. MAT idiomorphs were not present in equal proportion, with 144 strains containing the MAT1-1 idiomorph, 105 containing MAT1-2, and 11 strains presenting significant alignments to both MAT1-1 and MAT1-2. All strains with the MAT1-2 idiomorph also contained the MAT1-2-4 gene, including the 11 strains with alignments to both MAT1-1 and MAT 1–2. Strain names with the *** suffix denote the 8 strains demonstrating significant admixture between 2 or more clades. The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265.

Fig 5. GO enrichment.

Significantly enriched Biosynthetic Process GO terms associated with accessory gene families present exclusively in different A. fumigatus pan-genome categories (A–C) and clades (D, E). Enrichment analysis was conducted for each experimental set against a background of all GO terms for all gene families in the A. fumigatus pan-genome; p-values represent results of Fisher’s exact test of the top 10 most significantly enriched terms (for terms enriched at p < 0.05) (A) core, (B) accessory, (C) singleton, (D) Clade 1, (E) Clade 3. Clade 2 had no significantly enriched terms. The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265. GO, Gene Ontology.

Fig 6. Distribution of CAZymes and occurrence of cyp51A mutations across the A. fumigatus phylogeny.

(A) CAZyme profiles display patterns of both clade-specific gene gains and clade-specific gene absences across 28 CAZyme classes. CAZyme classes with differential abundance between the 3 clades were determined by Kruskal–Wallis tests at p < 0.001 after Bonferroni adjustment for multiple comparisons. Gene counts are normalized on a 0–1 scale (by CAZyme class) for visualization. (B) Identification of non-synonymous variants in the cyp51A gene across the phylogeny demonstrated structured occurrence of both known resistance variants (framed in pink) and variants with unknown functional impacts (framed in orange). Reference strain Af293 is highlighted in gray and with triangle. While the Leu98His variants (corresponding to the azole resistant TR34/Leu98His genotype) as well as azole resistant variant Gly138Cys were found exclusively in Clade 2, other characterized variants were scattered in low abundance throughout Clade 1, 2, or both. No characterized resistance variants were found in Clade 3. While some non-synonymous variants with unknown functional impacts occurred frequently in cyp51A and represented changes in a single branch leading to the reference strain Af293 (Glu255Asp, Thr248Asn), others were absent form this branch and absent in Clade 3 (Lys427Glu, Val172Met, Tyr46Phe), while others were low abundance and exclusive to Clade 1 (Ala9Thr, Lys427Arg, Ile242Val, Ala284Thr), or Clade 2 (Ser297Thr, Phe495Ile). The data underlying this figure can be found in DOI: 10.5281/zenodo.5775265.