Single representative genome sequences of 23 species and 10 families within the order Helotiales were used for the phylogenetic analysis. These included the three Monilinia genomes reported in this study and 20 genomes retrieved in the whole-genome shotgun (WGS) database in GenBank, which were selected from among those Helotiales species that have high levels of completeness and annotations of the protein sequences (accessed April 12, 2019). A complete list of the genomes used is given in Supplementary Table 1. Single-copy orthologous proteins within Sclerotiniaceae or within Helotiales were identified with the software package GET_HOMOLOGUES v3.2.1 (Contreras-Moreira and Vinuesa, 2013), which clusters homologous protein families using the bidirectional best-hit, COGtriangles (-G -t 0 -D) and OrthoMCL (-M -t 0 -D) algorithms. A core-cluster was obtained by running the “compare_clusters.pl” script from the clusters generated by each of the two algorithms.

A phylogenetic tree was then constructed using the orthologous protein families identified by GET_HOMOLOGUES, which were then concatenated and aligned using MegAlign Pro (Lasergene v. 15.0.1; DNASTAR Inc., Madison, WI, United States) with the MAFFT algorithm (BLOSUM30 scoring matrix) and default settings. The tree was inferred using the maximum likelihood method implemented in MEGA7 (Kumar et al., 2016), based on the JTT matrix-based method (Jones et al., 1992) set to default parameters. Statistical support for the branches was evaluated using bootstrap analysis (Felsenstein, 1985) with 1,000 replicates.

Orthology-based conserved syntenic blocks among the genomes of the three Monilinia strains and B. cinerea and S. sclerotiorum (chosen as closely related outgroup species in the Sclerotiniaceae family) were analyzed using the CoGe computational pipeline² (Lyons and Freeling, 2008).

The SynMap tool (Lyons et al., 2008; Castillo et al., 2018) was used to identify syntenic regions and large-scale changes in genome organization. Specifically, for each genome dot plot, comparisons were performed in pairs. Tandem gene duplicates were identified and filtered by using the blast2raw programme. Orthologous genes were identified using the BlastZ algorithm, and collinear series of syntenic genes were identified using DAGchainer (Haas et al., 2004), with the “relative gene order” option and the maximum distance between two matches parameter set to 5. Syntenic Path Assembly (SPA) (Lyons et al., 2011) in the SynMap tool was used to generate a whole genome assembly based on synteny. The rate of synonymous substitution (Ks) values was calculated as the temporal calculation of synteny for the detected syntelogous gene pairs using the Needleman–Wunsch algorithm (Needleman and Wunsch, 1970) in the CodeML software (Haug-Baltzell et al., 2015) of the PAML package (Yang, 1997), integrated in SynMap. Ks values were plotted as a histogram, and a color scheme was applied across the histogram. We further estimated the Ks values by classifying five evolutionary classes, from the younger to the older orthologous gene pairs: (1) Ks divergence ≤ 0.2; (2) 0.21 ≤ Ks ≤ 0.4; (3) 0.41 ≤ Ks ≤ 1.0; (4) 1.1 ≤ Ks ≤ 79; and (5) Ks ≥ 80. Microsynteny analysis of specific regions was performed using the SynFind (Tang et al., 2015) and GEvo (Lyons et al., 2008) tools in CoGe. In detail, we analyzed the syntenic relationships among specific genes using genes from one genome against the other analyzed genomes as the query. Then high-resolution analysis of neighboring genomic regions was performed using GEvo for comparison of multiple genomic regions. In addition, we used SynFind to estimate syntenic depth among genomes, i.e., the number of times a region is syntenic to target genome regions, and to suggest the orthologous and paralogous or co-orthologous genomic regions. SynFind analysis was performed using the BLAST-variant with a gene window size of 40 and setting 4 as the minimum number of genes.

TEs in the three Monilinia genomes, as well as in the B. cinerea and S. sclerotorium genomes, were identified through both de novo and homology searching methods. As a first step, de novo repeat detection was performed on the assembled genomes using the software package RepeatModeler, version 1.0.10³ (Smit and Hubley, 2015), which integrates the following programmes: rmblastn v. 2.2.28; RECON v. 1.08; Tandem Repeats Finder v. 4.07b; and RepeatScout v. 1.0.5 (Benson, 1999; Bao and Eddy, 2002; Price et al., 2005; Smit et al., 2015). The sequences identified were then classified using CENSOR⁴ (Jurka et al., 1996; Kohany et al., 2006) according to the “fungi” TEs database (Kapitonov and Jurka, 2008), to obtain genome-specific libraries. For the homology-based approach by searching in TE repositories, RepeatMasker (v 4.0.7) was used to annotate the newly identified repeat elements on the genome assembly of each strain using the resulting library merged with libraries included in RepeatMasker (RepBase-20181026 library; RepeatMaskerLib.h5, Df Dfam, v 3.2; Date: July 2, 2020; Families: 318,520) to classify the unknown TE families and subfamilies, and to generate the complete annotations.

