De novo assembly and comparative transcriptome analysis of Monilinia fructicola, Monilinia laxa and Monilinia fructigena, the causal agents of brown rot on stone fruits

Background Brown rots are important fungal diseases of stone and pome fruits. They are caused by several Monilinia species but M. fructicola, M. laxa and M. fructigena are the most common all over the world. Although they have been intensively studied, the availability of genomic and transcriptomic data in public databases is still scant. We sequenced, assembled and annotated the transcriptomes of the three pathogens using mRNA from germinating conidia and actively growing mycelia of two isolates of opposite mating types per each species for comparative transcriptome analyses. Results Illumina sequencing was used to generate about 70 million of paired-end reads per species, that were de novo assembled in 33,861 contigs for M. fructicola, 31,103 for M. laxa and 28,890 for M. fructigena. Approximately, 50% of the assembled contigs had significant hits when blasted against the NCBI non-redundant protein database and top-hits results were represented by Botrytis cinerea, Sclerotinia sclerotiorum and Sclerotinia borealis proteins. More than 90% of the obtained sequences were complete, the percentage of duplications was always less than 14% and fragmented and missing transcripts less than 5%. Orthologous transcripts were identified by tBLASTn analysis using the B. cinerea proteome as reference. Comparative transcriptome analyses revealed 65 transcripts over-expressed (FC ≥ 8 and FDR ≤ 0.05) or unique in M. fructicola, 30 in M. laxa and 31 in M. fructigena. Transcripts were involved in processes affecting fungal development, diversity and host-pathogen interactions, such as plant cell wall-degrading and detoxifying enzymes, zinc finger transcription factors, MFS transporters, cell surface proteins, key enzymes in biosynthesis and metabolism of antibiotics and toxins, and transposable elements. Conclusions This is the first large-scale reconstruction and annotation of the complete transcriptomes of M. fructicola, M. laxa and M. fructigena and the first comparative transcriptome analysis among the three pathogens revealing differentially expressed genes with potential important roles in metabolic and physiological processes related to fungal morphogenesis and development, diversity and pathogenesis which need further investigations. We believe that the data obtained represent a cornerstone for research aimed at improving knowledge on the population biology, physiology and plant-pathogen interactions of these important phytopathogenic fungi. Electronic supplementary material The online version of this article (10.1186/s12864-018-4817-4) contains supplementary material, which is available to authorized users.


Background
Brown rots are major diseases of stone and pome fruits that can cause severe yield losses during both field production and post-harvest processing [1]. Several species of Monilinia genus are responsible for the diseases, but Monilinia fructicola (G. Winter) Honey (MFRC), Monilinia laxa (Aderh. & Ruhland) Honey (MLAX) and Monilinia fructigena Honey (MFRG) are the most common species all over the world [1][2][3]. MLAX and MFRG have been the only ones present in Europe until 2001 when MFRC has been introduced, rapidly spread and become prevalent on the former indigenous species [4,5]. Some phenotypic differences in fitness parameters such as growth rate in response to temperature, fungicide sensitivity and virulence have been recorded among the three species [5][6][7] but the reasons underlying its successfulness have not yet well clarified.
Although the three fungi have been deeply studied, the availability of genomic and transcriptomic data in public databases is still scant. There are only few published studies on the genetic variation in Monilinia species using different molecular markers, i.e. internal transcribed spacer (ITS) [8,9], random amplified polymorphic DNA (RAPD) [6], simple sequence repeat (SSRs) and inter SSR (ISSR) [10,11], and the characterization of specific genes, such as genes responsible for fungicide resistance [12] or pathogenicity [13].
In this study, we used Illumina sequencing of mRNA from conidia and mycelia of MFRC, MLAX and MFRG to obtain a de novo Trinity-based assembly of their transcriptomes. Orthologous transcripts were identified and compared to detect those apparently unique or differentially expressed in each species.

