Comparative transcriptome analysis reveals significant differences in gene expression between pathogens of apple Glomerella leaf spot and apple bitter rot

Apple Glomerella leaf spot (GLS) and apple bitter rot (ABR) are two devastating foliar and fruit diseases on apples. The different symptoms of GLS and ABR could be related to different transcriptome patterns. Thus, the objectives of this study were to compare the transcriptome profiles of Colletotrichum gloeosporioides species complex isolates GC20190701, FL180903, and FL180906, the pathogen of GLS and ABR, and to evaluate the involvement of the genes on pathogenicity. A relatively large difference was discovered between the GLS-isolate GC20190701 and ABR-isolates FL180903, FL180906, and quite many differential expression genes associated with pathogenicity were revealed. The DEGs between the GLS- and ABR-isolate were significantly enriched in GO terms of secondary metabolites, however, the categories of degradation of various cell wall components did not. Many genes associated with secondary metabolism were revealed. A total of 17 Cytochrome P450s (CYP), 11 of which were up-regulated while six were down-regulated, and five up-regulated methyltransferase genes were discovered. The genes associated with the secretion of extracellular enzymes and melanin accumulation were up-regulated. Four genes associated with the degradation of the host cell wall, three genes involved in the degradation of cellulose, and one gene involved in the degradation of xylan were revealed and all up-regulated. In addition, genes involved in melanin syntheses, such as tyrosinase and glucosyltransferase, were highly up-regulated. The penetration ability, pathogenicity of GLS-isolate was greater than that of ABR-isolate, which might indicate that GLS-isolate originated from ABR-isolates by mutation. These results contributed to highlighting the importance to investigate such DEGs between GLS- and ABR-isolate in depth.