To estimate TE “age” and transposition history in the Monilinia species, we performed copy-divergence analysis of TEs based on their Kimura 2-parameter distances (K-values) (Kimura, 1980), as determined using calcDivergenceFromAlign.pl and createRepeatLandscape.pl (available on RepeatMasker util directory) on the alignment files. The proportions (%) and normalized numbers of TEs in each genome were calculated according to their classes.

With the aim to predict potentially active (autonomous) or degenerate (non-autonomous) TEs (Wicker et al., 2007; Santana et al., 2012), the presence of coding regions was assessed according to the following strategies. Overlap between TEs and annotated genes were detected using CLC Genomics Workbench and were submitted to InterProScan for conserved domain searching using the OmixBox 2.0.10 platform (Conesa et al., 2005). In addition, open reading frames (ORFs) were predicted using ORFfinder⁵ and analyzed using SmartBLAST⁶ against the UniProtKB/Swiss-prot database provided by the National Center for Biotechnology Information (NCBI⁷). The TEs selected were analyzed using the GEvo and FeatView tools on the CoGe platform to find their position features on the genomic sequences.

The whole-genome sequences of M. fructicola, M. fructigena, and M. laxa were subjected to repeat-induced point mutation (RIP) analysis using the RIPper programme (van Wyk et al., 2019), which measures changes in the frequencies of the targeted dinucleotides (substrates) and products that result from RIP. These analyses were performed using a 1,000-bp window size and 500-bp step size, with the default cut-off settings for the RIP parameters. Large RIP-affected regions (LRARs) were identified as genomic regions longer than 4,000 bp that included at least seven consecutive RIP-positive windows. For comparative purposes, all of these analyses were further performed on the genomes of B. cinerea and S. sclerotiorum. Changes in gene density for RIP-affected regions were investigated using the CLC Genomics Workbench. Dinucleotide frequencies were calculated using the RIPCAL tool (Hane and Oliver, 2008) with the default settings.

BLASTP searches (E-value < 1 × 10^–5) were used to identify genes encoding proteins involved in the RIP pathway in each Monilinia genome. These included homologs of the key enzymes described in Neurospora crassa, 5-cytosine methyltransferase RID (RIP deficient) and DIM-2 (defective in methylation), the RIP-associated cofactors DIM-3, -5, -7, -8, and -9, and HP-1 (heterochromatic protein) (Clutterbuck, 2011). Their expression in each Monilinia species was evaluated by exploiting RNA-seq data (De Miccolis Angelini et al., 2018).

Secreted proteins in the proteome of each strain were identified by selecting sequences with a signal peptide as predicted by SignalP 5.0 (Armenteros et al., 2019), and no multiple transmembrane helices or a single helix in the N-terminal 60 residues signal peptide, as predicted by TMHMM 2.0 (Krogh et al., 2001). Selected sequences were submitted to InterProScan annotation in OmicsBox 2.0.10. Effector protein candidates for M. fructicola, M. fructigena and M. laxa were predicted using the EffectorP 1.0 (Sperschneider et al., 2016) and EffectorP 2.0 (Sperschneider et al., 2018) prediction tools⁸. A combined EffectorP 1.0/2.0 classifier was finally used to achieve the highest prediction accuracy, to limit false-positive results. Candidate effector proteins were screened for functional annotation based on analysis of sequence homology and searches for conserved motifs, using the amino-acid sequences as an input and an E-value cut-off of 0.001. Carbohydrate-active enzymes (CAZymes) were annotated using dbCAN2 meta server (Zhang et al., 2018). Only CAZymes that were identified by at least two of the three tools, HMMER, DIAMOND and Hotpep, were selected for further analysis. Ortholog clusters of the predicted effectors and CAZymes were identified using OrthoVenn2 (Xu et al., 2019).

Putative BGCs for each Monilinia strain were predicted using AntiSMASH 5.2.0⁹ (Blin et al., 2019). The web-based software SMURF¹⁰ (Khaldi et al., 2010) was used to detect further potential gene clusters. The default parameters were used for both of these analyses. Furthermore, RNA-seq data (De Miccolis Angelini et al., 2018) were used to identify co-regulated gene clusters using FunGeneClusterS¹¹ (Vesth et al., 2016), with the parameter setting as follows: window-size = 3, gene skipping = 0 and Pearson correlation.

