Transcriptome analysis of fungicide-responsive gene expression profiles in two Penicillium italicum strains with different response to the sterol demethylation inhibitor (DMI) fungicide prochloraz

Background Penicillium italicum (blue mold) is one of citrus pathogens causing undesirable citrus fruit decay even at strictly-controlled low temperatures (< 10 °C) during shipping and storage. P. italicum isolates with considerably high resistance to sterol demethylation inhibitor (DMI) fungicides have emerged; however, mechanism(s) underlying such DMI-resistance remains unclear. In contrast to available elucidation on anti-DMI mechanism for P. digitatum (green mold), how P. italicum DMI-resistance develops has not yet been clarified. Results The present study prepared RNA-sequencing (RNA-seq) libraries for two P. italicum strains (highly resistant (Pi-R) versus highly sensitive (Pi-S) to DMI fungicides), with and without prochloraz treatment, to identify prochloraz-responsive genes facilitating DMI-resistance. After 6 h prochloraz-treatment, comparative transcriptome profiling showed more differentially expressed genes (DEGs) in Pi-R than Pi-S. Functional enrichments identified 15 DEGs in the prochloraz-induced Pi-R transcriptome, simultaneously up-regulated in P. italicum resistance. These included ATP-binding cassette (ABC) transporter-encoding genes, major facilitator superfamily (MFS) transporter-encoding genes, ergosterol (ERG) anabolism component genes ERG2, ERG6 and EGR11 (CYP51A), mitogen-activated protein kinase (MAPK) signaling-inducer genes Mkk1 and Hog1, and Ca2+/calmodulin-dependent kinase (CaMK) signaling-inducer genes CaMK1 and CaMK2. Fragments Per Kilobase per Million mapped reads (FPKM) analysis of Pi-R transcrtiptome showed that prochloraz induced mRNA increase of additional 4 unigenes, including the other two ERG11 isoforms CYP51B and CYP51C and the remaining kinase-encoding genes (i.e., Bck1 and Slt2) required for Slt2-MAPK signaling. The expression patterns of all the 19 prochloraz-responsive genes, obtained in our RNA-seq data sets, have been validated by quantitative real-time PCR (qRT-PCR). These lines of evidence in together draw a general portrait of anti-DMI mechanisms for P. italicum species. Intriguingly, some strategies adopted by the present Pi-R were not observed in the previously documented prochloraz-resistant P. digitatum transcrtiptomes. These included simultaneous induction of all major EGR11 isoforms (CYP51A/B/C), over-expression of ERG2 and ERG6 to modulate ergosterol anabolism, and concurrent mobilization of Slt2-MAPK and CaMK signaling processes to overcome fungicide-induced stresses. Conclusions The present findings provided transcriptomic evidence on P. italicum DMI-resistance mechanisms and revealed some diversity in anti-DMI strategies between P. italicum and P. digitatum species, contributing to our knowledge on P. italicum DMI-resistance mechanisms.


(Continued from previous page)
Conclusions: The present findings provided transcriptomic evidence on P. italicum DMI-resistance mechanisms and revealed some diversity in anti-DMI strategies between P. italicum and P. digitatum species, contributing to our knowledge on P. italicum DMI-resistance mechanisms.
Keywords: Transcriptome, Penicillium italicum, Demethylation inhibitor (DMI)-resistance, Prochloraz-responsive genes Background Penicillium digitatum (green mold) and P. italicum (blue mold) are well known as the predominant citrus pathogens causing postharvest diseases during fruits storing and transportation. The resulted economic losses are so great that aroused enormous attentions all over the world [1]. The sterol demethylation inhibitor (DMI) fungicides, such as imazalil and prochloraz, have been widely applied to control citrus molds [2][3][4][5][6]. However, resistance to these DMI fungicides has frequently occurred for the Penicillium molds in the past decade, especially for P. digitatum isolates with high DMIresistance [5,7], considerably reducing the efficacy of the fungicides. Up to date, we have got some understanding on the mechanism of azole fungicide resistance in P. digitatum [8][9][10][11][12][13]. However, little information is available to explain P. italicum resistance induced by the DMI fungicides. It would be theoretically important to address molecular background of P. italicum isolates causing their DMI resistance.
RNA sequencing (RNA-seq) technology has become a powerful tool to profile transcriptomic response to reveal azole-resistance mechanism for some pathogenic fungi including prochloraz-resistant P. digitatum [11], voriconazole-resistant A. fumigatus [81], tetraconazoleresistant Cercospora beticola [82], tebuconazole-resistant Fusarium culmorum [83], and fluconazole-resistant Candida glabrata [84]. Our earlier report has elucidated the mechanism of P. digitatum resistance to DMIfungicide prochloraz through RNA-seq analysis [11]. Nevertheless, the molecular mechanism(s) of P. italicum resistance to such fungicides are poorly understood. Now we have isolated two P. italicum strains exhibiting desirably contrasting response to common DMI fungicides including prochloraz, i.e. Pi-R (highly resistant to prochloraz with EC 50 = 30.2 ± 1.5 mg·L − 1 ) versus Pi-S (highly sensitive to prochloraz with EC 50 = 0.007 ± 0.002 mg·L − 1 ). The purpose of this work was to compare transcriptomic profiles between these two P. italicum strains with and without prochloraz treatment, to identify differentially expressed genes (DEGs) involved in the azole-class fungicide resistance, and to provide theoretical cues to explain P. italicum anti-azole mechanism.