Introduction
Apple Glomerella leaf spot (GLS) and apple bitter rot (ABR) are two devastating foliar and fruit diseases of apple, which are a serious threat to apple production worldwide [1][2][3]. GLS is an epidemic disease leading to early defoliation and small sunken lesions (1 to 3 mm) Open Access *Correspondence: lian305@qau.edu.cn † Bowen Jiang and Ting Cai contributed equally to this work. 1 College of Plant Health and Medicine, Qingdao Agricultural University, Qingdao 266109, China Full list of author information is available at the end of the article that rarely develop into a rot on fruit [1,4]. Apple cultivars descending from the Golden Delicious group, especially the cultivar Gala, are susceptible to GLS, and the disease may cause more than 80% defoliation and diseased fruit before harvest, reducing productivity in subsequent seasons [1,5]. Nevertheless, apple cultivars from Red Delicious, such as the cultivar Fuji, are highly resistant to GLS [6]. GLS was first reported in the southwestern USA by Taylor in 1971 [7] when the disease was named "necrotic leaf blotch". Later, GLS was found in Brazil but it was also given another name as "Gala leaf spot", and finally in the 1990s, Sutton and coworkers changed the disease name to the current one as "Glomerella leaf spot" [3,8]. In China, GLS was first reported in 2011 at Fengxian in Jiangsu province and has become a widespread disease in apple-producing areas, such as Shangdong, Hebei, Liaoning, and Gansu Provinces [1,9]. ABR is an ordinary fruit rotting disease of apples that occurs worldwide causing extensive fruit rot [10]. About 30 to 60% of fruits rot before harvest while serious occurring in commercial orchards [11]. All apple cultivars have been considered susceptible to ABR, especially those belonging to the late-harvest group, such as Cripps Pink and Granny Smith, are particularly susceptible [12], whereas few reports on the disease infecting apple leaves.
The causal agent of GLS was identified as Colletotrichum species [7,13]. More Colletotrichum species were described as pathogens since GLS is vigorously occurring in apple-producing regions over the years in several countries [14]. C. fructicola is one of the main pathogenic species causing GLS worldwide [14]. The Colletotrichum species have been reclassified through morphological and molecular methodologies [15][16][17] and organized into complexes. The GLS pathogen C. fructicola, C. aenigma, C. fioriniae, C. alienum, C. siamense, C. tropicale belong to the C. gloeosporioides species complex (CGSC) [18,19], while C. karstii belong to the C. boninense species complex [17], and both C. fioriniae and C. nymphaeae fall within the C. acutatum species complex [14,16]. The CGSC was the dominant pathogens of GLS, the complex was systematically described by Weir et al. based on phylogenetic analyses of up to eight genes in 2012, within which 22 species and one subspecies were included. All taxa accepted within this clade were morphologically similar to the broadly defined C. gloeosporioides, as it has been applied in the literature for around 50 years [15].
The causal agent of ABR have been identified as seven Colletotrichum species, including C. fructicola, C. gloeosporioides, C. alienum, C. nymphaeae, C. siamense, and C. orientalis [2,5,10,20]. The Colletotrichum species, such as C. fructicola, C. alienum, and C. siamense, are common causal agents of GLS and ABR, whereas the pathogenicity of isolates that isolated from the two diseases were various on apple. For example, C. fructicola was identified as GLS type and ABR type according to its pathogenicity [2].
Melanized appressoria formed by Colletotrichum species serves as the invasive structure into the host [21]. Appressoria produce a penetration peg that penetrates the plant cuticle and cell wall layers [21]. Appressorial melanization is a key process for successful invading, and three genes PKS1, SCD1, and THR1 are associated with the melanin biosynthesis and one regulatory gene CMR1 plays a crucial role in the penetration [22][23][24][25]. Mutants in which the genes encoding these enzymes were disrupted were defective in appressorial melanization and their abilities to penetrate host plants [21].
Pathogens of ABR or GLS exhibit a particular organ specialization being able to cause quite different symptoms [2,11,26]. However, the mechanisms for such differences are unclear until now although several reports tried to make clear the relations between GLS-and ABR-isolates. The genotypes causing GLS were may be originated from ABR-isolates under the selection force of susceptible apple Gala [27]. The difference in the ability to infect fruits and leaves could be because of extracellular enzymes produced by the fungi, which degrade the plant cuticle and cell wall components while invading [26]. The extracellular enzymes produced by ABRisolate, such as lipolytic enzymes, proteolytic enzymes, pectin lyase (PNL), polygalacturonase (PG), and laccase (LAC) had higher activity than that produced by GLSisolate in vitro [26]. However, those differences in extracellular enzymes of ABR-and GLS-isolate have not been observed in vivo [26].
In this context, the transcriptome profile comparison of GLS-and ABR-isolates were conducted through RNA-Seq. GLS-isolate GC20190701 and ABR-isolates FL180903, FL180906 were showed a basal gene expression pattern. The transcriptional variations among the three isolates were further identified. In addition, the potential roles for the identified DEGs between GLSand ABR-isolates were analyzed. These results will contribute to a better understanding of the mechanism of difference between GLS-and ABR-isolates and provide insight into the relations of GLS-and ABR-isolates.

Results
The pathogenicity of GLS-and ABR-isolates GC20190701, FL180903, FL180906 were diverse to gala apple The GC20190701 was isolated from leaves of Gala apple, which showed a typical symptom of apple Glomerella leaf spot, and the FL180903, FL180906 were isolated from Fuji apple fruits that showed a typical symptom of apple bitter rot. The pathogenicity of the three isolates was diverse on apple leaves and fruits. The conidia of GC20190701 could infect Gala leaves, and cause necrotic spots, but the conidia of isolates FL180903 and FL180906 could not complete the infection on Gala leaves without a wound (Fig. 1A). Whereas, the three isolates could infect Gala apple leaves by the conidia through micro-wound (Fig. 1B), the necrotic spot caused by GC20190701 and FL180906 were similar in size, which was larger significantly than that caused by FL180903. The three isolates could infect apple fruit through micro-wound by conidia, and cause rot lesions (Fig. 1C). The rot lesions caused by GC20190701 were larger than those caused by FL180903 or FL180906, while rot lesions were similar in size caused by FL180903 and FL180906.

Multi-locus phylogenetic analysis of GLS-and ABR-isolates GC20190701, FL180903, FL180906
To elucidate the phylogenetic position of the GLS-and ABR-isolates GC20190701, FL180903, FL180906, five genes, such as ITS, GAPDH, TUB2, ACT , CHS-1, were concatenated to form a supermatrix of 1966 bp, and the phylogenetic analysis of the concatenated data set was conducted using the neighbor-joining (NJ) method. The results showed that GC20190701 were clustered in the same clade with C. aenigma, however, FL180903 and FL180906 were clustered in the same clade with C. gloeosporioides (Fig. 2). The three isolates GC20190701, FL180903, and FL180906 were closely related species, which were clustered in Colletotrichum gloeosporioides species complex (Fig. 2).