The results from the phylogenetic analysis are shown in Figure 1 as a maximum likelihood phylogenetic tree that was obtained using 129 single-copy orthologous proteins across 23 Helotiales genomes. The taxonomic relationships inferred by this analysis are congruent with the current family level taxonomy and with the family clades within Helotiales defined by Johnston et al. (2019). A first group of families that included Mollisiaceae, Drepanopezizaceae, and Ploettnerulaceae were grouped in the “mollisioid” clade with 96% bootstrap support. The family Discinellaceae member of the “discinelloid” clade was closely related to Hyaloscyphaceae (“hyaloscyphoid” clade) and distant from Helotiaceae (“helotioid” clade) and Dermateaceae. As expected, the Monilinia species were closely related to B. cinerea and Sclerotinia species in the Sclerotiniaceae family, which together with the sister family Rutstroemiaceae and the more distant Chlorociboriaceae, are recognized as members of the “sclerotinioid” clade. A strong molecular phylogenetic structure was observed with bootstrap values >95% for most tree key branches, with a few branches showing lower bootstrap values. The relationship among Sclerotinia borealis and the closest Sclerotinia and Botrytis species within Sclerotiniaceae was poorly supported by the bootstrap analysis. The phylogenetic distance inside Sclerotiniaceae was also confirmed by maximum likelihood analysis based on a larger set (2,006) of single-copy ortholog proteins identified among the same species of the family, with Glarea lozoyensis used as an outgroup (data not shown).

Syntenic gene pairs among the Monilinia genomes were identified using SynMap. There were totals of 13,749, 12,424 and 10,800 coding genes for M. fructicola, M. laxa and M. fructigena, respectively. M. fructicola returned 9,015 to 7,765 syntenic genes when paired with M. laxa and M. fructigena, respectively, while 7,813 genes were syntenic between the M. laxa and M. fructigena genomes. The overall comparison highlighted 7,398 genes shared among these genomes (Supplementary Figure 1). The Sclerotiniaceae inter-genus comparison showed high syntenic levels: 7,930, 7,827 and 6,824 of the 13,703 coding genes of B. cinerea were syntenic with M. laxa, M. fructicola and M. fructigena, respectively, while 7,640, 7,572 and 6,623 of the 11,130 S. sclerotiorum genes were syntenic with M. laxa, M. fructicola and M. fructigena, respectively. Overall, 5,715 coding genes were shared by the five fungal strains tested (Supplementary Figure 1).

The investigation of collinear orthologous genes revealed a highly conserved gene order among the Monilinia genomes. In particular, the syntenic collinear blocks identified between M. fructicola and M. laxa were extensively conserved, with an average of 180 genes per 50 syntenic blocks. Then, 136 and 143 syntenic blocks with an average of about 50 genes per block were detected in the M. laxa/M. fructigena and M. fructicola/M. fructigena comparisons, respectively. As expected, greater fragmentation of syntenic blocks was detected by comparing the Monilinia genomes with B. cinerea and S. sclerotiorum. In the comparisons with B. cinerea, a range of 411–419 syntenic blocks were counted, while they were in the range 364–392 in comparisons with S. sclerotiorum (Figure 2, Table 1, and Supplementary Table 2).

The dot plots of whole genome synteny performed with SPA tools are shown in Figure 3, while detailed SPA mapping results are reported in Supplementary Table 2. A high syntenic level was observed among the Monilinia genomes, while the syntenic dot plots that included B. cinerea and S. sclerotiorum showed many discontinuities linked to highly fragmented syntenic blocks (Figure 3).

The dot plot investigation between the M. fructicola and M. laxa genomes showed discontinuities of diagonal lines in several genomic regions included in syntenic blocks (Supplementary Figure 2 and Supplementary Table 2). Microsynteny investigation showed that these genome sequences corresponded to AT rich (GC% < 20%) non-coding homologous regions that ranged from ∼70,000 to 110,000 bp. In these regions, high densities of non-autonomous/degenerated TEs were identified. Similar homologous regions were absent in M. fructigena. For each of these regions, the flanking genes in the M. fructicola and M. laxa genomes were syntenic with the genes located at the terminal positions of two M. fructigena scaffolds, which were adjacent to non-syntenic genes annotated as active TEs. These consisted of RNA-directed DNA polymerase (RdDP) linked to a non-long terminal repeat (non-LTR)/Tad1-2_ACap element, which showed the highest GC content (>62%) of the codon wobble positions, suggesting probable horizontal gene transfer, and of hypothetical proteins linked to DNA/Mariner transposons or LTR/retrotransposons (Supplementary Figure 2 and Supplementary Table 2).