Transcriptome sequencing and assembly
In the present study, Pi-R and Pi-S were treated with or without DMI-fungicide prochloraz to prepare four RNAseq samples, i.e., Pi-R-I, Pi-R-NI, Pi-S-I and Pi-S-NI .  After Illumina sequencing, the four transcriptomic  libraries contained 61,610,574, 70,012,472, 61,976,398 and 67,336,730 raw reads, respectively. By removing adaptor sequences and undesirable reads (ambiguous, low quality, and duplicated sequence reads), 58,744,798, 66,490,626, 59,134,840 and 64,262,170 clean reads were generated from the four libraries with Q30 > 90%, suggesting high quality for the present sequencing results. These clean reads were predominantly distributed in exon and intergenic regions (Additional file 4: Figure  S2). Using reference genome (PHI-1) as mapping template, clean reads were assembled into 47,195,871, 54, 176,219, 48,955,731 and 53,362,929 unigenes for the four libraries, respectively. All unigene expression levels in the four libraries were classified into five intervals, according to FPKM values (Table 1), and more than 50% of the total unigenes in each library were defined as highly expressed (i.e., FPKM interval ≥ 15).

Identification and analysis of differentially expressed genes (DEGs)
Based on the above FPKM values, hierarchical cluster (i.e., heat map) analysis was performed to visualize DEG profiles between Pi-R-I, Pi-R-NI, Pi-S-I and Pi-S-NI libraries (Fig. 1). Pi-R and Pi-S were gathered into two independent groups each containing two clusters (i.e., with and without prochloraz induction). Noticeably, prochloraz induced more dramatic change in gene expression profile between Pi-R-I and Pi-R-NI than between Pi-S-I and Pi-S-NI, suggesting the involvement of more DEGs in Pi-R response to prochloraz.
Further, the q-value 0.005 (i.e., corrected p-value 0.005) and an absolute value of log2(fold change) ≥ 1 were set as cut-off standard to identify DEGs between different libraries, including a) Pi-R-I vs Pi-R-NI, b) Pi-S-I vs Pi-S-NI, c) Pi-R-I vs Pi-S-I, and d) Pi-R-NI vs Pi-S-NI (Fig. 2). We  Table S6), representing different genetic background between the two P. italicum strains. Among these DEGs, we identified a considerable amount of common-accepted target protein genes associated with DMI resistance, including cytochrome P450 genes and drug efflux pump genes (ABC and MFS genes rather than  Table 2). Based on the volcano plot analysis, we applied Venn diagrams to profile the DEG distribution between Pi-R-I vs Pi-R-NI and Pi-S-I vs Pi-S-NI (Fig. 3a) and between Pi-R-I vs Pi-S-I and Pi-R-NI vs Pi-S-NI (Fig. 3b). As shown in Fig. 3b, the overlap part of circles Pi-R-I vs Pi-S-I and Pi-R-NI vs Pi-S-NI comprised 513 DEGs that might represent DEGs irrelevant to prochloraz induction. In contrast, only 110 DEGs were distributed in the overlap part of circles Pi-R-I vs Pi-R-NI and Pi-S-I vs Pi-S-NI (Fig. 3a), indicating a proportion of DEGs potentially involved in prochloraz response in both resistant and sensitive P. italicum strains.

