Species boundaries in the human pathogen Paracoccidioides

Accepted Manuscript Regular Articles Species boundaries in the human pathogen Paracoccidioides David A. Turissini, Oscar M. Gomez, Marcus M. Teixeira, Juan G. McEwen, Daniel R. Matute PII: DOI: Reference:

S1087-1845(17)30093-2 http://dx.doi.org/10.1016/j.fgb.2017.05.007 YFGBI 3060

To appear in:

Fungal Genetics and Biology

Received Date: Revised Date: Accepted Date:

31 January 2017 12 April 2017 31 May 2017

Please cite this article as: Turissini, D.A., Gomez, O.M., Teixeira, M.M., McEwen, J.G., Matute, D.R., Species boundaries in the human pathogen Paracoccidioides, Fungal Genetics and Biology (2017), doi: http://dx.doi.org/ 10.1016/j.fgb.2017.05.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

TITLE: Species boundaries in the human pathogen Paracoccidioides David A. Turissini1, Oscar M. Gomez2,3, Marcus M. Teixeira4, Juan G. McEwen2,5, and Daniel R. Matute1, ¶ 1

Biology Department, University of North Carolina, Chapel Hill, N.C., USA

2

Corporación para Investigaciones Biológicas (CIB), Medellín, Colombia.

3

Biology Institute, Universidad de Antioquia, Medellín, Colombia.

4

Northern Arizona Center for Valley Fever Research, Flagstaff, AZ, USA

5

School of Medicine, Universidad de Antioquia, Medellín, Colombia.

¶ Correspondence: [email protected] Biology Department, University of North Carolina, Chapel Hill, North Carolina 250 Bell Tower Drive, Genome Sciences Building Chapel Hill, NC 27510, USA RUNNING TITLE:

Speciation in the Paracoccidioides genus

KEY WORDS: Speciation, introgression, Paracoccidioides spp DATA PACKAGE: http://dx.doi.org/10.5061/dryad.j1f83

ABSTRACT

The use of molecular taxonomy for identifying recently diverged species has transformed the study of speciation in fungi. The pathogenic fungus Paracoccidioides spp has been hypothesized to be composed of five phylogenetic species, four of which compose the brasiliensis species complex. Nuclear gene genealogies support this divergence scenario, but mitochondrial loci do not; while all species from the brasiliensis complex are differentiated at nuclear coding loci, they are not at mitochondrial loci. We addressed the source of this incongruity using 11 previously published gene fragments, 10 newlysequenced nuclear non-coding loci, and 10 microsatellites. We hypothesized and further demonstrated that the mito-nuclear incongruence in the brasiliensis species complex results from interspecific hybridization and mitochondrial introgression, a common phenomenon in eukaryotes. Additional population genetic analyses revealed possible nuclear introgression but much less than that seen in the mitochondrion. Our results are consistent with a divergence scenario of secondary contact and subsequent mitochondrial introgression despite the continued persistence of species boundaries. We also suggest that yeast morphology slightlybut significantlydiffers across all five Paracoccidioides species and propose to elevate four of these phylogenetic species to formally described taxonomic species.

INTRODUCTION Identifying fungi species is a challenging task. Several approaches have proposed using morphology, reproductive isolation (e.g., Ajello 1967, Bao et al. 2004, Liti et al. 2006, Naumov 1996), or genetic divergence (e.g., Burt et al. 1997, Dettman et al. 2003 a,b, Hughes et al. 2009) to recognize species boundaries. Currently, the most widely used of these approaches, the use of genetic divergence as a proxy of isolation employs molecular traits to delimit species boundaries by identifying reciprocally monophyletic groups following the phylogenetic species concept (PSC; Taylor et al. 2003). The PSC has transformed the study of fungal speciation resulting in the identification of previously unrecognized microbial lineages (Behegaray and Caccone 2007). The application of the PSC has replaced species identification based solely on morphology (Taylor et al. 2003, Sites and Marshall 2003). Both of these approaches however, have shortcomings (reviewed in Dettman et al. 2003a,b; Sukumaran and Knowles 2017). An integrative approach is thus required to overcome the limitations of individual frameworks (Frantz et al. 2009, Sites and Marshall 2004, Cadena et al. 2017). For example, coupling molecular data with morphological characterizations and population genetic analyses can reveal genetic discontinuities. Species that differ morphologically are more likely to have completed the speciation process than species that are identical (Taylor et al. 2006; Uyeda et al. 2011; Restrepo et al. 2014). Similarly, identifying gene flow between nascent species can explain why the PSC expectations, such as reciprocal, monophyly are not always fulfilled, even in ‘good species’. In this integrative approach, the phylogenetic framework reveals the existence of potentially isolated groups. Phenotypic differences between the species can then be incorporated into taxonomical descriptions. Finally, population genetic approaches determine whether or not cryptic species have had the chance to interbreed, exchange genes, and possibly even fuse into a single lineage. In addition to revealing species boundaries, this combined approach can also illuminate processes by which species have diverged and uncover the evolutionary dynamics of morphological evolution (Taylor et al. 2000). The order Onygenales contains several of the most prevalent fungal pathogens and has been the subject of several successful applications of molecular taxonomy (Koufopanou et al. 1997; Fisher et al. 2000; 2001; Koufopanou et al. 2001; Kasuga et al. 2003; Brown et al. 2013). The order contains at least 30 genera and an indeterminate number of species (Matute et al. 2006a). Coccidioides, the virulent agent of Valley

Fever, was once thought to be a single species (Koufopanou et al. 1997; Fisher et al. 2000; 2001), but molecular markers revealed it to be a complex of two species. These cryptic lineages were initially identified using a PSC and were later determined to differ morphologically as well as physiologically, facilitating their formal taxonomical description (Fisher et al. 2002). These examples have revealed the importance of precisely identifying species boundaries: human pathogens, from potentially different species could represent different pathotypes. Hence, understanding species boundaries could benefit disease control. Cryptic species have also been found in other Onygenales genera: two in the saprophyte Uncinocarpus (Koufopanou et al. 1997), three in the saprophyte Auxarthron (Koufopanou et al. 2001), and an still undefined number in the pathogens Histoplasma (Kasuga et al. 2003), and Blastomyces (Brown et al. 2013). The genus Paracoccidioides, another human pathogen, has also been the subject of cryptic species research. The fungus is responsible for paracoccidioidomycosis (PCM), a disease endemic to South America that currently infects at least 10 million people (Brummer et al. 1993). Paracoccidioides was described in 1930 (Almeida 1930) and was considered to be a single species for 75 years. Nonetheless, several factors suggested that Paracoccidioides might be composed of several genetic groups. First, it has a large geographic range extending from Mexico to Argentina, which could allow for population structure, or even multiple allopatric species. Second, Paracoccidioides was shown to harbor extensive genetic variability when analyzed by molecular tools such as random amplified polymorphic DNA (Calcagno et al. 1998), restriction fragment length polymorphism (Nino-Vega et al. 2000), electrophoretic karyotyping (Montoya et al. 1997; Feitosa et al. 2003), microsatellite analyses (Nascimento et al. 2004; Matute et al. 2006b) and DNA sequence differences (Hebeler-Barbosa et al. 2003; Matute et al. 2007). Despite these clues suggesting that the genus Paracoccidioides encompassed multiple species, the possibility of cryptic speciation remained unexplored until recently. The use of molecular taxonomy led to the identification of five phylogenetic species in P. brasiliensis sensu lato, Paracoccidioides lutzii and four species in the brasiliensis species complex (Matute et al. 2006a; Carrero et al. 2008; Theodoro et al. 2012; Teixeira et al. 2013; 2014). Attempts to use universal barcoding probes (ITS:(Hebeler-Barbosa et al. 2003; Carrero et al. 2008; Imai et al. 2009); COII: (SalgadoSalazar et al. 2010) revealed almost complete monomorphism within Paracoccidioides which dampened interest in studying speciation within the genus. The application of the

PSC based on gene genealogies to Paracoccidioides revealed the existence of three phylogenetic species with different levels of divergence (Matute et al. 2006a). One of these species was thought to be primarily restricted to Colombia (PS3), while the other two had wider geographical ranges that spanned most of South America (PS2 and S1). A second study, reported a fourth, most diverged species, P. lutzii, that is not only genetically isolated but also morphologically distinct from the rest of the Paracoccidioides phylogenetic species (Teixeira et al. 2009). Finally, a close examination of molecular variation led to the identification of a fifth species restricted to Venezuela (PS4; (Teixeira et al. 2014)). Collectively, these studies suggest that Paracoccidioides species span a range of divergences from recent speciation (PS3 and PS4) to ancient divergence (P. lutzii vs. P. brasiliensis sensu lato). The identification of cryptic species in Paracoccidioides has not been without controversy. A study of mitochondrial loci found polymorphism in the Paracoccidioides genus that was not partitioned along species lines but instead corresponded to geographical location (Salgado-Salazar et al. 2010). Contrary to Matute et al. (Matute et al. 2006a) and Teixeira et al. (Teixeira et al. 2009; 2014), Salgado-Salazar et al. (Salgado-Salazar et al. 2010) argued that only two species were present in the Paracoccidioides complex, while the previously identified PS3, PS2 and S1 were proposed to represent population structure. Mitochondrial markers have higher mutation and substitution rates than nuclear markers which make them popular for detecting population structure and species boundaries (Jahnke et al. 1987; Smith & Anderson 1989; Guillamon et al. 1994; O'Donnell et al. 1998). Nonetheless, mitochondrial markers are also more prone than nuclear markers to cross species boundaries (Petit & Excoffier 2009) which might make them unsuitable for reconstructing species trees (Moore 1995). Combining nuclear and mitochondrial markers in an integrative framework can thus help reveal not only species boundaries but also genetic exchange between species (Moore 1995; Petit & Excoffier 2009; Dupuis et al. 2012). Cryptic species are not only of taxonomic interest. Identifying species boundaries is a prerequisite for determining the relative importance of hybridization in generating new diversity (Abbott et al. 2013). New traits can arise either through new mutations or through gene exchange between species via hybrids. Additionally, hybrids can help provide insight into the speciation process by providing a chance to study which genes can be exchanged and which genes cause reproductive isolation. Despite its importance, fungal hybridization remains understudied (for some exceptions see (Saari & Faeth

2012) and (Leducq et al. 2016)). In the case of Paracoccidioides, little is known regarding the origin and maintenance of genetic variation. Specifically, we do not understand how genetic variation is partitioned within and across genetic clusters. We aim to help fill this void. In this report, we reconcile the disagreement between the nuclear and mitochondrial based species designations in P. brasiliensis. We compiled published sequence data from both nuclear and mitochondrial studies and supplemented them with additional microsatellite data to show that all Paracoccidioides species are largely reciprocally monophyletic at nuclear loci. In contrast, mitochondrial genes have unusually low genetic divergence and are not monophyletic. Mitochondrial variation is also more similar among samples from the same geographic region than between samples from the same species but different regions. We demonstrate that nuclear and mitochondrial trees likely differ because of mitochondrial introgression, a common phenomenon across the tree of life (Petit & Excoffier 2009). Our results indicate that the four species of Paracoccidioides rarely exchange genes even though they are highly sympatric and have the chance to interbreed. Importantly, we show that Paracoccidioides is not fully clonal, that species can exchange alleles (though rarely do), and that several instances of genome admixture occurred at some point in the past. Even though interbreeding has occurred between Paracoccidioides species, the species boundaries have not collapsed as would be expected if there was no barrier to gene flow. Finally, we also show that the four species of the brasiliensis complex differ in yeast morphology and that each of the Paracoccidioides species can be molecularly and morphologically identified. We thus propose the taxonomical recognition of three new species of Paracoccidioides: P. americana for PS2, P. restrepiensis for PS3, and P. venezuelensis for PS4. For clarity, we refer to the species using the current names (S1, PS2, PS3 and PS4) throughout the manuscript and only use the proposed species names after their formal taxonomic description.

METHODS

Isolates

We gathered sixty-five P. brasiliensis sensu lato isolates representing six recognized endemic areas of PCM. Isolates were either collected or donated. The origin of the studied isolates is shown in Figure S1 and Matute et al. (2006a). All isolates were kept in glycerol at -80ºC.

Morphological differences

Conidia. The conidial area of three Paracoccidioides species had been previously measured (Teixeira et al. 2009). We assessed the heterogeneity of this trait between species. To calculate the area of each cell, we assumed that conidia were spherical and that the area of any given cell plane could be calculated as π × r 2. This dataset only included data for one strain of P. lutzii, one strain of PS2 and one strain of S1. For each species, we used data obtained from at least 30 cells (Table 1). Yeast area. We grew yeast cultures at 36ºC for thirteen strains of Paracoccidioides: B17, B10, B3 (S1), B26, B7 (PS2), C4, BAC, CNH (PS3), B18 (S1), V1, V3 (PS4), Pb66, Pb01 and EE (P. lutzii). This amounted to at least two strains per species which also ameliorates concerns of phenotypic variability within species and allowed us to compare both within and between species variance. Live cultures were inoculated and subcultured in a 50mL plastic Falcon tube (Beckton Dickinson Labware, Lincoln Park, NJ) with 5mL BHI + 1% glucose media. Subcultures were agitated (120rpm) for 48 hours. After that time, all subcultures were homogenized by gentle shaking, and 10 µl of the liquid were placed on a glass slide with a drop of Lactophenol Cotton blue. All slides were examined with a Nikon Eclipse E200 microscope (Nikon, Japan). Pictures were taken with an Optika Vision camera (OPTIKA Microscopes, Italy) and exported to the Micrometrics SE Premium software (ACCU-SCOPE INC., Commack, New York, USA) following the manufacturer’s specifications. The number of cells observed for each experiment are shown in Tables 2 and 3. We measured two traits: the diameter of the cell and the number of buds per yeast cell. Yeast area was calculated assuming each cell was a sphere in a similar manner to the calculation of conidial area (see above).