The rates of coding sequence evolution in syntenic regions were measured by synonymous changes (Ks) (Figure 3). Overall, one distinct peak was detected on the Ks histogram plots, which represent orthologous syntenic regions that follow taxa divergence during speciation. The Ks peak of Monilinia intra-genus comparisons was characterized by the prevalence of the blue-purple color, which suggested that there were no considerable changes in mutation rates at the genome-wide level. The mean Ks value between M. fructigena/M. laxa (Ks_mean = 0.273) was lower than both of those for M. laxa/M. fructicola (Ks_mean = 0.346) and M. fructigena/M. fructicola (Ks_mean = 0.442). Although the coding gene sequences between the three Monilinia genomes were very well conserved, these data support the hypothesis that M. fructigena is closer to M. laxa than to M. fructicola. The inter-genus comparison showed a prevalence of the green-cyan color, which emphasizes higher numbers of synonymous changes than older evolutionary changes. In particular, the highest Ks values in pairwise comparisons of the Monilinia genomes with B. cinerea and S. sclerotiorum were detected for M. fructigena (Ks_mean = 1.470, 1.023, respectively) and M. fructicola (Ks_mean = 1.356, 1.011, respectively), while the lowest were for M. laxa (Ks_mean = 1.058, 0.931, respectively) (Figure 3). These results suggest that the M. laxa genome is genetically closer to those of B. cinerea and S. sclerotiorum, as compared to the M. fructigena and M. fructicola genomes. An analysis of gene distribution in classes of Ks values (Figure 4) showed that homologous gene pairs with Ks ≤ 0.2 were more than 60% in the M. laxa/M. fructigena comparison and were as low as 28.3% and 21.0% in comparisons between M. fructicola and M. laxa or M. fructigena, respectively. This distribution of Ks values suggested a more ancient separation of M. fructicola from the other two Monilinia species. The frequency of homologous gene pairs with Ks ≤ 0.2 was lower than 1% in comparisons between each Monilinia species and B. cinerea or S. sclerotiorum.

The highest proportions of genes with a syntenic depth of 1 (indicative of a single orthologous genomic region) ranged from 80 to 95% in all of the pairwise comparisons among the tested Monilinia genomes, and from 57 to 73% in the comparisons with B. cinerea and S. sclerotiorum, respectively. A syntenic depth > 1 was more frequently detected (up to 17–18%) in the pairwise comparisons with M. fructigena used as the reference genome (Table 2).

For the Monilinia genome sequences, approximately 2.85 Mb (6.6%) of the 43.13 Mb of M. fructigena, 2.70 Mb (6.1%) of the 44.05 Mb of M. fructicola and 3.10 Mb (7.2%) of the 42.81 Mb of M. laxa consisted of TEs. Lower frequencies were detected in B. cinerea (2.19 Mb, 5.1% of 42.81 Mb) and S. sclerotiorum (2.3 Mb, 6.0% of 38.90 Mb) (Figure 5A and Supplementary Table 3). A strong heterogeneity was detected in the TE typology. Some TEs (e.g., LTR/Gypsy-2_PaPe-I or DNA/Mariner1_AO) occurred in large quantities and sizes in M. fructicola and M. laxa, but were absent or present as low numbers of small fragments (<500 bp) in M. fructigena, B. cinerea and S. sclerotiorum (Supplementary Table 4).

The Class I TEs was the most represented for all the species tested. In the Monilinia strains, the frequency ranged from 3.8% in M. fructicola to 5.3% in M. laxa and M. fructigena, and it was 3.7% and 3.3% of the genomes of B. cinerea and S. sclerotiorum, respectively. Six superfamilies of non-LTR retrotransposons were identified. The highest percentage was detected in the M. fructigena genome, due to the abundance of non-LTR/Tad1 (2.3%) compared to the other genomes (0.1–0.9%). LTR retrotransposons were mainly represented by the ubiquitous LTR/Gypsy superfamily followed by the LTR/Copia superfamily, and were predominant in M. laxa, M. fructicola, B. cinerea and S. sclerotiorum. LTR/Gypsy occupied 4.1% of the M. laxa genome and from 1.7% to 2.8% of the others. Twenty superfamilies of Class II DNA TEs were detected in all of the analyzed genomes, with very few exceptions. Overall, the M. fructicola genome showed the highest proportion of DNA TEs (2.4%), like S. sclerotiorum (2.7%). DNA/Mariner and DNA/hAT were the most represented superfamilies in M. fructicola, M. laxa, B. cinerea and S. sclerotiorum, while lower frequencies of these elements were detected in the M. fructigena genome (Figure 5A and Supplementary Table 3).

Monilinia fructicola showed the highest frequencies (TE counts/Mb genome) of Class I retrotransposons and Class II DNA TEs in the genome, although the total fraction (Mb) of the genome occupied by TEs was lower than in the other strains. This suggests a greater fragmentation of the annotated TEs in this genome.

