Mitogenome Variations in a Global Population of Aspergillus fumigatus

1. Introduction

Aspergillus fumigatus is a ubiquitous airborne saprophytic fungus, and a leading cause of human fungal infections [1]. Playing an important role in the carbon and nitrogen cycles, A. fumigatus’ ecological niche is primarily in organic debris and soil [2]. However, they can grow in a broad range of environmental conditions including indoor, outdoor, and aqueous and even extraterrestrial environments such as the International Space Station (ISS) [3,4]. Reproducing mainly asexually, A. fumigatus spreads via sporulation of minute conidia, roughly 2 – 3 µm in diameter, light enough for airborne dissemination and small enough to penetrate lung alveoli [5]. The combination of widespread environmental habitats and airborne spores leads to a high frequency of incidental inhalation every day for the majority of people throughout the world [6].

Despite the near constant inhalation of spores, a healthy host’s innate immune system is usually able to eliminate the foreign cells, such that clinical symptoms are mild and rarely seen [6]. However, in immunocompromised patients, such as transplant recipients, leukemia patients, and those with acquired immunodeficiency syndrome (AIDS), inhaled A. fumigatus spores often develop into invasive aspergillosis (IA), with mortality rates reaching as high as 90% [5]. Consequently, the World Health Organization (WHO) has designated A. fumigatus a Critical Priority fungal pathogen [7]. Antifungal treatment is essential for effectively clearing the fungal pathogen in acute aspergillosis cases [6]. However, shared fundamental metabolic pathways between fungal and human cells limit the availability of antifungal drug targets and treatment options. Furthermore, the emergence of drug-resistant strains in both environmental and clinical settings has contributed to a noticeable increase in mortality rates, surpassing non-resistant infections by more than 25% [8,9].

For effective growth and reproduction both in vitro and in vivo, like most other eukaryotes, A. fumigatus requires energy generated by mitochondria, through oxidative respiration. In addition, aside from energy generation, fungal mitochondria have also been implicated in various cellular processes, including temperature and oxidative stress responses, pathogenicity, antifungal resistance, iron homeostasis, and fungal dormancy [10,11,12,13,14,15]. In most fungal cells, there may be multiple mitochondria, with each containing one or more mitochondrial genome (mitogenome) that replicates independently of the nuclear genome and of the cell cycle [16,17]. Typically, each mitogenome is a single circular DNA molecule. The core gene content in fungal mitogenomes is relatively well conserved among species [18]. However, among species, the mitogenome sizes can vary significantly, from about 20kb to over 200kb [17,19,20]. Among strains within the same species, variations in mitogenome size and gene content have also been observed [17]. The inheritance of mitochondria in most eukaryotic organisms is uniparental, typically maternal, following a non-Mendelian pattern [16]. However, other modes of inheritance have been observed in fungi (for an in-depth review, see [17]). At present, the mode of mitochondrial inheritance in A. fumigatus and potential variations among strains in their mitogenomes of this fungal pathogen are unknown.

Since mitogenomes are small, relatively abundant, and well conserved yet often distinct among species, mitogenome analysis is a popular method for population genetics, physiological, and evolutionary studies [21,22,23,24]. Indeed, for most animals the mitochondrial cytochrome c oxidase 1 gene (COI) has been adopted as a DNA barcode, a standard for species identification in many animal taxa [25]. In fungi, however, COI is less commonly used due to the occurrence of introns within the gene and limited sequence divergence among species in certain fungal lineages [26]. Nevertheless, due to their small size and relatively high abundance within each cell, mitogenomic analysis could provide insights into fungal evolution, population structure, and reproduction at a fraction of the cost of conventional nuclear genomic sequencing [19].