GO and KEGG enrichments of prochloraz-responsive DEGs
The DEGs were classified into three GO categories by the Blast2GO (GOseq R package: http://www.geneontology. org), including biological process (BP), cellular component (CC), and molecular function (MF). The number of total GO terms and its distribution in the three categories for each comparison are listed in Table 3. In the comparison Pi-R-I vs Pi-R-NI ( Importantly, the up-regulated DEGs mapped to specific GO terms included a number of typical genes related to fungicide resistance. As summarized in Table  4, drug-pump genes (ABC1, ABC2, MFS1, MFS2, MFS3 and MFS4, mapped to GO:0016020 (membrane)), drugtarget P450 gene (CYP51A, mapped to GO:0055114 (oxidation-reduction process)), steroid biosynthesisrelated genes (ERG2 and ERG6, mapped to GO: 0006694 (steroid biosynthetic process)) and MAPK/calcium signaling-related genes (Mkk1, Hog1, CaMK1, CaMK2 and EF-hand1, mapped to GO:0016301 (kinase activity) and GO:0005509 (calcium ion binding)) were up-regulated in prochloraz-treated Pi-R, as compared to drug-untreated Pi-R or to drug-treated Pi-S. In contrast, most of these prochloraz-responsive DEGs, except for CYP51A, were down-regulated or unchanged in prochloraz-treated Pi-S, comparing to untreated Pi-S. GO enrichment also indicated lower transcript abundance of some of these prochloraz-responsive DEGs in Pi-R when compared with Pi-S under fungicide-free conditions (Table 4), including ABC2, MFS1, MFS2, MFS4, and CaMK2. The GO-term map distribution (i.e., hit and ranking records) of the prochlorazresponsive DEGs mentioned above was summarized in Additional file 9: Table S7.
Further, KEGG enrichment was applied to identify pathways associating the prochloraz-responsive DEGs with resistance mechanisms. In the present four comparisons, KEGG analysis enriched prochlorazresponsive DEGs into only two pathways, i.e., steroid biosynthesis (KEGG ID: pcs00100; q value = 0.013) and MAPK signaling pathway (KEGG ID: pcs04011; q value = 0.021) ( Table 5): the former pathway exclusively included up-regulated DEGs, i.e., CYP51A (PITC_ 083360) in the comparisons Pi-R (I/NI) and Pi-S (I/NI), ERG2 (PITC_020620) in the comparisons Pi-R (I/NI) and Pi-S (I/NI), and ERG6 (PITC_014340) in the comparisons Pi-R (I/NI) and I (Pi-R/Pi-S); while the latter pathway included 1) up-regulated DEGs (i.e., Mkk1 and Hog1) in Pi-R-involved comparisons, i.e., Pi-R (I/NI) and I (Pi-R/Pi-S) and 2) down-regulated DEG (i.e., Hog1) only in comparison Pi-S (I/NI). All the KEGGenriched DEGs, as components of metabolic and/or signal-transduction pathway(s), were well coincident with the results of GO enrichment. In other words, the present GO-enriched DEGs, if involved in specific biological pathway(s), were exclusively KEGG-included, and certainly, pathway-irrelevant genes, e.g., drugpump genes and drug-target genes, were KEGGexcluded, without exception.  Real-time quantitative PCR (qRT-PCR) validation of prochloraz-responsive DEGs According to the GO and KEGG enrichments combined, we selected 15 prochloraz-responsive DEGs to perform qRT-PCR validation. The 15 prochloraz-responsive DEGs, never reported before, potentially involved in P.