To investigate TE dynamics, the distribution of TEs in the genome of each species were assessed based on their K-values (Figure 5B). The TE landscape in M. fructigena showed a burst of TE insertions where the highest accumulation occurred more recently (90% of TEs with K-values < 4) compared with M. laxa and M. fructicola, which showed one major peak corresponding to K-values of 15–20. In M. laxa, after the main burst, all of the transposon categories were contracted to lower levels, whereas the Class II DNA TEs that were more consistently detected in M. fructicola. B. cinerea and S. sclerotiorum genomes showed a recent TE burst that was similar, although to a lesser extent, to that in M. fructigena. The two genomes showed both an older burst (K-value = 34–42) mostly involving LTR and DNA TEs in B. cinerea and non-LTR elements in S. sclerotiorum. Moreover, a younger burst of TE insertion events (K-value = 18–20) that mainly involving LTR elements was observed in S. sclerotiorum.

The prediction of potentially active or degenerated TEs revealed that non-autonomous TEs and TE remnants (<200 bp in size) massively overlapped with genomic regions annotated as coding genes. TEs < 500 bp, however, had no detectable TE-related domains, whereas differences in the levels of TEs ≥ 500 bp with or without functional conserved domains were observed among the genomes analyzed (Figure 5C). The complete list of the TE types, the relative functional domains and the annotated genes is given in Supplementary Table 4. Active TEs with ORFs coding for proteins required for transposition were rarely detected in M. laxa and M. fructicola, while they were more frequent in the other genomes, with the highest proportion in M. fructigena, followed by S. sclerotiorum and B. cinerea. None of the 938 and 1,141 TEs ≥ 500 bp (up to ∼10,000 bp in size) that were identified in M. fructicola and M. laxa, respectively, overlapped with annotated gene sequences. Only 9 TEs in M. fructicola and 39 in M. laxa contained TE-related domains that mainly belonged to LTR/Gypsy. On the contrary, 641 of 1,031 TEs ≥ 500 bp identified in M. fructigena overlapped with 739 genes, with 86% of these showing more than 50% overlap. Among 633 TEs with at least one TE-related domain identified in M. fructigena, 562 were in-gene TEs and 96 were not in-gene TEs. In a similar way, 290 of 719 TEs ≥ 500 bp in S. sclerotiorum overlapped with 325 genes, with 73% of them showing more than 50% overlap; among 483 TEs carrying at least one TE-related domain, 284 were in-gene TEs and 199 were not in-gene TEs. Twenty-one of 604 TEs ≥ 500 bp in B. cinerea overlapped with 22 genes, with the 80% of them showing more than 50% overlap; among 110 carrying at least one TE-related domain, 9 were in-gene TEs and 106 were not in-gene TEs. In the M. fructigena, S. sclerotiorum, and B. cinerea genomes, LTR/Gypsy elements overlapped with genes crucial for retrotransposition; i.e., structural (Gag), structural-catalytic polyprotein (Gag-Pol), derivate reverse transcriptase, integrase and retrotransposon nucleocapsid. The LTR/Copia overlapped frequently with genes coding for ribonuclease H, retrotransposon hobase, Gag, Gag-Pol, reverse transcriptase and integrase. The non-LTR/Tad1 in the M. fructigena genome largely overlapped genes coding for RNA-directed DNA polymerase from transposon X-element, while in S. sclerotiorum and B. cinerea, they mainly overlapped with hypothetical proteins.

Only small proportions of the three Monilinia genomes were RIP positive, with RIP product indices > 1.1, RIP substrate indices < 0.75 and RIP composite index ≥ 0. The extent to which RIP mutations were recorded varied slightly among the genomes, with ranking as follows: M. fructicola > M. laxa > M. fructigena. These were, in any case, higher compared to B. cinerea and S. sclerotiorum (Table 3). In M. fructicola, 2.3% of the entire genome assembly contained RIP mutations with proportions varying from 1.5% to 3.4% in 16 large scaffolds (VICG01000001-VICG01000016), which appear to represent the core chromosomes, and from 9.3% to 10.4% in the remaining small scaffolds (VICG01000017-VICG01000019) (Supplementary Table 5). RIP activity was also recorded in the scaffold VICG01000020, which is made up of repeat units of ribosomal DNA, and where it was limited to intergenic spacers. Fewer windows suggestive of RIP (1.0%) were detected in the M. laxa genome, which ranged from 0.2% to 1.3% in large scaffolds (VIGI01000001-VIGI01000017). Even fewer RIP regions were recorded in the M. fructigena genome, with a small number of affected windows that ranged from 0.0% to 0.6% in most of the assembled scaffolds, corresponding to 0.4% of the whole genome sequences. In the B. cinerea and S. sclerotiorum genomes, RIP-affected windows were ≤ 0.1% in the core chromosomes, and 0.5% to 0.6% in the two B. cinerea mini-chromosomes. RIP mutations were more obviously co-localized with TE-rich regions in M. fructicola and M. laxa than in M. fructigena (Figure 6). Reduced GC content (17.1% in M. fructicola, 18.8% in M. laxa, 31.9% in M. fructigena) was recorded in genomic regions suggestive of RIP, compared to that of the entire genome (40–42%). Gene density was 2.3 and 1.9 genes per 100 kb in RIP-affected windows of M. fructicola and M. laxa, while the average contents of genes in the genomes were 27.5 and 26.1 genes per 100 kb, respectively (Table 3). Genes in these regions (121 in M. fructicola, 29 in M. laxa) were largely without functional annotations (43–55%). In M. fructigena, the genomic gene density (24.4 genes per 100 kb, on average) increased up to 41.8 in RIP-affected regions (Table 3). The genes contained here (n = 204) were mostly involved in biological processes (e.g., endonucleolytic RNA phosphodiester bond hydrolysis [GO:0090502] and DNA integration [GO:0015074]) and molecular functions (e.g., nucleic acid binding [GO:0003676] and RNA-DNA hybrid ribonuclease activity [GO:0004523]), with prevalent nuclear localization (GO:0005634) (data not shown).

