Molecular Characterization of Pathogenic Members of the Genus Fonsecaea Using Multilocus Analysis

Members of the fungal genus Fonsecaea causing human chromoblastomycosis show substantial geographic structuring. Genetic identity of clinical and environmental strains suggests transmission from plant debris, while the evolutionary processes that have led to spatially separated populations have remained unexplained. Sequences of ITS, BT2, ACT1, Cdc42, Lac and HmgA were analyzed, either by direct sequencing or by cloning. Thirty-seven clinical and environmental Fonsecaea strains from Central and South America, Asia, Africa and Europe were sequenced and possible recombination events were calculated. Phylogenetic trees of Cdc42, Lac and HmgA were statistically supported, but ITS, BT2 and ACT1 trees were not. The Standardized Index of Association (I A S) did not detect recombination (I A S = 0.4778), neither did the Phi-test for separate genes. In Fonsecaea nubica non-synonymous mutations causing functional changes were observed in Lac gene, even though no selection pressures were detected with the neutrality test (Tajima D test, p.0.05). Genetic differentiation of populations for each gene showed separation of American, African and Asian populations. Strains of clinical vs. environmental origin showed genetic distances that were comparable or lower than found in geographic differentiation. In conclusion, here we demonstrated clonality of sibling species using multilocus data, geographic structuring of populations, and a low functional and structural selective constraint during evolution of the genus Fonsecaea. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


Introduction
The genus Fonsecaea comprises etiologic agents of human chromoblastomycosis, a chronic (sub)cutaneous infection eventually leading to cauliflower-like eruptions on the skin [1,2]. The fungus is present in human tissue in the form of muriform cells. The disease has been reported worldwide, but mostly in tropical and subtropical climate zones, with high incidence in endemic areas [3][4][5][6][7].
Inoculation of contaminated thorns or wooden splinters has been hypothesized to be a main route of infection [8,9]. Thus far the etiologic agents within Fonsecaea are limited to three closely related siblings composing a clearly delimited clade [10]: Fonsecaea pedrosoi, F. monophora and F. nubica. Environmental sampling to recover the species from their supposed natural habitat has been done [8,9]. F. pedrosoi and F. monophora were only rarely encountered. However, the majority of Fonsecaea-like strains concerned non-virulent species, which were not frequently isolated from on human infections [8]. Either the natural habitat of pathogenic Fonsecaea species has to be found somewhere else, or, alternatively, the species have some kind of advantage of being carried by a mammal host. The existence of evolutionary processes supporting the latter hypothesis may be revealed by comparing patterns of variability and distribution of potential etiologic agents.
The pathogenic strains form a well-supported clade in the Chaetothyriales [11], but specific delimitation within this clade is still a debated issue. Analysis of global genetic diversity using AFLP showed that five groups were distinguishable, which were considered to belong to three different species. Fonsecaea pedrosoi was relatively homogeneous and was found nearly exclusively in Central and South America, while F. monophora and F. nubica each comprised several AFLP groups and had worldwide distribution. Cases were found in a tropic climate zone around the equator, while the few clinical cases outside endemic areas were supposed to have been distributed by recent migration of the human host [11].
In the present study, we investigate patterns of variability of pathogenic Fonsecaea species using multilocus analysis of five functional genes with anonymous sequence and AFLP markers. The set of strains analyzed comprised clinical and environmental strains from three continents.

Ethical Standards
The present study has been fully reviewed and approved by Sun Yat-Sen University's Academic Committee. All subjects provided written informed consent and the procedures have been approved by the Sun Yat-sen University Medical Ethics Committee.

Fungal Strains and Culture Conditions
Seventeen strains of F. pedrosoi, 12 of F. monophora, 8 of F. nubica ( Table 1) and one of a neighbouring Cladophialophora species were obtained from the reference collection of the Centraalbureau voor Schimmelcultures Fungal Biodiversity Centre (Utrecht, the Netherlands), in addition to fresh strains recovered from patients, and environmental isolates. Stock cultures were maintained on slants of 2% malt extract agar (MEA) and oatmeal agar (OA) at 24uC.

