Clinical samples were selected from the Pathogenic Fungi Culture Collection of Pasteur Institute, Iran. Initially, the phenotypically distinct isolates were identified, followed by re-identification through the use of the sequencing related to the ITS region of rDNA. To this end, a total of 99 strains consisting of 95 dermatophyte clinical isolates, including T. rubrum (n = 19), T. interdigitale (n = 20), T. mentagrophytes (n = 6), T. tonsurans (n = 28), E. floccosum (n = 22), as well as standard strains provided from the Persian Type Culture Collection (PTCC) of T. Rubrum PTCC 5143, T. mentagrophytes PTCC 5054, T. tonsurans CBS 130924, and E. floccosum CBS 767.73 were used in the sequence analysis. Clinical features of these isolates are provided in Supplementary Table 1. All the strains were isolated from the skin and hair specimens. The names of the species were determined according to de Hoog et al. (2017).

All fungal strains were cultured on a mycobiotic agar (Merck, Germany) and incubated at 28°C for 4–7 days. DNA was extracted using the phenol-chloroform-isoamyl alcohol according to Makimura et al. (1999).

PCRs were carried out in 50 μl reaction volumes for TEF-1α, BT2, CaM, ACT, HSP70, ITS, and the D1/D2 fragments. Each mixture contained 25 μl of Premix (Ampliqon, Denmark), 3 μl of DNA template, and 0.8 μM of each primer, and the remaining volume was filled with water to reach a final volume of 50 μl. Negative controls (water instead of the fungal DNA) were added to each PCR. PCR amplification with different primer pairs was performed for all the isolates, following the seven loci adapted from the earlier genotyping studies (Table 1). The reaction mixture was initially denatured at 95°C for 5 min, followed by 30 cycles of 30 s at 94°C, annealing at changes in temperatures (from 55 to 67°C) for 40 s, 45 s at 72°C, and a terminal extension step of 72°C for 5 min. Five microliters of the PCR products was electrophoresed on 1% agarose gel in TAE buffer (Tris 40 mM, acetic acid 20 mM, and EDTA 1 mM), and photographed under ultraviolet irradiation.

All seven genes were successfully amplified and sequenced in forward and reverse for all 99 isolates. The PCR products were subjected to the ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA, United States). The forward and reverse sequences of each sample were subjected to ClustalW pairwise alignment using the MEGA7.0.21 software, and edited manually to improve the alignment accuracy (Kumar et al., 2001). Two-web databases of the CBS (for sequencing ITS regions) and BLASTn (for six other gene fragments) were used for identifying and comparing the fungi. All isolates were verified to the species level by sequencing the ITS regions of rDNA.

The bioinformatics data were analyzed and used to diagnose the interspecies and intraspecies nucleotide variation of seven loci in all of the 99 isolates in this study. The sequence data from each gene were aligned using the MEGA7.0.21 software; then, they were manually adjusted and concatenated (4,620 bp). Also, the sequences were analyzed by maximum parsimony (MP) for finding the informative sites of genes of the isolates that belonged to each species and PAUP version 4.0b109 to the Bayesian analysis for drawing the phylogenetic tree (Swofford, 2002). The best-fit model of molecular evolution was estimated in the jModelTest 2.1.10 (Darriba et al., 2012). The programs MrBayes version 3.2 (Ronquist et al., 2012) and RAxML version 8.2 (Miller et al., 2010; Stamatakis, 2014) were run on the CIPRES Science Gateway. Subsequently, two simultaneous analyses with eight Metropolis-coupled Markov chain Monte Carlo (MCMC) chains and incremental heating of 0.2 were run for 20 million generations which were sampled in every 1,000 generations. The convergence of parameter estimates was verified and effective sample sizes were obtained >200 for all parameters using the Tracer version 1.6 (Drummond and Rambaut, 2007). Maximum-likelihood (ML) analyses were performed using the RAxML. Optimization in RAxML was carried out using the GTRCAT option. Bootstrap values for maximum likelihood were 1,000 replicates with one search replicate per bootstrap replicate. This analysis was performed for 99 isolates, an individual gene, and a combination of seven genes. Fusarium solani was used as an outgroup to root the dendrogram.

First, species involving all three sequences, BT2, TEF-1α, and ITS, were deposited in the GenBank and identified for phylogenetic reconstruction. Then, these sequences became the datasets. The combined multilocus (BT2, TEF-1α, and ITS) dataset of the 99 isolates of this study along with the published sequences dataset was aligned using the MEGA7.0.21 software. ML analyses of the mentioned sequences were performed using the RAxML as described above. F. solani was used as an outgroup in a phylogenetic tree, constructed with Maximum likelihood using the RAxML v. 8.0.0 under the GTRCAT model and 1,000 bootstrap replications. Bootstrap support above 95% is shown above the branches.