Fig. 1
Lesions caused by the Colletotrichum aenigma isolate GC20190701 and C. gloeosporioides isolates FL180903, FL180906 on leaves and fruits of Gala apple inoculated with conidia at 10 days after inoculation. A GC20190701 infected Gala leaves causing necrotic spots, ABR-isolates FL180903 and FL180906 could not infect on Gala leaves by the conidia. B All the three isolates infected Gala apple leaves by the conidia through micro-wound. C The three isolates infected Gala apple fruit through micro-wound by conidia, and cause rot lesions

Transcriptome analysis of GLS-and ABR-isolates GC20190701, FL180903, FL180906
To comprehend the mechanism of pathogenicity of the three Colletotrichum isolates GC20190701, FL180903, and FL180906, transcriptome comparison of three Colletotrichum isolates were performed. The mycelia were harvested from PDA medium at 96 hpi, and RNA was extracted and purified for transcriptome sequence. A total of nine biological samples, three biological repeats for each isolate were sequenced using next-generation sequencing on the Illumina sequence platform. Totally 56.58 billion bp clean data was achieved, which were further spliced to 38,844 assembled unigenes ( Table 1, Table S1). The total length of unigenes was 51,182,491 bp, the average length was 1318 bp, and the longest one was 23,113 bp (Table 1).

Functional annotation and classification of unigenes
Gene annotation was performed to analyze the functions of the expressed genes. Totally (Table S2). There were 3655 unigenes (9.41%) were annotated in all the databases (Table S2).
The annotated unigenes were compared to known nucleotide sequences of microbe species. The best matched to the known nucleotide sequences were C. gloeosporioides CG-14 (42.53%) and C. fructicola Nara gc5 (29.28%). Only 5.86% of unigenes matched to other four Colletotrichum species (Fig. 3A).
Clusters of Gene Ontology (GO) classification were calculated by BLAST2GO. A total of 14,445 annotated unigenes (37.19% of 38,844) were assigned to at least one of the 48 GO terms (Fig. 3B, Table S3). The unigenes were assigned to biological process, cellular component, and molecular function, respectively. The unigenes in the molecular function category were abundant in catalytic activity (GO:0003824) and binding functions (GO:0005488) (Fig. 3B). The biological process category of the unigenes was predominantly associated with metabolic process (GO:0008152), single-organism process (GO:0044699), and cellular process (GO:0009987) (Fig. 3B). The cellular component category of the unigenes was predominantly associated with cell (GO:0005623), cell part (GO:0044464), and membrane (GO:0016020) (Fig. 3B). The fundamental biological processes of the isolates GC20190701, FL180903, FL180906 were identified according to the above findings.

Differential expression analysis of unigenes of isolates GC20190701, FL180903, FL180906
To analyze the differential expression of unigenes (DEGs) among the three isolates, the transcriptomes were compared in pairs, and the expression level of unigenes was calculated by the FPKM method. In the isolate GC20190701, a total of 8302 unigenes were differentially expressed compared to the isolate FL180906, and the up-regulated unigenes and down-regulated unigenes were 4788 and 3514, respectively (Fig. 4A). There were 9455 DEGs between the isolate GC20190701 and FL180903, and the up-regulated genes and down-regulated genes were 5185, and 3640 compared to FL180903 (Fig. 4A). Whereas only a total of 1115 DEGs between the isolate FL180906 and FL180903, the up-regulated unigenes were 378, and the down-regulated unigenes were 737 compared with FL180906 (Fig. 4A).
A total of 145 up-regulated unigenes and 447 downregulated unigenes overlapped between GC20190701 and FL180903 compared with FL180906 (Fig. 4B). The transcriptome expression profile of GC20190701 was quite different from that isolates FL180903 and FL180906, whereas the transcriptome expression profile of FL180903 was not that different from that FL180906.