Discussion
In the past decades, conventional synthetic DMIfungicides, such as prochloraz and imazalil, were widely applied to control Penicillium decay, but undesirably, a considerable number of resistant isolates including P. digitatum and P. italicum strains have developed [5][6][7]. The mechanisms underlying DMI-fungicide resistance have been elucidated for P. digitatum species by transcriptomic analysis [11]. However, how to develop DMIresistance in P. italicum species is still not clear, might due to rare opportunity to find highly DMI-resistant P. italicum strain(s). The EC 50 values of P. italicum isolates towards DMI-fungicide(s) (e.g., imazalil), published to date, were ≤ 0.92 ± 0.09 mg·L − 1 , no more than moderate resistance level [5,85]. Nevertheless, some recent investigations have suggested evolutional potential to develop high DMI-resistance in P. italicum species [85][86][87]. Now we have isolated a P. italicum strain (Pi-R) with extremely high resistance to some common DMIfungicides including prochloraz (Additional file 1: Figure  S1 and Additional file 2: Table S1). We believed that this strain could be useful to investigate DMI-fungicide resistance mechanism in P. italicum. Fungal resistance to azole-fungicides including a number of DMI-fungicides has been usually ascribed to overexpression of specific drug-efflux pumps such as ABC and MFS transporters [8, 42-50, 53, 54, 57, 58, 88]. Specially, ABC and MFS transporter-encoding genes, each containing multiple isoforms, were reported to be simultaneously up-regulated in the prochloraz-resistant P. digitatum [11]. The similar up-regulation of multiple  ABC and MFS gene members (i.e., ABC1-2 and MFS1-4) was also observed in the prochloraz-treated Pi-R rather than Pi-S (Table 4 and Fig. 6a), indicating the need of drug-efflux pumps to develop DMI-resistance in P.
italicum. The evidence that the knockout of mfs1 (i.e., MFS1-encoding gene) in Pi-R decreased the fungal resistance to prochloraz confirmed the involvement of particular drug-efflux pump in P. italicum anti-DMI mechanism (Additional file 12: Figure S4), and more biological support would be required for the other drugefflux pumps. Interestingly, regarding these drug-efflux pump-encoding genes, over-expressed in the prochlorazinduced P. italicum transcriptome, their orthologous genes in the prochloraz-induced P. digitatum transcriptome [11], as shown in Additional file 13: Table S9, were not responsive or negatively responsive to the DMIfungicide treatment. Such different isogene responsive profiles suggested a special drug-pump preference for different Penicillium species to develop DMI-fungicide resistance. Here we provide the first evidence that the two Penicillium species (i.e., green mold and blue mold) with close evolutionary association did differ in selecting drug-transporters to support their DMI resistance, and the mechanism(s) underlying need further research. More interestingly, regarding another class of drugpump-encoding genes, i.e., MATEs, none was found to be differentially expressed in the prochloraz-induced Pi-R and Pi-S transcriptomes ( Table 2). In contrast, Liu et al. [11] reported up-regulation of three MATE transporter-encoding genes (MATE1-3) in P. digitatum resistance to prochloraz. The orthologous gene of PdMATE3 has been identified in the present Pi-R and Pi-S genomes, however, no mRNA transcripts of this gene was detected in the prochloraz-treated P. italicum strains. This might suggest transcriptional irrelevance of MATE transporter with P. italicum resistance to prochloraz. Unlike ABC and MFS transporters, MATE (See figure on previous page.)