Samples
Two isolates of opposite mating type per each species, determined by PCR-based mating type assays [14], collected from naturally infected fruits sampled in orchards located in Italy, were used in this study ( Table 1). The isolates were grown under different conditions to obtain comprehensive transcriptomes: i) mycelium grown at 21 ± 1°C on cellophane disks overlaid on potato dextrose agar (infusion from 200 g peeled and sliced potatoes kept at 60°C for 1 h, 20 g dextrose, adjusted at pH 6.5, 20 g agar Oxoid No. 3, per litre) in the dark for 4 days; ii) mycelium grown as above but in the dark for 2 days and then exposed 12 h per day to a combination of 2 daylight (Osram L36 W/20) and 2 near-UV (Osram, L36/73) lamps for 2 days; iii) conidia (1 × 10 5 mL − 1 ) germinating in malt extract medium (20 g malt extract Oxoid, per litre) after 14 h at 24 ± 1°C in darkness under shaking (120 rpm). Mfrg344 C ¥ D = mycelium grown in the dark for 4 days; L = mycelium grown in the dark for 2 days and then exposed to light for 2 days; C = germinating conidia

RNA extraction, library preparation and sequencing
Total RNA was extracted from a total of 18 samples (6 per species) made up by 100 mg of mycelium or germinated conidia with the RNeasy Plant Mini Kit (Qiagen, Milan, Italy), following the manufacturer's protocol. cDNA libraries were prepared from 4 μg of total RNA using the TruSeq RNA Sample Preparation Kit v2 (Illumina, Inc., San Diego, CA, USA) and validated according to Illumina's low-throughput protocol. The protocol was adjusted to obtain an average library size of about 400 bp (insert length 130-290 bp), by reducing RNA fragmentation time to 2 min at 94°C. Sequencing was carried out on an Illumina HiScanSQ platform using TruSeq SBS kit v3 (Illumina, Inc.) to obtain paired-end reads, 92 nt in length. RNA and DNA quantity and quality were determined with a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA) and a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). After removing indexed adapters, reads from each library were filtered for quality score (QS ≥ 30) using CASAVA v1.8 software (Illumina, Inc.).

