DNA was extracted using the E.Z.N.A. fungal DNA extraction kit® according to manufacturer’s instructions (Omega bio-tek Ltd., USA). Four markers were selected to generate sequences for quality analysis: partial Actin (ACT), partial β-Tubulin (TUB2), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and internally transcribed spacer regions I and II (ITS) of the rDNA gene regions. The universal primer pair ITS4/5 was used in PCR to amplify the ITS region (496 bp) of the nuclear ITS1-5.8S-ITS2 rDNA (White et al. 1990); Bt2a/b primers allowed amplification of a 560 bp fragment of the TUB2 gene region (Glass and Donaldson 1995); GDF/GDR (Templeton et al. 1992) and ACT-512 F/ACT-783R (Carbone and Kohn 1999) primers were used to amplify 300 bp and 290 bp fragments, respectively. The standard calmodulin (CL1/CL2A) primers (O’Donnell et al. 2000) were not successful in generating any amplions for most isolates and this gene region was not subsequently used in the analysis. Each 25 μL PCR reaction contained 1 × PCR buffer; 1.5 mM MgCl₂, 0.2 mM dNTP, 2.5 U Taq DNA Polymerase (Invitrogen by Life Technologies Co., USA) and 50 pmoles of each primer (Integrated DNA Technologies, USA). Thermal cycling conditions applied were: an initial denaturation of 5 min at 94°C followed by 35 cycles of 1 min at 94°C, 1 min at 55°C, 1 min at 72°C with a final extension of 5 min at 72°C. PCR products were sequenced directly (Amplicon Express, WA, USA).

Only holotype and epitype sequences of C. gloeosporioides sensu stricto, C. asianum, C. fructicola, C. siamense and C. tropicale were used in the study (Phoulivong et al., 2010; Cannon et al. 2012) (Table 1). These selected species have been analyzed in previous independent phylogenetic studies, and thus provides an objective reference in phylogenetic relationships.

¹Strain:

ATCC: American type culture collection, 10801 University Boulevard, Manassas, Virginia, USA.

CBS: Culture collection of the Centraalbureau voor Schimmelcultures, Fungal Biodiversity Centre, Utrecht, The Netherlands.

ICMP: International Collection of Microorganisms from Plants, Landcare Research, Auckland, New Zealand.

IMI: Culture collection of CABI Europe UK Centre, Egham, UK.