Fig. 5 Distribution of up-and down-regulated genes in the top 30 enriched GO terms for Pi-R-I vs Pi-R-NI (a), Pi-S-I vs Pi-S-NI (b), Pi-R-I vs Pi-S-I (c), and Pi-R-NI vs Pi-S-NI (d).
In each panel, from top to bottom shows three GO categories, i.e., biological process, cellular component, and molecular function, comprising the top 30 GO terms in total, and the up-and down-regulated genes in each term is represented by red and green bars, respectively. X-and Y-axis indicate GO terms and corresponding number of DEGs, respectively   (Table S2- 5), and the slash '/' indicates the unfeasibility to generate log2(fold change) when the DEG is not included in the comparison transporters are more associated with bacterial drugresistance [59][60][61][62]. To date, the transcriptional response of MATE isogenes to DMI fungicides was only documented in prochloraz-resistant P. digitatum transcriptomes [11]. But such transcriptomic response of MATE transporter-encoding genes did not occur in the present P. italicum transcriptomes. Thus a potential debate on real function of MATE family member(s) in citrus Penicillium pathogens' resistance to DMIfungicides would be an interesting study topic. Over-expression of ERG11s, the P450-dependent sterol 14α-demethylase (CYP51)-encoding genes (e.g., CYP51A/ B/C), has been accepted as a primary strategy to develop fungal DMI-resistance [9,27,33]. Under prochloraz induction, the increasing fold of CYP51A in Pi-R (~9-fold) was more than that in Pi-S (~4.5-fold) (Fig. 6b),  All values obtained in the qRT-PCR analysis were expressed as the mean ± SD from 5 biological repeats each containing 3 technical replicates, and independent samples t-test (n = 5) was applied in the SPSS Statistics 17.0 context to assess the significance of differences between the means (*p < 0.05 and **p < 0.01) suggesting a positive correlation between CYP51A abundance and anti-DMI potential in P. italicum strains. On the other hand, the simultaneous up-regulation of ERG11 isoforms (CYP51A/B/C) was observed only in prochlorazinduced Pi-R, rather than Pi-S (Table 4 and Fig. 6b). This indicated that CYP51B and CYP51C could contribute to P. italicum prochloraz-resistance. However, according to the prochloraz-induced P. digitatum transcriptome [11], the orthologous genes of CYP51B and CYP51C in the prochloraz-resistant P. digitatum strain were not upregulated (Additional file 13: Table S9). The difference in choice of ERG11-encoding targets suggested diverse strategies for P. italicum and P. digitatum species to develop DMI-resistance. ERG2 and ERG6, the other two ERGs in ergosterol-biosynthesis pathway, as recommended to be potential multidrug targets [89], were up-regulated in the prochloraz-treated Pi-R (Tables 4, 5 and Fig. 6b). These ERGs could be contributors to the P. italicum DMIresistance. Previous reports have documented azoleinduced mRNA increase of ERG2 and ERG6 [39,81], correlating these ergosterol-biosynthesis genes with acute response of fungal membrane rigidity caused by azolefungicides. Such ERGs expression profile was observed in the DMI-resistant P. italicum, but not in the prochlorazresistant P. digitatum strain [11] (Additional file 13: Table  S9). These results suggested different strategies in ERGs response of the two Penicillium species to develop DMIresistance. Fungal multidrug resistance has been linked with MAPK signaling pathways including Slt2-MAPK [73,[90][91][92]. The Slt2-MAPK middle-stream components, two MAPKK-encoding genes ScMkk1 and ScMkk2, have been characterized in the model yeast S. cerevisiae to function in CWI signaling-mediated fungicide resistance [90,91]. In the present study, Mkk1 (PITC_088710), the homolog of ScMkk1 in P. italicum species, was GO-and KEGG-enriched as up-regulated DEG in the prochloraztreated Pi-R (Tables 4 and 5), indicating the involvement of Slt2-MAPK in the P. italicum DMI-resistance. Bck1 (PITC_061930) and Slt2 (PITC_008290), the MAPKKKand MAPK-encoding genes located upstream and downstream of Mkk1, respectively, were both up-regulated in prochloraz-treated Pi-R ( Fig. 6c and Additional file 10: Table S8). These results supported Slt2-MAPK contribution to the P. italicum prochloraz-resistance. Disruption of slt2 led to increased sensitivity to DMI-fungicide imazalil in wheat pathogen M. graminicola [93], also suggesting Slt2-MAPK function in fungal DMIresistance. However, slt2-deleted P. digitatum exhibited no obvious change in DMI-resistance [73], indicating that Slt2-MAPK might be irrelevant to the fungal anti-DMI strategy. Actually, according to prochloraz-resistant P. digitatum transcriptomes [11], the orthologs of PiSlt2-MAPK genes were not included in up-regulated DEGs (Additional file 13: Table S9). These contrasting results indicated an interesting difference in MAPK choice for different fungal pathogens to resist DMI fungicide(s). Fungal resistance to azole fungicides additionally required some cellular Ca 2+ -signaling processes via CaMKs, specially, via CaMK1, CaMK2 and their homologues [78][79][80]94]. Here we also identified two CaMK homologues CaMK1 (PITC_087700) and CaMK2 (PITC_025800), both up-regulated only in prochloraztreated Pi-R (Table 4 and Fig. 6d), thus suggesting the contribution of CaMKs to fungal DMI-resistance. CaMKs mediated fungal responses to cell-wall and oxidative stresses induced by fungicides [77][78][79][80]. Thus the involvement of CaMKs in P. italicum DMI-resistance would be a strategy to cope with some azole-induced stresses. Unlike Pi-R, the prochloraz-resistant P. digitatum strain showed no up-regulation of CaMK homologues (Additional file 13: Table S9). The difference in CaMK response to DMI-fungicides in the two Penicillium species needs further research.

Conclusions
In conclusion, the present work for the first time provided transcriptomic analysis of prochloraz-responsive gene expression profiles for two P. italicum strains with contrasting response to prochloraz, revealing potential mechanisms underlying P. italicum resistance against DMI fungicides. The strategies that P. italicum species adopt to overcome the azole stresses, based on DEG enrichment analysis, are summarized as 1) up-regulation of specific isogenes encoding ABC and MFS transporters, 2) simultaneous induction of all major EGR11 isoforms including CYP51A/B/C, 3) over-expression of ERG2 and ERG6 to modulate ergosterol anabolism, and 4) concurrent mobilization of Slt2-MAPK and CaMK signaling processes to adapt azole-induced stresses. Some differences in the choice of anti-DMI strategy between P. italicum and P. digitatum species were also discussed.