Statistical tests. All statistical tests were conducted in R (Team 2014). We tested whether there was heterogeneity in any of the three measured traits (one for conidia, two for yeast) using linear models. We also assessed whether any between species pairwise comparison was significant. These comparisons hold the potential to reveal whether any single trait (or a combination of them) differentiates among the phylogenetic Paracoccidioides species and is thus diagnostic. For conidia measurements, we found that the residuals of the linear model were not normally distributed for any of the three traits (Shapiro-Wilk normality test: W = 0.9002, P < 1 × 10-15). We thus used nonparametric tests for the conidia measurements. We used a Kruskal-Wallis rank sum test with multiple post-hoc comparisons (function kruskalmc, R package ‘pgirmess’, (Giraudoux 2011)) in which the traits were the response and the species was the only comparison level. P-values were automatically adjusted for multiple comparisons. We fit linear models for the yeast traits (one per trait), in which the trait was the response, the species and the isolate were fixed effects, and the isolate was nested within the species. Pairwise comparisons between species were carried out with a Tukey Honest Significant Difference posthoc test with a correction for multiple comparisons (R, library ‘multcomp’ (Hothorn et al. 2008)). For interspecific differences to be useful in identifying species it is necessary to show that an individual sample can be placed in the correct species. We used measurements for two lines per species and tried to assign individuals from a third line to the correct species (B16 for S1, C12 for PS3, EE for P. lutzii, B15 for PS2, and V6 for PS4). To do so, we calculated the Mahalanobis distances of each to-be-assigned individual to the phenotypic distributions for each of the five species. The distance is based on both the mean and variance of the predictor variables (i.e., area of individual cells and number of buds per yeast cell), plus the covariance matrix of all the variables (Mahalanobis 1936). The Mahalanobis distance between a to-be-assigned individual (i) and the distributions for each of the five species was calculated as: MDij = (Fi – μLB)T × SLB-1 × (Fi – μLB) where T denotes matrix transposition, S is the covariance matrix of a given dataset (the genotypes of the individuals used to characterize each species), Fi is the vector of

phenotypic observations in a to-be-classified individuals I, and μLB is the vector of average phenotypic observations for the individuals used to characterize each species. We calculated the Mahalanobis distance between each of the yet-to-be-assigned individuals and each of the five species and then assigned the individual to the species with the shortest distance. Distances are shown in Table S4 and could not be converted to p-values due to the non-normality of their distribution. The proportion of individuals assigned to the correct species was used as a proxy of efficacy of the combination of morphological traits to identify species identity.

Molecular polymorphism data

DNA extraction. We extracted DNA from the 65 Paracoccidioides isolates (see above). Mycelia from each Paracoccidioides strain were grown in modified liquid synthetic medium (McVeigh and Morton) for 15 days at 20ºC in a gyratory shaker set at 120 rpm (Gallenkamp, Leicester, England). Subsequently, mycelia were washed with sterile distilled water and macerated thoroughly in liquid nitrogen using a cooled pestle and mortar (van Burik et al. 2009). Macerated mycelia (0.1 g) were immediately transferred into a micro-centrifuge tube (1.5 ml). Additionally, total DNA was extracted from each yeast culture using either of two protocols, either glass beads (van Burik et al. 2009) or maceration of frozen cells (Morais et al. 2000).

DNA Sequencing. We amplified ten regions from five different supercontigs of the Paracoccidioides genome (Table 4). These regions were identified by mining EST libraries and searching for repetitive regions (i.e., microsatellites) as described in Matute et al. (Matute et al. 2006b; See below "Microsatellite marker development and sequencing"). None of these regions coded for proteins as listed in Table 4. We used OLIGO 4.0 (National Biosystems, Plymouth, Minn.) to design the PC and sequencing primers (Operon Technologies Inc., Alameda, CA, USA). The primers and their respective annealing temperatures are listed in Table S1. All regions were amplified using PCR following previously published protocols (Matute et al. 2006b). Sequencing of PCR fragments was done using the cyclic reaction termination method (Sanger sequencing) using the BigDye Terminator Cycle Sequencing Kits (Applied Biosystems, Foster City, CA, USA). To minimize false polymorphism, we sequenced both DNA strands, aligned them, and examined them with Sequence Navigator v. 1.0.1 (Applied

Biosystems, Foster City, CA, USA). Since Paracoccidioides isolates are haploid, no heterozygous sites were expected (Almeida et al. 2007).

Public Data. We took advantage of the fact that previous work had generated population-wide datasets for the Paracoccidioides genes. We used previously published sequences for seven nuclear loci (Matute et al. 2006a; 2007; Almeida et al. 2009; Theodoro et al. 2013) and five mitochondrial loci from four Paracoccidioides species (Salgado-Salazar et al. 2010). The accession numbers are shown in Table S1.

Phylogenetic Reconstruction

We collected sequence data for 12 loci that were publically available: five mitochondrial, and seven nuclear protein coding (see 'Public data’' section). It must be noted that since mitochondrial loci show broadly similar patterns (Salgado-Salazar et al. 2010) and mitochondrial recombination is rarer than nuclear recombination (Ladoukakis & Eyre-Walker 2004; Fritsch et al. 2014), all mitochondrial loci were considered as a single locus. Additionally, we sequenced ten noncoding regions associated with single repeat sequences (see above, “Molecular polymorphism data” section). Our aim was twofold. First, we assessed whether phylogenetic species of Paracoccidioides were reciprocally monophyletic at nuclear and mitochondrial loci. Second, we investigated whether nuclear and mitochondrial topologies were consistent with each other. Alignment. Sequences from each locus for the 66 individuals (65 Paracoccidioides brasiliensis sensu lato and one P. lutzii, the outgroup) were individually aligned using Clustal Omega with default parameters (Sievers et al. 2011). Briefly, Clustal Omega uses the hHalign algorithm (Söding et al. 2005). The gap opening penalty was 6 bits, while the gap extension penalty was 1 bit. The subsequent alignments were visually inspected to ensure indels were aligned properly and to confirm nucleotide homology using Mesquite version 3.04 (Maddison & Maddison 2001). We separately concatenated nuclear coding loci (66 taxa, 2,520 characters), nuclear noncoding loci (65 taxa, 6,741 characters) and mitochondrial loci (52 taxa, 2,731 characters) since we were interested in comparing the phylogenetic signals of each data set.

Topologies. Loci were concatenated into three datasets: mitochondrial, nuclear coding, and nuclear noncoding. Analysis files were formatted as nexus files and converted with BEAUti (Drummond et al. 2012) specifying a Gamma site model and a random local clock. For tree priors, we used a pure birth (Yule 1925) tree with Poisson rate changes and, the best-fitting model according to jMODELTEST (Posada 2008) was JC69. The chain length for all datasets was 1 × 10 8 steps. Each file was stored as a .xml file and imported into BEAST (Drummond et al. 2012; Bouckaert et al. 2014). All analyses were run five times to confirm that the different runs converged on the same likelihood and tree topology. Topologies were rooted with P. lutzii, the sister species to the Paracoccidioides brasiliensis species complex. We tried burnins of 1 × 106, 5 × 106, and 1 × 107 steps, and the results remained unchanged. Estimates of branch support were generated with TreeAnnotator (Rambaut & Drummond 2013) and visualized with FigTree v1.4.2 (Rambaut 2009). We also built maximum likelihood trees using RAxML (Stamatakis 2006; 2014). RAxML was run with parameters m=JC69 and –no-bfgs for 1,000 bootstrap replicates and was rooted with P. lutzii.

Quantifying phylogenetic discordance. To quantify the extent of phylogenetic discordance among nuclear and mitochondrial genes, we used the approximately unbiased (AU) test implemented in CONSEL (Shimodaira & Hasegawa 2001). The goal of this analysis was to determine whether both sets of loci supported the hypothesized species tree (Matute et al. 2006a) or whether they supported alternative topologies (Salgado-Salazar et al. 2010). Since we aimed to detected large inconsistencies in the phylogenetic relationships, we used consensus sequences for each of the four species of the brasiliensis complex. Using all of the individual isolates would result in 64! possible topologies mainly driven by different within species relationships. We tested the support for four topologies (Figure 1). We first produced maximum likelihood topologies for the concatenated nuclear and mitochondrial alignments using RAxML (Stamatakis 2014). We used the JC69 molecular model for both datasets as stated above. Next, we estimated the level of support for each topology using the approximately unbiased (AU) test which implements LRT tests in a similar manner as described for STRUCTURE (see below). We also determined how many loci supported each topology by generating individual gene genealogies. Each gene genealogy was built with RAxML using the same parameters as in the concatenated sequence analyses.

We also performed a second set of AU tests. We tested the hypotheses of hybridization between S1 and P. lutzii and between PS2 and P. lutzii (Theodoro et al. 2012). To this end, we built topologies that included the two most divergent species of the brasiliensis complex, S1 and PS2, P. lutzii and a highly divergent species, Histoplasma capsulatum Nam1 (Wang et al. 2009). Histoplasma capsulatum Nam1 homologous sequences were obtained through reciprocal BLAST to the publicly available Histoplasma genome (Bioproject accession number: PRJNA12654). The three possible topologies between these species are shown in Figure S2. We used this streamlined phylogenetic tree to detect possible gene exchange between P. lutzii and the two species of the brasiliensis complex that are known to occur frequently in Brazil. The support for each topology was calculated as described immediately above. Population genetic analyses

Nuclear data: We calculated the number of segregating sites, the nucleotide diversity, and Tajima’s D statistic to test for neutral evolution at each locus using MEGA 6.06 (Tamura et al. 2013). For each phylogenetic species, we calculated nucleotide polymorphism (π) and the average pairwise divergence between groups (Dxy). In the case of coding regions (N= 6 loci), we calculated diversity and divergence using the third position of the codon. We restricted our analysis to potentially neutral sites (πS). For noncoding regions (N= 10), we used all sites (πNC). For intraspecific diversity and interspecific distance, we calculated the variance of the estimates by bootstrapping the dataset and generating 1,000 replicates. Microsatellite markers: Coding loci evolve slowly which might make them uninformative when assessing whether species are exchanging genes. We used fast evolving microsatellites to assess the possibility of gene flow between species. Isolates were genotyped using ten microsatellite loci (simple sequence repeats, SSR markers; Table S1). Five of these markers were previously published, and we developed an additional five. Public data. A previous study developed five microsatellites that could differentiate the three species of the Paracoccidioides species complex (Matute et al. 2006b). The accession numbers for those sequences are: DQ313963-DQ314154, DQ313835-