In recent years, there has been a notable increase in the number of studies focusing on mitochondrial genomics. However, the application of mitochondria-based analyses to investigate fungal population structure, particularly within the Aspergillus species, is not well represented in the literature. Furthermore, investigations into the intraspecies structure or whole-genome level analysis of mitogenomic single nucleotide polymorphisms (mtSNPs) in A. fumigatus are currently nonexistent. Previous studies examining interspecies population dynamics and phylogenetics in the Aspergillus genus have predominantly relied on single gene such as cytochrome b or concatenated alignments of multiple mitochondrial genes, most of which were conducted more than a decade ago. Notably, Kozlowski and Stepien [27] pioneered the use of mitochondrial restriction enzyme analysis to establish interspecies phylogenetic relationships within the Aspergillus genus. Subsequent studies by Wang et al. [28,29] successfully investigated the phylogenetic relationships between A. fumigatus and other pathogenic Aspergillus species using cytochrome b gene sequences. More recently, Joardar et al. [19] published the first complete mitogenomes of nine Aspergillus and Penicillium species, including employing 14 concatenated mitochondrial proteins for inter-species and inter-genus phylogenetic analyses. Although relatively less common than single gene focused assays, genome-wide mtSNPs for phylogenetic and population analyses has been successfully demonstrated in studies on crops, cattle, humans, and selected fungi [30,31,32,33].

With the increasing abundance and availability of whole genome sequencing (WGS) reads in A. fumigatus, we conducted a large-scale intraspecific analysis of mtSNPs in a global A. fumigatus population. Specifically, raw sequencing reads of about 2000 isolates of A. fumigatus were available (as of May 2023) in the National Centre for Biology Information’s (NCBI) Sequence Read Archive (SRA), each including both the nuclear and mitochondrial sequence reads. While the nuclear genomes of these strains have been analysed and published in a diversity of papers and journals [34,35,36,37,38], the accompanying mitochondrial DNA (mtDNA) sequences have received limited attention. Furthermore, due to the higher abundance of mtDNA molecules than nuclear chromosomal DNA in the cell and the nature of gene sequencing, the sequence read depths of mtDNA are generally high, enabling accurate SNP calls [19]. Our study aimed to identify mtDNA sequence variation, examine the potential geographic patterns of the observed variation, and investigate evidence of recombination in shaping mitogenome diversity in this clinically significant human fungal pathogen.

2. Materials and Methods

2.1. Data Acquisition and Processing

NCBI’s SRA, which includes raw data from published experiments, the European Bioinformatics Institute (EBI), and the DNA Database of Japan (DDBJ), served as the data source and repository for this study. A query was performed on the SRA to identify all accessions associated with “Aspergillus fumigatus” and Whole Genome Sequencing (WGS) and Whole Genome Amplicon (WGA) assay types in May 2023. The raw sequence data for each accession was obtained using the sra-toolkit (v3.0.0), and subsequently subjected to quality checking and trimming using Fastp (v0.23.1) [39].

To ensure data integrity, reads were filtered based on quality (requiring at least 60% of bases with a phred score of 30 or higher), length (minimum of 30 base pairs), and complexity. The filtered reads were aligned to a reference A. fumigatus mitochondrial genome (strain Af293, GenBank JQ346808.1) using Burrows-Wheeler Aligner (BWA-MEM v0.7.17-r1188) [40].

The alignments were further processed using samtools (v1.17) [41] to ensure quality, length, and alignment depth. Strains with alignments with an average sequence depth below 10x were considered of poor quality and excluded from further analysis. Duplicate reads were removed using MarkDuplicates in the Picard tool (v2.26.3) and genetic variants were called using GATK tools, including HaplotypeCaller (with the `-ERC GVCF` flag), GenomeDBImport, and GenotypeGVCFs (v4.2.5.0) [42]. The resultant multisample VCF file was subjected to additional filtering using bcftools (1v.16) [41] to remove insertions/deletions.

Country and continent of origin for each strain was initially determined from associated metadata collected from NCBI. In cases where metadata was incomplete, a manual literature search was conducted to confirm the origin information. Accessions with unconfirmed or artificial origins were labelled as having "Unknown" geographic origins. All isolates retrieved and analyzed, including their relevant metadata are shown in Supplementary Table S1

2.2. Data Analysis

2.2.1. Phylogeny

To investigate the relationships among strains and mitogenomes, we conduct phylogenetic analysis using concatenated SNPs. The SNPs were aligned using the program vcf2phylip (v2.0) [43] and the alignment was used to construct a neighbour joining phylogenetic tree using the R programming language (v4.2.3) of the fastreeR package (v1.2.0) [44]. The tree was rooted using the ape package (v5.7) [45] with an A. fischeri isolate as the outgroup. The Interactive Tree Of Life (iTOL) web viewer (v5) [46] was used to visualize the resulting tree and incorporate sample metadata.

2.2.2. Population Genetic Analysis