Methods
Strains, cultivation and treatments P. italicum strains Pi-R and Pi-S, used in this study, were isolated from rotten citrus fruits in local packinghouses (Yunnan Province) and orchards (Hainan Province), respectively. Pi-R exhibits dramatically high prochloraz-resistance (EC 50 = 30.2 ± 1.5 mg·L − 1 ), while Pi-S is prochloraz-sensitive (EC 50 = 0.007 ± 0.002 mg·L − 1 ) (Additional file 1: Figure S1). Meanwhile, Pi-R and Pi-S did show contrasting response to the other two DMI fungicides (imazalil and triadimefon) and the two benzimidazole-class fungicides (carbendazim and benomyl), according to their EC 50 values; however, the two fungal strains were both highly sensitive to the phenylpyrrole fungicide fludioxonil (Additional file 2: Table  S1). Fungal strains were routinely cultivated on potato dextrose agar (PDA) medium for 5 days to prepare respective conidial suspension (10 7 spores mL − 1 ) as previously described [11]. Afterwards, for each fungal strain, 200 μL of conidial suspension (approximately 2 × 10 6 spores) was incubated with 200 mL potato dextrose broth (PDB) medium for 2 days at 28°C and 180 rpm shaking, and the resulting mycelia were treated with or without prochloraz. In detail, prochloraz was added to the PDB medium at final concentration that was in agreement with EC 50 value for each P. italicum strain (i.e., 30.2 mg·L − 1 for Pi-R and 0.007 mg·L − 1 for Pi-S). Prochloraz was pre-dissolved in 100 μL DMSO, and the same volume of DMSO was added to 200 mL PDB medium to prepare control samples. The prochlorazinduced and no-induced (control) samples were cultured at the same conditions (at 28°C and 180 rpm shaking) for 6 h before RNA extraction. The present study collected 4 samples in total for the following RNA manipulations, i.e. prochloraz-induced and no-induced Pi-R (designated as Pi-R-I and Pi-R-NI, respectively), and prochloraz-induced and no-induced Pi-S (designated as Pi-S-I and Pi-S-NI, respectively).

RNA extraction, RNA-seq library construction and Illumina sequencing
Total RNA was extracted for the four fungal samples with three biological replicates, according to the method described before [11]. The integrity of RNA samples was first examined by 1% (w/v) agarose gel electrophoresis analysis and then confirmed by Bioanalyzer 2100 (Agilent Technologies, CA, USA). Each RNA sample was further checked by Nano-Photometer (Implen, CA, USA) to ensure its purity, and quantified using Qubit RNA Assay Kit (Life Technologies, CA, USA). Sequencing libraries were constructed according to the protocol of NEBNext Ultra™ RNA Library Prep Kit for Illumina sequencing (NEB, USA), using the amount of 3 μg RNA for each sample. Briefly, Poly (A) mRNA was purified and enriched from total RNA by oligo-dT paramagnetic beads, and fragmented to provide templates for the first-and second-strand cDNA synthesis. The cDNA products were end-repaired to be blunt fragments and adenylated at their 3′ ends. NEBNext adaptors (sequencing adapters) were ligated with the 3′-adenylated DNA fragments. Then 3 μL USER Enzyme (NEB, USA) was applied to produce doublestranded cDNA templates for polymerase chain reaction (PCR) amplification with Pfu DNA polymerase and specific primers. The obtained index-coded samples were clustered by cBot Cluster Generation System for sequencing using Illumina HiSeq X platform, and the resulting paired-end reads (raw reads) in~150 bp length were deposited at NCBI database under accession number PRJNA421419.