Gene ontology and KEGG enrichment analysis of DEGs of isolates GC20190701, FL180903, FL180906
Gene ontology (GO) enrichment analyses were conducted using a cut-off of P < 0.05 to determine the functional roles of DEGs of the three isolates. The DEGs of isolate GC20190701 and FL180906 were assigned in 3565 enriched GO terms. There were 184 significant enriched GO terms (p < 0.05), which the numbers related to biological process, molecular function and cellar component were 89, 81, and 14, respectively (Table S5) Table S5). The DEGs of isolate FL180903 and FL180906 were assigned in 909 enriched GO terms. Totally 109 significant enriched GO terms (p < 0.05) were identified, the numbers related to biological process, molecular function, and cellar component were 43, 61, and 5, respectively (Table S6).
KEGG enrichment analyses were performed to identify the basal level biological pathways of the three isolates. All DEGs between GC20190701 and FL180906 were enriched into 110 KEGG pathways (Table S7). A total of 14 pathways were significantly enriched with P values < =0.05 (Fig. 5B, Table S7). The DEGs between FL180903 and FL180906 were enriched into 60 KEGG pathways, and eight of which were significantly enriched (Table S8).

Functional analysis of the genes that differentially expressed in GLS-pathogen and ABR-pathogen
To compare the function of DEGs of GLS-and ABRpathogen, we first removed the DEGs from GC20190701 which were co-differential expressed with FL180903 compared with the expression of FL180906 since the FL180903 and FL180906 were common pathogens of ABR. Second, we selected the highly differential expression genes with a |log2FoldChange| > 3, the p value< 0.01, and the adjusted p-value < 0.01 from the transcriptome of GC20190701. A total of 1124 DEGs were selected, and 649 genes were up-regulated, 475 genes were downregulated compared with the expression of FL180906 (Table S9).
The functions of the 1124 DEGs were analyzed by searching the Swissprot database. Among the DEGs, 42 unigenes associated with the biosynthesis of secondary metabolites including cytotoxin, mycotoxin, and antibiotics were revealed, of which 30 unigenes were up-regulated whereas 12 were down-regulated (Table 2). There were 17 cytochrome P450, five methyltransferase protein, three short-chain dehydrogenase reductases, three FAD binding proteins, two multicopper oxidases, two efflux pump antibiotic resistance, two NmrA-like family proteins were involved in the 42 DEGs, besides that one of each Mgt family protein, ABC-2 type transporter, 4-hydroxyacetophenone monooxygenase, Aminotransferase classes I and II family protein, dynamin GTPase, integral membrane protein, oxidoreductase protein, thioredoxin were also included in the 42 DEGs (Table 2).
Four genes associated with the degradation of the host cell wall were identified, which were up-regulated. One Rhamnogalacturonate lyase (DN26894_c0_g1), two glucosidases (DN21088_c0_g1, DN32627_c0_g1), involve in degradation of cellulose, and one acetylxylan esterase (DN30086_c0_g1) involves in degradation of xylan.
Two tyrosinase (DN33240_c0_g2, DN34063_c0_g2), and one glucosyltransferase (DN25823_c0_g1) were upregulated, which were involved in melanin synthesis. Three CFEM domain-containing proteins, which were associated with pathogenicity were also identified. Two of them (DN36403_c6_g2, DN21054_c0_g1) were up-regulated, and one (DN27480_c0_g1) were down-regulated. Three carboxypeptidases (DN36348_c7_g7, DN34219_ c0_g1, DN28578_c1_g1) were identified, which were associated with pathogenicity, and both of them were up-regulated.

Validation of RNA-Seq data by quantitative real-time RT-PCR
We validated the RNA-Seq data by quantitative RT-PCR for six representative genes that showed strong up-regulation or down-regulation in GC2190701 compared with FL180906. The genes used for validation, their log2 fold change, and the primer sequences are presented in Table S10. For quantitative RT-PCR, we prepared new samples following the same procedures that were used to prepare samples for RNA-Seq. The expression patterns of the selected six genes all agreed with the RNA-Seq results (Fig. 6), suggesting that the RNA-Seq results were reliable in this study.