DQ313898, and DQ314091-DQ314154. We used those sequences to develop primers that would amplify the loci in the four phylogenetic species of the brasiliensis complex (Table S1). Microsatellite marker development and sequencing. We followed the same protocol described in Matute et al. (Matute et al. 2006b). Briefly, we mined EST libraries for repetitive sequences and developed primers to amplify the region containing the repetitive elements and ~200bp on each side of the tract. We identified ten suitable markers distributed on eight supercontigs of the Paracoccidioides genome (Table 4). PCR primers and conditions are listed in Table S1. PCR products were visualized on a 1.5% agarose gel to confirm amplicon size. Amplicons of the expected size were lyophilized and sent to Macrogen, Korea (South Korea) using express mail. PCR amplicons were purified and sequenced using Sanger sequencing and BigDye terminator chemistry on an ABI 3730 instrument (Applied Biosystems). Both ends of the amplicon were sequenced to correct for sequencing error. All electropherograms were visually inspected for sequencing quality and a single haplotype per strain was generated using the two sequencing products (forward and reverse). Alignments of the repetitive region and the flanking sequences were done as described above (See‘ 'Phylogenetic reconstruction: Alignment’'). All loci were blasted against the B26 (Pb03; PS2) genome (http://genome.jgi.doe.gov/pages/blast.jsf?db=Parbr1; (Desjardins et al. 2011)) to determine their genomic location. Allele counts are shown in Table S2. All sequences were deposited in GenBank under the accession numbers KX719990KX720560. Population Structure. To assess between-population differentiation and the contribution of individual alleles to population structure, we used Principal Component Analysis (PCA), a general method for summarizing high dimensional data (Novembre & Stephens 2008). PCA transforms a set of possibly correlated variables into a reduced set of orthogonal variables. Sampled individuals are then projected in a two dimensional space where the axes are the new uncorrelated variables, or principal components. Related individuals aggregate in the bidimensional space, thus revealing population relationships or within-population stratification. We used the R package ‘adegenet’ to run a PCA analysis for the combined microsatellite data across species (Jombart 2008; Jombart et al. 2009). The allele frequencies were simply the counts of each allelic class as

Paracoccidioides is haploid (Almeida et al. 2007). We plotted the first five principal components. 95% confidence ellipses were drawn using the dataEllipse() function from the ‘car’ package in R (Fox & Weisberg 2011). To estimate the number of populations that best explained the genetic variance in the isolates we studied, we used the Bayesian model-based clustering program STRUCTURE v2.3 (Pritchard et al. 2000). To ascertain the most likely clustering in the P. brasiliensis species complex we inspected two different summary statistics: LnL (Falush et al. 2003) and ∆K (Evanno et al. 2005). LnL compares the log-likelihood of two different clustering schemes using a likelihood ratio test (LRT). LRT follows a χ2 distribution and compares hierarchically nested models. We used it to determine which clustering model best fit the observed genetic variance. Since adding an extra population constitutes an additional parameter, models with different K values can be compared with LRTs. We did sequential LRTs to determine when adding a new population improved the model likelihood. We also used the ∆K method (Evanno et al. 2005). This method infers the most likely number of clusters by evaluating the rate of change in the data log probability between successive K values for each re-sampled dataset. We ran STRUCTURE with K, the number of populations, set to 1 through 8. Each run involved 100,000 MCMC steps with a burn-in of 100,000 and used the following parameters: NOADMIX=0, LINKAGE=0, INFERALPHA=1, ALPHA=1.0, UNIFPRIORALPHA=1, ALPHAMAX=10.0, and FREQSCORR=0. Divergence time

No fossil record is available for the dimorphic fungi of the Onygenales order. To avoid relying on deeply divergent nodes to calibrate our tree, we used two approaches that give approximate divergence ages. The first approach used indirect and sequential calibration (e.g., Santos et al. 2009) assuming that P. lutzii and the brasiliensis complex diverged 33 million years ago (sensu Teixeira et al. 2009). This estimation assumes constant population sizes, and a clock-like accumulation of substitutions that follows the function: A = 2µT, where A is the number of synonymous substitutions, µ is the mutation rate and T is the divergence time. (The factor of 2 represents twice the split time in a haploid population.) The assumed substitution rate for nuclear DNA was 1 × 10-9 substitution/site/year, while the assumed rate for mitochondrial substitution was 6.80 × 10−7 substitution/site/year (following Torriani et al. 2014). However, it is worth noting that

nuclear and mitochondrial mutation rates have not been measured in Paracoccidioides as has been done in yeast (Lynch et al. 2008) and fission yeast (Farlow et al. 2015). Second, we established a possible range of ages for each of the speciation events assuming different rates of substitution across the tree (i.e., different models of molecular clock). We estimated the branch lengths of each of the speciation events in Paracoccidioides using coding sequences and five types of molecular clock on BEAST: strict clock, random local clock, fixed local clock, uncorrelated relaxed clock with an exponential distribution, uncorrelated relaxed clock with a lognormal distribution and uncorrelated relaxed clock with a gamma distribution. For the uncorrelated relaxed clock models, we fitted models with and without a continuous quantile parametrization. Each of these estimations was done using BEAST’s specifications described above (See above “Phylogenetic Reconstruction” section). These branch lengths were then converted to relative ages and are expressed in relative units of divergence where S1 and P. lutzii have a divergence of 1. Both sets of estimations constitute an indirect source of evidence, but the information provides a relative scale of the divergence between Paracoccidioides species.

Assessment of gene flow between species

Direction of gene flow. We looked for evidence of introgression in the nuclear genome and its possible directions using the model implemented in TreeMix (Pickrell & Pritchard 2012). TreeMix estimates the most likely evolutionary history of a group of populations by estimating the levels of genetic drift at a set of markers. The analysis is done in two steps. First, it estimates the relationships between sampled populations and estimates the most likely dendrogram. Second, it compares the covariance structure modeled by this dendrogram to the actual observed covariance between populations. If a pair of populations is more closely related than expected by the strictly bifurcating tree, then the model suggests an admixture event in the history of those populations and determines the most suitable topology and the direction of migration events. Since we had allele frequency data for four populations, we could estimate between 0 and 3 migration events (m). TreeMix was run allowing for m=0, 1, 2, or 3. We split S1 into two populations based on the STRUCTURE results (S1 Brazil and S1 southern South America). We determined the most likely demographic scenario in terms of migration between species in two ways. First we compared the log likelihood values of

each demographic event using a Likelihood Ratio Test (LRT) with two degrees of freedom to account for changes in the number of migration edges and the direction of gene flow. Second, we calculated the Akaike Information Criterion for each demographic scenario (Akaike 1974). To determine what migration scenario (i.e., number of migrations) best explained the data, we compared the AIC values for each migration model by ascribing weights to each model using weighted AIC (wAIC) using the following equation (Burnham & Anderson 2003):

RESULTS

Paracoccidioides species differ in their yeast morphology

We studied three different morphological traits: conidial area, yeast cell area, and the number of buds produced by yeast mother cells. The results for each of the traits are described as follows:

Conidia. We compared the mean conidial area between pairs of the three out of five Paracoccidioides species for which we had data (P. lutzii, S1 and PS2). A Kruskal-Wallis test revealed significant heterogeneity among species (Kruskal-Wallis χ2 = 167.3063, df = 2, P < 1× 10-15, Figure 2). Table 1 shows the mean conidial area and standard deviation for each of the three sampled species. Additionally, it shows the three possible posthoc pairwise comparisons (Kruskal-Wallis multiple comparisons corrected for multiple comparisons). We found that conidia could be used to distinguish between P. lutzii and the two species from the brasiliensis complex. Nonetheless, PS2 and S1 did not significantly differ in conidial size indicating that this trait does not contain sufficient taxonomic information to separate species within the brasiliensis complex.

Yeast cell area. We measured the yeast stage cell area of the five phylogenetic species of Paracoccidioides (at least two strains per species). Additionally, given that B18 showed conflicting phylogenetic information and may be either misassigned or of hybrid origin (see below, Figure 3), we categorized that isolate as an independent genotype (i.e., the ‘species’ effect on linear model had six levels). We detected significant heterogeneity in yeast cell area across genotypes (Species effect, Linear model: F 5,373 = 15.5863, P = 6.19 ×10-14, Figure 2) indicating phenotypic differentiation among species. We found no significant intraspecific variation (Isolate effect, Linear model F1,373 = 2.19, P = 0.1398). Table 2 and Figure 2 show the mean yeast area and the standard deviation for each of the five species and B18. Pairwise comparisons revealed that P. lutzii has a different cell size from PS2 but not from S1, PS3, or PS4 (Table 2). Within the brasiliensis complex, cell area differs in five out of the six possible comparisons and the mean value is a diagnostic trait for PS2 and S1 because it differentiates these two clades from all other species. PS3 also differs from S1 and PS2 in cell size but not from PS4 (Table 2). This last result is not surprising given that they are sister species. Our

results indicate that yeast area is an effective diagnostic trait for most, but not all, species of Paracoccidioides.

Number of buds in yeast cells. The third and final trait we studied was the number of buds produced by each of the species in the brasiliensis complex. We counted the number of buds sprouting from the mother cell in the yeast life stage. We detected significant heterogeneity in the number of buds across genotypes (Species effect, Linear model: F5,318 = 6.8865, P = 4.039 × 10-6, Figure 2) indicating phenotypic differentiation among species. Again, we found no significant intraspecific variation (Isolate effect, Linear model: F1,318 = 1.55, P = 0.2139). Table 3 shows the mean number of buds per yeast for each of the five species and B18. We found that PS4 shows a significantly lower number of daughter cells when compared to any of the other Paracoccidioides species (Table 3). Notably, the number of buds per yeast cell can differentiate between PS3 and PS4, the only species comparison for which yeast cell size could not be used as a differentiating trait (Tables 2 and 3). Interestingly, B18, the isolate that had conflicting phylogenetic positions and has been previously associated with S1, had a yeast cell size akin to S1 and different from all the other species (Table 2). In terms of the number of buds produced, B18 is similar to PS2 and PS3 yet different from S1 and PS4 (Table 3). Thus, morphological features produce conflicting assignments for B18. The magnitude of the phenotypic differences is subtle. For that reason, we studied to what extent differences in yeast cell area and the number of yeast cell buds had predictive power in species assignment. We collected data for 20 cells from each the five species, and individuals were assigned to the species with which they had the lowest Mahalanobis distance (Table S4). We found that individuals from all species can be assigned to the correct species with a confidence of over 90%. Individuals from P. lutzii could be assigned correctly 100% of the time. Our results indicate that even though the phenotypic characterization is currently modest (~2 lines per species), the combination of these two yeast traits is an effective tool for discriminating between Paracoccidioides species.

The species within the brasiliensis complex are reciprocally monophyletic

We generated rooted phylogenetic trees for Paracoccidioides species using available sequence data. Bayesian and maximum likelihood reconstructions gave similar results so we only discuss the former. The results for the Bayesian analyses are presented in Figure 3 and show two different patterns. First, we found that the four phylogenetic species from the brasiliensis species complex are reciprocally monophyletic at nuclear loci (coding, Figure 3A and Figure S3; noncoding, Figure 3B and Figure S4), which is consistent with the observation that when these genealogies are unrooted they effectively differentiate the four genetically isolated groups (S1: green, PS2: red; PS3: blue, and PS4: yellow). The only isolate that showed conflicting positions between the two types of nuclear loci (coding vs. noncoding) was B18, an isolate from Brazil (Restrepo-Moreno & Schneidau 1967; Robledo et al. 1968). Notably, B18 appears as the outgroup of the dyad PS3/PS4 in the tree built with nuclear coding loci (Figure 3). In the phylogeny built with noncoding nuclear loci, B18 is part of a polytomy with PS3, PS4, and S1. These results conflict with the original reports from Matute et al. (Matute et al. 2006a), Teixeira et al. (Teixeira et al. 2009), and Teixeira et al. (Teixeira et al. 2014) where the same isolate was classified as S1. Thus, our sampled loci cannot unequivocally determine the identity of B18 but it seems possible it is the only isolate with recent hybrid origin. The mitochondrial DNA genealogical pattern differed from that for nuclear loci. First, two of the four species, S1 and PS3, were not monophyletic (Figure 3C and Figure S5). Most notably, the reciprocal relationships of the four species differed from those obtained from nuclear markers. Nuclear markers showed that PS3 and PS4 are sister species, S1 is the sister of the dyad, and PS2 is the earliest branching species in the species complex. In the mitochondrial phylogeny, PS2 and PS4 are sister groups which is clearly inconsistent with the nuclear phylogenies. Since the geographic ranges of PS2 and PS4 overlap in Venezuela, the mitochondrial clustering seems to be dictated by geography rather than by species identity. This pattern suggests gene exchange between PS2 and PS4. Notably, both PS2 and PS4 are also reciprocally monophyletic in the mitochondrial topologies and are not interspersed among each other. Additionally, two S1 isolates (A3 and U1) from southern South America also cluster together and are the outgroup of the Venezuelan dyad of PS4 and PS2. These results make S1 a paraphyletic group unlike in the nuclear phylogenies. Similarly, six isolates from S1 (henceforth referred to as the B21 mitochondrial haplotype) appear to be closely related to the vast majority of PS3 isolates (20 out of

21), even though nuclear alleles reveal no such close relationship (Figure 3A and 3B). The S1 isolates from the B21 mitochondrial haplotype come from Brazil, Paraguay, Argentina and Peru and form a monophyletic sister group to PS3. Other inconsistencies exist between the mitochondrial and nuclear phylogenies. In the mitochondrial tree, C10, a PS3 isolate is nested within S1 and has the same mitochondrial haplotype as eight Brazilian isolates. However, in the nuclear trees, C10 is found with the other PS3 isolates. Consistent with the nuclear results, B18 is difficult to place in the tree and is nested within PS3.

Inconsistency between nuclear and mitochondrial genomes indicates shared ancestry among recently diverged species of Paracoccidioides

We next quantified the level of inconsistency between nuclear and mitochondrial markers using AU tests for genealogical concordance. The AU test provides information about the set of topologies that are plausible at a certain confidence level. We tested four possible topologies (Figure 1). We excluded B18 from this analysis since its position might confound species relationships. Nuclear markers overwhelmingly support PS3 and PS4 as sister species, S1 as the sister group of the PS3/PS4 dyad, and PS2 as the most early-diverging species of the Paracoccidioides species complex. For mitochondrial markers, the best supported topology had PS2 and PS4 as sister clades (supported by 13 out of 15 genealogies; Table 5). Similarly, PS3 and S1 appear as sister species (supported by 12 out of 15 genealogies; Table 5). These results are consistent with the results from the phylogenetic reconstructions and indicate that the evolutionary histories of nuclear and mitochondrial markers differ. Specifically, they show that there is shared ancestry between PS4 and PS2 and between S1 and PS3 for the mitochondrial genome but nor for the nuclear genome. We also used AU tests to explicitly test the possibility of introgression between S1 and P. lutzii, or PS2 and P. lutzii (as proposed by Theodoro et al. 2012). We evaluated three possible topologies (Figure S1). We find that the topologies with P. lutzii as sister species of either S1 or PS2 shows no support in any of the iterations in either mitochondrial or nuclear tests (Table 5). These results indicate there is no support for the hypothesis that either pair (P. lutzii and S1, or P. lutzii and PS2) has any shared ancestry, thus contesting the proposal they have produced hybrids in the recent past.

Population genetic analyses

We employed a population genetics framework to quantitatively assess if the patterns of molecular variation were consistent with mitochondrial introgression. First, we calculated levels of within-species population variation. The levels of genetic variability (nucleotide diversity, π) of these species fall within the range of genetic diversity at nuclear loci in other species of fungi (Table 6, Lower boundary = 0.000843 (Cryptococcus neoformans; (Simwami et al. 2011; Leffler et al. 2012)), Upper boundary = 0.0046 (Saccharomyces paradoxus; (Elyashiv et al. 2010; Leffler et al. 2012)). Assuming a constant mutation rate in the Paracoccidioides genus, our results indicate that the four species of Paracoccidioides have similar effective population sizes (Kruskal-Wallis on π of nuclear loci only: 2= 7.184, df = 3, P = 0.066). Two pairwise comparisons are almost significant: PS2 vs. PS4 (observed difference = 15.732) and S1 vs. PS4 (observed difference = 14.420). The critical Kruskal-Wallis difference for significance was 17.419 and the P-value associated for both values was larger than 0.05 in both cases. Second, we calculated levels of divergence between the six possible species pairs at both nuclear and mitochondrial loci. The results for these estimates are shown in Table 4. In the absence of mitochondrial gene flow, mitochondrial loci should be on average more divergent than nuclear DNA given observations from other organisms that mitochondrial mutation rates are larger (Bruns & Szaro 1992). The ratio of mitochondrial to nuclear divergence should be equivalent across comparisons if no introgression has occurred. If there has been no mitochondrial introgression, one would expect the average genetic distance in mitochondrial loci to be consistently proportional to that of nuclear loci in all pairwise interspecies comparisons. However, this is not the case. While the ratio of mitochondrial to nuclear divergence in PS3 and PS4 is close to 8, the ratio of PS2 and S1 and of PS2 and PS4 is lower than 2 (Table 7). In general, mitochondrial DNA does not consistently reflect the levels of divergence seen for nuclear markers (range of mitochondrial divergence: [0.01211-0.03094]; range of nuclear divergence: [0.00306-0.009445]. These calculations however assume constant mutation rates in the Paracoccidioides genus. Finally, we assessed molecular divergence in a comparative framework. Since PS3 and PS4 are sister species, one would expect PS2, an outgroup, to be equally differentiated from both species. Nuclear differentiation in both species pairs (PS3 vs.

PS2 and PS4 vs. PS2) is close to 1%. That is not the case for the mitochondrial genome. While PS3 and PS2 have a mitochondrial differentiation of 3.1%, PS4 and PS2 have a differentiation of 1.2% (Table 7). Notably, PS4 and PS2 also have incongruent phylogenetic trees (see above) and share a large portion of their geographic range. This genetic pattern is incompatible with retention of ancestral polymorphism pre-dating speciation because one would expect similar retentions of polymorphism in both PS3 and PS4. Coupled with the discordant nuclear and mitochondrial phylogenetic trees, the population genetic results strongly suggest horizontal gene exchange of mitochondria occurred between PS2 and PS4.

Paracoccidioides species show a continuum of speciation times

We calculated approximate dates for the speciation events in the Paracoccidioides genus using two different approaches. First, we assumed constant substitution rates in the Paracoccidioides genus. We relied on previous calculations that indicated that the divergence between P. lutzii and the brasiliensis complex occurred approximately 33 million years ago (mya). We find that species of Paracoccidioides show a large range of divergence (Table 8, row 1) from a fraction of a million years (PS3-PS4) to dozens of millions of years (P. lutzii-brasiliensis complex). This type of calculation also allows us to roughly quantify when mitochondrial exchanges might have occurred. Assuming that Paracoccidioides shows a mitochondrial substitution rate similar to other ascomycete fungi (Torriani et al. 2014), the differentiation in mitochondrial loci between PS2 and PS4 indicates that the mitochondrial homogenization occurred around 0.1mya, and indicates that coalescence of these mitochondrial haplotypes is a much more recent event than the nuclear divergence between PS2 and PS4 (~0.71 mya, Table 8, row 1). We also calculated divergence ages relaxing our molecular clock assumptions. Results for the estimated divergence for nuclear coding loci using multiple clock models are shown in Table 8 (rows 2-10). These results are consistent with the constant substitution rate calculations and further support that species within the Paracoccidioides genus show a large range of speciation times. Assessment of population structure

We used two additional approaches to further study genetic variation in Paracoccidioides. First, we did a Principal Component (PC) Analysis (PCA) to analyze allelic frequencies in microsatellite data. The first five principal components of the PCA explained about 50% of the genetic variance. The results suggested three genetic patterns (Figure 4). First, PC1 (variance explained = 15.6%) differentiates between three of the four species of the brasiliensis species complex: S1, PS3 and PS4. The differentiation is not complete as there is some overlap between isolates from S1 and PS4. Second, PC2 (variance explained = 9.3%, Figure 4) differentiates PS2 from the other three species. Third, PC3 explains variation primarily within PS2 and to a lesser extent within the three other species (8.3%). No obvious patterns explain PC4 (7.7%) or PC5 (7.3%). We also estimated the most likely number of clusters in the P. brasiliensis species complex using STRUCTURE. Allele frequencies best explained by five populations (Figure 5; K = 5, logL= -726.5; LRT (K=4 vs. K=5): LRT = 33.0; df = 1; P = 8.34 × 10-8). Population assignments for higher values of K are shown in Figure S6. Notably, in the most likely clustering scenario, each of the species of Paracoccidioides represents one of the identified genetic clusters confirming the results from the nuclear analyses and the PCA. The clustering, however, does not precisely adhere to the species boundaries. S1 seems to be substructured and has two well-defined populations corresponding to geographical locations (S1_Brazil and S1_southern South America). These two subpopulations have disconcordant nuclear genealogies (e.g., B22 and A8 have different locations in the nuclear coding and nuclear noncoding trees, Figure 2) and do not satisfy the requirements to be considered two different species. Nonetheless, we treat them as two differentiated populations within S1. In terms of clustering, most individuals from PS3, PS4, and S1 are largely straightforward to include in their respective groups. Some exceptions (Bayesian probability < 60%) do exist (e.g., B18). The exceptions can be explained by gene flow, ancestral polymorphism, or simply marker convergence. Our STRUCTURE results show that the observed allele frequencies are best explained by four species of Paracoccidioides brasiliensis one of which shows considerable population structure. We also find uncertainty in the clustering of individuals from all species. Given these results, we formally explored the possibility of gene flow between these species. Paracoccidioides species exchange genes infrequently

Since we detected mitochondrial introgression, and clustering algorithms were unclear about the identity of some isolates, we evaluated the possibility of gene flow at nuclear loci. Since the nuclear regions we evaluated show little intraspecific variance, we used microsatellites, a rapidly evolving class of loci. We assessed the extent to which S1, PS2, PS3, and PS4 exchange genes in natural populations. We ran Treemix with 0, 1, 2, and 3 migration events. Two migration edges were the most likely migration scenario when assessed by both LRT (one migration edge vs. two migration edges: LRT= -10.799; P = 8.003 × 10-13) and wAIC (wi = 0.682; Figure 6). Adding an additional migration event did not improve the log likelihood (two migrations edges vs. three migration edges: LRT= -10.522; P =0.757). wAIC also indicated that two migrations were the most likely demographic scenario. The results for the most-likely topology are shown in Figure 8 and reveal that there have been two independent migrations from PS2 to S1. The drift parameter, which is proportional to the age of the migration, revealed that the genetic exchange was unlikely to be recent. These results strongly suggest that the four species of Paracoccidioides from the brasiliensis complex are differentiated and in some instances have had the chance to interbreed. DISCUSSION

In this report, we show that Paracoccidioides is composed of several species that differ genetically and morphologically. More importantly, even though the Paracoccidioides species coexist in their geographic ranges (and in some cases within the same host) and thus have had the chance to interbreed, the level of gene flow between them is low and is mostlyyet not exclusivelyrestricted to the mitochondria. Our results indicate that despite low levels of gene flow, the five genetic groups of Paracoccidioides (P. lutzii and the four species in the P. brasiliensis complex) are indeed well separated species. Our results expand on previous descriptions of phylogenetic species within Paracoccidioides and address the genealogical relationships between them. Our phylogenetic and population genetics results reveal that nuclear markers, both coding and noncoding, are consistent with a speciation scenario in which PS3 and PS4 are sister species which are, in turn, closely related to S1. PS2 is the earliest diverged species within the Paracoccidioides brasiliensis species complex and is an outgroup of

the PS3/PS4/S1 triad. S1 was considered a paraphyletic clade (i.e., it did not include all of the descendants of a common ancestor; (Matute et al. 2006a)) and could not be classified as a species. Our study reveals that S1 is a monophyletic group and thus satisfies the requirements of a phylogenetic species. Our results also resolve a disagreement regarding how many species exist in Paracoccidioides. Nuclear genealogies identified five phylogenetic species in the genus (Matute et al. 2006a; Teixeira et al. 2009; Theodoro et al. 2012; Teixeira et al. 2014), while mitochondrial genealogies identified just two (Salgado-Salazar et al. 2010). The main source of controversy is the number of species in the brasiliensis species complex: while nuclear gene genealogies suggested the existence of four, mitochondrial genealogies suggested only one. We jointly analyzed nuclear and mitochondrial markers and found that they differ in two ways. First, there have been at least three independent mitochondrial gene exchange events between Paracoccidioides species (PS3 and S1, PS2 and PS4, and possibly S1 and PS2/PS4). The most pronounced gene exchange event occurred between PS2 and PS4 which appear as sister clades in mitochondrial phylogenies. This relationship is not consistent with the one observed with nuclear markers and might be the result of secondary contact and hybridization as both PS2 and PS4 coexist in Venezuela. Moreover, PS2 and PS4 are not interspersed with each other and instead form reciprocally monophyletic groups in mitochondrial phylogenies, indicating that the admixture event was ancient. Mutations have had time to accumulate resulting in reciprocal monophyly. A similar pattern is observed between PS3 and S1. The B21 mitochondrial type from S1 is most similar to the predominant mitochondrial haplotype in PS3. Yet, the two species do not appear intermixed and instead form two monophyletic groups which suggests an old admixture event. A third, and perhaps even older, hybridization event occurred between the ancestor of the PS2/PS4 dyad and two S1 isolates from southern South America (S1 and A1). A fourth hybridization event characterized by mitochondrial replacement may have occurred in the PS3 isolate, C10. Even though it is nested within PS3 in phylogenies built from both sets of nuclear loci, its mitochondrial haplotype is identical to eight isolates from S1. Since the C10 haplotype is identical to one observed within S1, the hybridization event must have happened recently enough that the C10 mitochondria has not accumulated any new mutations. It is possible that the identical haplotypes across species are not explained by introgression but by alternative hypotheses. Our phylogenetic (tree reconstructions, AU

tests) and population genetic (molecular divergence) analyses suggest gene-tree discordance between mitochondrial and nuclear loci, but they do not differentiate between incomplete lineage sorting (ILS) and hybridization. The discordance between nuclear and mitochondrial topologies, however, is more likely to be caused by introgression than retention of ancestral polymorphism. First, given the expected levels of divergence between species and the reciprocal monophyly at nuclear loci, ancestral variation is not expected to still be segregating in either the mitochondrial or nuclear genomes (following Hudson & Coyne 2009 and Gao et al. 2015). Second, the clustering of mitochondrial haplotypes is dictated by geography. While there is no reason to expect geographic clustering with an ILS scenario, shared ancestry confined to geographic regions is precisely the expectation of introgression. Our tests and results based on consensus sequences have two limitations. First, they do not address conflicting gene genealogies within species. Second, levels of gene flow will be underestimated since only introgressions present at high frequencies in the recipient species will be detected. A third possibility is that the similarity at mitochondrial genomes is not caused by shared ancestry (either ILS or introgression) but instead is caused by parallelism (i.e., convergent evolution caused by the same genetic changes, sensu Rosenblum et al. 2014). More powerful tests such as genome-wide phylogenetic trees based on individual sequences might be more suitable for addressing intraspecific genetic variation and measuring interspecific gene flow. Muñoz et al. (2016) have generated a population level genome wide dataset that should be useful to answer this question. The prevalence of introgression and the porosity of species boundaries in Paracoccidioides will be addressed in future contributions. More generally, our observations indicate that mitochondrial markers are poorly suited for resolving phylogenetic relationships in the Paracoccidioides genus. These limitations have also been observed in other organisms (Moore 1995; Nichols 2001). Mitochondrial DNA evolves faster than nuclear DNA (Elena et al. 1998; Song et al. 2005), but mitochondrial gene genealogies in Paracoccidioides do not reflect this expected greater divergence. This result seems to be inconsistent with a simple allopatric divergence between species as the coalescent time is expected to be longer for mitochondrial than nuclear loci (Elena et al. 1998). Instead Paracoccidioides species appear to have exchanged genes in the recent past resulting in homogenized mitochondrial genotypes. Notably, nuclear microsatellite markers also indicate that gene flow might have occurred but from PS2 into S1, a different species pair than that seen

from the mitochondrial alleles. Overall, the combination of these results strongly suggests that Paracoccidioides species have exchanged genes in the past either by sexual or somatic processes (i.e., regular anastomosis with exchange of mitochondria but no nuclei; Sanders 2006; e.g., Yoo et al. 1987, Madhosingh 1992, Roca et al. 2004). Analyses of nuclear microsatellites also revealed gene flow between Paracoccidioides species. Genetic material from PS2 can be found in both subpopulations of S1 but not in the other two Paracoccidioides species. This result is unexpected because we hypothesized that nuclear and mitochondrial introgressions should co-occur. However, our results do not support such a scenario. There are two potential explanations for the discordance between nuclear and mitochondrial phylogenies. First, nuclear and mitochondrial introgressions might have occurred at different times. This is likely to be the case as the mitochondrial introgressions (PS4 and PS2, PS3 and S1, PS2/PS4 and S1) seem to be ancient. Given the high mutation rate of microsatellites, they might be better suited for understanding relatively recent instances of gene flow. The second possibility is that the introgressions happened at the same time but that selection purged mitochondrial and nuclear allelic combinations differently in different species pairs. This would be likely if there was postzygotic isolation between Paracoccidioides species and if the genetic basis of hybrid breakdown differs among different species pairs. Mitonuclear interactions may be involved between some species pairs but not between others. Both hypotheses are speculative and should be systematically assessed. Our results also show that even though interbreeding and genetic exchange have occurred between Paracoccidioides species, it is not common enough to merge the species into a single genetic entity. Given extensive sympatry between Paracococcidioides species, observed genetic differentiation argues for the existence of barriers to gene flow. Nonetheless, if Paracoccidioides is strictly clonal, then once lineages have diverged they will not interbreed even when fully sympatric. Three lines of evidence argue against this possibility. First, previous studies (Matute et al. 2006a; Teixeira et al. 2009) have found a molecular signal of recombination. It remains unknown whether recombination is meiotic, mitotic, or somatic/parasexual. These possibilities are not mutually exclusive either (Engelstaedter 2017). Regardless, Paracoccidioides does not seem to be strictly asexual. Second, Paracoccidioides lines have intact mating type loci indicating that they are under the influence of purifying selection and may be functional (Torres et al. 2010; Teixeira et al. 2013). Similarly,

Torres et al. (Torres et al. 2010) showed it is possible to produce sexual structures (perithecia) in crosses between Paracoccidioides isolates. None of the perithecia were fertile, and they carried no ascospores. However, only interspecific crosses produced progeny (PS3 × S1, PS3 × P. lutzii). It remains unclear whether the incomplete development of perithecia represents truncated sexual development in Paracoccidioides, or whether it arises from hybrid sterility. Finally, mitochondrial haplotypes, and to a lesser extent nuclear loci, have been exchanged between Paracoccidioides species, indicative of historic interbreeding. We also show that admixture in the Paracoccidioides genus is not restricted only to ancient instances but might still be ongoing. The isolate B18 shows a mixture of molecular and morphological traits from S1 and PS3 which suggests, but does not confirm, that it might be a hybrid. Notably, this makes studying hybridization in the Paracoccidioides genus possible. Furthermore, since some, but not all, Paracoccidioides species are resistant to antifungals, gene flow between species might have important medical implications (Johann et al. 2010).

Species description

We show that even though the four species of the Paracoccidioides brasiliensis species complex have had the chance to interbreed which has ultimately led to mitochondrial replacements, they are highly differentiated at nuclear loci. Since the early days of speciation research, genetic groups that can coexist in the same geographic area and persist as largely (if not completely) isolated gene pools have been considered to be ‘good species’ (Mayr 1963; Coyne & Orr 2004). PS4 and PS2 overlap in Venezuela, and the fact that there is no detectable gene flow at the nuclear level in spite of the mitochondrial similarity indicates that there has been no recent gene flow between these two species. Gene flow would be possible since they have been collected from the same locations. The same is true for S1 and PS3. While PS3 was initially only thought to be restricted to Colombia, the species has also been found in Venezuela (Roberto et al. 2016) and possibly Brazil (Gonzalez & Hernandez 2016). All five Paracoccidioides species (P. lutzii and the four species of the brasiliensis complex) also seem to differ morphologically. Yeast morphological traits are largely, but not completely, predictive of the species identity. Taken together, our results indicate that the phylogenetic species of Paracoccidioides should be elevated to taxonomical species status. We thus propose to

conventionally adopt the names P. americana for PS2, P. restrepiensis for PS3, and P. venezuelensis for PS4. We also propose to restrict the use of P. brasiliensis for S1.

Paracoccidioides americana sp. nov. (Matute and McEwen 2016)

Mycobank MB816860. Formerly known as PS2. Mycelial colonies on McVeigh and Morton are white to cream and cottony with a radius of under 10 mm after 20 days of incubation at 25°C. Yeast colonies are variable in size, cream-colored (white to brown) on Sabouraud’s agar. Hyaline conidia are usually thick-walled arthroconidia, but thinwalled empty disjunctor cells, single cell conidia, and terminal/ lateral aleuroconidia are also present. Yeast size (i.e., area) differentiates P. americana sp. nov. from the other species from the brasiliensis complex (median = 48.7695; SD = 9.7437). Median number of buds per yeast is similar to that of other species in the complex and it is not a diagnostic trait (median = 2; SD = 1.6800). No sexual stage is known. Paracoccidioides americana can be discriminated from other cryptic species of the P. brasiliensis complex by PCR amplification of GP43: forward – 5’ CCAGGAGGCGTGCAGGTGTCCC 3’ and reverse – 5’ GCCCCCTCCGTCTTCCATGTCC 3’ or α-Tubulin 5’ CTTAACGGTGCCTTCTTTGCGG 3’ and reverse – 5’ CTTAACGGTGCCTTCTTTGC 3’. Other diagnostic SNPS are listed in (Matute et al. 2006a). Additionally, P. americana can be identified by microsatellite typing (Matute et al. 2006b), RFLP patterns (Nino-Vega et al. 2000; Roberto et al. 2016), and genomic divergence (Desjardins et al. 2011; Muñoz et al. 2014). HOLOTYPE: Isolated in Botucatu, São Paulo State, Brazil from a chronic PCM patient. Ex-type culture is preserved in the Corporación para Investigaciones Biologicas (CIB_Pb03, = B26). Etymology: americana; refers to the distribution of the species in the South American continent (Matute et al. 2006a).

Paracoccidioides restrepiensis sp. nov. (Matute and McEwen 2016)

Mycobank MB816861. Formerly known as PS3. Mycelial colonies on McVeigh and Morton are white to cream and cottony with a radius of under 10 mm after 20 days of incubation at 25°C. Yeast colonies are variable in size, cream-colored (white to brown) on Sabouraud’s agar. Hyaline conidia are usually thick-walled arthroconidia, but thin-

walled empty disjunctor cells, single cell conidia, and terminal/ lateral aleuroconidia are also present. Large variance in the production of conidia in Yeast-Extract, Glucose salts and Water agar (Bustamante-Simon et al. 1985). Yeast size (i.e., area) differentiates P. restrepiensis sp. nov. from the other species from the brasiliensis complex (median = 91.1005; SD = 34.1232), with the exception of P. venezuelensis. The median number of buds per yeast is similar to that of other species in the complex but higher than the median number in P. venezuelensis (median = 2; SD = 1.8565). No sexual stage is known. Paracoccidioides restrepiensis can be discriminated from other cryptic species of the P. brasiliensis complex by PCR amplification of GP43: forward – 5’ CCAGGAGGCGTGCAGGTGTCCC 3’ and reverse – 5’ GCCCCCTCCGTCTTCCATGTCC 3’ or α-Tubulin 5’ CTTAACGGTGCCTTCTTTGCGG 3’ and reverse – 5’ GTGAAAGTATTGTTGCCCAGCG 3’. Other diagnostic SNPS are listed in (Matute et al. 2006a). Additionally, P. restrepiensis can be identified by microsatellite typing (Matute et al. 2006b), RFLP patterns (Nino-Vega et al. 2000; Roberto et al. 2016), or genomic divergence (Desjardins et al. 2011; Muñoz et al. 2014). HOLOTYPE: Isolated in the Department (State) of Antioquia from a chronic PCM patient in Colombia. Ex-type culture is preserved in the American Type Culture Collection (ATCC 60855). Etymology: restrepiensis; a tribute to Dr. Angela Restrepo who has been working with P. brasiliensis for over 50 years. She has published over 250 papers on the fungus and has trained a generation of mycologists in the field. Paracoccidioides venezuelensis sp. nov. (Matute and McEwen 2016)

Mycobank MB816862. Formerly known as PS4. Mycelial colonies on McVeigh and Morton are white to cream and cottony with a radius of under 10 mm after 20 days of incubation at 25°C. Yeast colonies are variable in size, cream-colored (white to brown) on Sabouraud’s agar. Yeast size (i.e., area) differentiates P. venezuelensis sp. nov. from the other species from the brasiliensis complex (median = 92.2886; SD = 33.3980) with the exception of P. restrepiensis. Median number of buds per yeast is lower than the median of other Paracoccidioides species and constitutes a diagnostic trait (median = 1; SD = 0.5793). No sexual stage is known. Yeast colonies are variable in size, creamcolored (white to brown) on Sabouraud’s agar. Hyaline conidia are usually thick-walled arthroconidia, but thin-walled empty disjunctor cells, single cell conidia, and terminal/

lateral aleuroconidia are also present. Large variance in the production of conidia in Yeast-Extract, Glucose salts and Water agar. Paracoccidioides venezuelensis can be discriminated from other cryptic species of the P. brasiliensis complex by PCR amplification of GP43: forward – 5’ CCAGGAGGCGTGCAGGTGTCCC 3’ and reverse – 5’ GCCCCCTCCGTCTTCCATGTCC 3’ (Matute et al. 2006a), microsatellite typing (Matute et al. 2006b), RFLP patterns (Nino-Vega et al. 2000). HOLOTYPE: Isolated in the State of Miranda, Venezuela from a sample of soil by M. Albornoz (De Albornoz 1971). Ex-type culture is preserved in the Corporación para Investigaciones Biologicas (Pb300). Etymology: venezuelensis; refers to the distribution of the species in Venezuela (Teixeira et al. 2014).

Paracoccidioides brasiliensis sensu stricto

Mycelial colonies on McVeigh and Morton are white to cream and cottony with a radius of under 10 mm after 20 days of incubation at 25°C. Yeast colonies are variable in size, cream-colored (white to brown) on Sabouraud’s agar. Hyaline conidia are usually thickwalled arthroconidia, but thin-walled empty disjunctor cells, single cell conidia, and terminal/ lateral aleuroconidia are also present. Yeast size (i.e., area) differentiates P. brasiliensis sensu stricto from the other species from the brasiliensis complex (median = 68.0023; SD = 37.7582). Median number of buds per yeast is similar to that of other species in the complex and it is not a diagnostic trait (median = 3; SD = 2.27). No sexual stage is known. Paracoccidioides brasiliensis sensu stricto can be discriminated from other cryptic species of the P. brasiliensis complex by PCR amplification of GP43: forward – 5’ CCAGGAGGCGTGCAGGTGTCCC 3’ and reverse – 5’ GCCCCCTCCGTCTTCCATGTCC 3’ or α-Tubulin 5’ CTTAACGGTGCCTTCTTTGCGG 3’ and reverse – 5’ CTTAACGGTGCCTTCTTTGC 3’. Other diagnostic SNPS are listed in (Matute et al. 2006a). Additionally, P. brasiliensis can be identified by microsatellite typing (Matute et al. 2006b), RFLP patterns (Nino-Vega et al. 2000; Roberto et al. 2016), and genomic divergence (Desjardins et al. 2011; Muñoz et al. 2014). HOLOTYPE: Isolated in Brazil from a chronic PCM patient in 1929. Isolate was donated by C. Fava Netto. Ex-type culture is preserved in the Corporación para Investigaciones Biologicas (Pb18 = B17).

Future directions

The study of speciation in Paracoccidioides may have important consequences for our understanding of the disease it causes. Paracoccidioidomycosis is an important disease in Latin America, and yet, very little is known regarding how the genetic background of the pathogen affects clinical symptoms. This is especially important because the ecological niche and specific mode of infection has remained elusive despite coordinated efforts to find natural reservoirs and ecological conditions in which Paracoccidioides might be found in nature (Restrepo et al. 2001; Bagagli et al. 2008; Arantes et al. 2016). Additionally, Paracoccidioides isolates differ in their ability to infect mice (Molinari-Madlum et al. 1999). More generally, a full characterization of the morphological, physiological, and epidemiological characteristics of each of the Paracoccidioides species will reveal to what extent the morphological differences reported here are affected by phenotypic plasticity. Whether, this variation corresponds to the different Paracoccidioides species remains unexplored. Virulence factors have also been shown to differ between species (Matute et al. 2008; Tamayo et al. 2016) which is a provocative observation but deserves formal testing. Our study is not the first one to propose the existence of different species within the genus Paracoccidioides. As early as 1935 a second species of Paracoccidioides was reported: P. cerebriformis (Moore 1935, Artagaveytia-Allende 1949). This isolate clearly differed in colony morphology which led to its recognition as a potentially distinct group. However, molecular genotyping revealed that P. cerebriformis actually belongs to the genus Aspergillus. This of course exemplifies the dangers of coining taxonomical binomial names without enough taxonomic evidence but also the self-correcting process of taxonomy, as the binomial has not been used since 2005 (Cavalcanti et al. 2005). The latest example of a newly coined binomial in Paracoccidioides is P. lutzii (Teixeira et al. 2014). Paracoccidioides lutzii differs from P. cerebriformis because it was diagnosed by genetic differentiation and not merely morphology. Our argument is similar to the one proposed for P. lutzii but for species that are more closely related. Experiments involving crosses and studying the magnitude of genome differentiation and whether such lines of evidence corroborate our hypotheses are of primary need. Our study opens the door for the study of speciation and divergence in Paracoccidioides in two ways. First, we formally recognize the taxonomic stature of three new species of Paracoccidioides encouraging the study of potential interspecific

differences in morphology, physiology, and virulence across the whole genus. Currently, all of these traits have only been examined in a handful of isolates. A full phenotypic characterization of each species is sorely needed. Such efforts are currently under way to assess the magnitude of phenotypic differentiation between P. lutzii and other species of Paracoccidioides (e.g., Medeiros et al. 2016) Second, we demonstrate that in spite of their divergence, species of Paracoccidioides have exchanged genetic material such that mitochondrial divergence, which is expected to be greater than nuclear divergence, has been reduced. Future studies will reveal the extent of introgression in Paracoccidioides and establish the role of geography in such genetic exchange.

DATA AVAILABILITY

All DNA sequences were deposited in GenBank under the accession numbers KX719990-KX720560. Additional data and analytical code have been deposited in Dryad (http://dx.doi.org/10.5061/dryad.j1f83). ACKNOWLEDGEMENTS

We would like to thank Angela Restrepo, Alejandro Rojas, and Silvia Restrepo for scientific discussions and comments at every stage of the manuscript.

TABLES TABLE 1. Pairwise comparisons of the surface area per conidium between three different Paracoccidioides species. All comparisons were done with a Kruskal-Wallis test and were adjusted for multiple comparisons. Numbers on the upper triangular matrix are the observed differences between species. N Lines refers to the number of scored lines. NCells refers to the total number of scored cells. Numbers on the lower triangular matrix are the expected differences required for significance for each species pair. Significant comparisons are marked in bold.

Kruskal-Wallis pairwise comparisons Species

NLines

NCells

Mean

SD (μm2)

P. lutzii

P. brasiliensis

P. americana

2

(μm ) P. lutzii

1

167

10.2899

5.1548

*

160.3007

166.9146

P. brasiliensis (Pb18)

1

166

4.7970

1.3956

33.3520

*

6.6139

P. americana (Pb03)

1

112

4.6345

1.5702

38.4580

38.1600

*

TABLE 2. Pairwise comparisons of the area per yeast cell between the five Paracoccidioides species. All comparisons were done with a Tukey Honest different test with multiple comparison correction. Numbers on the upper triangular matrix are t value of the pairwise species comparisons. NLines refers to the number of scored lines. NCells refers to the total number of scored cells. Numbers on the lower triangular matrix are the associated p-values for the t-test (Pr>|t|). Significant comparisons are marked in bold. Posthoc pairwise comparisons (Tukey test)

P. lutzii

NLin

NCel

Median

SD

es

ls

(μm2)

(μm2)

2

138

93.2275

49.899

P. lutzii

P.

B18

P.

P.

P.

brasilien

american

restrepie venezuel

sis

a

nsis

ensis

*

5.074

0.0689

-7.987

-2.929

-0.653

<0.001

*

-0.663

2.855

-1.927

-3.362

2.726

0.9850

*

-2.796

0.735

1.947

<0.001

0.0488

0.0574

*

4.721

5.832

0.0396

0.3764

0.9762

<0.001

*

1.666

0.9860

0.0104

0.3646

<0.001

0.5444

0 P. brasiliensis (S1)

2

62

68.0023

37.758 2

B18

2

23

82.0335

37.328 8

P. americana (PS2)

2

50

52.4522

19.394 2

P. restrepiensis (PS3)

2

67

79.4851

32.751 6

P. venezuelensis (PS4)

2

44

92.2886

33.398 1

*

TABLE 3. Pairwise comparisons of the mean number of buds per yeast between the five Paracoccidioides species. All comparisons were done with a Tukey Honest Significant Different test with multiple comparison correction. Numbers on the upper triangular matrix are t value of the pairwise species comparisons. Numbers on the lower triangular matrix are the associated p-values for the t-test (Pr>|t|). Significant comparisons are marked in bold.

`

Posthoc pairwise comparisons (Tukey test) NLin

NCel

Median

SD

es

ls

(cells)

(cells)

P. lutzii

P. brasiliensis

B18

P. americana

P. restrepiensis

P. venezue lensis

P. lutzii

2

83

1

2.51701

*

1.399

-2.533

1.371

0.610

-3.519

P. brasiliensis (S1)

2

62

3

2.2646

0.7221

*

-1.483

-0.056

0.764

4.521

B18

2

23

3

1.9681

0.1143

0.6691

*

-1.395

-2.056

-4.871

P. americana (PS2)

2

50

2

3.8704

0.7389

1.0000

0.7244

*

-0.777

-4.363

P. restrepiensis

2

67

3

1.8772

0.9899

0.9724

0.3073

0.9703

*

-3.899

2

44

1

0.5793

0.0063

< 0.001

< 0.001

< 0.001

0.00158

*

(PS3) P. venezuelensis (PS4)

TABLE 4. Molecular polymorphism at five mitochondrial and sixteen nuclear genes across the Paracoccidioides genus. The genomic location of each allele was determined with respect to the genome of Pb03. N is the number of genotypes isolates. S represents the number of segregating sites across the genus. The majority of genes show a negative value for Tajima’s D (Tajima 1989) indicating an excess of low frequency polymorphisms relative to expectation

locus

Location

N

S

Tajima’s D

cob2

mitochondrial

52

22

2.1494

atp6

mitochondrial

54

13

-0.4374

cox3

mitochondrial

54

14

-1.0145

rns

mitochondrial

55

12

0.6929

rnl

mitochondrial

56

72

-1.3741

gp43

Nuclear, coding: supercont_11:757587-758105

65

32

-1.0032

prp8

Nuclear, coding: supercont_2:212295-214013

22

9

1.2853

chs2

Nuclear, coding: supercont_25:32391-33004

65

12

-1.1740

arf

Nuclear, coding: supercont_1:1477749-1478155

65

6

-0.5177

cdc42

Nuclear, coding: supercont_7:474536-475621

17

11

-2.0875

fks

Nuclear, coding: supercont_6:1419296-1419889

65

5

-1.5883

orn1

Nuclear, non-coding: supercont_1:2156966-2157496

64

29

-1.2627

11b12b

Nuclear, non-coding: supercont_11:512815-512941

63

14

-0.9517

15b16b

Nuclear, non-coding: supercont_3:1487381-1487740

63

19

-1.7571

AB

Nuclear, non-coding: supercont_18:111672-112076

63

18

-1.8475

KL

Nuclear, non-coding: supercont_7:1050097-1050587

63

31

-2.1789

MN

Nuclear, non-coding: supercont_6:371103-371516

63

49

-2.5097

R56

Nuclear, non-coding: supercont_11:281070-281281

64

10

-1.0991

tub

Nuclear, non-coding: supercont_17:267063-267325

65

10

0.1193

III-IV

Nuclear, non-coding: supercont_9:1436424-1436795

64

41

-2.1239

XI-XII

Nuclear, non-coding: supercont_4:2836783-2837186

64

23

-1.7017

TABLE 5. Analyses of genealogical concordance between mitochondrial and nuclear loci. Comparison of four phylogenetic relationships among the four species of the Paracoccidioides brasiliensis species complex for mitochondrial (mtDNA) and nuclear regions. The four topologies (A-D) used to study introgression within the brasiliensis complex are shown in Figure 1. The three topologies to assess introgression between P. lutzii and species of the brasiliensis complex are shown in Figure S2.

brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group brasiliensis group P. lutzii-S1 P. lutzii-S1 P. lutzii-S1 P. lutzii-S1 P. lutzii-S1 P. lutzii-S1

Dataset

Tree comparison

|ΔL|

mtDNA mtDNA mtDNA mtDNA mtDNA mtDNA nuclear DNA nuclear DNA nuclear DNA nuclear DNA nuclear DNA nuclear DNA mtDNA mtDNA mtDNA nuclear DNA nuclear DNA nuclear DNA

Tree A-Tree B Tree A-Tree C Tree A-Tree D Tree B-Tree C Tree B-Tree D Tree C-Tree D Tree A-Tree B Tree A-Tree C Tree A-Tree D Tree B-Tree C Tree B-Tree D Tree C-Tree D Tree E-Tree F Tree E-Tree G Tree F-Tree G Tree E-Tree F Tree E-Tree G Tree F-Tree G

45.6 21.20 73.11 94.55 43.18 8.1 40.14 72.10 95.61 3.61 7.82 9.11 34.12 57.12 2.11 76.21 45.17 4.12

p value from Au test < 1 × 10-3 0.02 < 1 × 10-3 < 1 × 10-3 < 1 × 10-3 0.32 < 1 × 10-3 < 1 × 10-3 < 1 × 10-3 0.45 0.63 0.19 < 1 × 10-3 < 1 × 10-3 0.79 < 1 × 10-3 < 1 × 10-3 0.57

TABLE 6. Nucleotide diversity in the four species of the Paracoccidioides brasiliensis species complex. For coding genes, nucleotide diversity, was calculated based on the third position of the codon (πS). For noncoding genes, we used all sites (πNC). Both indices (πS and πNC) are a proxy of intraspecific diversity. The number of loci per category, the mean and the standard error (in parentheses) are shown. Numbers shown are true values, not percentages. Raw values are shown in Table S3.

Locus type

S1

PS2

PS3

mitochondrial

0.01881

0.00585

0.0071

coding

(0.0203)

(0.0062)

(0.0077)

mitochondrial

0.007005

0.0113

7.65 × 10-4

noncoding

(0.0052)

(0.0160)

(0.0011)

0 (0)

nuclear coding

0.0020

0.0021

0.0025

8.75 × 10-4

(0.0018)

(0.0017)

(0.0044)

(0.0018)

0.0045

0.00684

0.0047

0.0027

(0.0024)

(0.0056)

(0.0054)

(0.0045)

nuclear noncoding

PS4

0 (0)

TABLE 7. Genetic distance between the six possible pairs of species in the Paracoccidioides brasiliensis species complex. We calculated the median and standard deviation for nuclear (Mediannucl, SDnucl) and mitochondrial (Medianmitoc, SDmitoc) absolute divergence (Dxy).. Nnuc and Nmit are the number of nuclear and mitochondrial loci used for each calculation.

Nnuc Mediannucl

SDnucl

Ratio

Comparison

Nmit

Medianmitoc

SDmitoc

S1-PS2

6

0.01687

0.021179

5

0.00864

0.005249

1.9526

S1-PS3

6

0.01939

0.025823

5

0.00514

0.004201

3.7724

PS2-PS3

6

0.03094

0.035970

5

0.009445

0.040380

3.2758

PS2-PS4

6

0.01211

0.042972

5

0.009335

0.008185

1.2973

S1-PS4

6

0.02945

0.028531

5

0.00519

0.005708

5.6744

PS3-PS4

6

0.02451

0.019969

5

0.00306

0.007974

8.0100

Medianmitoc/Mediannucl

TABLE 8. Estimates of speciation time among branches of the Paracoccidioides genus. Row 1 shows divergence calculations in million years ago (mya) based on a constant substitution rate. Columns 2-10 show different models of molecular evolution implemented in BEAST. All branch lengths were normalized with respect to the divergence between the brasiliensis complex and P. lutzii (last column). CQP: Continuous quantile parametrization.

1 2 3 4 5 6 7 8 9 10

Substitution model

PS3-PS4 0.1270

(PS3, PS4)-S1 0.5000

(PS3, PS4),S1PS2 0.7140

brasiliensis complex-lutzii 33.000

Constant substitution rate (MYA) Strict clock Random local clock Fixed local clock Uncorrelated relaxed clock- Exponential Uncorrelated relaxed clock- Lognormal Uncorrelated relaxed clock- Gamma Uncorrelated relaxed clock- Exponential CQP Uncorrelated relaxed clock- Lognormal CQP Uncorrelated relaxed clock- Gamma CQP

0.0324 0.0200 0.0333 0.0331

0.1142 0.0736 0.1047 0.0990

0.1812 0.1172 0.1967 0.1656

1 1 1 1

0.0256

0.0997

0.1641

1

0.0382

0.0760

0.1354

1

0.0465

0.1176

0.1595

1

0.0576

0.1148

0.1942

1

0.0323

0.0965

0.1645

1

FIGURE LEGENDS

FIGURE 1. Four possible topologies used to quantify the level of discordance between nuclear and mitochondrial loci. Topology A shows the hypothesized species tree derived from nuclear loci. Trees B-D show P2 and PS4 as sisters, consistent with shared genetic variation dictated by geography. The topologies were used to perform the AU test and the results are shown in Table 5.

A P. lutzii

B PS2

S1

PS3

PS4

P. lutzii

PS3

S1

PS2

PS4

P. lutzii

C P. lutzii

PS3

S1

PS2

PS4

S1

PS3

PS2

PS4

D

FIGURE 2. Species of the Paracoccidioides brasiliensis complex differ in their mean values for three morphological traits. Each panel shows boxplots with the distribution of the three traits with the median marked by the thickest line and the 25% and 75% percentiles marked by the edges of the box. A. Conidial area. At least 20 individuals were measured per species. Data were taken from Teixeira et al. ((Teixeira et al. 2009)) which measured only three of the five Paracoccidioides species. B. Yeast cell area. We measured two lines for each species and at least 20 cells per line. C. Number of budding cells. We measured two lines for each species and at least 23 number of individuals per line. Boxes are color-coded using the same notation as Figure 1 with P. lutzii drawn as a white box and B18 shown as a grey box.

FIGURE 3. Rooted phylograms for the species of the Paracoccidioides brasiliensis complex using coding three types of molecular data. A. Nuclear concatenated coding loci. B. Noncoding nuclear concatenated loci. C. Mitochondrial concatenated loci. We built Bayesian-based topologies depicting the phylogenetic relationships among Paracoccidioides isolates. The topologies are rooted with P. lutzii (lut). The number above each branch indicates the support of a given branch (only branches with a posterior probability higher than 0.75) were retained. Terminal branches are color coded: PS3, blue; PS2: red; PS4: purple; S1: orange. The Bayesian posterior probability supporting each node is shown next to the relevant node and is color coded (see legend). Paracoccidioides lutzii (lut) and the isolate B18 are not color coded. Cladograms of each data type are shown in Figures S1-S3.

A.

C.

B.

p o s te r io r 1

0

S 1_B 21 S 1_P 2 S 1_A 1 S 1_P 1 S 1_A 6 S 1_P E 1 S 1_B 8 S 1_A 7 S 1_B 17 S 1_B 24 S 1_B 9 S 1_B 10 S 1_A 5 S 1_A 3 S 1_B 12 S 1_B 22 S 1_B 16 S 1_A 2 S 1_A 4 S 1_B 6 S 1_B 3 S 1_B 2 S 1_B 4 S 1_B 5 S 1_B 20 S 1_B 14 S 1_B 1 S 1_B 19 S 1_A 8 S 1_B 25 S 1_B 11 S 1_U 1 P S 4_V 6 P S 4_V 4 P S 4_V 3 P S 4_V 5 P S 4_V 1 P S 3_C 20 P S 3_C 8 P S 3_C 12 P S 3_C 15 P S 3_C 7 P S 3_C 5 P S 3_C 13 P S 3_C 10 P S 3_C 4 P S 3_C 16 P S 3_C 9 P S 3_C 17 P S 3_C 18 P S 3_C 6 P S 3_C 1 P S 3_C 11 P S 3_C 2 P S 3_C 21 P S 3_C 14 P S 3_C 19 P S 3_C 3 S 1_B 18 P S 2_B 23 P S 2_B 26 P S 2_V 2 P S 7_B 7 P S 2_B 15 P S 2_B 13 lut

p o s te r io r 100%

0 .0 1 %

P S 2_B 23 S1_PE10 S1_P2 S1_P1 S 1_B 21 S 1_A 6 S 1_A 1 S 1_U 1 S 1_B 11 S 1_B 25 S 1_A 4 S1_B6 S 1_A 7 S 1_A 3 S 1_B 17 S 1_B 24 S 1_B 10 S 1_B 22 S 1_B 16 S 1_A 5 S 1_B 12 S1_B9 S 1_A 2 S1_B8 S1_B4 S1_B5 S 1_B 20 S1_B2 S1_B1 S 1_A 8 S 1_B 19 S1_B3 S 1_B 14 P S 3_C 17 PS3_C5 PS3_C2 P S 3_C 11 PS3_C9 PS3_C1 P S 3_C 19 P S 3_C 18 P S 3_C 14 P S 3_C 13 P S 3_C 12 PS3_C3 P S 3_C 15 PS3_C4 P S 3_C 16 P S 3_C 21 P S 3_C 20 P S 3_C 10 PS3_C7 PS3_C8 PS3_C6 PS4_V3 PS4_V6 PS4_V4 PS4_V1 S 1_B 18 PS4_V5 P S 2_B 26 PS2_V2 P S 2_B 13 P S 2_B 15 PS2_B7 lut

p o s te r io r 100%

0 .0 2 %

P S 2_B 15 P S 2_B 26 P S 2_B 23 P S 2_B 7 P S 2_V 2 P S 4_V 1 P S 4_V 5 P S 4_V 4 P S 4_V 6 S 1_A 3 S 1_U 1 P S 3_C 1 P S 3_C 11 P S 3_C 4 P S 3_C 5 P S 3_C 16 P S 3_C 9 P S 3_C 12 P S 3_C 6 P S 3_C 13 P S 3_C 2 P S 3_C 7 P S 3_C 20 P S 3_C 18 P S 3_C 3 P S 3_C 21 P S 3_C 19 S 1_B 18 S 1_A 1 S 1_A 6 S 1_P 1 S 1_P E 1 S 1_B 21 P S 3_C 10 S 1_B 12 S 1_B 20 S 1_B 5 S 1_B 6 S 1_B 10 S 1_B 14 S 1_B 2 S 1_B 19 S 1_B 16 S 1_B 17 S 1_A 2 S 1_A 5 S 1_A 7 S 1_A 8 S 1_B 25 S 1_B 9 S 1_B 4 lutz

FIGURE 4. Principal component analysis for the Paracoccidioides brasiliensis species complex. The genetic variance in the genus Paracoccidioides as determined by microsatellite markers largely separates the four cryptic species in the genus. Individuals of S1 are shown in blue, PS2 in orange, PS3 in green, and PS4 in purple. 95% confidence ellipses are drawn for each species. A. Principal components (PC) 1 and 2 are shown and the percentage of variance explained by each eigenvalue is shown within parentheses on each axis. B. PC3 and PC4. No principal component explained more than 16% of the variance. Isolate B18 is marked with a green circle.

FIGURE 5. STRUCTURE analysis for the Paracoccidioides brasiliensis species complex. Each isolate is represented by a vertical bar. The length of the colored segments represents the probability of belonging to a given cluster. A. K = 2. B. K = 3. C. K = 4. D. K = 5. The logL value and the p-value of the LRT test is shown on the title of each panel. ΔK values overwhelming support the existence of five clusters corresponding to the four species of the brasiliensis species complex, with S1 having two well-structured populations. Results for other K assignments are shown in Figure S1. Isolate B18 is marked with an asterisk.

FIGURE 6. Treemix results show genetic differentiation between species from the Paracoccidioides brasiliensis species complex. A. Best fitting genealogy without admixture for Paracoccidioides populations calculated from the variance-covariance matrix of allele frequencies. at nine microsatellites. We estimated the likelihood of three demographic scenarios with 1 to 3 migration events (m). B. m = 1. C. m = 2. D. m = 3. To determine the best fitting demographic scenario we used LRTs to compare nested models (one migration vs. two migrations, two migrations vs. three migrations). We also used wAIC to find the best supported migration scenario (See text). The title of each panel contains the number of migrations, the likelihood, and the associated p-value when comparing that demographic scenario with the previous one. The best fitting demographic history involved two migration events from PS2 to the two subpopulations of S1 (panel B).

REFERENCES Abbott R, Albach D, Ansell S et al. (2013) Hybridization and speciation. Journal of evolutionary biology, 26, 229–246. Ajello L, Cheng S-L (1967) Sexual reproduction in Histoplasma capsulatum. Mycologia, 59, 689. Akaike H (1974) A new look at the statistical model identification. IEEE transactions on automatic control, 19, 716–723. Almeida AJ, Cunha C, Carmona JA et al. (2009) Cdc42p controls yeast-cell shape and virulence of Paracoccidioides brasiliensis. Fungal Genetics and Biology, 46, 919– 926. Almeida AJ, Matute DR, Carmona JA et al. (2007) Genome size and ploidy of Paracoccidioides brasiliensis reveals a haploid DNA content: Flow cytometry and GP43 sequence analysis. Fungal Genetics and Biology, 44, 25–31. Almeida FD (1930) Differences between the etiological agents of coccidioidal granuloma in the United States and in Brazil. New genus for the Brazilian fungus. Compte rendu des seances de la Societe de biologie 105: 315–316. Arantes TD, Theodoro RC, Teixeira M de M, Bosco S de MG, Bagagli E (2016) Environmental mapping of Paracoccidioides spp. in Brazil reveals new clues into genetic diversity, biogeography and wild host association. PLoS neglected tropical diseases, 10, e0004606. Artagaveytia-Allende RC (1949) Some biological characteristics of the pathogenic fungi named Paracoccidioides brasiliensis and Paracoccidioides cerebriformis. Journal of dental research, 28: 242-247. Bagagli E, Theodoro RC, Bosco SMG, McEwen JG (2008) Paracoccidioides brasiliensis: phylogenetic and ecological aspects. Mycopathologia, 165, 197–207. Bao D, Kinugasa S, Kitamoto Y (2004) The biological species of oyster mushrooms (Pleurotus spp.) from Asia based on mating compatibility tests. Journal of wood science, 50, 162–168. Beheregaray LB, Caccone A (2007) Cryptic biodiversity in a changing world. Journal of Biology, 6, 9. Bouckaert R, Heled J, Kühnert D et al. (2014) BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLoS computational biology, 10, e1003537. Brown EM, McTaggart LR, Zhang SX et al. (2013) Phylogenetic analysis reveals a cryptic species Blastomyces gilchristii, sp. nov. within the human pathogenic fungus Blastomyces dermatitidis PLoS one, 8, e59237. Brummer E, Castaneda E, Restrepo A (1993) Paracoccidioidomycosis: an update. Clinical microbiology reviews, 6, 89–117. Bruns TD, Szaro TM (1992) Rate and mode differences between nuclear and mitochondrial small-subunit rRNA genes in mushrooms. Molecular biology and evolution, 9, 836–855. Burnham KP, Anderson DR (2003) Model selection and multimodel inference: a practical information-theoretic approach. New York: Springer-Verlag. Burt A, Dechairo BM, Koenig GL et al. (1997) Molecular markers reveal differentiation among isolates of Coccidioides immitis from California, Arizona and Texas. Molecular ecology, 6, 781–786. Bustamante-Simon B, McEwen JG, Tabares AM, Arango M, Restrepo-Moreno A (1985) Characteristics of the conidia produced by the mycelial form of Paracoccidioides brasiliensis. Sabouraudia, 23, 407–414. Cadena CD, Jiménez I, Zapata F. 2017 Atlantean evolution in Darwin's finches - issues and perspectives in species delimitation using phenotypic data. bioRxiv.

Calcagno AM, Niño-Vega G, San-Blas F, San-Blas G (1998) Geographic discrimination of Paracoccidioides brasiliensis strains by randomly amplified polymorphic DNA analysis. Journal of clinical microbiology, 36, 1733–1736. Cavalcanti SD, Levi JE, Dantas KC, Martins JE. (2005) Analysis of the genetic polymorphism of Paracoccidioides brasiliensis and Paracoccidioides cerebriformis" Moore" by random amplified polymorphic DNA (RAPD) and 28S ribosomal DNA sequencing: Paracoccidioides cerebriformis revisited. Revista do Instituto de Medicina Tropical de São Paulo. 47:119-23. Carrero LL, Niño-Vega G, Teixeira MM et al. (2008) New Paracoccidioides brasiliensis isolate reveals unexpected genomic variability in this human pathogen. Fungal Genetics and Biology, 45, 605–612. Coyne JA, Orr HA (2004) Speciation. Sinauer Associates Incorporated. De Albornoz MB (1971) Isolation of Paracoccidioides brasiliensis from rural soil in Venezuela. Sabouraudia, 9, 248–253. Desjardins CA, Champion MD, Holder JW et al. (2011) Comparative genomic analysis of human fungal pathogens causing paracoccidioidomycosis. PLoS genetics, 7, e1002345. Dettman JR, Jacobson DJ, Taylor JW (2003a) A multilocus genealogical approach to phylogenetic species recognition in the model eukaryote Neurospora. Evolution. 57, 2703. Dettman JR, Jacobson DJ, Turner E, Pringle A, Taylor JW (2003b) Reproductive isolation and phylogenetic divergence in Neurospora: comparing methods of species recognition in a model eukaryote. Evolution. 57, 2721–2741. Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular biology and evolution, 29, 1969–1973. Dupuis JR, Roe AD, Sperling FAH (2012) Multi‐locus species delimitation in closely related animals and fungi: one marker is not enough. Molecular ecology, 21, 4422– 4436. Elena J, Himla S, Paula J et al. (1998) Mitochondrial mutation rate revisited: hot spots and polymorphism. Nature genetics, 18, 109–110. Elyashiv E, Bullaughey K, Sattath S et al. (2010) Shifts in the intensity of purifying selection: an analysis of genome-wide polymorphism data from two closely related yeast species. Genome research, 20, 1558–1573. Engelstaedter J. (2016) Asexual but not clonal: evolutionary processes in populations with automictic reproduction. bioRxiv. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular ecology, 14, 2611– 2620. Falush D, Stephens M, Pritchard JK (2003) Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics, 164, 1567–1587. Farlow A, Long H, Arnoux S et al. (2015) The spontaneous mutation rate in the fission yeast Schizosaccharomyces pombe. Genetics, 201, 737–744. Feitosa LDS, Cisalpino PS, Machado dos Santos MR et al. (2003) Chromosomal polymorphism, syntenic relationships, and ploidy in the pathogenic fungus Paracoccidioides brasiliensis. Fungal Genetics and Biology, 39, 60–69. Fisher MC, Koenig GL, White TJ, Taylor JW (2002) Molecular and phenotypic description of Coccidioides posadasii sp. nov., previously recognized as the nonCalifornia population of Coccidioides immitis. Mycologia, 94, 73–84. Fisher MC, Koenig GL, White TJ et al. (2001) Biogeographic range expansion into South

America by Coccidioides immitis mirrors New World patterns of human migration. Proceedings of the National Academy of Sciences, 98, 4558–4562. Fisher MC, Koenig G, White TJ, Taylor JW (2000) A test for concordance between the multilocus genealogies of genes and microsatellites in the pathogenic fungus Coccidioides immitis. Molecular biology and evolution, 17, 1164–1174. Fox J, Weisberg S (2011) An R Companion to Applied Regression. SAGE Publications. Frantz AC, Cellina S, Krier A, Schley L, Burke T (2009) Using spatial Bayesian methods to determine the genetic structure of a continuously distributed population: clusters or isolation by distance? Journal of Applied Ecology, 46, 493–505. Fritsch ES, Chabbert CD, Klaus B, Steinmetz LM (2014) A Genome-Wide map of mitochondrial dna recombination in yeast. Genetics, 198, 755–771. Gao Z, Przeworski M, Sella G. (2015) Footprints of ancient‐balanced polymorphisms in genetic variation data from closely related species. Evolution. 69:431-46. Giraudoux P (2011) pgirmess: data analysis in ecology. R package version 1.5. 8. Gonzalez A, Hernandez O (2016) New insights into a complex fungal pathogen: the case of Paracoccidioides spp. Yeast 33, 113–128. Guillamon JM, Barrio E, Huerta T, Querol A (1994) Rapid characterization of four species of the Saccharomyces sensu stricto complex according to mitochondrial DNA patterns. International Journal of Systematic and Evolutionary Microbiology, 44, 708–714. Hebeler-Barbosa F, Morais FV, Montenegro MR et al. (2003) Comparison of the sequences of the internal transcribed spacer regions and PbGP43 genes of Paracoccidioides brasiliensis from patients and armadillos (Dasypus novemcinctus). Journal of clinical microbiology, 41, 5735–5737. Hothorn T, Bretz F, Westfall P (2008) Simultaneous Inference in General Parametric Models. Biometrical journal, 50, 346–363. Hudson RR, Coyne JA (2009) Mathematical consequences of the genealogical species concept.. Evolution, 56, 1557-1565. Hughes KW, Petersen RH, Lickey EB (2009) Using heterozygosity to estimate a percentage DNA sequence similarity for environmental species’ delimitation across basidiomycete fungi. New Phytologist, 182, 795–798. Imai T, Sano A, Mikami Y et al. (2009) A new PCR primer for the identification of Paracoccidioides brasiliensis based on rRNA sequences coding the internal transcribed spacers (ITS) and 5·8S regions. Medical Mycology. Jack W Sites J, Marshall JC (2004) Operational criteria for delimiting species. Annual Review of Ecology, Evolution, and Systematics, 35: 199-227. Jahnke KD, Bahnweg G, Worrall JJ (1987) Species delimitation in the Armillaria mellea complex by analysis of nuclear and mitochondrial DNAs. Transactions of the British Mycological Society, 88, 572–575. Johann S, Sá NP, Lima LA et al. (2010) Antifungal activity of schinol and a new biphenyl compound isolated from Schinus terebinthifolius against the pathogenic fungus Paracoccidioides brasiliensis. Annals of Clinical Microbiology and Antimicrobials, 9, 30. Jombart T (2008) adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics (Oxford, England), 24, 1403–1405. Jombart T, Pontier D, Dufour AB (2009) Genetic markers in the playground of multivariate analysis. Heredity, 102, 330–341. Kasuga T, White TJ, Koenig G et al. (2003) Phylogeography of the fungal pathogen Histoplasma capsulatum. Molecular ecology, 12, 3383–3401. Koufopanou V, Burt A, Taylor JW (1997) Concordance of gene genealogies reveals reproductive isolation in the pathogenic fungus Coccidioides immitis. Proceedings of

the National Academy of Sciences, 94, 5478–5482. Koufopanou V, Burt A, Szaro T, Taylor JW (2001) Gene genealogies, cryptic species, and molecular evolution in the human pathogen Coccidioides immitis and relatives (Ascomycota, Onygenales). Molecular biology and evolution, 18, 1246–1258. Ladoukakis ED, Eyre-Walker A (2004) Direct evidence of recombination in human mitochondrial DNA. Heredity 93: 321-321. Leducq J-B, Lou Nielly-Thibault, Charron G et al. (2016) Speciation driven by hybridization and chromosomal plasticity in a wild yeast. Nature Microbiology, 1, 15003. Leffler EM, Bullaughey K, Matute DR et al. (2012) Revisiting an old riddle: what determines genetic diversity levels within species? PLoS Biology, 10, e1001388. Liti G, Barton DBH, Louis EJ (2006) Sequence diversity, reproductive isolation and species concepts in Saccharomyces. Genetics, 174, 839–850. Lynch M, Sung W, Morris K et al. (2008) A genome-wide view of the spectrum of spontaneous mutations in yeast. Proceedings of the National Academy of Sciences of the United States of America, 105, 9272–9277. Madhosingh C. (1992) Interspecific hybrids between Fusarium oxysporum lycopersici and Fusarium graminearum by mycelial anastomoses. Journal of Phytopathology. 136:113-23. Restrepo AM (2009) The ecology of Paracoccidioides brasiliensis: a puzzle still unsolved. Sabouraudia: Journal of Medical and Veterinary Mycology. Maddison WP, Maddison DR (2001) Mesquite: a modular system for evolutionary analysis. (2001). Mahalanobis PC (1936) On the generalized distance in statistics. Proceedings of the National Institute of Sciences (Calcutta), 2, 49–55. Matute DR, McEwen JG, Puccia R et al. (2006a) Cryptic speciation and recombination in the fungus Paracoccidioides brasiliensis as revealed by gene genealogies. Molecular biology and evolution, 23, 65–73. Matute DR, Quesada-Ocampo LM, Rauscher JT, McEwen JG (2008) Evidence for positive selection in putative virulence factors within the Paracoccidioides brasiliensis species complex. PLoS neglected tropical diseases, 2, e296. Matute DR, Sepúlveda VE, Quesada LM et al. (2006b) Microsatellite analysis of three phylogenetic species of Paracoccidioides brasiliensis. Journal of clinical microbiology, 44, 2153–2157. Matute DR, Torres IP, Salgado-Salazar C, Restrepo A, McEwen JG (2007) Background selection at the chitin synthase II (chs2) locus in Paracoccidioides brasiliensis species complex. Fungal Genetics and Biology, 44, 357–367. Mayr E (1963) Animal Species and Evolution. Harvard University Press, Cambridge, MA and London, England. Medeiros SI, Falcomer FC, Correa AA, Oliveira SA, Souza JM, Raimundo CJ, Grace MK, Antônio IC, Ribeiro AM, Burguel PH, Felipe MS, Tavares AH, Bocca AL. (2016) Distinct patterns of yeast cell morphology and host responses induced by representative strains of Paracoccidioides brasiliensis (Pb18) and Paracoccidioides lutzii (Pb01). Medical mycology. 54: 177-188 Molinari-Madlum EE, Felipe MS, Soares CM (1999) Virulence of Paracoccidioides brasiliensis isolates can be correlated to groups defined by random amplified polymorphic DNA analysis. Medical Mycology, 37, 269–276. Montoya AE, Moreno MN, Restrepo A, McEwen JG (1997) Electrophoretic Karyotype of Clinical Isolates of Paracoccidioides brasiliensis. Fungal Genetics and Biology, 21, 223–227.

Moore, M., (1935). A new species of the Paracoccidioides Almeida (1930): P. Cerebriformis Moore,(1935). Moore WS (1995) Inferring Phylogenies from mtDNA Variation: Mitochondrial-Gene Trees Versus Nuclear-Gene Trees. Evolution, 49, 718. Morais FV, Barros TF, Fukada MK, Cisalpino PS, Puccia R (2000) Polymorphism in the gene coding for the immunodominant antigen gp43 from the pathogenic fungus Paracoccidioides brasiliensis. Journal of clinical microbiology, 38, 3960–3966. Muñoz JF, Gallo JE, Misas E et al. (2014) Genome update of the dimorphic human pathogenic fungi causing paracoccidioidomycosis. PLoS neglected tropical diseases, 8, e3348. Rosenblum EB, Parent CE, Brandt EE. (2014).The molecular basis of phenotypic convergence. Annual Review of Ecology, Evolution, and Systematics. 45: 203-26. Muñoz JF, Farrer RA, Desjardins CA, Gallo JE, Sykes S, Sakthikumar S, Misas E, Whiston EA, Bagagli E, Soares CM, Teixeira MD, Taylor JW, Clay OK, McEwen JG, Cuomo CA. (2016) Genome diversity, recombination, and virulence across the major lineages of Paracoccidioides. mSphere. 1: e00213-16. Nascimento É, Martinez R, Lopes AR et al. (2004) Detection and selection of microsatellites in the genome of Paracoccidioides brasiliensis as molecular markers for clinical and epidemiological studies. Journal of clinical microbiology, 42, 5007– 5014. Naumov GL (1996) Genetic identification of biological species in the Saccharomyces sensu stricto complex. Journal of Industrial Microbiology & Biotechnology. Nichols R (2001) Gene trees and species trees are not the same. Trends in Ecology & Evolution, 16, 358–364. Nino-Vega GA, Calcagno AM, San-Blas G et al. (2000) RFLP analysis reveals marked geographical isolation between strains of Paracoccidioides brasiliensis. Medical Mycology, 38, 437–441. Novembre J, Stephens M (2008) Interpreting principal component analyses of spatial population genetic variation. Nature genetics, 40, 646–649. O'Donnell K, Kistler HC, Cigelnik E, Ploetz RC (1998) Multiple evolutionary origins of the fungus causing Panama disease of banana: concordant evidence from nuclear and mitochondrial gene genealogies. Proceedings of the National Academy of Sciences, 95, 2044–2049. Petit RJ, Excoffier L (2009) Gene flow and species delimitation. Trends in Ecology & Evolution, 24, 386–393. Pickrell JK, Pritchard JK (2012) Inference of population splits and mixtures from genome-wide allele frequency data. PLoS genetics, 8, e1002967. Posada D (2008) jModelTest: phylogenetic model averaging. Molecular biology and evolution, 25, 1253–1256. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 155, 945–959. Rambaut A (2009) FigTree. Rambaut A, Drummond AJ (2013) TreeAnnotator v1. 7.0. Available as part of the BEAST package. Restrepo A, McEwen JG, Castaneda E (2001) The habitat of Paracoccidioides brasiliensis: how far from solving the riddle? Medical Mycology, 39, 233–241. Restrepo S, Tabima JF, Mideros MF, Grünwald NJ, Matute DR (2014) Speciation in Fungal and Oomycete plant pathogens. dx.doi.org, 52, 289–316. Restrepo-Moreno A, Schneidau JD (1967) Nature of the skin-reactive principle in culture filtrates prepared from Paracoccidioides brasiliensis. Journal of bacteriology, 93, 1741–1748.

Roberto TN, Rodrigues AM, Hahn RC, de Camargo ZP (2016) Identifying Paracoccidioides phylogenetic species by PCR-RFLP of the alpha-tubulin gene. Medical Mycology, 54, 240–247. Robledo VM, Ospina CS, Restrepo IM (1968) Distribution of paracoccidioidin sensitivity in Colombia. Am J Trop Med Hyg 17: 25-37. Roca MG, Davide LC, Davide LM, Mendes-Costa MC, Schwan RF, Wheals AE. 2004. Conidial anastomosis fusion between Colletotrichum species. Mycological research. 108:1320-6. Saari S, Faeth SH (2012) Hybridization of Neotyphodium endophytes enhances competitive ability of the host grass. New Phytologist, 195, 231–236. Salgado-Salazar C, Jones LR, Restrepo A, McEwen JG (2010) The human fungal pathogen Paracoccidioides brasiliensis (Onygenales: Ajellomycetaceae) is a complex of two species: phylogenetic evidence from five mitochondrial markers. Cladistics, 26, 613–624. Sanders IR. 2006. Rapid disease emergence through horizontal gene transfer between eukaryotes. Trends in Ecology & Evolution. 21:656-8. Santos JC, Coloma LA, Summers K, Caldwell JP, Ree R, Cannatella DC. (2009) Amazonian amphibian diversity is primarily derived from late Miocene Andean lineages. PLoS Biol. 10:e1000056. Shimodaira H, Hasegawa M (2001) CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics (Oxford, England), 17, 1246–1247. Sievers F, Wilm A, Dineen D et al. (2011) Fast, scalable generation of high‐quality protein multiple sequence alignments using Clustal Omega. Molecular Systems Biology, 7, 539–539. Simwami SP, Khayhan K, Henk DA et al. (2011) Low diversity Cryptococcus neoformans variety grubii Multilocus sequence types from thailand are consistent with an ancestral African origin. PLOS Pathogens, 7, e1001343. Sites JW Jr, Marshall JC (2003) Delimiting species: a Renaissance issue in systematic biology. Trends in Ecology & Evolution, 18, 462–470. Smith ML, Anderson JB (1989) Restriction fragment length polymorphisms in mitochondrial DNAs of Armillaria: identification of North American biological species. Mycological Research, 93, 247–256. Song S, Pursell ZF, Copeland WC et al. (2005) DNA precursor asymmetries in mammalian tissue mitochondria and possible contribution to mutagenesis through reduced replication fidelity. Proceedings of the National Academy of Sciences, 102, 4990–4995. Söding J, Biegert A, Lupas AN (2005) The HHpred interactive server for protein homology detection and structure prediction. Nucleic acids research, 33, W244– W248. Stamatakis A (2006) RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics, 22, 2688–2690. Stamatakis A (2014) RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics, 30, 1312–1313. Sukumaran J, Knowles LL. (2017) Multispecies coalescent delimits structure, not species. Proceedings of the National Academy of Sciences. 114:1607-12.

Tajima F (1989) The effect of change in population size on DNA polymorphism. Genetics, 123, 597–601. Tamayo D, Muñoz JF, Lopez Á et al. (2016) Identification and analysis of the role of superoxide dismutases isoforms in the pathogenesis of Paracoccidioides spp. PLoS neglected tropical diseases, 10, e0004481. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular

Evolutionary Genetics Analysis version 6.0. Molecular biology and evolution, 30, 2725–2729. Taylor JW, Jacobson DJ, Fisher MC (2003) The evolution of asexual fungi: reproduction, speciation and classification. Annual Review of Phytopathology 37: 197–246. Taylor JW, Jacobson DJ, Kroken S et al. (2000) Phylogenetic species recognition and species concepts in Fungi. Fungal Genetics and Biology, 31, 21–32. Taylor JW, Turner E, Townsend JP, Dettman JR, Jacobson D (2006) Eukaryotic microbes, species recognition and the geographic limits of species: examples from the kingdom Fungi. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 361, 1947–1963. Team RC (2014) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2013. Teixeira M de M, Theodoro RC, Derengowski LDS et al. (2013) Molecular and morphological data support the existence of a sexual cycle in species of the genus Paracoccidioides. Eukaryotic Cell, 12, 380–389. Teixeira MM, Theodoro RC, de Carvalho MJA et al. (2009) Phylogenetic analysis reveals a high level of speciation in the Paracoccidioides genus. Molecular Phylogenetics and Evolution, 52, 273–283. Teixeira MM, Theodoro RC, Niño-Vega G, Bagagli E, Felipe MSS (2014) Paracoccidioides species complex: ecology, phylogeny, sexual reproduction, and virulence. PLOS Pathogens, 10, e1004397. Teixeira MM, Theodoro RC, de Oliveira FF, Machado GC, Hahn RC, Bagagli E, SanBlas G, Felipe MS. (2014) Paracoccidioides lutzii sp. nov.: biological and clinical implications. Medical mycology. 52:19-28. Theodoro RC, de Melo Teixeira M, Felipe MSS et al. (2012) Genus Paracoccidioides: species recognition and biogeographic aspects. PLoS one, 7, e37694. Theodoro RC, Scheel CM, Brandt ME, Kasuga T, Bagagli E (2013) PRP8 intein in cryptic species of Histoplasma capsulatum: evolution and phylogeny. Infection, genetics and evolution : journal of molecular epidemiology and evolutionary genetics in infectious diseases, 18, 174–182. Torres I, García AM, Hernández O et al. (2010) Presence and expression of the mating type locus in Paracoccidioides brasiliensis isolates. Fungal Genetics and Biology, 47, 373–380. Torriani SF, Penselin D, Knogge W, Felder M, Taudien S, Platzer M, McDonald BA, Brunner PC. (2014) Comparative analysis of mitochondrial genomes from closely related Rhynchosporium species reveals extensive intron invasion. Fungal Genetics and Biology. 62: 34-42. Uyeda JC, Hansen TF, Arnold SJ, Pienaar J (2011) The million-year wait for macroevolutionary bursts. Proceedings of the National Academy of Sciences of the United States of America, 108, 15908–15913. van Burik JAH, Schreckhise RW, White TC, Bowden RA, Myerson D (2009) Comparison of six extraction techniques for isolation of DNA from filamentous fungi. Medical Mycology 36: 299-303. Wang H, Xu Z, Gao L, Hao B. (2009) A fungal phylogeny based on 82 complete genomes using the composition vector method. BMC evolutionary biology. 10:195. Yoo YB, You CH, Park YH, Lee YH, Chang KY, Peberdy JF. (1987) Interspecific protoplast fusion and sexuality in Pleurotus. The Korean Journal of Mycology.;15:135-41. Yule GU (1925) A mathematical theory of evolution, based on the conclusions of Dr. JC Willis, FRS. Philosophical Transactions of the Royal Society of London. Series B, Containing Papers of a Biological Character 213: 21–87.

HIGHLIGHTS

1. Paracoccidioides brasiliensis is composed by four species. 2. Hybridization among Paracoccidioides species has occurred but is rare.

3. Mitochondrial introgression occurred in the Paracoccidioides genus in the recent past.

4. Paracoccidioides species show slight differences in their yeast morphology.