MFU: fungarium of Mae Fah Luang University, Thailand (cultures in BCC (BIOTEC Culture Collection, Thailand).

STE-U: Culture collection of the Department of Plant Pathology, University of Stellenbosch, South Africa.

In addition to the sequences derived in this study (Table 2), query sequences were also extracted from a previous study by Rampersad et al. (2013) and used in the final data set. Representative sequences were submitted to GenBank (GenBank: JQ218143, JQ218144, JQ218145, JQ218146). A total of 143 sequences for all markers was used in the final data set - 35 taxa for GAPDH; 34 taxa for ACT; 37 taxa for TUB2; 37 taxa for ITS. Alignments were carried out using the MAFFT server (http://mafft.cbrc.jp/alignment/server/) under default parameters. BioEdit software (http://www.mbio.ncsu.edu/BioEdit/bioedit.html). was used to edit the alignments of each marker data set. A schematic outlining the approach used to explore sequences is illustrated in Figure 1.

Homogeneity tests can indicate different rates of substitution characteristic of certain lineages across different genes (Hedges and Kumar 2003). Several factors can influence substitution rates including the functionality of proteins and RNA, generation time, metabolic rate, population size, and life histories (Thomas et al. 2006; Smith and Donoghue 2008). Base composition for each gene according to test species was calculated in MEGA 5 (Tamura et al. 2011). The disparity index (ID) corresponds to the observed difference in base composition between two sequences compared to the compositional difference that would be expected under homogeneity of the evolutionary process. I _D equals 0 when the homogeneity assumption is satisfied. Values greater than 0 indicate that the larger differences in base composition bias than expected under homogeneity. Disparity Index (I _D) per site was calculated for all sequence pairs. All positions containing gaps were omitted from the data set (complete deletion option). The significance of a given I _D value was then assessed using a Monte Carlo approach.

Base substitution saturation decreases the amount of phylogenetic information contained in a sequence data set and can disrupt analysis involving deep phylogeny (Xia and Xie 2001). Homoplasy due to multiple substitutions was tested with the index of substitution saturation (ISS) (Xia et al. 2003a), which assumes a critical index of substitution saturation (ISSc) that defines a threshold for significant saturation in the data. The level of substitution saturation was assessed by using the substitution saturation test of the program package DAMBE v. 5.3.46 software package (Xia et al. 2003b) which determines an “index of substitution saturation”, based on the notion of entropy in information theory.

Entropy is a measure of sequence variation and can be used to quantify the level of phylogenetic information available in sequence alignments (Shannon 1948; Krüger et al. 2012). A position has zero entropy (invariant alignment column) if it has the same character state in every sequence in the alignment, thereby making this nucleotide position phylogenetically uninformative. The maximum Shannon entropy is dependent on the number of discrete variable in the data set. We investigated the entropy landscape of each alignment by constructing entropy plots for each marker in BioEdit.

Multiple recombination detection methods were applied as implemented in RDP3 to detect recombination events in the aligned sequences for each gene. Using RDP3, four methods were applied: RDP (Martin and Rybicki 2000), GENECONV (Padidam et al. 1999), MaxChi (Maynard 1992) and 3Seq (Boni et al. 2007). Sequences were treated as linear and the threshold P-value was set at 0.05, using Bonferroni correction. All other settings were set to default. Methods like GENECONV and MAXCHI do not include phylogenetic information in their recombination detection algorithms and phyogenetic trees generated did not include any recombination event information.

The genetic distances were calculated under the Kimura 2 parameter (K2P) model using MEGA 5 for all sequences of the data set for each marker. Frequency histograms were generated to examine the intra- and possible interspecific variation among genetic distances. Using the genetic distances from pair-wise comparisons, the percentage comparisons with genetic distances greater than 5% are reported. These cut-off values are arbitrary but may indicate less inclusive clades. The distribution and size of the modes depend on the degree and range of sequence divergence and the number of intra-specific and interspecific taxa in the data set.

Errors can occur during amplification when two different sequences are used as template DNA. The PCR product may include a combination of these original sequences (Jumpponen 2011). Such chimeric sequences may be misinterpreted as novel which can artificially inflate estimates of diversity and interfere with phylogenetic inference and species discrimination if undetected (Hugenholtz and Huber 2003). ITS sequences were checked for putative chimeras using the UNITE PlutoF Chimera checker (Edgar et al. 2011) and the Chimera Test developed in the Fungal Metagenomics Project at the University of Alaska (https://biotech.inbre.alaska.edu/fungal_portal/?program=chimera_test).

Alignments were carried out using the online version of the sequence alignment program MAFFT version 6 ((http://mafft.cbrc.jp/alignment/server/) to ensure that the 5.8S region of the ITS array forms an anchor region in the middle of the alignment. This check will determine if the ITS sequences were composed of stochastic, artifactual nucleotide data.

Multiple sequence alignments of each gene were made and manually adjusted where necessary with BioEdit. Individual gene trees were reconstructed using the Bayesian Markov chain Monte Carlo (MCMC) approach as implemented in Mr. Bayes 3.2.1 (Huelsenbeck and Ronquist 2001), and Neighbour joining (NJ) and Maximum Likelihood (ML) algorithms implemented in PAUP* version 4.0b8 (D. L. Swofford, Sinauer Associates, Sunderland, MA, Swofford 2003). For each of the markers, tree topologies were compared and evaluated for disagreement.

The Mr. Bayes analyses were performed under the GTR + Γ + I model. To ensure accuracy, four parallel runs were carried out with one cold and three heated MCMC chains per run for 1,100,000 generations; trees were sampled every 200 generations and the first 100,000 trees were discarded as these represented the trees of the burn-in phase of the analysis. Convergence was assessed by examining the stationarity of the ln-likelihood and effective sample size which was >500 for all analyses. The remaining trees were summarized as a Maximum Clade Credibility (MCC) tree (i.e. tree topology with the highest clade posterior probabilities across all nodes).

Bootstrap-supported (1,000 replicates) maximum likelihood (ML) and neighbour-joining (NJ) trees were estimated under their respective best-fit models of nucleotide substitution as determined using ModelTest v.3.7 (Posada 2008) with corrected Akaike information criteria (AICc). The models used were TUB2: HKY85; ACT: F81; ITS: GTR and GAPDH: GTR. Separate sequences of the ITS1 and ITS2 loci were used in a second and third data set to generate Bayesian trees under the GTR + Γ + I model, and ML and NJ phylogenetic trees based on the best-fit evolutionary model which were HKY85 for ITS1 and GTR for ITS2.

We included sequence data from representative epitype isolates which were selected as authentic cultures by the CBS (Centraalbureau voor Schimmelcultures institute, Utretcht, Netherlands), IMI (CABI Europe – UK, Bakeham Lane, Egham, Surrey, UK) and/or as cultures listed as authentic by Cannon et al. (2012) (Table 1). We also sought to reduce the level of incongruence between individual gene trees by using the approach of conditional editing of sequences (Salichos and Rokas 2013) where (i) sites containing gaps and any aberrant sequences that produced a “bad” alignment were removed from the data set, (ii) quickly evolving taxa that are attracted by the distant outgroup or taxa that resulted in long branch attraction were removed, (iii) a limited number of taxa in each gene tree (query and epitype), (iv) only established or known ex-epitype and epitype sequences and (v) genes that recover specific internodes and where species placement has proven consistent among independent studies were used. FigTree was used to annotate trees and export graphics.

Tests to determine base substitution saturation were conducted for all markers separately, for combined ITS array data set, and for the stem and loop regions of ITS1 and ITS2 regions separately (Table 4). Results indicated that for each data set, ISS values were consistently lower than the ISSc values which indicated little saturation in base substitution. The test also revealed little difference between paired (stem) and unpaired positions (loops). In the stem and loop regions of all relevant ITS data sets, the ISS value was always lower significantly (P <0.0001) than the observed ISSc. The complete results of the saturation tests are summarised in Table 3. The findings of the saturation test suggested that there was no need to partition the stem and loop regions of the ITS1 and ITS2 data sets for phylogenetic analyses as these regions had little saturation.Estimation of saturation can also be determined based on the slope of the regression line in the saturation plot; data sets without any saturation have a slope =1, data sets with high saturation have a slope =0. Generally, base saturation by transitions was higher than base saturation by transversions except for the ITS markers which had near equivalent saturation for both transition- and transversion-type base substitution (Figure 2). Slopes for transversion-type and transition-type base substitutions were closer to 1 for ACT, GAPDH, TUB2 and ITS, in the order of least to most saturated.

¹Only the symmetrical tree topology is presented.

Shannon entropy distributions within the GAPDH, ACT, TUB2 and ITS DNA sequence alignments were examined. Regions of high entropy were the most dispersed in the GAPDH marker alignment. Entropy plots (Figure 3) revealed only two high entropy clusters extending greater than eight sites. GAPDH had the highest number of peaks over 0.25 (18 peaks with three over 0.5 maximum), followed by ACT with eight peaks over 0.25 (only one over 0.5 maximum), ITS had seven peaks (three over 0.5 maximum) and TUB2 had five peaks (all of which were over a 0.5 maximum. Similar results were obtained for the translated amino acid alignment for reading frames 1, 2 and 3 of the TUB2 gene, data not shown). With respect to the ACT marker, the plot indicated a distinct cluster of sites towards the end of the alignment (alignment position 110 to 141) with entropy values greater than 0.25. For the ITS alignments, most of the entropy values occurred in two distinct clusters and which corresponded to the more variable ITS1 and ITS2 loci that flank the 5.8S region.

Aligned sequences of each gene were analyzed with recombination detection software, using the three methods of analyses i.e. RDP, GENECONV and MaxChi, with the aim of identifying possible recombinant isolates and to locate the recombined regions in these recombinant sequences. Recombination was detected in the aligned ACT gene sequences. No significant recombination events above the P-value threshold of 0.05were detected for MAXCHI analyses. 3Seq detected 2 isolates with significant (P <0.05) recombination events; PAW-Cg-107 and PAW-Cg-41; however, as none of the other methods detected these sequences, it is not certain whether these isolates were actual recombinants. There was no evidence of recombination in any of the other gene sequences.

The genetic distances for all sequences of the isolates in the data set for each marker were calculated under the Kimura 2-parameter (K2P) model using MEGA5 (Table 5). Frequency histograms were generated to examine the intra- and possible inter-specific variation based on pair-wise genetic distances. Summary statistics for genetic distance data are presented in Table 3. There is an apparent difference in the ability of a given markers to distinguish between intra-specific and interspecific sequences according to marker. The ACT and GAPDH markers were able to produce a discernible modal gap in the frequency histograms which indicate separation of intra- and inter-specific isolates (Figure 4). There is some uncertainty for TUB2 whether a true gap exists because the separation is not an easily discernible bimodal distribution as for the other markers. No discernible gap was observed for the ITS marker.

¹Summary statistics of genetic distances calculated using the Kimura 2-paramter (K2P) model and include: Maximum genetic distance (Max GD), Minimum genetic distance (Min GD), Median genetic distance (Med GD) with standard deviation (Std Dev GD) and the number of pair-wise comparisons of genetic distances greater than or equal to 5% (≥5%).

²The number of pair-wise comparisons of genetic distances greater than or equal to 2% (≥2%) was also determined to be zero for ITS.

The phylogenetic hypotheses for the Bayesian MCC trees are described here (Figures 5, 6, 7, 8, 9 and 10). Individual loci were analyzed separately to assess the topological congruence among the different sequence data sets. The ability of each data set to efficiently resolve terminal lineages with robust branch support for relationships within the C. gloeosporioides sensu lato species complex was determined. The topologies for the Bayesian MCC, ML and NJ trees were compared for each marker and were found to be congruent.

The GAPDH Bayesian MCC tree (Figure 5) also revealed a distinct C. fructicola clade with one subclade consisting of sister taxa PAW-Cg-110 and PAW-Cg-8 similar to the ACT tree. C. tropicale formed a separate clade along with nine isolates but forming a polytomic phylogeny. The C. asianum clade consisted of PAW-Cg100 and PAW-Cg-101 with one epitype sequence. For the ACT and GAPDH analyses, four isolates (PAW-Cg-104, PAW-Cg-107, PAW-Cg-108, PAW-Cg-113) were commonly positioned in both trees.

For the ACT tree (Figure 6), there was a distinct C. fructicola clade but with only moderate clade support, and with two subclades. C. siamense, C. gloeosporioides sensu stricto, C. asianum and C. tropicale epitype sequences formed a polytomic phylogeny with seven isolates in the data set. PAW-Cg-107 and PAW-Cg-108 appeared to be sister taxa and formed a separate strongly-supported clade with PAW-Cg-104. Similarly PAW-Cg-109 and PAW-Cg-121 also appeared to be sister taxa forming a well-supported clade. PAW-Cg-118, PAW-Cg-107 and PAW-Cg-108 appear to be faster evolving or more genetically divergent isolates than the others given the comparatively longer branch lengths.

The TUB2 Bayesian MCC tree (Figure 7) indicated a strongly-supported separate C. tropicale clade, consisting of 12 isolates which was identical to the C. tropicale clade of the GAPDH tree. There may be an apparent C. siamense clade consisting of 15 isolates but these formed a basal polytomic phylogeny.

The ITS Bayesian MCC tree (Figure 8) displayed with a distinct C. gloeosporioides sensu stricto clade consisting only of the three epitype sequences used in the data set. Similarly, a distinct C. asianum clade was discernible with strong branch support. C. siamense, C. fructicola and C. tropicale epitype sequences formed a basal polytomic phylogeny with the majority of the isolates in the data set.

The separate ITS1 and ITS2 loci were considered. For the ITS1 region tree (Figures 9 and 10), C. fructicola and C. tropicale formed a polytomic phylogeny with the majority of the isolates. C. siamense formed a separate polytomic clade with seven isolates; the three C. gloeosporioides sensu stricto and the C. asianum epitypes formed two separate clades with higher branch support than for the ITS2 region tree. For the ITS2 region tree, C. fructicola and C. tropicale together with the majority of isolates formed a polytomic phylogeny; the three C. gloeosporioides sensu stricto epitypes, one C. siamense epitype and one query isolate and the two C. asianum epitypes formed separate moderately-supported clades.

Overall, analyses indicated that the ACT marker resolved only the C. fructicola clade; the GAPDH marker resolved the C. fructicola, C. tropicale and C. asianum clades; the TUB2 marker resolved the C. tropicale and C. siamense clades; and the ITS marker (entire) resolved the C. gloeosporioides sensu stricto and the C. asianum clades; the ITS1 region resolved the C. gloeosporioides sensu stricto and the C. asianum clades and the ITS2 region resolved the C. gloeosporioides sensu stricto, C. siamense and the C. asianum clades.