Discussion
Apple Glomerella leaf spot (GLS) and apple bitter rot (ABR) are two serious plant diseases that threaten the production of apples worldwide. The causal agents of GLS and ABR both belong to the genus Colletotrichum, and there are several common species causing both GLS and ABR. However, the common causal agents exhibit a particular organ specialization, which causes quite different symptoms on apples. Although several reports tried to explore the mechanism for such difference [2,26,27], the causes for such differences remain unknown until now. In context, we performed a transcriptome comparison of three GLS-and ABR-isolates, GC20190701, FL180903, and FL180903. A relatively large difference was discovered between the GLS-and ABR-isolate and quite many differences in expression genes associated with pathogenicity were revealed. Using the multi-locus phylogenetic analysis, we made it clear that the isolate GC20190701 belonged to C. aenigma, whereas FL180903 and FL180906 belonged to C. gloeosporioides, however, the three isolates were all fall in C. gloeosporioides species complex (CGSC), and taxa accepted within this clade were morphologically similar each other [15]. We acquired high-quality transcriptome data by sequencing nine samples, and the transcripts were assigned to the genus Colletotrichum ( Fig. 2A), therefore, the transcriptome data meet the requirement for further comparing analysis.
The transcriptome expression patterns were divergent in the three isolates GC20190701, FL180903, and FL180906. A total of 1115 unigenes were differentially expressed between FL180903 and FL180906, both of which cause ABR on fruit (Fig. 4A). Nevertheless, the number of DEGs was 9455 and 8302 in GC20190701 compared with FL180903 and FL180906, respectively (Fig. 4A). As a causal agent of GLS and ABR, the isolates not only cause the various symptoms but also Fig. 4 The number of differentially expressed genes (DEGs) among the Colletotrichum aenigma isolate GC20190701 and C. gloeosporioides isolates FL180903 and FL180906. A The DEGs were compared in pairs of the three isolates. In the isolate GC20190701, a total of 8302 unigenes were differentially expressed compared to the isolate FL180906, and there were 9455 DEGs between the isolate GC20190701 and FL180903. Whereas only a total of 1115 DEGs between the isolate FL180906 and FL180903. B Venn diagram displaying the distribution of DEGs (genes with > 2-fold change in expression) in GC20190701 and FL180903 compared with FL180906 possess significant differences in the transcriptome profiles.
The enrichment of GO terms and KEGG analysis of DEGs were conducted. The DEGs between the GC20190701 and FL180906 were significantly enriched in secondary metabolites, however, the categories of degradation of various cell wall components did not significantly enrich. The secondary metabolites were significantly divergent between the GLS-isolate and ABR-isolate, while the degradation enzyme may not change that much. However, the enrichment of GO terms of DEGs between isolate FL180903 and  FL180906 had identical trends with that DEGs between GC20190701 and FL180906 although the number of unigenes involved in the GO terms was less than that of GC20190701 and FL180906. We combined several measures in the analysis procedure to reveal rigorous DEGs between GLS-isolate and ABR-isolate. First, we removed the DEGs common in GC20190701 and FL180903 compared with FL180906 because we considered that the common DEGs were not the causes of different symptoms since both the isolates FL180903 and FL180906 cause ABR on fruit. Second, only DEGs showing more than 8-fold differences, the p value< 0.01, and the adjusted p-value < 0.01 were kept for further functional analysis. A total of 1124 DEGs were selected, which were considered as the DEGs of GLS-and ABR-isolate (Table S9). Looking through the DEGs, secondary metabolites, such as cytotoxin, mycotoxin, and antibiotics, were various between GC20190701 GLS-and ABR-isolate while extracellular enzymes were mostly up-regulated.
Killing the host tissues using the produced toxins, such as sesquiterpenoids is the main strategy for the full pathogenicity of necrotrophic fungal pathogens [28]. Several putative sesquiterpene synthases (STS) were revealed during plant infection in C. graminicola and C. higginsianum [29]. Cytochrome P450s (CYP) play essential roles in the fungal biosynthesis of secondary metabolites and detoxification of toxic compounds [28,30,31]. In this study, a total of 17 Cytochrome P450s (CYP) were revealed in the DEGs, 11 of which were up-regulated while six were down-regulated, which demonstrated that the secondary metabolites might be diversiform between GLS-and ABR-isolate. The functional roles of each CYP in the GLS-and ABR-isolate should be further explored. The methyltransferases catalyze an important pathway in the metabolism of many drugs and toxic compounds [32]. In the maize pathogen Cochliobolus heterostrophus, the Lae1-like methyltransferases act as a regulator of T-toxin production and thus impacts virulence to the host [33]. Five Methyltransferase proteins were revealed from the 1124 DEGs, which were all up-regulated. Those methyltransferases may involve in the pathways that increased the toxin level in GLS-isolate and therefore result in large necrotic spots and premature defoliation on leaves.  As hemibiotrophic fungi, CGSC infection procedures include penetration, growth inside living host cells (biotrophy), and tissue destruction (necrotrophy) [29]. To complete the pre-penetration, an invasion structure having mechanical pressure and some cuticle-and cell wall degrading enzymes, such as cutinases, pectinases, hemicellulases and cellulases are necessary [34]. Besides that, many enzymes, such as amylases, lipases, and proteases are also secreted to degrade plasma-membrane components and provide nutrients to help the fungus spread in plant tissue [35]. We retrieved four genes associated with the degradation of the host cell wall, which was up-regulated in GLS-isolate compared with ABR-isolate. Three genes were involved in the degradation of cellulose, and one gene was involved in the degradation of xylan (Table 1). Those results demonstrated that GLS-isolate secreted more cuticle-and cell wall degrading enzyme than that ABR-isolate, which might be the reason that GLS-isolate could infect Gala leaves and develop fast but ABR-isolate could not. However, Velho et al. proved that ABR-isolate had higher activity of pectin lyase (PNL), polygalacturonase (PG) and laccase (LAC) than GLS-isolate in culture broth [26]; nonetheless, they could not find a significant difference between the two isolates for those enzymes in infected apple leaves [26].
CGSC penetrates host epidermal cells through the melanized appressoria by producing a penetration peg, and melanization of appressoria is crucial for appressorial function [36,37]. In the current study, we identified genes involved in melanin synthesis, such as tyrosinase and glucosyltransferase, which were highly up-regulated ( Table 2). The melanin accumulation level might be higher in GLS-isolate than that in ABR-isolate causing a greater invasive ability on host plant tissues. In addition, we revealed several genes associated with pathogenicities, such as two CFEM domain-containing proteins and three carboxypeptidases that were upregulated. Combined with the above founding, we concluded that the penetration ability, pathogenicity of GLS-isolate was greater than that of ABR-isolate.