To investigate the potential influence of geographic isolation on mitogenomic population structure, population genetic analyses were performed using the R statistical software and the R package poppr (v.2.9.4) [47]. The first step involved calculating basic diversity statistics and examining the distribution of unique multi-locus genotypes (MLGs) based on mitogenome SNPs present in the entire population. From the initial multisample VCF, four distinct datasets were created to facilitate subsequent analyses. These datasets included: (1) the full dataset without the A. fischeri root, (2) a clone-corrected (CC) full pooled dataset, (3) a dataset excluding countries with fewer than nine samples and samples with unknown or non-geographic origins, and (4) a country wise CC dataset where only one representative for each unique MLG type from each country was included.

Population structure was explored by performing a Discriminant Analysis of Principal Components (DAPC) on both the full and country wise clone corrected datasets. For both datasets, Kernel (K) values were chosen from where the Bayesian Information Criterion (BIC) began to decrease minimally, based on 10 repetitions using the find.clusters function from the adegenet R package (v2.1.10) [48]. DAPC was then applied to the dataset clustered using the chosen K values using the dapc function from the adegenet R package [48]. For the full dataset, a K of three was chosen and both three and four were chosen for the country wise clone corrected dataset for comparison.

The impact of geography on genetic distribution was further explored through an Analysis of MOlecular VAriance (AMOVA). The country filtered and country wise CC datasets were subjected to AMOVA at both the country and continental levels. Additionally, pairwise Nei's Fst values were calculated between pairs of countries to assess the genetic differentiation among them. To evaluate the correlation between geographic and genetic distances, a Mantel test was performed using the non-CC filtered dataset. This analysis involved comparing the geographic distances (as measured by latitude and longitude) with the corresponding genetic distances (measured as pairwise Nei’s Fst in this implementation) between countries. Significance testing was conducted using 999 permutation bootstrap iterations, with the null hypothesis assuming no genetic differentiation among the geographic populations.

2.2.3. Recombination and Linkage Disequilibrium

Three tests were conducted to assess the presence of recombination in the mitogenome of A. fumigatus. The first test employed the common multilocus Index of Association (I_A) [49,50] to investigate allele associations across different loci. To determine the statistical significance of the observed I_A, we generated 999 artificially recombined datasets through random allele permutation while maintaining allelic proportions at each locus, and evaluated each for I_A. A significantly higher observed I_A compared to the randomly recombined permutations would indicate a lack of random recombination in the dataset, leading to the rejection of the null hypothesis of panmixia (i.e., no association of alleles between different loci). To facilitate comparisons between studies, we standardized I_A by dividing it by the number of loci analyzed, yielding the measure rD. I_A and rD calculations were performed on the clone-corrected dataset, where only one representative of each unique MLG was included. These calculations were conducted using the R package poppr [47].

To further investigate the potential recombination patterns in the mitogenome, the second test investigated allelic associations between all pairs of SNP loci across the mitogenome, using the clone-corrected pooled dataset. For each loci pair, both the coefficient (D) of linkage disequilibrium (LD) and the correlation coefficient (r²) were computed, with statistical significance assessed through a χ² test against the null hypothesis of panmixia (i.e., no association of alleles between different loci). Considering the extensive number of locus pairs, p-values from the pairwise χ² tests were adjusted using the Bonforroni method. Both LD calculations and χ² analyses were performed in R.

The third test complimented the LD analysis as it focused on determining the proportion of phylogenetically compatible pairs of loci using the Four Gamete Test between pairs of SNP loci. Phylogenetic compatibility of two loci is established when it is possible to account for all observed genotypes in the sample without inferring homoplasy (reversals, parallelisms, or convergences) or recombination [23]. For biallelic SNP loci observed in this study, a maximum of four possible haploid genotypes per pair of loci can exist. If no more than three of the four genotypes are observed in the sample, all combinations can be explained without homoplasy or recombination, indicating phylogenetic compatibility. However, when all four possible genotypes are present, the presence of all four cannot be explained without invoking homoplasy or recombination, signifying phylogenetic incompatibility [23]. The analysis was performed in R, and bootstrap procedures were manually implemented with the assistance of the shufflepop function in the poppr package [47]. Visualizations for both LD and the Four Gamete test were generated using the corrplot R package [51].

