Nonrandom Distribution of Azole Resistance across the Global Population of Aspergillus fumigatus

Of the 4,049 A. fumigatus isolates, a total of 385 alleles were recorded across the nine microsatellite loci analyzed. All loci were polymorphic, with the number of alleles per locus ranging from 25 to 101. Average Nei's gene diversity for all loci was 0.87, varying between 0.76 and 0.97 for individual loci. There were 2,293 represented multilocus genotypes (MLGs).

Analysis of Bruvo’s distances, to investigate the genetic relationship between isolates, showed no significant bootstrap support for any major nodes (>70%) but some support for nodes connecting a small number of isolates (see Fig. S4 in the supplemental material). A predominantly distinct subset of isolates did form a small, well-supported clade, and these isolates appeared to be fundamentally different from the majority of the other isolates.

Hierarchical clustering of Bruvo’s distances by Ward’s method identified two broadly divergent clades, with 1,003 isolates being assigned to clade A and 1,781 being assigned to clade B (Fig. 1a). Principal-coordinate analysis (PCoA) on Bruvo’s distances further confirmed the existence of these two clusters when visualized, with moderate overlap on the first two PCoA axes, which accounted for the majority of the genetic diversity in the data (46.78% and 34.24% of the total variance, respectively), being seen (Fig. 1b). A 2-dimensional kernel density overlaid on the PCoA scatter plot, used to highlight compactness within the two clusters, suggested the existence of additional, more loosely defined subclusters within clades A and B (Fig. S5). Discriminant analysis of principal components (DAPC) cross validation further confirmed that the two clades were indeed genetically distinct, with the rate of success in assigning isolates to their correct clade using only their allelic profiles being 97%, a rate far exceeding random expectations (Fig. 1c; Fig. S6). Allele richness, which here indicates variations in genetic diversity, differed between the two clades, with group A having 268 alleles and group B having 345 alleles.

