A total of 18 Colletotrichum isolates were included in this study. Seven of these isolates (H-8, 05-161, 05-200, H-3, 05-155, H-24, H-19) were retrieved from fern fronds previously collected between USA and Costa Rica (MacKenzie et al., 2009). Colletotrichum filicis CBS 101611 (Damm et al., 2012; Crous et al., 2021) is the ex-type the only genome available, and for these reasons was selected as reference genome in our dataset too. Isolates and data collected from each isolate are listed in Table 1 .

All isolates sequenced in this work were grown on PDA plates for seven days before DNA extraction. The fungal mycelium was scraped from the surface of a PDA plate using a sterile scalpel and transferred to a sterile 1.5 mL tube. Genomic DNA was extracted from the isolates using a modified CTAB method (Prodi et al., 2011).

Mycelium was ground by adding 700 μL of 3% CTAB solution using the Motor-Driven Tissue grinder G50 (Coyote Bioscience Company, Beijing, China). The tubes were placed in a water bath at 65°C for 30 minutes and then centrifugated at 13,225 x g for 10 minutes. 500 μL of the supernatant was transferred in a new sterile 1.5 mL tube and 500 μL of isoamyl alcohol-chloroform (1:24) was added and mixed. After centrifugation at 13,225 x g for 10 minutes, the supernatant was transferred into a new sterile tube. DNA was precipitated by adding an equal volume of ice-cold iso-propanol. The samples were incubated overnight at -20°C. The DNA was centrifuged at 13,225 x g at 4°C for 10 minutes and the isopropanol was removed. Two washes of the pellet were then performed by the addition of 500 μL of 70% ethanol and a subsequent centrifugation at 13,225 x g at 4°C for 5 minutes. The DNA was resuspended in 50 μL of sterile nuclease-free water, quantified, and assessed for quality using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, DE, USA). DNA was then stored at -20°C before sequencing.