DNA Extraction and Identification
DNA extraction and quality test were performed as previously reported [12,13]. DNA concentrations were measured with nanodrop DNA concentration detector at 260 nm (Thermo Scientific, U.S.A.).

Degenerate Primer Design, Cloning and Specific Primer Design for Cdc42, Lac and HmgA
Degenerate and specific primers of Cdc42 refer to the study of Xie et al. [14]. The degenerate primers of HmgA and Lac were designed using a complete alignment of the amino acid sequences of species listed in Table 2. Multiple sequence alignments were generated with the software Clustal W [15] using the amino acid substitution matrix BLOSUM62 [16,17]. Highly conserved areas were chosen for degenerate primer design. Degenerate forward and reverse primers were designed with minimal degenerate degree using Primer 5.0 software ( Table 2).
DNA of type strains of the genus Fonsecaea were used as the PCR amplification template. Optimal amplification condition was optimized by temperature gradient PCR amplification. Specific amplicons were purified using gel extraction kit (Qiagen, Germany), cloned using a cloning kit (Promega, Madison, WI, U.S.A.) and confirmed by direct PCR amplification with the primer set M13fw (59-GTA AAA CGA CGG CCA GT-39) and M13rv (59-GGA AAC AGC TAT GAC CAT G-39) according to the manufacturer's instructions. PCR amplicons were then purified with Sephadex G-50 fine (GE Healthcare Bio-Sciences, Uppsala, Sweden) and sequencing was done on an ABI 3730XL automatic sequencer (Applied Biosystems, Foster City, CA, U.S.A.). Sequence data were edited using the SeqMan of Lasergene software (DNAStar, Madison, WI, U.S.A.). The resulting sequences were aligned using BioNumerics software v. 4.61 (Applied Maths, Kortrijk, Belgium). The specificity of these sequences for three genes was confirmed by BLASTx search on GenBank (http://blast.ncbi.nlm.nih.gov). Specific primers for HmgA and Lac were obtained by comparison with the degenerate primer of HmgA and Lac respectively. The resulting specific primers cdc42-SF1s, cdc42-SR1s, Lac-Is, Lac-IAs, HmgA-F2s and HmgA-R12s (Table 2) were subsequently tested with the aim to establish amplification conditions, and then were used to test the 38 strains listed in Table 1.

Multilocus Gene Amplification and Sequencing
PCR amplification and sequencing of ITS, BT2, ACT1 was done according our earlier study [18]. PCR amplification of Cdc42, Lac and HmgA was performed with cdc42-SF1s and cdc42-SR1s, LacIS and LacIAS, and HmgA-F2s, HmgA-R12s and HmgA-R22s, respectively. PCR was performed in a 50 ml volume of a reaction mixture containing 14 ml Go Taq master mix (Promega) containing dNTPs, MgCl 2 , reaction buffer, 2 ml of each primer (10 pmol) and 1 ml DNA. Amplification was performed in an ABI PRISM 2720 (Applied Biosystems) thermocycler as follows: 95uC for 5 min, followed by 35 cycles consisting of 95uC for 45 sec, 49.5uC for 30 sec and 72uC for 1.5 min, and a delay at 72uC for 7 min. Annealing temperature was changed to 52uC for Lac. A seminested PCR was performed to amplify the HmgA gene, the first run with primer HmgA-F2s and HmgA-R22, as follows: 95uC for 5 min, followed by 35 cycles consisting of 95uC for 45 sec, 49.5uC for 30 sec and 72uC for 1.5 min, and a delay at 72uC for 7 min. One ml amplicon of the first run were used as templates for the second run with primer HmgA-F2s and HmgA-R12s under the same reaction conditions. Sequencing of PCR amplicons was done on an ABI 3730XL automatic sequencer (Applied Biosystems, U.S.A.). Sequence data were edited using the SeqMan of Lasergene software (DNAStar Inc., Madison, U.S.A.).