The frequency of the TR₃₄/L98 and TR₄₆/Y121F/T289A resistance genotypes was not evenly distributed across the two clades both with and without a clone correction step (Pearson's chi-square test, χ² = 205.39, degrees of freedom [df] = 2, and P < 2.2e−16 and χ² = 152.67, df = 2, and P < 2.2e−16, respectively; Fig. S7). Clade A, which had 639 isolates with an annotated CYP51A allele, consisted of 44.4% TR₃₄/L98H, 19.4% TR₄₆/Y121F/T289A, and 36.2% wild-type (WT) isolates. Clade B, which had 473 isolates with an annotated CYP51A allele, consisted of 12.1% TR₃₄/L98H, 13.7% TR₄₆/Y121F/T289A, and 74.2% WT isolates (Fig. 1d).

Allelic richness (the total number of MLGs) was lower in the TR₃₄/L98H (n = 175) and TR₄₆/Y121F/T289A (n = 74) groups than in the WT group (n = 439). As richness is influenced by sample size, the number of effective multilocus genotypes (eMLGs) was also analyzed; both resistant groups—TR₃₄/L98H (n = 110) and TR₄₆/Y121F/T289A (n = 74)—retained a lower number of eMLGs than the WT group (n = 164), despite rarefaction to the smallest sample size.

According to all three diversity indices (Fig. 2a to c), the Shannon-Wiener diversity index (H; df = 2, P < 0.001), Simpson’s index (λ; df = 2, P < 0.001), and the standardized Stoddart and Taylor's index (G′; df = 2, P < 0.001), the resistant populations were significantly less diverse than the WT population (Fig. 2a to c). This held true with the introduction of a clone correction step, although, as expected, the diversity of populations with the TR₃₄/L98H and TR₄₆/Y121F/T289A genotypes did increase due to the reduction in the frequency of clonal genotypes (Fig. S8). There was a significant difference between the TR₃₄/L98H and TR₄₆/Y121F/T289A genotypes for H (df = 2, P < 0.001) but no significant difference for λ (df = 2, P = 0.268) and G′ (df = 2, P = 0.948). Isolates with the TR₃₄/L98H and TR₄₆/Y121F/T289A genotypes had a higher standardized index of association than the wild-type isolates, demonstrating a relationship among azole-resistant isolates more clonal than that found among wild-type isolates (df = 2, P < 0.001) (Fig. 2d).

There was no significant correlation between pairwise genetic and geographic distances for the TR₃₄/L98H (r² = 0.0299; P = 0.444), TR₄₆/Y121F/T289A (r² = 0.0476, P = 0.434), or WT (r² = −0.0762, P = 0.729) genotype (Fig. S9). Therefore, there were no spatially structured patterns of genetic variation (isolation by distance [IBD]) for either the resistant or the WT isolates. Resistant MLGs were globally distributed, spanning multiple countries and multiple continents (Fig. 3). For example, isolates harboring TR₃₄/L98H with identical microsatellite profiles (MLGs) were collected in Tanzania (Africa), Romania (Europe), and India (Asia). Additionally, TR₄₆/Y121F/T289A-containing isolates were also widely distributed, with one MLG being found in four European countries (France, Germany, Ireland, and The Netherlands) and another closely related MLG being found in both Germany (Europe) and Colombia (South America). Resistant isolates with matching MLGs were also distributed among clinical and environmental sources (Fig. 4). There were six occasions where TR₃₄/L98H-containing isolates with matching MLGs were found in both environmental and clinical sources and five occasions where TR₄₆/Y121F/T289A isolates were found in both sources.

There were 144 short tandem repeat of A. fumigatus (STRAf) genotype duplicates found in both environmental and clinical locations. Of the 999 replicated simulations, zero genotype duplicates were generated, showing that the chance that two STRAf genotypes would randomly occur in both a clinical location and an environmental location without a relationship is very low (P = 0.01) (Fig. S10).

Short tandem repeat (STR) genotypes for the 4,049 isolates used in this study were obtained from the short tandem repeat of A. fumigatus (STRAf) barcoding databank for clinical and environmental A. fumigatus strains, which is based at the Centre of Expertise in Mycology, Radboudumc/Canisius-Wilhelmina Ziekenhuis (CWZ), in Nijmegen, The Netherlands (45). The collection spans 26 countries across four continents (see Fig. S1 in the supplemental material) and includes metadata on source environment (environmentally or clinically sourced) and the presence or absence of azole resistance genotypes TR₃₄/L98H and TR₄₆/Y121F/T289A, where available. All isolates in the database were barcoded using a panel of nine short tandem repeat (STR) loci (namely, STRAf 2A, 2B, 2C, 3A, 3B, 3C, 4A, 4B, and 4C), as previously described (46) and validated (47). All STRAf genotypes were generated at CWZ. Of the 4,049 isolates, a subset of 1,112 had been genotyped for either TR₃₄/L98H, TR₄₆/Y121F/T289A, or wild-type (WT) CYP51A.

Due to the predominantly clonal reproductive strategy of A. fumigatus, we used a multivariate, distance-based clustering approach, as opposed to methods that rely on assumptions such as linkage equilibrium and panmixia, to investigate the presence of the population structure (48, 49). First, considering that pockets of oversampled clones can complicate the identification of structure (50), a clone correction step was performed using the package poppr (v2.8.0) (51) for R software (52). Clone correction was based on geographic origin, source environment, and resistance genotype, in which a single isolate per multilocus genotype (MLG) was represented for each of the three strata. A total of 2,784 isolates were retained following the clone correction step. Before analysis, a genotype accumulation curve was generated to confirm that the STRAf data set had the correct number of loci (power) to discriminate between individuals in a population (Fig. S2) (50).

Genetic distances between individual isolates were determined using Bruvo’s distance in the R package poppr, which utilizes a stepwise mutation model for microsatellite loci (53). Genetic distances were calculated with 1,000-sample bootstrap support. Pairwise distances were hierarchically clustered using Ward’s method with the R function hclust and visualized on an ultrametric dendrogram. Population definitions (K) were selected using the R function cutree, where, on inspection, the tree was parsimoniously divided into two broad clusters (K = 2). A principal-coordinate analysis (PCoA) (54) of genetic distances between isolates was calculated using the R package ade4 (v1.7-11) (55), after ensuring that the distances were Euclidean using the Cailliez transformation (ade4) (56, 57). The first two principal axes of PCoA were retained and visualized using a scatter plot where the two identified clusters were overlaid.

To identify whether populations defined by Ward’s clustering were fundamentally different, a modified discriminant analysis of principal components (DAPC) (58) in which the PCoA axes were used as input variables was used. The resulting analysis identified linear combinations of the PCoA axes exhibiting the strongest discrimination between clusters. DAPC was carried out using the package adegenet (v2.1.1), using cross validation to select the optimal number of variables retained in the initial dimension reduction step.

An interactive, web-based tool, AfumID, using information generated from this analysis was created using the R package Shiny (v1.0.5) and can be accessed at http://afumid.shinyapps.io/afumID. AfumID allows users to input microsatellite genotypes, together with metadata, to receive a brief characterization of their isolate relative to a subset of the STRAf database (Fig. S3). AfumID outputs include a notification as to whether the STRAf profile has been previously found and to which of the two populations the input isolate belongs. Users can also browse the CWZ database and visualize their isolate’s position on the PCoA scatter plot.

The frequency of resistance genotypes across the predicted populations (K = 2) was calculated, and the significance of the association was determined using Pearson’s chi-square test in R. The association was visualized using an extended association plot in the R package vcd (v1.4-4) (59). The association of a resistance allele to a specific population was deemed significant when P was <0.05.

A subset of 1,112 isolates with an annotated CYP51A allele (TR₃₄/L98H, TR₄₆/Y121F/T289A, or WT) was retained to determine the genetic diversity of azole-resistant isolates, which was performed without previously generated population associations. Allelic richness (total number of MLGs), the number of effective multilocus genotypes (eMLGs), Simpson’s index (λ), the Shannon-Wiener diversity index (H), the standardized Stoddart and Taylor's index (G′), and the standardized index of association ( ${\overline{r}}_d$ ) were all calculated using the R package poppr. Rarefied MLG, given here as eMLG, is the number of MLGs for each group at the lowest sample size across all groups. Simpson’s index is a measure of MLG diversity that takes into account the number of different MLGs present, as well as the MLG abundance (60). Shannon-Wiener diversity indices measures genotypic diversity by combining both MLG richness and MLG evenness (61). The standardized Stoddart and Taylor's index measures genotypic diversity while taking into account sample sizes (62, 63). The standardized index of association has been widely used as a tool to detect clonal reproduction within populations (64–66) and is used here to measure the relatedness of resistant and wild-type A. fumigatus isolates. All indices were calculated with (n = 750) and without (n = 1112) a clone correction step, to determine the importance of clones and identify patterns of clonal expansion. Confidence intervals were generated using 1,000 bootstrap samples. An analysis of variance was used to compare the diversity indices of different resistance genotypes, and Tukey’s test was used to identify which genotypes differed significantly.

To examine the global distribution of drug-resistant clones, a subset of resistant isolates with identical MLGs and with at least two or more representative individuals was used to generate minimum-spanning networks (MSN) based on Bruvo’s distances, using the R package poppr. Each node of the MSN represents a single MLG that consists of two or more individuals. Nodes or MLGs that are more closely related have darker and thicker edges, whereas nodes that are more distantly related have lighter and thinner edges.

To test for spatial patterns of genetic variation via isolation by distance (IBD), a Mantel test was performed to elucidate an association between geographic distance and the genetic divergence of azole-resistant isolates (67). Pairwise geographic distances between countries were based on the great-circle distance according to the haversine formula, which was computed using the R package geosphere (v1.5-7) (68). A pairwise fixation index (F_ST) between countries was calculated using the R package hierfstat (v0.04-22) (69) for the TR₃₄/L98H, TR₄₆/Y121F/T289A, and WT genotypes. A Mantel test was performed on the matrices of the represented pairwise F_ST and geographic distances using the R package ade4 and 999 permutations.

A nonparametric bootstrapped (1,000-sample) unweighted pair group method using average linkages (UPGMA) dendrogram was produced using Bruvo’s distances to highlight the relatedness of a subset of resistant isolates with identical MLGs sourced from clinical and environmental locales. We used a nonparametric, Monte Carlo test to assess whether the occurrence of isolates with identical MLGs in both environmental and clinical settings deviated significantly from random expectations. The number of observed genotype duplicates across both settings was used as a test statistic. To generate the expected distribution of the statistic, we simulated environmental and clinical genotypes from the groups’ allele frequencies using the function hybridize in the adegenet package. A total of 999 independent random replicates were used for the test. The P value was calculated as the proportion of simulations greater than or equal to the original test statistic.