Since ITS-genotyping is only determined for the T. interdigitale/T. mentagrophytes species complex and not for other species and genes, we used the nomenclature from the studies by Heidemann et al. (2010) and Taghipour et al. (2019) for the ITS-genotyping of the T. interdigitale/T. mentagrophytes complex. We chose the borders of the ITS region of T. interdigitale as deposited under the accession number MK312735 for homogeneity, and did not consider the four nucleotides (GGTT) at the 5′−end of our deposited sequences in the GenBank. For other species and genes in this study, the numeration was used, starting from 1.

Because an MLST database for dermatophytes does not exist, the ST numbers started from 1. The ST was put in numerical tags as a contract, which was based on the case suggested by Bougnoux et al. (2002).

According to the MLST standards, particular allele sequences for each of the housekeeping genes were labeled as a genotype, and each different combination of genotypes was determined as a distinct ST. STs were also analyzed using the eBURST package¹ to determine the possible relationships between isolates based on the single allele differences (Feil et al., 2004; Spratt et al., 2004). The discriminatory power was calculated according to Hunter’s formula (Hunter, 1990). Pairwise genetic distances (p-distances) between species were calculated using the MEGA7.0.21. The possibility of selective pressure at each of the loci was computed by the ratio of non-synonymous nucleotide substitutions (dN/dS) calculated by Nei and Gojobori (1986). Fisher’s exact test was used to determine the associations between the ST and anatomical origin, age, gender, and genotype using the statistical SPSS package (version 19). Values with P < 0.05 were considered significant.

Susceptibility of 95 clinical isolates and four reference strains to griseofulvin, lanoconazole, butenafine (Sigma-Aldrich), and terbinafine (Dr. Reddy’s Laboratories) were determined according to the CLSI M38-A2 broth microdilution method (Wayne, 2008). MIC results were visually read after 4 days of incubation at 35°C. The MIC was defined as the point at which the growth of dermatophyte was inhibited by 80% or more reduction for eight antifungals after a visual inspection in comparison with the control. MIC range, geometric mean, MIC₅₀, and MIC₉₀ were then calculated for the isolates that belonged to each species. All tests were performed in duplicate. Candida krusei (ATCC 6258) and Candida parapsilosis (ATCC 22019) were used as quality controls.

Trichophyton interdigitale and T. mentagrophytes were often isolated from tinea pedis, and E. floccosum was the most common cause of tinea cruris. T. rubrum and T. tonsurans are distributed among different areas (Supplementary Table 1).

PCR of the ITS, D1/D2, BT2, TEF-1α, HSP70, ACT, and CaM genes of the 99 isolates and the standard strains showed bands, sizes of which ranged from 350 to 800 bp. The nucleotide substitutions of different species by the gene regions are represented in Table 2. The gene regions without a nucleotide change are not included in Table 2. The HSP70 sequences did not show any intraspecies variation in all species (Figure 1). All the newly generated sequences were deposited in the GenBank (Supplementary Table 2).

The characteristics of the multilocus obtained from the sequences of the PCR products from T. interdigitale loci were compared to those of the corresponding loci in T. mentagrophytes, T. rubrum, T. tonsurans, and E. floccosum (Table 2). The sequences were obtained from the combination of seven genes which contained 4,620 nucleotides. The HSP70 region of E. floccosum was found to be the shortest gene with only 346 bp, while BT2 was the longest gene of T. rubrum with 794 bp. The number of nucleotide substitutions from the dermatophytes genome sequence was studied, which were not conserved at the amino acid level, and the most amino acid changes were found in T. interdigitale (three) at the locus TEF-1α. For example, the HSP70 sequences failed to reveal intraspecies variation in all species. Also, the ITS, CaM, and D1/D2 sequences did not show intraspecies variation in T. rubrum, adding that the ITS, BT2, ACT, CaM, and D1/D2 sequences did not show intraspecies variation in T. tonsurans. The species with no intraspecies variations (with only one genotype) are excluded from Table 2. In the combination of three loci (ITS, BT2, and TEF-1α), there was a relationship between STs and taxa (except for T. tonsurans that ST1 and 2 were in the same taxa and T. interdigitale that ST1 was distributed among the other STs) (Figure 2), whereas, in the combination of seven genes, this relationship was seen in 21 of the 27 singletons.

A phylogenetic tree for a combined seven loci is shown in Figure 1. Based on the phylogenetic analysis, all isolated strains studies were placed in four supported clades (pp = 1). E. floccosum clade takes a basal position to the Trichophyton clades (Figure 1).