Few LRARs were identified in the genomes of M. fructicola (n = 9) and M. laxa (n = 2), with none identified in the M. fructigena, B. cinerea, and S. sclerotiorum genomes. The LRARs presented low GC contents and did not contain genes, with very few exceptions. For six of the LRARs predicted in M. fructicola and one of those predicted in M. laxa, fine-scale analysis of their flanking sequences revealed that they were included in large RIP-affected regions (14,000–103,000 bp in size) with marked variation in RIP index values. All of these regions were rich in TEs (mainly DNA/Mariner, DNA/hAT, LTR/Gipsy), and represented synteny breakpoint regions (Supplementary Figure 3).

Di-nucleotide frequency analysis of the RIP-affected regions identified in the Monilinia genomes showed the expected pattern of classical CpA→TpA type RIP mutations and the occurrence of CpG→TpA and CpC→TpT double mutations. Moreover, there was depletion of GpC, and to a lesser extent of TpC, that resulted from changes in the second nucleotides, as already reported for B. cinerea and S. sclerotiorum (Clutterbuck, 2011) (Supplementary Figure 4).

The Monilinia genomes contained genes that encoded homologs (with conserved domains) of the RID and DIM-2 enzymes, which are essential for RIP and RIP-associated methylation in N. crassa, and the associated co-factors HP-1 and DIM-3, -5, -8, -9, but not DIM-7 (Table 4). RNA-seq data provided evidence of transcription for all of these genes detected (data not shown).

In M. fructicola, there were 759 secreted proteins, with 753 in M. laxa and 532 in M. fructigena, corresponding to 5% to 6% of the total annotated proteins. These data agree with those reported for S. sclerotiorum (Guyon et al., 2014), B. cinerea (Amselem et al., 2011; van Kan, 2017) and other Botrytis species (Valero-Jiménez et al., 2019). The prediction results of EffectorP 1.0 and EffectorP 2.0 and the number of proteins of each Monilinia strain labeled as effectors by both prediction tools are given in Table 5. The final numbers of the predicted effectors ranged from 56 (M. fructigena) to 68 (M. fructicola), which represented ∼0.5% of the total proteome of each strain, and from 8.1% (M. laxa) to 10.5% (M. fructicola) of the secretome, in agreement with the findings of Marcet-Houben et al. (2021). Orthologous clusters were identified among putative effectors identified in each Monilinia genome (Supplementary Figure 5). A set of 15 effector candidates had orthologs in all three genomes, while 12 were common to M. fructicola and M. laxa, 9 were common to M. laxa and M. fructigena, and 5 were common to M. fructicola and M. fructigena. Several effectors of M. fructicola (n = 36), M. laxa (n = 25), and M. fructigena (n = 27) did not have orthologs among the candidate effectors predicted in the other genomes. A list of all of the orthologous clusters and their conserved functional domains is given in Supplementary Table 6.

Among nine candidate effectors common to the five genomes, one was a Nep1-like protein (NLP) with a necrosis-inducing protein NPP1 domain (PF05630). The other ones comprised, for instance: a protein carrying a peptidase inhibitor I9 and peptidase S8 domains (G3DSA:3.30.70.80); a Rare lipoprotein A (RlpA) with a conserved region with a double-psi β-barrel (DPBB) fold (PF03330); a protein with a CVNH (cyanovirin-N homology) domain (PF08881); and a member of the PTHR40845 panther superfamily associated to putative effectors of fungal plant pathogens.