Reads mapping to the reference genome
Raw reads stored in fastq format after the Illumina sequencing were first processed through in-house perl scripts to thoroughly remove low quality reads, as described before [11], to generate clean reads with high quality that were assessed by parameters Q20, Q30 and GC content. The clean reads were mapped to P. italicum PHI-1 reference genome (GenBank accession number: JQGA01000000) [95] using TopHat version 2.0.11 [96]. Prior to the reads mapping, the gene annotation files were downloaded from the website (http://genome.jgi. doe.gov/Pendi1 /Pendi1.home.html), and the reference genome was indexed by Bowtie version 2.0.6.

Gene expression analysis and functional enrichments
Gene expression level (transcript and/or unigene abundance) in each sample was estimated using HTSeq v0.6.1, according to FPKM analysis [97], which processed based on uniquely mapped reads and thus eliminating the experimental bias due to sequencing discrepancies. Then, the read counts were adjusted by edgeR program package through one scaling normalized factor, as previously described [98], to prepare for differential gene expression analysis. To identify differentially expressed genes (DEGs) between samples, fold-changes of expression level for each gene, defined as the ratio of the RPKM values, were calculated by DEGSeq R package (1.12.0), and P-values were statistically corrected to assess the significance for the differences in transcript abundance according to Benjamini & Hochberg method [99]. In the present study, DEGs were selected as transcripts and/or unigenes differentially expressed with at least 2-fold change (i.e., the absolute value of log2 Fold change ≥1.0) and corrected P-value ≤0.005 between two groups of comparison. The identified DEGs were hierarchically clustered by Cluster 3.0 [100], and then subjected to heat-map analysis by Plotly (Montreal, Quebec) software and Venn diagram analysis at the website http://bioinfogp.cnb.csic.es/tools/venny/ index.html. Further, DEGs were functionally enriched to Gene Ontology (GO) database (http://www. geneontology.org) using the GOseq R package based on Wallenius' non-central hyper-geometric distribution [101], and also enriched to Kyoto Encyclopedia of Genes and Genomes (KEGG) public database (http://www.genome.jp/kegg/), using the KOBAS software to test the significance of enriched DEGs in particular metabolic and signal transduction pathways [102].

Quantitative real-time PCR (qRT-PCR) validation
Nineteen unigenes related to azole-drug resistance in the present P. italicum transcriptoms were selected for qRT-PCR validation, including 6 drug transporter genes, 5 ergosterol biosynthesis-related genes, 4 MAPK signaling pathway-related genes, and 4 Ca 2+ signal transducer-related genes. Total RNA was extracted from the same P. italicum samples used for RNA-seq, according to Fungi RNA Kit user guide (OMEGA, USA). First-strand cDNA synthesis was performed with PrimeScript™RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China), according to the manufacturer's instructions, and the qRT-PCR was performed using a BIO-RAD CFX96 qPCR system with SYBR Green I fluorescent dye detection as previously described [11]. The specific primers used in this study are shown in Additional file 3: Table  S2, and the two housekeeping genes, i.e., β-actin and glyceraldehyde-3-phosphate dehydrogenase (GAPD), were respectively applied as internal reference to calculate the relative mRNA abundance for the selected unigenes, according to the 2 -ΔΔCt method [103] with 5 biological repeats each containing three technical replicates. Relative ratios for the expression of each selected unigene were further calculated in the 4 comparison groups, including Pi-R (I/NI) (i.e., Pi-R-I relative to Pi-R-NI), Pi-S (I/NI) (i.e., Pi-S-I relative to Pi-S-NI), I (Pi-R/Pi-S) (i.e., Pi-R-I relative to Pi-S-I), and NI (Pi-R/Pi-S) (i.e., Pi-R-NI relative to Pi-S-NI). All values obtained in the qRT-PCR analysis were expressed as the mean ± SD (standard deviation of the mean), and based on 5 independent experiments (i.e., 5 biological repeats), independent samples t-test (n = 5) was applied in the SPSS Statistics 17.0 context to assess the significance of differences between the means (*p < 0.05 and **p < 0.01).
Construction of mfs1-knockout mutant (Δmfs1) from Pi-R The mutant Δmfs1 was derived from its parental strain Pi-R by introduction of mfs1-knockout cassette into the DMI-resistant P. italicum protoplasts that were prepared according to the method of Zhao et al. [104]. Doublejoint PCR was performed to construct the desirable knockout cassette, also according to the method of Zhao et al. [104]. In detail, the hygromycin phosphotransferase (hph) resistance gene (approximately 2.1 kb) from the plasmid pTFCM (generously provided by Dr. Daohong Jiang, Huazhong Agricultural University, China) was inserted between the upstream and downstream flanking sequences of mfs1 gene (approximately 2.0 kb in total) from the Pi-R genomic DNA. The knockout cassette was introduced into wild-type Pi-R protoplasts as described by Zhao et al. [104], and thus the Pimfs1 gene was replaced by the hph resistance gene.