2.2.4. Mitochondria to Nuclear Genome Comparison

To assess the relative abundance of the mitogenome and the nuclear genome within a cell, we used two methods. In the first, we obtained the average read depth of a single-copy gene, cyp51A (GeneID: 3509526), to represent the nuclear genome and compared it with the average read depth of the mitogenome. The cyp51A read depth for each strain was obtained by aligning the adapter-trimmed and quality-filtered reads from each sample to the reference genome Af293 using the same methodology as described earlier. The entire mitogenome sequence was used to obtain the average read depth of the mitogenome for each strain. The read depth for each nucleotide position for both the cyp51A gene and the mitogenome in each strain was determined using samtools depth (v1.17) [41], incorporating the unmapped positions through the inclusion of the ‘-a’ flag. This approach ensured that even the unmapped positions were accounted for in the read depth calculation for the specific gene fragments. The resulting read depths were then averaged across all nucleotide positions, allowing for a more accurate estimation of the average depth by considering both mapped and unmapped sites. For each sample, we calculated the mitogenome to nuclear genome ratio by dividing the average mitogenome read depth by the cyp51A read depth.

In the second method, we estimated the nuclear genome read depth based on whole-genome sequencing coverage. To estimate coverage, the total count of bases in each sequence file that remained post quality filtering and trimming was divided by the aggregate genome size of the A. fumigatus reference strain, Af293 (GenBank: GCA_000002655.1), minus that of the mitogenome. The obtained coverage was then compared with the mitogenome coverage described in method #1 above. The obtained ratios were then analyzed for their potential phylogenetic and geographic patterns.

3. Results

In this study, we retrieved the raw DNA sequence reads of all the strains of A. fumigatus that had their whole-genome sequences deposited in the SRA by May 2023, and one Aspergillus fischeri sample used as an outgroup for later phylogenetic analysis. The retrieved 1939 strains came from 22 countries representing six continents, as well as from outer space (2 isolates) and from within the International Space Station (2 isolates) (Table 1). At the continental level, 53.43% of the isolates were from Europe, followed by 29.55% from North America, 9.33% from Asia, 1.24% from Oceania, 0.46% from Africa, and 0.1% from South America. About 5.9% of the isolates didn’t have any associated geographic information (Table 1). At the country level, USA had the highest number of isolates sequenced, followed by the Netherlands, Germany, United Kingdom, France, Japan, and Ireland, with the remaining 16 countries each contributed less than 30 isolates (Table 1).

3.1. Distribution of Multilocus Genotypes Based on Mitogenome SNPs

Using the mitogenome sequence of Af293 as reference and the retrieved sequence reads from the 1939 isolates, we identified a total of 39 SNP sites in this global collection of isolates. These 39 SNPs enabled identification of 79 mitogenome multilocus genotypes. Similar to the sample size pattern, at the continental level, the largest number of MLGs were found in Europe, followed by North America, Asia, Oceania, Africa, and South America. At the country level, the highest number of MLGs were from the USA, followed by Germany, France, the UK, Japan, Ireland, the Netherlands, and Spain. The remaining 14 countries each had less than 10 mitogenome MLGs in the SRA.

Among the 79 MLGs, 39 were shared by at least two countries and 30 were shared by at least two continents. The remaining 40 were each found in only one country so far and 35 of these were only represented by one isolate each (Figure 1). The two most frequently shared MLGs were broadly distributed, representing 420 isolates from 11 countries and 418 isolates from 18 countries, respectively (Figure 1). The remaining MLGs were relatively less frequent. However, there were 19 other MLGs each shared among at least three countries by 10 or more isolates.

3.2. Phylogenetic Relationships among Strains and MLGs Based on Mitogenome SNPs

The relationships among strains and MLGs based on mitogenome SNPs are shown in Figure 2. This neighbour-joining phylogenetic tree was constructed from the concatenated 39 SNPs. Overall, a few country- and/or continent-based small clusters of strains and MLGs were found (Figure 2). For example, small clusters containing two or more closely related MLGs were found within several countries such as the US, Germany, Japan, Spain, UK, etc. However, the overall phylogeny revealed that all major clusters of MLGs contained strains and MLGs from different countries and/or different continents. Taken together, the results are consistent with frequent gene flow of mitogenomes among geographic regions but that there are unique mitochondrial genotype(s) within many geographic regions.