Phylogenetic Reconstruction and DNA Polymorphism
The Cipres Portal (http://www.phylo.org) was used to construct maximum likelihood trees with RAxML v. 7.2.6 for ITS, BT2, ACT1, Cdc42, Lac and HmgA. Maximum likelihood searches for the best scoring tree were made after a bootstrap estimate of the proportion of invariable sites automatically determined the number of bootstrapping runs. RAxML will then automatically determine the point at which enough bootstrapping replicates have been produced [19]. Bootstrap values equal to or greater than 80% were considered significant. After repeated construction for all six markers, the combined single file was used to calculate the standardized Index of Association, I A S [20] using the LIAN 3.5 webserver (http://pubmlst.org). The test options were set to Monte Carlo with 1,000 iterations/random resamplings. The same alignments were used to show split-decomposition trees using SPLITSTREE 4 v. 4.8. The same software package was used to apply Phi test (pairwise homoplasy index) to distinguish recurrent mutations (or homoplasies) from recombination in generating genotypic diversity.
DNA polymorphism analyses were carried out using D NA SP 5.10.00 software. A subset of Fonsecaea strains and genotypes was used to calculate haplotype and nucleotide diversity, as well as Tajima's D neutrality test that is based on the number of pairwise differences and the number of segregating sites in a sample of sequences and the number of parsimonious informative sites [21]. The same software package was used to calculate F ST , showing the genetic differentiation among populations, for ITS, BT2, ACT1, Cdc42, Lac and HmgA.

AFLP Genotyping Assay
AFLP genotyping data were taken from our previous study, where a detailed description of the methodology is provided [11].

Laccase and Homogentisate 1,2-dioxygenase Enzyme Activity Assays
All strains representing F. pedrosoi, F. monophora and F.nubica indicated in Table 1 were tested for laccase and homogentisate 1, 2-dioxygenase enzyme activities. Tests were repeated three times for each strain. Laccase was tested according to Mander et al. [22]. Solid MM with a pH of 5 supplemented with 5 mM 2, 2-azino-di-(3-ethylbenzthiazolinsulfonate) (ABTS) which is oxidized by laccase and results in colored compounds. Cultures were pre-incubated at 25uC for 7 days. Subcultures were cut with a cork borer 2 mm diam and placed at the centre of the plate with three replicates. Diameters of colored metabolite halos were measured from day 1 to day 7. For the homogentisate 1,2-dioxygenase enzyme activity test, we followed Ye & Szaniszlo [23]. Solid MM was supplemented separately with 5 mM L-phenylalanine (Sigma, U.S.A.) and 5 mM Ltyrosine (Sigma) which served as artificial substrates to evaluate the homogentisate 1,2-dioxygenase enzyme activity. Culture conditions were the same as in the laccase test. After two weeks of culture, colony diameters were measured.

Statistics
Metabolite diameters were analyzed by one way ANOVA using Prism 5.0 software, followed by Tukey's HSD Post-hoc test. Mean diameters are the result of triplicate experiments. The error bars indicate standard error of the mean; p,0.05 was considered to indicate a significant difference.

Primer Development for Cdc42, Lac and HmgA
Highly conserved domains were found in Cdc42, HmgA and Lac genes after comparison of sequences downloaded from GenBank (Table 3) and these were used for degenerate primer design. PCRs with degenerate primer pairs Cdc42-F and Cdc42-R, Lac-Ds and Lac-Das, HmgA-F2 and HmgA-R12/HmgA-R22 yielded multiple bands. After cloning and alignment analysis, the specific primers cdc42-SF1s and cdc42-SR1s, Lac-IS and Lac-IAS, and HmgA-F2s, HmgA-R12s/HmgA-R22s were obtained ( Table 2). The sets of specific primers each yielded single PCR products of about 0.85 kb, 1 kb and 0.9 kb, respectively (data not shown). The introns were taken out when used for further analysis. The primer sets proved to amplify all Fonsecaea agents of chromoblastomycosis successfully. To establish an outgroup, degenerate primers were used to amplify the target gene, and multiple bands were cloned and sequenced. BLAST searches using translated amino acid sequences in GenBank showed that the amplified fragments of Cdc42, Lac and HmgA had high homology with published target genes [24][25][26]. The conserved domain search revealed that Cdc42 contains a Raslike GTPase superfamily (aa 1-120 ) which involved a GTP/Mg 2+ binding site (aa 45 -100 ) and switch I and II regions (aa 20-25 , aa  ) [27]. Lac contained a Cu-oxidase superfamily which typically exists in the laccase family [27]. HmgA contained the HgmA superfamily (aa 1-204 ), a hexamer arrangement consisting of a dimer of trimers with which the active site iron ion is coordinated [27].

Multilocus Recombination Analyses
Standardized Index of Association I A S [28] performed using LIAN 3.5 [20]

Strain Polymorphism Statistics
The calculated parsimonious informative sites, monomorphic sites, segregating sites and the total number of mutations are summarized in Table 4. In total 37 strains were used for all six genes. The haplotype diversity for ITS (0.761), BT2 (0.743), ACT1  Table 4).
The values of F ST lie between 0 (panmictic) and 1 (total separation). The tested F ST values based on six genes by comparing geographic origins of the strains (Table 5A) show  Table 5B). The comparisons of geographic origins between continents low values were found (Table 5) suggesting separation of populations. Synonymous and non-synonymous changes of the genus Fonsecaea in amino acid sequence in six genes are listed in Table 4. In total 787 amino acid codons were used for the comparison, and 8 1 st base, 2 2 nd base and81 3 rd base mutations were found within the three species. All 1 st base mutations caused non-synonymous changes, but the 2 nd base and 3 rd base mutations caused synonymous changes. A further analysis showed that the non-synonymous changes in ACT1 and BT2 both did not occur in the functional domain (ACT1aa 135 , BT2aa 81 ), while non-synonymous changes in Lac and HmgA both occurred in functional domains (Lacaa 159 , HmgAaa 38 , aa 88 , aa 164 , aa 175 ). Most non-synonymous changes were observed in F. nubica, where all strains isolated to date originate from chromoblastomycosis patients ( Table 6).