Sequencing read quality and trimming
Reads were analysed for quality statistics, nucleotide distribution and redundancy using FastX-tools (http:// hannonlab.cshl.edu/fastx_toolkit), and trimmed with Trimmomatic 0.36 (http://www.usadellab.org/cms/ index.php?page=trimmomatic) setting the parameters as follows: i) LEADING and TRAILING = 3, removing bases from the two ends of the reads if below a threshold quality of 3; ii) SLIDING WINDOW = 4:2, cutting the reads when the average quality within the window composed of 4 bases falls below a threshold equal to 2; iii) MINLEN = 50, removing the reads shorter than 50 bp [15].

Transcriptome de novo assembly
Trinity software v.2.1.1 (https://github.com/trinityrnaseq/ trinityrnaseq/wiki) was used for the de novo assembly of the transcriptomes using together sequencing data from the six libraries (2 isolates and 3 growing conditions) from each species (Table 1). Default assembly parameters of Trinity were used, with the addition of the "-jaccard_clip" function because a high gene density with overlapping of UnTranslated Region (UTR) was expected [16].

Functional annotation
The annotation of the putative transcripts obtained from each assembly was performed using local BLAST+ 2.3.0 (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/) [17] and Blast2GO PRO v4.0.7 [18]. BLASTx analysis was carried out by searching against the NCBI non-redundant protein database (downloaded 22 November 2016) and setting E-value cut off at 10 − 3 . Blast2GO PRO was used to predict Gene Ontology (GO) terms, to assign the assembled sequences to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and to analyse protein domains using the InterProScan tool. Blast2GO annotation search was conducted with 10 − 6 as the E-value hit filter, 55 as the annotation cut-off and 5 as the GO weight; no HSP-hit coverage cut-off was considered. Moreover, all GO terms retrieved via InterProScan analysis were added and used to validate GO annotations; further enhancement of GO terms was conducted with Annotation Expander (ANNEX).
The Trinity script analyze_blastPlus_topHit_coverage.pl (https://github.com/trinityrnaseq/trinityrnaseq/wiki/Counting-Full-Length-Trinity-Transcripts) was launched to determine the number of full-length or nearly full-length transcripts using BLAST+ with an E-value cut off 10 − 20 and a protein database built from the proteome of B. cinerea (ASM83294v1; http://fungi.ensembl.org/Botrytis_cinerea/Info/Index, downloaded 10 April 2016), used as the closest related organism.

Comparative transcriptome analysis
tBLASTn search was used to identify orthologous transcripts in the three Monilinia species. In detail, Botrytis cinerea B05.10 proteins (ASM83294v1) integrated with the mitochondrial proteins (ftp://ftp.broadinstitute.org/ pub/annotation/fungi/botrytis_cinerea/genomes/botry-tis_cinerea_b05.10_mito, downloaded 10 April 2016) were queried against the assembled transcriptomes of MFRC, MLAX and MFRG, to identify homologous sequences. tBLASTn was run using CLC Genomics Workbench with a threshold E-value< 10 − 3 . Reads were then mapped on the selected transcripts and the unmapped reads re-aligned on the complete Trinity transcriptomes, to retrieve Monilinia transcripts with not homologs in the B. cinerea proteome; all transcripts with counted reads > 50 and significant matches to fungal proteins in BLASTx search were included for further analysis. Homology between the putative Monilinia orthologs was assessed through BLASTn pairwise alignments in all possible combinations (MFRC vs MLAX; MFRC vs MFRG; and MFRG vs MLAX) with the threshold parameter of E-value< 10 − 10 .
Reads from each isolate per species were mapped on the identified set of orthologous transcripts using CLC Genomics Workbench with the previously reported alignment parameter setting. The values of gene expression were measured in reads per Kb of transcript per million mapped reads (RPKM). The data from the three growing conditions for each isolate were pooled and used as two biological replicates per each species. Fold Change (FC) was calculated comparing the RPKM values in all pairwise combinations between the three fungal species. False Discovery Rate (FDR) was determined using the edgeR Bioconductor package [20] and transcripts with FC > |8| and FDR ≤ 0.05 in both the comparisons of each species with the other two were considered as differentially expressed transcripts (DETs) and submitted to functional analysis. For an accurate comparison among the species, we removed from the identified DET sets: i) incorrectly associated transcripts revealed by pairwise BLASTn alignments, including those with doubtful RPKM values due to multiple isoforms; and ii) chimeric transcripts identified by BLASTx, i.e. assembled sequences derived from two or more adjacent transcribed genes [21].
WEGO was used to perform functional classification of Trinity unigenes based on the GO annotation and compare the overall distribution of gene functions in the three species using the Pearson Chi-Square test [22]. All unigenes from each transcriptome were also submitted to search against the EggNog (Evolutionary genealogy of genes: Non-supervised Orthologous Groups) database, integrated in the Blast2GO pipeline, to predict and classify gene functions based on sequence similarity within clusters of orthologous groups (OGs) [23].

Sequencing and transcriptome de novo assembly
Illumina mRNA sequencing from germinating conidia and actively growing mycelia of two isolates of opposite mating type per each species generated a total of 19.5 Gb, more than 6 Gb per species corresponding to about 70 million of paired-end reads (Additional file 1: Table S1).
After trimming of low quality reads, less than 1% of input paired-end reads were discarded, and the remaining reads used for de novo assembly. The most relevant data are in Table 2. Overall, Trinity assembly generated 33,861 contigs for MFRC, 31,103 for MLAX and 28,890 for MFRG corresponding to putative transcripts, including isoforms. The number of unigenes assembled by Trinity exceeded the expected number of protein coding genes, likely due to fragmented transcripts, i.e. unigenes representing the same transcript that could not be assembled because containing a gap.

Quality assessment of de novo assembly
More than 80% of the paired-end reads used for the assembly successfully mapped to the respective assembled transcriptome used as reference (Additional file 1: Data obtained from transcriptome completeness analysis by BUSCO, based only on conserved fungal orthologs, showed that over 90% of the assembled transcripts for the three Monilinia species were complete and that the percentage of duplicated transcripts was always less than 14%. Less than 10% of transcripts were fragmented or missing (Table 3).

Functional annotation
More than 50% of the assembled contigs had a significant hit in BLASTx search ( Table 4). The top-hits by species distribution analysis showed the highest similarities with B. cinerea (strains BcDW1, T4 or B05.10), Sclerotinia sclerotiorum (1980 UF-70) and Sclerotinia borealis (F-4157) (Fig. 1a). The E-value distribution of the top blast hits showed about 50% of hits with an E-value equal to zero and about 70% with E-values ranging from 0 to 1e-61 ( Fig. 1c). Sequences with no significant hits were mostly short fragments (Fig. 1b)  In the three transcriptomes, more than 80% of the unigenes with a significant hit were associated to at least one GO term, and the GO terms were functionally classified and plotted. The three transcriptomes were very similar in their profiles of unigene distribution in functional categories (Fig. 2). Many unigenes were assigned to cellular and metabolic processes and to biological processes related to morphogenesis, localization, pigmentation, development and growth, reproduction, response to stimulus, multi-organism and multicellular organismal processes. In the molecular function category, binding and catalytic activities represented the majority followed by transporters, structural molecules and enzymes, transcription and translation regulation. In the cellular component category, cell, cell part, organelle and macromolecular complex were prevalently represented.

Comparative transcriptome analysis
Ten-thousand four-hundred and fifty-two putative orthologous transcripts were identified among the three Monilinia species. In detail, 10,045 orthologs were identified by tBLASTn analysis using the proteome of B. cinerea as reference, whereas 407 transcripts were specific for at least one of the Monilinia species and had not homologous sequences in B. cinerea proteins.
Heat-map representation of expression profiles of the DETs from the two isolates per species under the different growing conditions confirmed the diversity among the species revealed by mapping pooled reads of each isolate and showed differences between germinating conidia and mycelium samples that are not analysed in the present paper ( Fig. 5).
Orthologues of the BcNRPS1, coded T8874, were fragmented in MLAX and over-expressed in MFRG (FC = 8.1-27.9) like T6684 (FC = 54.0-72.0), which encodes a protein sharing homology with a FAD-binding-domain containing protein of S. borealis (E-value = 0; id = 89%) and bifunctional solanapyrone synthases from Phialophora attae, Diplodia seriata and Penicillium chrysogenum (best E-value = 1e-137; id up to 50%). Three transcripts related to secondary metabolism (T10344, T10360 and T10358) were assembled only in MFRG. T10344 encodes a conserved hypothetical protein containing a phenylacetate-coenzyme A (CoA) ligase domain, with homologs in B. cinerea, S. borealis and several Penicillium and Aspergillus species (E-value = 0; id up to 85%). Two taurine catabolism dioxygenases TauD/TdfA (T10360 and T10358), unique in MFRG and having homologs in B. cinerea, produce sulphite and aminoacetyldehyde from the sulphur-containing amino acid taurine that is used by several bacteria and fungi for growth under sulphate starvation [25]. T10360 exhibits also significant homology (best E-value = 0; id up to 90%) with proteins of Penicillium, Aspergillus and other fungi containing a clavaminic acid synthetase (CAS)-like domain.

Nucleic acid modification and metabolism
A zinc ion binding protein (T9754) containing a WLM (WSS1-like metalloprotease) domain was over-expressed in MFRC (FC = 10.8-21.4); it is related to the DNA damage response proteins WSS1 involved in sister chromatid separation and segregation and shares homology with B. cinerea and S. borealis proteins (e-values = 0; id≥77%). In addition, a group II intron reverse transcriptase/maturase (T10103), with homologs in Neurospora intermedia (E-value = 1e-40; id = 32%) and other fungi but not Sclerotiniaceae, was uniquely assembled in MFRC.
A mitochondrial DNA polymerase type B protein (T91) was uniquely assembled in MLAX, as well as a ribonuclease III (RNAse III) protein (T10215).

Discussion
The complete transcriptomes of MFRC, MLAX and MFRG were reconstructed and annotated. Comparative analyses among orthologous transcripts revealed 65 transcripts over-expressed (FC≥8 and FDR≤0.05) or unique for MFRC, 30 in MLAX and 31 in MFRG.
DETs encoding hydrolytic or carbohydrate-active enzymes were identified; nine were GHs, each belonging to a different family according to the Carbohydrate-Active Enzymes database (CAZy) classification (http://www.cazy.org/Glycoside-Hydrolases.html), which represent the largest and most diverse family of biopolymer-degrading enzymes [26].
The GHs over-expressed in MFRC included a β-1,3-glucan binding protein that in Phanerochaete chrysosporium is involved in the degradation of laminarin [27]; an α-L-arabinofuranosidase which is homolog to the Magnaporthe oryzae MoAbfB required for full virulence and inducing host defence responses [28]; and a chitin-binding type 1 protein. Effector proteins with chitin binding activity are directly involved in pathogenic processes in Cladosporium fulvum, M. oryzae and Mycosphaerella graminicola [29] and along with chitinases are up-regulated during S. sclerotiorum infection on Brassica napus [30]. Pizzuolo et al. [6] compared isolates of MFRC, MLAX and MFRG for pectolytic activity and their isoenzyme patterns and reported no substantial differences in the production of pectolytic enzymes (i.e. pectin methylesterase, polygalacturonase and pectin lyase) useful to differentiate the three species. In our study, two enzymes involved in plant cell wall breakdown were highly over-expressed in MFRC: a pectate lyase and a galactose oxidase, which is a secreted enzyme well characterized in Fusarium spp. catalysing the oxidation of D-galactose and other primary alcohols to aldehydes with concomitant reduction of molecular oxygen to hydrogen peroxide [31]. Moreover, two GAL4 TFs, over-expressed in MFRC, are positive regulators for galactose-induced genes in S. cerevisiae and for genes encoding enzymes degrading plant cell-wall polysaccharides in pathogenic fungi. The GAL4 TF encoded by the AbPf2 (Alternaria brassicicola pathogenicity factor 2) gene is essential for the pathogenicity of the fungus [32]. Several α/β hydrolases, including proteins with carboxylesterase and lipase activities, and a putative amidohydrolase were over-expressed or uniquely assembled in MFRC. Amidohydrolases catalyse the hydrolysis of a wide range of substrates bearing amide or ester functional groups and include ß-lactamases, e.g. penicillin amidohydrolase involved in antibiotic metabolism [33].
GHs over-expressed in MFRG include a putative α-amylase (GH13) acting on α-glycosidic bonds and involved in lignocellulose saccharification processes [36] and a ß-mannosidase B-like protein (GH2), enzymes displaying different substrate specificity in Aspergillus spp. that likely reflect the diversity of their functions [37]. Two additional cell-wall degrading enzymes were identified, a carbohydrate esterase family 12 protein, acting on pectin, and a cellobiose dehydrogenase which plays a key role during pathogenesis in Sclerotium species [38] and in lignocellulose degradation by wood-rotting fungi [39].
A protease inhibitor containing a PEBP domain was over-expressed in MFRC. The yeast PEBP homologue TFS1 is supposed to be a bridge between cell signalling and intermediate metabolism in Saccharomyces cerevisiae [40]. Protease inhibitors are fungal effectors involved in plant-pathogen interactions; e.g., the inhibitor encoded by the C. fulvum avirulence gene Avr2 is secreted during tomato infection and acts on papain-like cysteine proteases involved in plant defence [41]. We also identified in MFRC an RHS repeat protein. This class of proteins includes secreted bacterial insecticidal and nematocidal toxins and intercellular signalling proteins mediating bacterial-host or bacterial-bacterial interactions (e.g., [42]) A different class of the same family includes the Ss-RHS1 secretory protein from S. sclerotiorum that is related to sclerotial development, appressoria differentiation and virulence [43]. A trypsin-like serine protease was over-expressed in MFRC while a secreted pepsin-like aspartic proteases was over-expressed in MLAX. Trypsin-related enzymes are supposed to play a specific role in pathogenicity since they are the major proteases produced by plant-pathogenic Verticillium spp. whereas there are no reports of trypsins being secreted by saprophytes [44]. Secreted aspartic proteases, commonly known as acid proteases, are endopeptidases involved in functions related to nutrition and pathogenesis [45].
Among DETs related to fungal morphogenesis and development, we found over-expressed in MFRC a developmental-specific protein homolog to the S. sclerotiorum sclerotium-specific protein (Ssp1) detected in all stages of sclerotium and apothecium development [46]. A phenotypic comparison among the three Monilinia species described in a recent study by Villarino et al. [7] showed that sclerotia and stromata differentiation, affecting fungal survival under unfavourable conditions, was the most important factor distinguishing MFRC from the other two species. Moreover, MFRC is described as able to produce apothecia from pseudosclerotia in mummified fruit under both natural and in vitro conditions [47]; further functional studies are needed to determine the role of Ssp1 homologous gene in pseudosclerotia development and evolution of Monilinia spp.. MFRC is also characterized by a faster colony growth and more abundant sporulation than MLAX and MFRG [5,7]. In this study, MFRC showed over-expression of a fatty acid oxygenase with homology to PpoA proteins which are cyclooxygenase-like enzymes responsible for oxylipin production and regulate the development of conidiophores and cleistothecia and the production of toxic metabolites and degradative enzymes affecting fungal development and plant-pathogen interactions [48]. Other DETs were over-expressed in MFRC: a MSS4-like protein playing a role in cell polarization, polar tip extension and hyphal growth in Neurospora crassa [49]; a fatty acid elongase (FEN1) involved in sphingolipid biosynthesis and in modulation of resistance to amphothericin B in yeasts [50]; and the SED1 protein which is believed involved in resistance to lytic enzymes. Besides, two MFRC cell surface proteins are homologous to the Blumeria graminis f.sp. hordei gEgh16 and Egh16H, M. grisea GAS1 and MAS3, and Colletotrichum CAS1, all virulence factors essential for the establishment of plant-pathogen interactions [51].
DETs over-expressed in MLAX include a GPI anchored protein, important for cell wall synthesis and integrity in N. crassa [52] and implicated in virulence and in planta proliferation in Fusarium graminearum [53]; a regulator of G protein signalling which play key roles in upstream regulation of fundamental processes in fungi, including vegetative growth, sporulation, mycotoxin and pigment production, pathogenicity and mating [54]; two ankyrin repeat proteins, mediating protein-protein interactions; and a putative ATPase of the AAA family, a large and diverse group of enzymes inducing conformational changes in a wide range of proteins associated with cell-cycle regulation, proteolysis, organelle biogenesis and intracellular transport [55].
A Ca2+/calmodulin-dependent (CaMK) and a serine/ threonine protein kinases were over-expressed in MFRC while a protein kinase-like over-expressed in MLAX shares homology with phosphotransferases responsible for antibiotic resistance in Aspergillus and Penicillium spp. [56].
A Zn finger C2H2 TF was over-expressed in MFRG. Such TFs are involved in pathogenicity, catabolite repression, acetamide regulation, differentiation of fruiting body, as well as stress responses and multidrug resistance in yeasts and human fungal pathogens. In B. cinerea they are involved in phytotoxin biosynthesis, secondary metabolism, carbohydrate metabolism, transport, virulence and detoxification mechanisms [57,58]. MFRG also showed over-expression of TPR-like helical protein involved in stress and hormone signalling in plants [59].
The full-length transcript of an HSP mitochondrial precursor was reconstructed in MFRC. HSPs are involved in various biological processes and play crucial roles in morphogenetic changes, stress adaptation and antifungal resistance [60].
Several transmembrane transporters were over-expressed or exclusive for MFRC: MFS transporters that are involved, like the ATP-binding cassette (ABC) transporters, in multidrug resistance to natural and synthetic toxicants and playing a key role in plant pathogens as a shield against plant defence compounds during the pathogenesis and in fungicide resistance [61]; two flavin-binding monooxygenase-like proteins with CzcO conserved domain associated with the cation diffusion facilitator CzcD, which are involved in metal tolerance or resistance [62]; and a putative signal sequence protein of the Tat secretion pathway that serves to actively translocate fully folded proteins across membranes [63].
The most common SMs in fungi are polyketides, terpenoids and shikimic acid-derived compounds. PKSs and NRPSs are key enzymes in the biosynthesis of SMs containing highly conserved functional domains. The corresponding genes in fungal genomes are frequently co-localized and co-expressed with other genes coding for other enzymes involved in the same biosynthetic pathway [64]. Bifunctional solanapyrone synthases, over-expressed in both MFRC and MFRG, are involved in the biosynthesis of the polyketide-derived phytotoxins solanapyrones in Alternaria solani, Ascochyta rabiei and other fungi. Their phytotoxicity is well documented but their contribution to the pathogenicity has been questioned, while antifungal activity against microorganic competitors has been reported during saprotrophic but not parasitic growth [65]. β-lactamase transcripts were reconstructed in MFRC. Bacterial β-lactamases are known to confer resistance to β-lactam antibiotics, while the fungal homologs have been supposed to be responsible for the degradation of plant or microbial lactam compounds; for instance, a β-lactamase-containing gene cluster (FDB1) of Fusarium verticillioides confers resistance to lactam phytoanticipins [66]. A phenylacetate-CoA ligase was reconstructed in MFRG. Enzymes of this family are involved in the degradation pathway of aromatic compounds and in the biosynthesis of the β-lactam antibiotic penicillin G in P. chrysogenum [67]. In MFRC we also found a highly expressed putative cinnamoyl-CoA reductase, enzymes that in plants are involved in the phenylpropanoid biosynthetic pathway and have homologues in some bacteria and fungi although lacking the full pathway [68]. Several DETs are involved in the biosynthesis of clavulanic acid, a potent inhibitor of bacterial class A serine β-lactamases, with weak antibiotic activity. The clavaminic acid synthetase (CAS) and clavulanic acid dehydrogenase (CAD) are involved in its biosynthesis [69]. The PhyH protein identified in MFRC shares homology with a CAS-like protein of the human fungal pathogen Glarea lozoyensis; a putative CAD protein was over-expressed in MLAX; and the TauD unique in MFRG exhibits significant homology to CAS-domain-containing proteins of Penicillium, Aspergillus and other fungi. Moreover, a PhyH protein coded by the thnG gene from Streptomyces cattleya is involved in the biosynthesis of the β-lactam antibiotic thienamycin [70].
Several DETs were involved in mycotoxin metabolism. In MFRC a putative enoyl-hydratase isomerase containing a transferase domain showed similarity with the fumigaclavine B O-acetyltransferase involved in the biosynthesis of fumigaclavine C, an ergot alkaloid produced by fungi of the Trichocomaceae family [71] and with the trichothecene 3-O-acetyltransferases from Fusarium spp., which modify the trichothecene mycotoxin deoxynivalenol (DON) and reduce its toxicity [72], while a cytochrome P450 monooxygenase is related to isotrichodermin C-15 hydroxylases involved in the biosynthesis of trichothecenes [73]. DETs identified in MLAX included a benzoate 4-monooxygenase cytochrome P450 protein homologous to fungal trichodiene oxygenases, like the Fusarium sporotrichioides Tri4, involved in the trichothecene biosynthesis [74], and a cytochrome P450 monooxygenase showing similarity to enzymes essential for the synthesis of the polyketide-derived mycotoxin sterigmatocystin [75]. DETs identified in MFRG included a cytochrome P450 monooxygenase sharing similarity with enzymes that in Fusarium spp. and Aspergillus spp. are involved in the biosynthesis of the trichothecene and gliotoxin [76].
Differences in detoxification mechanisms were observed. The glutathione S-transferase assembled in MFRG is likely involved in detoxification of xenobiotics. Two nitrilases/cyanide hydratases were highly expressed in MLAX. Phytopathogenic fungi use nitrilases/cyanide hydratase to detoxify cyanide or cyanide by-products generated by plants, like amygdalin and prunasin produced by stone fruits which contribute to defence against herbivores and fungal pathogens [77].
Group II introns were detected in MFRC and MFRG. In bacterial genome and eukaryotic organelles, they are supposed to be an ancient class of ribozymes and mobile retroelements using intron-encoded reverse transcriptase, maturase and DNA endonuclease activities for site-specific insertion into homologous intron-less alleles [78]. Several transcripts related to LTR retroelements of the Ty3/Gypsy and the Ty1/Copia families were generally identified as unique in MFRG, revealing an active retrotransposition in the pathogen genome.
An RNAse III protein was uniquely assembled in MLAX. Eukaryotic RNase III or RNase III like enzymes, such as Dicer, are involved in RNAi (RNA interference) and miRNA (micro-RNA) gene silencing [79].

Conclusions
In this study, we provide the first large-scale reconstruction and annotation of the transcriptomes of the phytopathogenic fungi MFRC, MLAX and MFRG which are the main causes of heavy yield losses on pome and stone fruits all around the world. The assembled transcriptomes included about 30,000 transcripts per species that were mostly complete, with low to moderate levels of duplication, fragmentation and estimated missing transcripts. About 50% of transcripts were functionally annotated. The transcriptomes of the three Monilinia species did not show any significant differences in the GO and OG cluster functional categories. Consequently, orthologous transcripts were identified in the three species on the base of the B. cinerea proteome and, additionally, more than 400 transcripts specific for at least one of the Monilinia species with no homologs in the reference were detected. A comparative analysis among orthologs in the three species based on transcript abundance revealed DETs over-expressed or exclusive for each of the three species. They were mostly associated with biopolymer-degrading enzymatic activities, detoxifying activity, secondary metabolism, such as biosynthesis and metabolism of xenobiotics, i.e. antibiotics and toxins, as well as processes affecting fungal morphogenesis and development, diversity and pathogen interactions with the host plants and the microbiotes.
Biopolymer-degrading enzymes, which are frequently pathogenicity factors, were more numerous in MFRC and MLAX than in MFRG. They included effectors, such as a chitin-binding protein, and enzymes involved in antibiotic metabolism in MFRC, and a methylglyoxal detoxifying enzyme in MLAX. Secreted proteins related to pathogenicity were detected, like a pepsin-like protein in MLAX and a trypsin-related enzyme in MFRC which displays also a putative protease inhibitor acting as effector during pathogenesis. These results are consistent with previous findings displaying less aggressiveness components of MFRG compared to MFRC and MLAX, in terms of lesions and incubation and latency period on artificially inoculated nectarine fruit at postharvest [7]. Moreover, MFRC was characterised by a very high expression of the developmental-specific protein Ssp1, cell surface proteins acting as virulence factors, enzymes responsible for oxylipin production and regulation of spore differentiation and plant-pathogen interactions, a MSS4-like protein playing a role in apical hyphal growth, numerous MFS transporters related to pathogenesis and multidrug resistance, enzymes involved in the biosynthesis of antimicrobials (solanapyrone, clavulanic acid) and mycotoxins, and β-lactamases involved in antimicrobial resistance. MLAX was characterised by enzymes responsible for antibiotic resistance, CAD enzymes involved in clavulanic acid biosynthesis, a P450 protein putatively involved in biosynthesis of mycotoxins, and cyanide hydratases detoxifying cyanogenic compounds which are typically produced by stone fruit trees. MFRG was characterised by a glutathione S-transferase involved in detoxification of xenobiotics, solanapyrone synthases, and enzymes involved in the biosynthesis of β-lactam antibiotics, clavulanic acid and gliotoxin-like mycotoxins. The fungus displayed numerous DETs related to transposable elements indicating an intense transposition activity in its genome.
The genes differentially expressed in the three pathogens play relevant roles in morphogenesis and development, diversity and pathogenesis and are worthwhile of further investigations since they might explain their different fitness. The data obtained represent new insights in the transcriptome analyses of the three species of Monilinia and provide a new and comprehensive genetic resource that will contribute to get deeper knowledge on the population biology, physiology and plant-pathogen interactions of these important phytopathogenic fungi which are essential for improving sustainable crop protection strategies.