3.3. Genetic Clusters of the Global Mitogenome MLGs

We analyzed whether the 79 MLGs could be grouped into distinct genetic clusters. To accomplish this, we used the Discriminant Analysis of Principal Components (DAPC) of pairwise distances among MLGs. The DAPC analysis showed a clustering pattern very similar to that observed in the cladograms, with three groups appearing as the most parsimonious. However, there were four cases of incongruence between the DAPC grouping of three genetic clusters and the NJ tree, highlighted by different colors shown in Figure 3.

The Bayesian Information Criterion (BIC) plot exhibited a sharp drop from a k value of one to two, followed by a less steep decline from two to three, and three to four (Figure 4.A). Subsequent increases in k values showed minimal differences and less pronounced drops (Figure 4.A). Plotting the first and second loadings for each grouping value between two and five revealed distinct clusters for a k value of three (Figure 4.B). Increasing the k value from three to four resulted in cluster #3 being split into two clusters, #3 and #4 (Figure 4.C-D). However, according to the NJ tree, with k=4, one MLG, MLG 74, belonging to cluster #3 showed a closer phylogenetic relationship with cluster #4 MLGs than with the remaining cluster #3 MLGs. Similarly, two MLGs, MLGs 32 and 33, were categorized as cluster #4 despite a closer phylogenetic relationship with clusters #3 and #1.

3.4. Geographic Structure

To investigate the geographic structure of mitogenome SNPs and estimate the potential quantitative relationships between geographic distance and genetic distance among our samples, we conducted AMOVA at the continental and country levels and the Mantel analyses at the country level.

AMOVA was conducted using three different datasets: (1) the non-clone-corrected dataset including all countries with sample sizes greater than eight, (2) the same dataset after clone correction on a country wise basis, and (3) the country wise clone corrected dataset but included only countries with more than nine MLGs. The AMOVA test of the first dataset revealed a statistically significant genetic differentiation among countries, but the differentiation was statistically insignificant among continents. However, after clone correction, no statistically significant genetic differentiations were observed at either the country or the continental levels for both the second and third dataset. Taken together, the results suggested that localized clonal expansion of mitogenome MLGs were responsible for the observed statistically significant genetic differentiations in the total samples among countries.