A pectate lyase (PF03211) was identified in the three Monilinia genomes and in B. cinerea, while a hydrophobic surface binding protein A (HsbA; PF12296) had orthologs in the three Monilinia genomes and in S. sclerotiorum. A family 2 peroxidase (PF01328) and a FK506-binding protein (FKBP; PF00254) had orthologs in all three of the Monilinia genomes. The remaining orthologs identified in more than one genome included a protein containing a MD-2-related lipid-recognition (ML) domain (PF02221), which has been implicated in the recognition of pathogen-related products like lipopolysaccharides, a putative Egh16-like virulence factor (PF11327), a protein with a cerato-platanin domain (PF07249), and a fungal hydrophobin (PF06766).

Among candidate effectors predicted only in the Monilinia genomes, one contained a phosphatidylethanolamine-binding protein (PEBP; PF01161), and one carried the tuberculosis necrotizing toxin (TNT) domain from Mycobacterium tuberculosis (PF14021), which induces necrosis of infected cells. Effectors predicted in M. fructicola but not in M. laxa and M. fructigena included, e.g., a FAS1 (fasciclin-like) domain protein (SSF82153) with orthologs in B. cinerea and S. sclerotinia; a cerato-platanin-like protein (PF07249); and a protein with chitin-binding, type 1 (PF00187) and fungal cellulose binding (PF00734) domains. A candidate effector identified only in M. laxa and M. fructigena was associated to the heat-labile enterotoxin family (PF01375), while an Egh16-like virulence factor (PF11327) was predicted only in M. fructigena, as well as a putative GPI-anchored protein (PTHR39599), a [PSI+] induction protein 2 (PTHR40018), and cupredoxin (SSF49503, G3DSA:2.60.40.420).

CAZyme prediction for the genomes revealed 183 enzymes in M. fructicola, 200 in M. laxa, 134 in M. fructigena, 271 in B. cinerea and 224 in S. sclerotinia, corresponding to 24.1%, 26.6%, 25.2%, 25.0%, and 29.5% of secreted proteins, respectively. These data agree with previous reports (e.g., Amselem et al., 2011; Allan et al., 2019; Marcet-Houben et al., 2021; Vilanova et al., 2021). Overall, 87 CAZyme families were identified in the predicted secretomes, including 51 glycoside hydrolases (GHs), 11 auxiliary activities (AAs), 11 carbohydrate esterases (CEs), 8 glycosyltransferases (GTs), and 3 polysaccharide lyases (PLs), as well as 9 carbohydrate-binding modules (CBMs), which can potentiate CAZyme activities (Supplementary Table 7). In all of the genomes analyzed, the largest numbers of predicted CAZymes were in the GH28 family of polygalacturonases, followed by AA3 in the Monilinia genomes, and by GH16 in B. cinerea and S. sclerotiorum.

The distributions of the CAZyme classes were similar in the five genomes, with very little difference for GHs, AAs and GTs (Figure 7). CAZymes were grouped into ortholog clusters (Supplementary Figure 6 and Supplementary Table 7). A chitinase of family GH18 and two glycosyltransferases (GT15, GT31) detected in the Monilinia genomes did not have orthologs in either B. cinerea or S. sclerotiorum. Members of the GH28 and GT71 families and a β-(1,3)-glucanosyltransglycosylase of the GH72 family were predicted only in M. laxa. CAZymes unique for M. laxa and M. fructicola included, for instance: an α-1,4-galactosaminogalactan hydrolase (GH135) homologous to cell surface spherulin 4-like protein; a glycosyl hydrolase family 43, involved in degradation of plant cell-wall arabinans; and choline dehydrogenase or related flavoprotein/GMC oxidoreductase of the subfamily AA3-2, involved in lignocellulose degradation. Two AA3-2 CAZymes were common to M. laxa, M. fructicola and B. cinerea, with no orthologs in M. fructigena and S. sclerotiorum, as also for mannosyltransferase (GT32), α-mannosidase (GH47), β-glucosidase (GH3), and a carbohydrate esterase family 1 (CE1). An AA3-2 member and a polygalacturonase protein (GH28) were also common to M. laxa, M. fructicola and S. sclerotiorum, but not to M. fructigena and B. cinerea. CAZymes common to M. laxa and B. cinerea but not to M. fructigena included, for instance, a cell-wall glycosyl hydrolase protein (GH105), while those common to M. laxa, B. cinerea and S. sclerotinia included several GHs, a mannosyltransferase, two laccases of the AA1-3 subfamily and an alginate lyase protein (PL7). CAZymes common to M. fructicola, B. cinerea and/or S. sclerotiorum with no orthologs in M. fructigena and M. laxa included endo-1,4-β-mannosidase (GH5), α-1,3-glucanase (GH71), a member of the GH13 family, a mannosyltransferase (GT32), and a gdsl lipase acylhydrolase family protein (CE16). A mannosidase of the GH92 family (FCW CAZy) and AA enzymes were uniquely found in M. fructigena, while a fungal cellulose binding domain (CE15 + CBM1) was unique for M. fructigena and M. laxa.