The library preparation and genome sequencing of the Colletotrichum genomes were both performed by an external service provider, Personal Genomics, located at Via Roveggia 43b, 37136, Verona, Italy. The sequencing was carried out using the Illumina NovaSeq 6000 platform with a 150 bp paired-end configuration. The quality of the reads was evaluated using FastQC v0.12.1 (Babraham Bioinformatics, Cambridge, UK). Sequence adapters and low-quality reads were trimmed with Trimmomatic v0.33 (Bolger et al., 2014). Pair-end reads were merged with FLASH v1.2.11 (Magoč and Salzberg, 2011). Merged and unmerged reads were then assembled using SPAdes v3.15.1 (Bankevich et al., 2012). Scaffolds with low coverage were removed as possible contaminations. The completeness of the assembly was assessed using BUSCO v3.1 (Simão et al., 2015) while statistics were evaluated with QUAST v5.0.2 (Gurevich et al., 2013). Gene prediction was performed using AUGUSTUS v3.5.0 (Stanke et al., 2006). We use data produced by previous work (Baroncelli et al., 2024) to build the gene model for the closely related species Colletotrichum lupini; RNAseq data publicly available (SRX2782478) were assembled using rnaSPAdes v3.15.1 (Bushmanova et al., 2019), assembled transcripts were used to train the genome assembly (GCA_030913515.1) through AUGUSTUS v3.5.0 (Stanke et al., 2008). To extract protein-coding sequences from the GFF3 output, we used a pre-configured Perl script, getAnnoFasta.pl, provided with the AUGUSTUS v3.5.0 package. Finally, a custom script, (ProtRename.py available at https://github.com/RiccardoBaroncelli) was employed to rename proteomes file, streamlining downstream analysis. The details of the bioinformatic approach employed are illustrated in the workflow shown in Figure 1A .

For each genome, the final assembly was imported in Geneious Prime v2024.0.2 (https://www.geneious.com) where a local database was created. BLAST analyses were performed on the local databases using as query reference sequences for each locus. Seven loci were retrieved and used for phylogenetic analysis: the internal transcribed spacer (ITS) region, a partial sequence of the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene, the glutamine synthetase (GS) gene, the partial sequence of the beta-tubulin 2 (TUB2) gene, the histone-3 (HIS-3), the chitin synthetase gene (CHS-1) and actin (ACT) (Damm et al., 2012). The sequences of the seven genes of each analyzed genomes were extracted using the ‘extract reads’ plug-in. These sequences were then aligned using MAFFT v7.525 (Katoh and Standley, 2013) and manually adjusted, where necessary, with Geneious Prime v.2024.0.2. The multiple sequence alignments were exported to MEGA11 (Tamura et al., 2021), where the optimal substitution model for each individual dataset was calculated. The multi-locus alignment was performed using Geneious Prime v2024.0.2. Phylogenetic analyses were conducted using Randomized Axelerated Maximum Likelihood (ML) (Stamatakis, 2006), Bayesian Inference (BI) (Ronquist and Huelsenbeck, 2003), and Approximate Maximum Likelihood (Price et al., 2010) for each gene and for the concatenated genes. Maximum likelihood analyses were constructed with the RAxML v8.2.7 software (Stamatakis, 2014) using the GTR CAT model with 1,000 bootstrap replicates. A Bayesian phylogenetic analysis was conducted using a Markov Chain Monte Carlo (MCMC) algorithm in MrBayes v3.6.2 (Ronquist et al., 2012), employing the “nst=2” and a rate variation “=equal” parameters. Four Markov Chain Monte Carlo (MCMC) chains were run from random starting trees for 1,000,000 generations, with sampling occurring every 1,000 generations. The initial 25% of the trees generated were discarded as burn-in, and posterior probabilities (PP) were calculated based on the remaining trees (AvgStdDev = 0.006403). For the approximate maximum likelihood analysis, FastTree v2.1.11 (Price et al., 2010) was used with standard settings as implemented in Geneious Prime v2024.0.2. Phylogenetic trees were compared as described by Shimodaira and Hasegawa (1999), Steenwyk et al. (2023) and Tsang et al. (2017).

To prepare the proteomes, an AUGUSTUS built-in Perl script (getAnnoFasta.pl) was used to extract the CDS and amino acid sequences from the AUGUSTUS output, and a renaming script was applied to standardize the names of the amino acid files for clarity. The phylogenomic analysis performed was based on clustering proteins into orthologous groups using Orthofinder v3.0 (Emms and Kelly, 2019), for the identification of single-copy orthologues. Single-copy orthologues sequences were aligned using MAFFT v7.525 (Katoh and Standley, 2013). To enhance alignment quality, we used Gblocks v0.91 (Castresana, 2000) to eliminate poorly aligned regions and gaps that could introduce noise into the phylogenetic analysis. Gblocks v0.91 automates the trimming process, retaining only conserved regions that are highly informative for phylogenetic reconstruction. After generating and trimming alignments for each aminoacid sequence alignment, we concatenated the individual alignments into a single dataset. To evaluate the quality of these alignments, we employed AMAS (Borowiec, 2016), which computes key alignment statistics, including alignment length and the number of conserved sites. These metrics were used to assess the robustness and informativeness of each alignment. The phylogenomic tree was constructed using MrBayes, FastTree, RAxML and MEGA, as described in the “MLST Alignment and Phylogenetic analysis” section, which inferred the tree through either Bayesian inference or approximate maximum likelihood methods.

The details of the bioinformatic approach employed are illustrated in the workflow shown in Figure 1B . This comprehensive pipeline integrated state-of-the-art tools with systematic quality control measures to ensure the generation of high-quality phylogenomic data. Key steps included robust gene prediction, precise orthogroup clustering, and accurate phylogenetic tree inference, all aimed at elucidating evolutionary relationships with high confidence.

Genome completeness scores, evaluated using BUSCO, indicate a high level of assembly quality across all analyzed genomes and species, with values ranging from 91.5% to 98.4% (calculated as the percentage of Sordariomycetes Single-Copy Orthologs identified in the genomes [C:] Figure 2 ). This high completeness reflects the retention of most single-copy orthologs. All genomes display minimal levels of missing or fragmented BUSCOs ([F:] from 0.3% to 7% and [M:] from 0.1% to 1.5% respectively), emphasizing the robustness of the assemblies. Duplication rates [D:] are generally low, staying below 1% for most isolates. However, a few isolates, notably 05-200, H-3, and H-24, exhibit slightly higher rates of duplication and fragmentation. These deviations may point to subtle strain-specific differences in genome structure, such as gene family expansions or assembly challenges, possibly influenced by repetitive sequences. The assembly statistics were computed using QUAST to assess the quality of the assemblies. The assembly sizes range from approximately 47 Mb to 58 Mb ( Figure 2 ), which aligns with the expected genome size for Colletotrichum species. Additionally, the GC content remains stable at approximately 51–52%, consistent with previous studies. However, it is important to note that the actual genome size may be larger, and the GC content may be lower, as the genomes were sequenced exclusively using short reads. Consequently, repetitive elements may have collapsed during the assembly process (Talhinhas et al., 2017; Kooij and Pellicer, 2020; Baroncelli et al., 2024). These metrics not only underscore the typical genomic features of the genus but also suggest that the sequenced isolates reflect a diverse representation of Colletotrichum’s genetic landscape.

The phylogenetic analyses of single-gene trees ( Figures 3A–G ) reveal diverse topologies and distinct clustering patterns among Colletotrichum species, emphasizing the complexity of their evolutionary relationships.

The ACT gene tree ( Figure 3A ) shows that isolate H-3 clusters with C. tamarilloi, whereas all other isolates group with C. filicis. In contrast, the CHS-1 gene tree ( Figure 3B ) associates isolate H-8 with C. tamarilloi, reflecting variability in gene-specific relationships. Meanwhile, the GAPDH gene tree ( Figure 3C ) presents a unique scenario, with isolates 05-155 and 05-200 forming a distinct cluster separate from the remaining isolates, suggesting potential divergence or lineage-specific variation. The GS gene tree ( Figure 3D ) highlights further variability. Here, only isolate 05-155 clusters with C. filicis, whereas the remaining isolates form a distinct group closely related to C. paranaense. The HIS-3 tree ( Figure 3E ) displays more consistency, with all isolates clustering near and grouping with C. filicis. A similar pattern is observed in the ITS gene tree ( Figure 3F ), where the isolates align closely with C. filicis, but with C. abscissum also appearing within the group, suggesting shared ancestral or conserved sequences. The TUB gene tree ( Figure 3G ) identifies two distinct groups: one composed of isolates from Florida, USA (05-155, 05-161, 05-200) clustering with C. tamarilloi, and another including isolates from Costa Rica (H-3, H-8, H-19, H-24) associated with C. filicis. These geographic and genetic splits may reflect potential population-level differentiation within the species complex.

None of the individual gene trees fully resolves the observed taxonomic ambiguities. However, trees derived from loci such as HIS-3, ITS, and TUB exhibit a closer alignment with the combined tree. An analysis of the correlation between tree resolution and the percentage of variable sites within alignments reveals a nuanced pattern. The ITS region, with the lowest percentage of variable sites (1.9%), demonstrates higher taxonomic concordance, while HIS-3 and TUB, with moderate variability (7.6% and 6.1%, respectively), also produce relatively well-resolved trees. Intermediate variability is observed in ACT and CHS-1 (6.9% and 4.2%, respectively), where Colletotrichum sp. mostly cluster with C. filicis, apart from one divergent isolate. Conversely, GAPDH and GS alignments, which have the highest variability (14.9% and 9.3%), generate the most divergent tree topologies. These findings suggest that, in this case, loci with lower variability are more effective in resolving taxonomic ambiguities compared to those with higher variability. This evidence may suggest that the alignment tools employed in the analysis may face challenges in accurately handling loci with high variability, potentially affecting the resolution of taxonomic relationships.

In the concatenated tree based on seven genes ( Figure 3H ), all 17 isolates form a single clade (Clade 1). Within this clade, isolate 05-155 clusters with C. filicis CBS 101611, while the remaining isolates are grouped together in a sister clade. The node connecting the sister clade to the group comprising C. filicis isolate CBS101611 and isolate 05-155 consistently exhibits high support across all three phylogenetic analyses. Support values are notably robust, with scores of 1.00, 0.98, and 83.2 from MrBayes, FastTree, and RAxML, respectively. However, several nodes show low bootstrap support, and the topologies inferred by the different phylogenetic methods are not always consistent, indicating limited statistical confidence in certain phylogenetic relationships. For example, the topologies generated by FastTree and RAxML reveal slight discrepancies compared to the MrBayes tree. Specifically, the nodes connecting C. tamarilloi, C. costaricense, and C. abscissum to the C. filicis cluster lack bootstrap support in the FastTree and RAxML trees. In contrast, the MrBayes tree provides low support values for these nodes (0.55 and 0.74), except for one node, which shows moderately high support (0.91). These limitations make the interpretation of the concatenated tree robust for definitive Colletotrichum species identification but challenging for resolving intraspecific evolutionary relationships, highlighting the constraints of multilocus sequence typing (MLST) in addressing certain phylogenetic ambiguities.

These findings underscore the importance of integrating single-gene and concatenated analyses to gain a more comprehensive understanding of Colletotrichum species’ evolutionary relationships, while accommodating potential discrepancies among loci.

The phylogenomic analysis was performed using a consensus alignment derived from 7,830 single-copy orthologous sequences across 18 strains. After trimming, the alignment spanned 3,900,074 sites, providing a comprehensive dataset for the study. Notably, only 46 undetermined characters were present, reflecting a near-complete alignment with minimal missing data. The dataset exhibited a high degree of conservation, with 3,711,450 sites conserved across all taxa. Additionally, 188,624 sites were identified as variable, including 28,310 singleton sites and 60,314 parsimony-informative sites, which provided critical phylogenetic signals for reconstructing evolutionary relationships. These metrics underscore the robustness and informativeness of the alignment, ensuring that the phylogenomic tree is based on a statistically solid foundation that captures both conserved regions and evolutionary divergence among the species studied.

The phylogenomic tree ( Figure 4 ) presents a detailed depiction of the evolutionary relationships among the nine species classified within Clade 1 of the acutatum species complex, as previously described by Damm et al. (2012) and by Crous et al. (2021), based on their shared orthologous gene content. Consistent with the MLST analyses, all isolates sequenced in this study and associated with Rumohra adiantiformis cluster closely with C. filicis CBS101611. The tree’s branching structure is robustly supported by high-confidence values generated through multiple analytical approaches, including Bayesian posterior probabilities, FastTree, RAxML, and MEGA bootstrap scores, ensuring reliability in the inferred phylogenetic relationships.

Orthogroup statistics ( Figure 4 ) further illuminate patterns of gene conservation and genome-specific adaptations. Across the 18 species analyzed, 260,819 genes were identified, with 256,333 genes (98.3%) assigned to orthogroups, underscoring a high degree of conservation across these genomes. Only 1.7% of the genes remained unassigned, reflecting a relatively small proportion of unique genes. A total of 19,218 orthogroups were identified, capturing the diversity and evolutionary significance of the gene content in these species. Within these orthogroups, 8,552 were found to contain genes present in all genomes, likely representing a core set of conserved genes essential for fundamental fungal functions. A smaller subset of 79 orthogroups, comprising 494 genes (0.2%), specific to individual genomes. These genome-specific orthogroups may indicate adaptations unique to particular species, potentially linked to ecological niches or pathogenicity-related traits. The high proportion of genes assigned to orthogroups (98.3%) highlights the evolutionary conservation among these fungal genomes. This conservation is also reflected in the strong support values observed in the phylogenomic tree. Although limited in number, the presence of genome-specific orthogroups offers a promising avenue for future research, particularly for uncovering unique functional traits or pathogenicity-related genes that may differentiate certain species or strains within this lineage.

Interestingly, the phylogenomic analysis reveals notable differences in topology compared to the MLST approach, particularly in the clustering of Rumohra adiantiformis-associated isolates and their geographic structuring ( Figure 5 ). The phylogenomic tree robustly places all R. adiantiformis-associated isolates together, forming a sister cluster to C. filicis CBS101611. This clustering is supported consistently across various methodologies, including MrBayes, FastTree, RAxML, and MEGA. By analyzing all single-copy genes conserved across genomes, the phylogenomic approach captures deeper evolutionary signals, resolving ambiguities that MLST cannot. By leveraging this extensive genomic dataset, the phylogenomic approach incorporates regions of the genome under evolutionary pressure, resulting in more precise and insightful interpretations. Additionally, the phylogenomic tree provides clearer resolution of geographic patterns. All seven isolates sequenced in this study cluster together but separately from C. filicis CBS101611.

Interestingly, the isolates associated with Rumohra adiantiformis formed two distinct clades, exhibiting a notable geographic pattern. One lineage consisted of isolates from the United States (05-161, 05-200, 05-155), while the other comprised isolates from Costa Rica (H-24, H-19, H-8, H-31) (MacKenzie et al., 2009). This geographic structuring suggests potential genetic differentiation between populations from these regions. However, this pattern was not evident in the MLST-based phylogenetic tree, where these isolates remained intermixed without clear separation ( Figure 5 ).

The ability of the phylogenomic approach to detect these patterns underscores its utility in understanding population structure and evolutionary dynamics. The clustering of isolates into distinct geographic lineages offers insights into potential adaptations to local environments or host preferences. The differences observed between MLST and phylogenomic topologies emphasize the limitations of MLST in resolving complex evolutionary relationships. While MLST relies on a limited set of genes and often overlooks subtle genomic signals, the phylogenomic approach analyzes thousands of single-copy orthologues, providing enhanced clarity and reliability in discerning species relationships. These findings demonstrate the superiority of the phylogenomic method in uncovering both evolutionary history and geographic differentiation, offering a robust framework for studying Colletotrichum species.