Laccase and Homogentisate 1,2-dioxygenase Enzyme Activity Assay
All strains tested yielded positive laccase activity. Colored metabolites were observed in all three species, but statistical analysis showed that F. nubica had higher enzyme activity than other species (F. nubica vs. F. pedrosoi, p,0.001, F. nubica vs. F. monophora, p.0.05, F. monophora vs. F. pedrosoi, p,0.01) (Fig. 3). The homogentisate 1,2-dioxygenase enzyme activity assay revealed that all strains are able to assimilate L-phenylalanine and Ltyrosine as sole carbon sources; no difference was observed within the three species (data not shown).

Discussion
In the evolution of black fungi (order Chaetothyriales) [30], we witness a functional change from a rock-inhabiting life style prevalent in ancestral Coniosporium (Knufia) and relatives to an increased ability to infect humans and other vertebrates in derived clades. Agents of chromoblastomycosis are particularly interesting because they exhibit a pathogenic phase in tissue, the muriform cell, which shows morphogenetic resemblance isodiametrically enlarging cell clumps of rock-inhabiting Coniosporium (Knufia) species. A functional change in the Cdc42 gene, involved in cellular polarity has been hypothesized [31]. The change of life style seems to have been quite successful in the F. pedrosoi clade, judging from the fact that three related species are nearly exclusively found on humans [32]. Nevertheless the shift was not seen to be reflected in the cytoskeleton-associated Cdc42 gene when compared over the order Chaetothyriales [33].
In the present study, six genes were compared in humanpathogenic Fonsecaea species. ITS was used as a standard for phylogenetic construction. ACT1, BT2 and Cdc42 play a role in cell cycle progression and actin cytoskeleton construction, and are involved in morphogenetic switching, leading to large spherical cells with subsequent cellular division giving rise to the infective muriform cell [34]. Lac and HmgA are well-documented virulence factors of black fungi, and participate in the synthesis of melanin. DHN melanin is negatively charged, hydrophobic and of high molecular weight, and arises by the oxidative polymerization of phenolic and/or indolic precursors [35]. Melanin enhances virulence in black fungi of the order Chaetothyriales [36][37][38][39][40][41][42]. We developed primers to amplify Cdc42, Lac and HmgA which proved to be specific for Fonsecaea. The sequenced genes were aligned and confirmed to be Cdc42, Lac and HmgA using BLAST oine search in GenBank. The genes contained the gene-specific conserved domains when searched with translated amino acid sequences [27].
The phylogenetic trees reconstructed with Cdc42, Lac and HmgA (Fig. 1) yielded high bootstrap support for the three sibling Fonsecaea species, while the ITS, ACT1 and BT2 trees were not supported. The lack of support was probably caused by incomplete lineage sorting, several mutations not having reached fixation. Based on the Standardized Index of Association (I A   S ) and Phi-test using six genes, no recombination events were detected among the three sibling species. This phenomenon is frequently observed in opportunistic members of Chaetothyriales, where clonality seems to be prevalent [43]. The neutrality test with Tajima's D yielded no significant results, suggesting that no positive selection was detected in the sequenced genes indicating a low functional and structural selective constraint during evolution.
Relatively low haplotype diversity was observed within the six genes analyzed. A total of 91 fixed synonymous and nonsynonymous changes were observed in coding regions. The nonsynonymous changes in the cytoskeleton genes ACT1 and BT2 are not responsible for morphogenetic changes [7,44] among the three species because the mutations occurred outside functional domains. The non-synonymous changes in Lac and HmgA both occurred in functional domains (Lacaa 159 , HmgAaa 38 , aa 88 , aa 164 , and aa 175 ) ( Table 6), but did not cause obvious functional changes when catalysis of substrates was tested in vitro. A possible explanation might be that the non-synonymous mutations did not cause any changes in the three-dimensional structure of the molecule. A systematic alignment of 223 plant and fungi laccase sequences showed that there are four signature sequence regions (L1-4) and 12 housekeeping amino acids [45], while the detected non-synonymous mutations (Lacaa 159 ) in this study occurred Table 5. between L2 and L3 and do not belong to a conserved region. DNA sequence alignment of HmgA showed that HmgAaa 38 , aa 88 , aa 164 , and aa 175 are not located in conserved regions either.
Therefore we conclude that the non-synonymous changes within two genes are not linked to functional or structural selective constraints within the genus Fonsecaea. Subsequent studies may  Table 6. Synonymous and non-synonymous changes in DNA and amino acid sequence in ACT1, BT2, Cdc42, Lac and HmgA genes of Fonsecaea spp.  reveal whether such changes have occurred in the analyzed genes in ancestral clades, where dramatic changes in life style are supposed to have taken place. Several studies reported on the molecular epidemiology of the sibling species Fonsecaea [45]. Ribosomal and mitochondrial DNA typing has been used to map the geographic origins of strains [46,47]. The molecular epidemiology of this genus showed substantial geographic structuring in all species with differences between American, African and Asian populations similar to what has been found by Kawasaki et al. [46] in mtDNA profiles. In conclusion, we demonstrated clonality of sibling species using multilocus data, geographic structuring of populations, and a detected low functional and structural selective constraint during evolution of the genus Fonsecaea.