Conclusions
The transcriptome profile between GLS-and ABR-isolate were relatively large differences, and genes involved in the secondary metabolism and extracellular enzymes were divergent. More and higher secondary metabolites were produced in GLS-isolate, the secretion of extracellular enzymes and melanin accumulation were increased, and the genes associated with pathogenicity were also up-regulated. Therefore, the pathogenicity of GLS-isolate was higher than that of ABR-isolate, which might indicate that GLS-isolate originated from ABR-isolates by mutation into a more virulent strain. The consistent deduction was concluded in previous studies [27], and the mutation was caused mainly by the high production of susceptible apple Gala acting as a selection force [2].

Fungi isolates
GLS-isolate GC20190701 and ABR-isolates FL180903, FL180906 were originated from leaves and fruits with GLS-and ABR symptoms, respectively. The GLS leaves and ABR fruit were collected from Gala and Fuji at the experimental station of Qingdao Agricultural University, Jiaozhou City of China, respectively. The causal agents were isolated from diseased leaves or fruits by a single spore isolation method as described previously [38]. The fungus was grown on potato dextrose agar (PDA) medium at 25 °C in an incubator (MGC-400HPY; Shanghai Bluepard Instruments Co. Ltd., Shanghai, China).

Pathogenicity test assays
Conidia of the three isolates GC20190701, FL180903, FL180906 were inoculated to apple leaves and fruits to test the pathogenicity. The three isolates were incubated on a PDA medium for 3 days at 25 °C, subsequently, the mycelium was scraped from the surface of the PDA medium to induce the production of conidia. The produced conidia by each isolate were scoured into the water and adjusted to 1 × 10 4 per mL using a hemocytometer. Twenty microliter conidia were dropped to Gala leaves or fruit by a micropipettor, and the inoculated leaves or fruit were sealed in a container, which keep moisture using a wetted tissue. The sealed containers with inoculated leaves and fruits were incubated at 25 °C in an incubator, and the symptoms caused by the three isolates were examined at 10 dpi. The pathogenicity test experiment was carried out three times at different times, and five apple leaves and five fruits were used for each treatment.