A total of 18 putative BGCs (4 NRPS, 3 NRPS-like, 10 T1PKS, 1 Terpene) was found in the Monilinia genomes (Supplementary Table 8): 17 BGCs on 10 scaffolds of the M. fructicola genome or on 12 scaffolds of the M. laxa genome, and 12 BGCs on 11 scaffolds of the M. fructigena genome. Generally, only one BGC per scaffold was detected, although up to 4 BGCs were on the M. fructicola scaffold VICG01000003.1 and on the M. laxa scaffold VIGI01000002.1.

The cluster coded as BGC-1 contained a NRPS Sid2-like gene, putatively involved in siderophore production because of its similarity to the ferrichrome biosynthetic gene cluster of Aspergillus oryzae (64–75% amino acid sequence identity). The core gene of BGC-2 showed 32 to 35% identity to the NRPS of the Aspergillus nidulans BGC that produces the bioactive lipopeptides aspercryptins. BGC-1 and BGC-2 were found in the three Monilinia genomes, and the genes included showed co-regulation in M. fructigena. BGC-3 was common to M. fructicola and M. laxa, but not to M. fructigena, and contained a core gene homologous (39% identity) to the NRPS involved in biosynthesis of destruxin A (a bioactive cyclohexadepsipeptide) in Metarhizium robertsii. Two NRPS-like clusters common to the three Monilinia genomes showed homology to BGCs involved in the biosynthesis of mycotoxins, including aspirochlorine in A. oryzae (BGC-4; 42–44% identity) and fusaric acid in Fusarium verticillioides (BGC-5; 34% identity). BGC-4 included 13 to 15 co-regulated genes in each of the three Monilinia genomes. Two further NRPS-like BGCs that were common to the three Monilinia genomes were putatively involved in the biosynthesis of antibiotic compounds, as they had similarity to the myxalamid BGC from Stigmatella aurantica (BGC-6; 29–36% identity) and the fusarielin H from Fusarium graminearum (BGC-7; 25–32% identity).

Among T1PKS-type BGCs, two clusters were putatively involved in DHN-melanin biosynthesis: BGC-8, containing a homolog to BcPKS12 of B. cinerea (88% identity), and BGC-9, containing a homolog to BcPKS13 of B. cinerea (87–88% identity). In B. cinerea, similar to other Leotiomycetes, the gene encoding PKS12 was physically linked with a gene encoding the transcription factor sclerotial melanin regulator 1 (SMR1) (Schumacher, 2016). In the present analysis, this pks12-smr1 module appeared to be conserved within a cluster of co-regulated genes in M. fructicola and M. laxa, whereas in M. fructigena, the PKS12-homolog was not detected in a cluster due to a genome rearrangement that occurred in the region surrounding the pks12 gene that resulted in a loss of synteny with M. fructicola and M. laxa (Figure 8). In more detail, the smr1 homolog in the M. fructigena genome was found about 60 kb upstream from pks12, probably due to inversion/translocation events and insertions of active TEs (i.e., Tad1-2_ACap, Gypsy-41_BG; Boty_LTR) in the genomic region in between smr1 and pks12 and in the downstream region.

Three T1PKS-like BGCs were identified in all of the three Monilinia genomes, homologous to BGCs involved in the biosynthesis of toxicants, like the xanthone dimer with antibacterial activity neosartorin in Aspergillus fumigatus (BGC-10; 67% identity), the insecticidal furanocoumarin neosporin A in N. crassa (BGC-11; 41–47% identity), and the maleidride with herbicidal activity cornexitin (BGC-12; 34–42% identity). BGC-10 clustered genes shown to be co-regulated in M. fructicola, and BGC-11 and BGC-12 clustered genes shown to be co-regulated in M. laxa. BGC-13 and BGC-14, which are co-regulated in M. fructicola and M. laxa and absent in M. fructigena, were putatively involved in the biosynthesis of phytotoxins. BGC-13 showed homology to the botcinic acid BGC from B. cinerea (82–84% identity). BGC-14 was found in two consecutive copies on the same scaffold and showed homology to the solanapyrone D BGC from Alternaria solani (49–51% identity). Two other T1PKS-like BGCs were identified in M. fructicola and M. laxa but not in M. fructigena: BGC-15 homologous to the brefeldin A BGC from Penicillium brefeldianum (49–51% identity), and BGC-16 homologous to the grayanic acid BGC from the lichen Cladonia grayi (33–37% identity). The PKS in BGC-17, which was identified only in M. fructigena, was homologous to the fumonisin BGC from Fusarium oxysporum (42% identity).

A single terpene biosynthetic gene cluster (BGC-18) was identified in the Monilinia genomes. It showed homology to the squalestatin S1 BGC in Aspergillus sp. (57–68% identity), which is an inhibitor of mammalian and fungal squalene synthases.