To further identify which pairs of national A. fumigatus samples were genetically differentiated, we obtained Fst values between pairs of national samples using both the non-clone-corrected data and the clone-corrected data. Similar to the analyses above, only countries with an initial sample size exceeding eight for the original dataset (i.e., dataset #1 above) and exceeding eight MLGs for clone-corrected data (i.e., dataset #3 above) were analyzed. The preliminary pairwise comparison revealed a considerable and diverse distribution of Fst values between countries (Figure 5). Among the 91 paired national samples in the non-clonal-corrected dataset, 51 pairs showed significant (p < 0.01, 999 permutations) differentiation. The two pairings that resulted in the highest Fst values both included Cote d’Ivoire. The two highest Fst values, 0.7538 and 0.7308, were observed between Cote d’Ivore and the Netherlands and New Zealand, respectively. However, after clone correction, there was a substantial reduction in variability among country pairs. In addition, six of the initial 14 countries assessed failed to meet our sample size criteria for comparison, and thus only eight countries were compared (Figure 5). As shown in Figure 5, clone-correction led to the reduction in pairwise Fst values across the board and none of the observed differences between countries were statistically significant in the clone-corrected samples. Consistent with the observations of frequent gene flows, Mantel test revealed slightly positive but no statistically significant correlation between geographic distance and population genetic differentiation between pairs of national populations (the Mantel r statistic was 0.109 and a p-value of 0.271 based on 999 bootstrap iterations).

3.5. Recombination and Phylogenetic Incompatability

Three analyses were performed to investigate the possibility of recombination in the mitogenome in A. fumigatus. Here, only the global sample of 79 clone-corrected MLGs were included in the analyses. In the first analysis of multilocus index of association (I_A) among the 39 SNP loci, our sample showed significant deviations from random mating and recombination, with the overall I_A and standardized rD values of 2.50 and 0.07, respectively. Both values were significantly higher than what would be expected in randomly recombining populations (p < 0.01). In the second test, phylogenetic compatibility analysis of the 39 SNP loci also demonstrated a significant departure from random mating. Approximately 80% of the loci pairs displayed phylogenetic compatibility, which is significantly higher than what would be observed under the hypothesis of random recombination (p < 0.01). In addition, out of the 741 SNP loci pairs, 73 were found to exhibit significant allelic associations, indicative of linkage disequilibrium (p < 0.01, Bonferroni correction). Taken together, these results indicated no evidence of random recombination in the mitogenome of the global A. fumigatus population.

However, signatures of recombination were found in the mitogenomes of this global sample. Specifically, 145 of the 741 pairs of SNP loci showed phylogenetic incompatibility (Figure 6), consistent with recombination between these SNP sites. In addition, upon comparison of the LD and four gamete test results per loci pair, the majority (86%) of the phylogenetically incompatible pairs failed to reject the null hypothesis of recombination (Figure 6). Taken together, our results suggest that while non-random mating and linkage disequilibrium dominate, there is unambiguous evidence for recombination in the mitogenomes of the global population of A. fumigatus.

3.6. Reletive Mitogenome to Nuclear Genome Copy Number Ratio

The average sequencing depth per nucleotide across the entire cyp51A gene varied widely among the 1939 isolates, from a few reads to over 400 reads (Figure 7.A). As expected, the sequencing depths estimated based on cyp51A and the whole-genome coverage are highly correlated with each other (r² = 0.95, p<0.001; Figure 7.B). Similarly, the average sequencing depth per nucleotide across the mitogenome varied widely among the 1939 isolates, from a few reads to over 15000 reads (Figure 7.C). Both the cyp51A read depth and the mitogenome read depth data showed a heavy right skew, indicated by the high values of γ (shape parameter), 1.79 and 2.83, respectively (Figure 7.A and C). Interestingly, the ratio between mitogenome sequencing depth and cyp51A sequencing depth also showed a very broad distribution, from less than 1 to over 100, with one ratio close to 600.

We checked the metadata and citations of strains with either very low or very high ratios. Our analyses identified that ten of the strains with either very high (>200) or very low (<1) ratios were either transcriptome data or synthetic metagenome data but were deposited as whole-genome sequence data. These ten strains were removed from subsequent analyses of mitogenome to nuclear genome ratios. Interestingly, though the raw ratios had a significantly skewed distribution (γ = 1.10), the square-rooted ratios had a normal distribution (γ = 0.18), with a mean ratio of 20.70 among the 1929 isolates (95% CI = (18.06, 23.34); Figure 7.D).

To investigate the potential patterns of mitogenome to nuclear genome ratios among strains and genetic and geographic populations, we compared the ratio patterns across both the DAPC grouping and among countries with sample size (>8). Kruskal Wallis testing revealed statistically significant heterogeneity in ratio distribution within both the DAPC groups and the countries examined (p < 0.05, p < 0.01). Subsequent post hoc pairwise Wilcoxon rank sum testing uncovered significant (p < 0.05) differences between Group #3 and both Group #1 and Group #2. Specifically, the Group #3 strains had an overall higher mitogenome to nuclear genome ratio than those in Groups #1 and #2, 29 compared to 21 and 18, respectively. In contrast, no significant differences were observed between Groups #1 and #2.

Among the 14 compared countries (each country with sample size greater than 8 isolates), 31 pairs (out of 91 total pairs) showed statistically significant (p < 0.01) differences in their mitogenome to nuclear genome ratios. Canada had the highest average ratio of 46.2 mitochondrial to nuclear genome depth ratios, more than double the overall average of 20.7. Conversely, Spain had the lowest average ratio, 11.4, just over half of the overall average. Canada had a significantly higher ratio than 11 of the 13 compared countries. In contrast, Spain had a significantly lower ratio than 6 of the 13 compared countries. Interestingly, China stood out as the only country that did not display any significant differences when compared to other 13 countries, primarily due to the large variations among strains within China. Furthermore, China and Canada were not significantly different, but Canada and New Zealand were (at a 0.05 significance level) despite New Zealand having a higher average ratio than China (22.0 and 20.0, respectively). Interestingly, there is evidence of a statistically significant negative correlation between the mitochondrial to nuclear genome ratio and climate. Comparing the samples from the 14 countries to the average annual temperature for 2021 resulted in a Pearson correlation coefficient of -0.61 (p < 0.05).

4. Discussion

In this study, we conducted a broad analysis of mitogenome SNPs in a global collection of 1939 isolates of A. fumigatus, a critical priority human fungal pathogen. Our analyses identified 39 SNPs and 79 multi-locus genotypes among the analyzed isolates. In addition, we found a range of mitogenome to nuclear genome copy number ratios, from about 11 to over 100, with a mean of ~20. Our examination of concatenated SNP sequences separated the strains into at least three genetic groups. However, though multiple genotypes are often overrepresented in certain geographic regions, there was limited geography-based clustering at either the country or the continent level for the 79 mitochondrial genotypes. Interestingly, evidence for mitogenome recombination was found between several pairs of SNP loci. Below, we discuss the relevance and potential implications of our results.

4.1. Geographic Structuring

The analyzed global population included 1939 A. fumigatus isolates representing 22 countries in 6 continents, including two strains from the International Space Station. Among the 79 MLGs from these 1939 isolates, 19 were each distributed in at least three countries, indicative of continental and global clonal expansion. Indeed, all the sequenced samples from Africa, Oceana, South America, and the International Space Station shared MLGs with samples from other geographic regions and none had private mitochondrial MLGs. Interestingly, the two samples from the International Space Station were genetically divergent from each other, similar to the extensive differences found in their nuclear genomes [3]. On the other hand, 40 of the 79 MLGs were found only in one country each, and AMOVA analysis of the total samples revealed statistically significant genetic differentiations among national samples. However, clone-corrected samples revealed limited genetic differentiations among national samples. Together, the results suggested that localized clonal expansion of certain MLGs played a major role in the observed genetic differentiations among countries in the total sample. For instance, MLGs 11 and 79 were the most abundant, representing 420 isolates from 18 countries and 418 isolates from 11 countries, respectively (Figure 1). However, over half and one third of MLG 11 and MLG 79 isolates were found in the Netherlands and the USA, respectively. Localized clonal expansion can significantly bias allele frequencies, contributing to the observed genetic difference among geographic populations [52].

Our phylogenetic and DAPC analyses showed our samples can be optimally clustered into three major clades. This finding differed from previous studies conducted by Pringle et al. [53], Klaassen et al. [54], and Ashu et al. [55] who found two, five, and eight clades, respectively, in their global samples. These three studies used different markers as ours: they all relied on nuclear genomic data, based on either multilocus sequence typing or microsatellite markers. Interestingly, while both Pringle et al. [53] and Laassen et al. [54] conclude no significant correlation between clades and geography, Ashu et al. [55] found limited but statistically significant geographic contributions to genetic variation with a statistically significant geographic contributions in the non-clone corrected AMOVA, but a non-significant Mantel test. The difference in population structure among the studies likely reflect the different samples used and the difference in inheritance patterns between nuclear and mitochondrial genomes, especially where mitochondria largely inherit in a non-Mendelian pattern during sexual mating in filamentous ascomycetes [56,57].

We note that despite the large sample size (1939 isolates) and broad geographic representation (22 countries from across six continents), the analyzed samples were still limited in their geographic distributions. Specifically, the 22 countries represented only about 12% of the total countries in the world. A previous study by Ashu et al. [55] revealed the Cameroonian population of A. fumigatus to be genetically highly distinct from other national populations. It’s likely that the mitogenomes from the Cameroonian populations will be genetically distinct as well. However, there is yet no genome sequence information from Cameroonian isolates of A. fumigatus. In addition, among the 22 sampled countries, the distributions of isolates were also highly skewed. Specifically, over half of the isolates (1036) originated from western Europe, and an additional 30% (563) were from the United States. Furthermore, we were unable to find the geographic origins for 114 of the isolates. Together, we believe that better data submission and including more strains from under-sampled regions and strains from not-yet sampled regions will reveal additional novel mitogenome diversity in this species.

4.2. Recombination

We observed signatures of recombination in the mitogenome of this global population of A. fumigatus. In A. fumigatus, the mitochondrial inheritance has not been critically examined in sexual crosses in the lab. However, sexual mating through hyphal fusion between strains of different mating types should result in significant cytoplasmic mixing, including the potential for mitochondrial fusion and mitogenome recombination. The observation of cleistothecia development along junction zones between mating partners is consistent with hyphal fusion being a dominant mode of mating in A. fumigatus [58]. Our observed mitogenome recombination in natural populations of A. fumigatus contrasts with those observed in the majority of eukaryotic organisms where mitogenomes are inherited uniparentally in sexual crosses without recombination [57]. However, a growing number of fungal species have shown evidence of mitogenome recombination in their natural populations. These included several mushroom species such as Agaricus bisporus [24], Thelephora ganbajun [59], Cantharellus cibarius species complex [60], and several human pathogenic fungi, such as Cryptococcus neoformans [22], Cryptococcus gattii [23], Candida albicans [61], and Candida auris [62]. Interestingly, studies on other members of the Aspergillus genus have provided potential evidence of intra- and inter-species mitochondrial recombination. For instance, investigations on Aspergillus nidulans using genetic crosses, gel electrophoresis, Southern hybridization, and electron microscopy suggested genetic exchange between mitogenomes [63,64]. Similarly, novel mitogenomes were observed in the progeny of a sexual cross in Aspergillus niger, a related opportunistic human pathogen, caused by the homing of mobile mitochondrial group I introns [65]. In both cases, restriction fragment length differences were observed between parental and progeny mitogenomes, consistent with horizontal transfer of mobile introns. However, in our A. fumigatus natural populations, evidence for recombination was observed between SNPs located across the mitogenome.

4.3. Relative Copy Number Ratios of Mitogenome to Nuclear Genome

Aside from mitogenome SNPs, we also investigated the relative copy number differences between mitochondrial and nuclear genomes among the 1939 strains of A. fumigatus. Wide variations of the ratio, from about 11 to over 100, were found among strains and among geographic populations. Based on the redundancy in their assembled sequence, a previous study estimated a clinical isolate of A. fumigatus Af293 having a mitochondrial to nuclear DNA ratio of ~12:1 [66]. In a more recent study based on real-time PCR of single copy genes in the nuclear and mitochondrial genomes, Neubauer et al. [67] estimated that the mitochondrial to nuclear genome ratio in a wild-type strain AfS35 was about 60(±15) :1. At present, the underlying mechanisms for such variations are unknown. Based on the literature, all the 1929 genome-sequenced strains were grown in favorable laboratory conditions before their mycelia were harvested and DNA extracted for sequencing. While some variations may exist in growth conditions among studies in preparing their strains for DNA extraction and genome sequencing, the observed differences among geographic populations suggest that the identified range of ratios might reflect their adaptations to the natural environments from which the strains originally resided. Specifically, we hypothesize that the relatively high mitochondrial to nuclear genome ratio for Canadian isolates may reflect the low temperature these isolates were adapted to and the need to have a high number of mitochondrial genes to generate more ATP to fulfill their metabolic needs in such environments. Controlled experimental evolution under different temperatures are needed to test this hypothesis. Furthermore, due to potential differences in their capacity to generate ATP among strains with different mitochondrial to nuclear genome ratios, we speculate that such strains may exhibit different properties such as antifungal drug resistance and pathogenesis. Targeted experimental studies are needed to determine the biological significance of the observed variations in mitochondrial to nuclear genome ratios among strains.

5. Conclusions and Perspectives

Our study provides novel insights into mitogenome variations in a large global collection isolates of A. fumigatus, a critical human fungal pathogen. Our analyses revealed evidence for both global and local clonal expansions of mitogenomes. The over-representation of relatively few mitochondrial genotypes in the global sample suggested potential selective advantages of these mitochondrial MLGs and/or their associated nuclear genomes in both natural environments and patients. In addition, evidence of mitogenome recombination was found among several mtSNP loci, indicating a departure from typical uniparental mitochondrial inheritance pattern in eukaryotes. Similarly, the investigation of relative copy number differences between mitochondria and nuclear genomes revealed a wide range of ratios. At present, the potential biological significance of the observed variations in mitogenome to nuclear genome ratios is not known. However, it’s tempting to speculate that the variable ratios likely reflect the dynamic nature of mitochondrial genome replication and maintenance which are regulated by various environmental and metabolic cues. Our study underscores the need for further investigation into the functional implications of the observed mitogenome variations in A. fumigatus. Collectively, these findings help lay the foundation for future studies into the genetic mechanisms underlying the diversity and evolution of this critical fungal pathogen.