RNA extraction, library construction, and sequencing
The three isolates GC20190701, FL180903, FL180906 were conducted transcriptome analysis. A total of nine isolate cultures, three replicas for each isolate, were subjected to RNA extraction. The isolates were incubated on a PDA medium for 96 h at 25 °C in an incubator, and mycelium was harvested from the surface of the medium for RNA extraction. Total RNA of the nine samples was extracted using the Trizol Reagent Kit (Sangon Biotech, Shanghai, China) and treated with RNase-free DNase I (TaKaRa) following the manufacturer's protocols. Promega PolyATtract mRNA Isolation Systems were used to purified poly(A) messenger RNA (mRNA) from the total RNAs, biotinylated beads with oligo (dT) were used to enrich mRNAs. The mRNAs were fragmented into short fragments of about 300 bp in length using Magnesium RNA Fragmentation Module (New England BioLabs). Subsequently, using the short fragments as templates, the first-strand cDNA was synthesized applying random hexamer primers, and then the second-strand cDNA was synthesized by adding the buffer, dNTPs, RNase H, and DNA polymerase I. Next, the short fragments were connected with sequencing adapters with respect to the result of agarose gel electrophoresis, and suitable fragments of about 450 bp in length were selected as templates for amplification with PCR to enrich the cDNA fragments. Next, the PCR products were purified with a QIAquick PCR Purification Kit and elution in EB buffer, the obtained products were considered as the final cDNA library for sequencing. Finally, the Quality of cDNA libraries was analyzed by Agilent 2100 Bioanalyzer, and the libraries were conducted paired-end sequencing using next-generation sequencing (NGS) on Illumina HiSeq ™ 4000 (Illumina, CA, USA). Preparation and sequencing of the cDNA library were implemented by Shanghai Personalbio Technology Co.,Ltd. (Shanghai, China).

De novo assembly of sequencing reads
To cleaned-up, the raw sequencing reads, the following criteria were applied to remove the low-quality sequences: adapter and sequence less than 50 bp, and low-quality reads with more than 50% of bases with quality lower than Q20 level. The high-quality clean reads were recovered after filtering and used for transcriptome de novo assembly with Trinity software (v2.8.4). The Trinity software first combined reads with a 30 bp length of the overlap to form longer fragments without N, and these N-free assembled reads were assembled to generated transcripts. Subsequently, the generated transcripts were clustered and the longest one of each transcript was considered as a unigene. Reads counts of each unigene were calculated by mapping the transcripts back to unigenes.

Functional annotation and differentially expressed genes
To functional analyze the transcriptomes, the assembled unigenes were firstly aligned by BLASTX to protein databases, such as NR (NCBI non-redundant protein sequences), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genome), eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups), Swiss-Prot, and Pfam, to retrieve proteins with the highest sequence similarity with the unigenes along with their protein functional annotations. Then, the gene expression level was calculated with RSEM (v1.1.12) using the transcripts as a reference sequence database. The clean reads of each sample were aligned to the transcripts database, the reads of each transcript were calculated for each sample. To compare the difference of gene expression among different samples, the FPKM (Fragments per kilobase of transcript per million mapped reads) method was used for normalization [41].

Gene ontology and KEGG pathway enrichment analysis
To examine the biological functions and pathways of the identified DEGs, we firstly annotated the DEGs with GO database (http:// www. geneo ntolo gy. org/) using a hypergeometric test [42]. GO terms that are significantly enriched in DEGs compared to the genome background were retrieved by GO functional enrichment analysis. Briefly, the GO functional enrichment analysis firstly maps all DEGs to GO terms in the database, calculating gene numbers for every term, then using the ultra-geometric test to find significantly enriched GO terms (P-value < 0.05) in DEGs compared to the genome background.
The significantly enriched metabolic pathways and signal transduction pathways in DEGs compared with the whole genome background identified by KEGG pathway enrichment analysis [43]. Using FDR = 0.05 and P-value < 0.05 as the threshold, pathways were defined as those with significant enrichment for DEGs.