Individual CaM tree failed to separate the clade of T. mentagrophytes/T. interdigitale complex from the clade of T. tonsurans (Figure 1), while individual ITS, D1/D2, BT2, TEF-1α, HSP70, and ACT trees provided altogether the satisfactory resolution in four clades (Figure 1). The results of individual trees showed that the only individual ITS, BT2, and TEF-1α tree can distinguish between T. mentagrophytes and T. interdigitale in distinct clusters and T. mentagrophytes strains are placed in a basal position to the T. interdigitale (Figure 1).

In the data obtained from the MP analysis of seven loci, no informative site was observed to be HSP70, while the most informative locus was TEF-1α. Overall, these results showed that the most informative sites were found in T. interdigitale isolates (Table 3).

Deposited sequence of 22 strains in the GenBank was used for the phylogenetic reconstruction. Twenty-two standard strains were selected in a way that all three BT2, TEF-1α, and ITS sequences were deposited in the GenBank. These 22 sequences along with 99 isolates of this study were subjected to ClustalW pairwise alignment using the MEGA7.0.21 software (Figure 2).

Five genera introduced by the study of de Hoog et al. (2017) (Epidermophyton, Microsporum, Nannizzia, Paraphyton, and Trichophyton) were used in this classification, and other genera were not included in this classification due to the absence of the TEF-1α, BT2, and ITS sequences in the GenBank. The highest ratio of bootstrap-support (100%) was observed in groups Epidermophyton, Microsporum, Nannizzia, and Paraphyton, and the lowest belonged to Trichophyton.

Pairwise distances of the combination of the three loci (based on overall average) were examined, indicating a high degree of overlap (difference 0.001%) between the species T. tonsurans and T. equinum (results are not shown). Also, P-distances of the combination of the three genes showed homology of about 99.994% between T. rubrum and T. violaceum species. On the other hand, P-distances of the combination of the three genes showed significant differences between the species N. nana and N. gypsea with E. floccosum (results are not shown).

A total of 19 out of 21 T. interdigitale isolates were classified in genotype II, and five out of six T. mentagrophytes isolates were classified in genotype II^∗. A total of two new genotypes were identified in this study; one genotype was obtained in two T. interdigitale isolates (XXVII) and the other in one isolate of T. mentagrophytes (XXVIII).

Genotypes (except for the ITS-genotyping of T. mentagrophytes/T. interdigitale complex) and STs numerically start from 1 in the contract due to the absence of the MLST database for dermatophyte species. The positions of the single-nucleotide polymorphisms (SNPs), the different genotypes identified at each of the seven loci, and the determination of the genotype, as well as their frequency, are shown in Table 2. Besides, these numbers correlate with the genotype numbers listed in Supplementary Table 2 which represents how ST is determined and shows nucleotide diversity at the level of dermatophytes species. Regarding the genetic variability, the highest number of genotypes was found within the T. interdigitale (five genotypes) species (Supplementary Table 2 and Table 2). The phylogenetic analysis of species with only one genotype is excluded from Table 2. The analysis of 99 isolates yielded 34 STs which are illustrated in Supplementary Table 2. All in all, 27 STs were singletons (represented by a single isolate).

The eBURST analysis was performed for five species. The discriminatory power of the MLST for T. interdigitale, T. mentagrophytes, T. rubrum, T. tonsurans, and E. floccosum was estimated to be 0.881, 0.9333, 0.4474, 0.3202, and 0.5217, respectively. The 32 SNPs among the seven loci led to 16 non-synonymous changes in encoding amino acids. Considering the species studied, dS occurred more frequently than dN; thus, the ratio of dN/dS was less than one for five loci (protein-coding): E. floccosum: −1.19, T. interdigitale: −1.18, T. mentagrophytes: −0.98 T. rubrum: −0.78, and T. tonsurans: −0.28. Pairwise distances of the combination of the seven loci (based on an overall average) to resolve the boundaries between the Trichophyton species are shown in Supplementary Table 3, whereas the sequence homology between the T. interdigitale and T. mentagrophytes species was noted to be greater than 99.999%. A comparison of these numbers indicates that these two species belong to a complex. The second-most similarities were observed between the species of T. mentagrophytes/T. interdigitale complex with T. tonsurans, and the most distant species were found to be T. rubrum.

By Fisher’s exact test, no significant association was found between ST and the anatomical source. Also, phylogenetic analysis of each locus did not show any statistically significant variation across clades, anatomical sources, and demographic information.

The results of antifungal susceptibility testing for the 99 isolates of dermatophytes against four antifungals are shown in Table 4. Overall, the lowest MIC were observed for terbinafine > butenafine > lanoconazole > griseofulvin, respectively. Terbinafine and griseofulvin had the lowest and highest geometric mean MICs which were 0.01 and 1.23 μg/ml for T. interdigitale and E. floccosum, respectively.