Skip to main content

Advertisement

Springer Nature is making Coronavirus research free. View research | View latest news | Sign up for updates

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

Abstract

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.

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 DMI-resistance [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.

The mechanism of fungal DMI-resistance involves strategies targeting ergosterol-biosynthesis enzymes. The site mutations in CYP51s (ERG11-encoding proteins) can alter drug-target interactions and increase DMI-resistance for various fungal pathogens, as reported in the model yeast Saccharomyces cerevisiae [14,15,16], the clinical pathogens Candida albicans [17,18,19,20] and Aspergillus fumigatus [21,22,23], and the plant pathogens Mycosphaerella graminicola [24, 25], Monilinia fructicola [26] and P. digitatum [27]. Fungal resistance to DMIs can also be ascribed to over-expression of CYP51s, especially by some enhancer elements [9, 27,28,29,30,31,32,33]. In addition to CYP51s, recently, other genes encoding fungal ergosterol biosynthesis-related enzymes have been proposed to be potential targets, including ERG2 (encoding C− 8 sterol isomerase) [34,35,36] and ERG6 (encoding C− 24 sterol methyltransferase) [37,38,39,40]. The importance of both ERG2 and ERG6 to cycloheximide resistance for S. cerevisiae has also been genetically emphasized [41].

Fungal DMI-resistance has also been ascribed to specific drug-transporter proteins that can reduce fungicide accumulation in fungal cells, including ATP-binding cassette (ABC) transporter family proteins, major facilitator superfamily (MFS) proteins, and multidrug and toxic compound extrusion (MATE) family proteins. ABC transporters have been functionally characterized in many fungal pathogens including green mold and verified to be up-regulated in their fungicide resistance [42,43,44,45,46,47,48,49,50,51,52,53,54]. MFS proteins constitute another class of broad-spectrum transporters to develop fungal DMI-resistance, including CaMDRl in C. albicans [55], MgMfsl in wheat pathogen Mycosphaerella graminicola [56], and PdMFS1 and PdMFS2 in P. digitatum strains [57, 58]. Unlike ABC and MFS transporters, MATE proteins function predominantly in bacterial drug-resistance [59,60,61]. To date, the MATE contribution to fungal drug-resistance was only reported in the ectomycorrhizal fungus Tricholoma vaccinum [62] and the citrus pathogenic fungus P. digitatum [11].

Fungicide resistance is further associated with particular protein kinase signaling and calcium (Ca2+) signaling. The mitogen-activated protein (MAP) kinase signaling pathways, ubiquitously found in eukaryotes (from yeasts to various pathogenic fungi), comprise a set of cascaded protein kinases, MAP kinase kinase kinase (MAPKKK), MAP kinase kinase (MAPKK) and MAP kinase (MAPK), acting in series to modulate target protein activities [63, 64]. Three major MAPK signaling pathways, Fus3/Kss1, Hog1, and Slt2, have been revealed in model yeasts [65,66,67] and filamentous fungi, including the citrus pathogens Alternaria alternata [68,69,70,71] and P. digitatum [72, 73], regulating pheromone/invasion processes, high osmolarity glycerol anabolism, and stress-induced cell wall remodeling, respectively. Hog1-MAPK (PdOs2)-mediated CWI signaling are involved in P. digitatum resistance to the fungicides iprodione and fludioxonil [72]. Hog1 homolog BcSak1 was identified in Botrytis cinerea and functionally required for iprodione resistance [74, 75]. FgOs2 also participated in Fusarium graminearum resistance to fludioxonil [76]. The latest evidence has suggested an essential role of PdSlt2 MAPK in regulating gene expression to develop azole-fungicide resistance [73]. Ca2+ signaling via Ca2+/calmodulin (CaM)-dependent kinases (CaMKs), usually linked with particular MAPK pathway(s), extensively participates in fungal responses to environmental stresses. The over-expression of CaMK2 (also named Cmk2) in the yeast S. cerevisiae facilitated its resistance to some azole-fungicides (e.g., dithiothreitol and miconazole) [77]. Recent studies also implied the essential role of CaMKs in protecting fungal cell wall integrity against oxidative and/or heat stresses [78,79,80].

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], tetraconazole-resistant 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 DMI-fungicide 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 EC50 = 30.2 ± 1.5 mg·L− 1) versus Pi-S (highly sensitive to prochloraz with EC50 = 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.

Results

Transcriptome sequencing and assembly

In the present study, Pi-R and Pi-S were treated with or without DMI-fungicide prochloraz to prepare four RNA-seq 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).

Table 1 FPKM intervals to assess unigene expression level for four P. italicum RNA-seq libraries

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.

Fig. 1
figure1

Hierarchical cluster analysis of differentially expressed genes (DEGs). Blue to red colors represent gene expression levels (i.e., FPKM values from −1 to 1)

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 identified 1) 1052 DEGs between Pi-R-I and Pi-R-NI (614 up-regulated and 438 down-regulated) (Fig. 2a and Additional file 5: Table S3), representing the drug-responsive genes in prochloraz-resistant strain; 2) 298 DEGs between Pi-S-I and Pi-S-NI (63 up-regulated and 235 down-regulated) (Fig. 2b and Additional file 6: Table S4), representing the drug-responsive genes in prochloraz-sensitive strain; 3) 1482 DEGs between Pi-R-I and Pi-S-I (811 up-regulated and 671 down-regulated) (Fig. 2c and Additional file 7: Table S5), representing difference in drug-induced gene expression between fungicide-resistant and -sensitive P. italicum strains; and 4) 958 DEGs between Pi-R-NI and Pi-S-NI (422 up-regulated and 536 down-regulated) (Fig. 2d and Additional file 8: 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 MATE genes) (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.

Fig. 2
figure2

Volcano plot of DEGs in the comparison between Pi-R-I and Pi-R-NI (a), Pi-S-I and Pi-S-NI (b), Pi-R-I and Pi-S-I (c), and Pi-R-NI and Pi-S-NI (d). X-axis indicates log2(fold change) of DEGs between each two samples. Y-axis indicates the -log10(q value) (i.e., corrected p value and abbreviated as qval.) of gene expression variations, and the qval. Was applied to assess statistical significance of the change of unigene expression. The up-regulated, down-regulated, and unchanged unigenes are dotted in red, green, and blue, respectively

Table 2 Analysis of target protein genes associated with azole resistance among identified DEGs
Fig. 3
figure3

Venn diagram of DEGs shared in DEG groups Pi-R-I vs Pi-R-NI and Pi-S-I vs Pi-S-NI (a) and DEG groups Pi-R-I vs Pi-S-I and Pi-R-NI vs Pi-S-NI (b). Yellow circle stands for number of DEGs between Pi-R-I and Pi-S-I (a) and between Pi-R-I and Pi-R-NI (b). Purple circle represents number of DEGs between Pi-R-NI and Pi-S-NI (a) and between Pi-S-I and Pi-S-NI (b). The overlapping region comprises the DEGs shared in the two DEG groups Pi-R-I vs Pi-R-NI and Pi-S-I vs Pi-S-NI (a) and another two DEG groups Pi-R-I vs Pi-S-I and Pi-R-NI vs Pi-S-NI (b)

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 (Fig. 4a), 770 DEGs were enriched into 2005 GO terms without significant enrichment. In the comparison Pi-S-I vs Pi-S-NI (Fig. 4b), 225 DEGs were enriched into 1025 GO terms with 11 significant enrichments (q value ≤0.05), and the top 5 terms significantly enriched were oxidoreductase activity (GO:0016491; q value 1.55E-07), oxidation-reduction process (GO:0055114; q value 3.05E-06), single-organism metabolic process (GO:0044710; q value 2.67E-04), catalytic activity (GO:0003824; q value 1.21E-03), and single-organism process (GO:0044699; q value 4.68E-03). In the comparison Pi-R-I vs Pi-S-I (Fig. 4c), 1086 DEGs were enriched into 2298 GO terms without significant enrichment. In the comparison Pi-R-NI vs Pi-S-NI (Fig. 4d), 711 DEGs were enriched into 1684 GO terms with 11 significant enrichments (q value ≤0.05), and the top 5 terms significantly enriched were oxidoreductase activity (GO:0016491; q value 1.73E-06), oxidation-reduction process (GO:0055114; q value 1.73E-06), hydrolase activity (hydrolyzing O-glycosyl compounds; GO:0004553; q value 1.60E-04), hydrolase activity (acting on glycosyl bonds; GO:0016798; q value 3.50E-04), and transmembrane transport (GO:0055085; q value 8.85E-04). Figure 5 reports the distribution of up- and down-regulated unigenes in the top 30 enriched GO terms for the 4 comparisons mentioned above. Interestingly, the DEGs enriched in the top 30 GO terms were found mostly up-regulated in the comparisons Pi-R-I vs Pi-R-NI and Pi-R-I vs Pi-S-I (Figs. 5a and c) and generally down-regulated in the comparisons Pi-S-I vs Pi-S-NI and Pi-R-NI vs Pi-S-NI (Figs. 5b and d).

Table 3 Summary of GO term distribution
Fig. 4
figure4

Gene ontology (GO) classifications of DEGs 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). For each comparison, GO enrichment classified DEGs into three categories (types) (i.e., biological process, cellular component, and molecular function), as shown in green, orange, and purple bars, respectively. Each GO category (type) displays 30 terms (listed on Y-axis) significantly or most enriched for DEGs in the given comparisons, and X-axis indicates the number of DEGs involved in particular GO term

Fig. 5
figure5

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

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)), drug-target P450 gene (CYP51A, mapped to GO:0055114 (oxidation-reduction process)), steroid biosynthesis-related 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 prochloraz-responsive DEGs mentioned above was summarized in Additional file 9: Table S7.

Table 4 Summary of GO-enriched DEGs associated with prochloraz resistance

Further, KEGG enrichment was applied to identify pathways associating the prochloraz-responsive DEGs with resistance mechanisms. In the present four comparisons, KEGG analysis enriched prochloraz-responsive 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 KEGG-enriched 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., drug-pump genes and drug-target genes, were KEGG-excluded, without exception.

Table 5 Summary of KEGG enrichments for prochloraz-responsive DEGs in the present 4 comparisons

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. italicum response to DMI fungicides, included 1) drug-pump genes: ABC1 (PITC_032590), ABC2 (PITC_006400), MFS1 (PITC_098100), MFS2 (PITC_012240), MFS3 (PITC_056240), and MFS4 (PITC_091150); 2) ergosterol biosynthesis-related genes: CYP51A (PITC_083360), ERG2 (PITC_027000), and ERG6 (PITC_014340); 3) MAPK signaling-related genes: Mkk1 (PITC_088710) and Hog1 (PITC_062470); 4) Ca2+ signal transducer-related genes: CaMK1 (PITC_087700), CaMK2 (PITC_025800), EF-hand1 (PITC_033750), and EF-hand2 (PITC_036760). Additionally, FPKM-based unigene expression quantification combined with local Blast-based annotation revealed differential expression patterns for particular prochloraz-responsive unigenes in the present 4 comparisons, including CYP51B (PITC_064600), CYP51C (PITC_028940), Bck1 (PITC_061930) and Slt2 (PITC_008290) (Additional file 10: Table S8). Considering 1) functional clustering of CYP51A/B/C (i.e., isofroms of drug-target gene CYP51) and 2) cascaded association of Bck1 (encoding MAPKKK), Mkk1 (encoding MAPKK) and Slt2 (encoding MAPK) in Slt2-MAPK pathway, we also performed qRT-PCR validation for the 4 prochloraz-responsive unigenes that were not included in the present DEG list for comparison (i.e., not included in Additional files 5-8: Tables S3–6). As shown in Fig. 6, the qRT-PCR expression patterns of the total 19 prochloraz-responsive DEGs (including 4 FPKM-defined DEGs) were all in agreement with the obtained RNA-seq results. Further, the qRT-PCR results using internal reference gene β-actin were confirmed by another dataset of qRT-PCR analysis based on a different housekeeping gene GAPD (Additional file 11: Figure S3).

Fig. 6
figure6

qPCR validation of 19 prochloraz-responsive DEGs including drug transporter genes (a), ergosterol biosynthesis-related genes (b), MAPK signaling pathway genes (c), and Ca2+ signal transduction genes (d). The housekeeping gene β-actin was used as internal reference to calculate the relative mRNA abundance for the selected unigenes. Relative ratios for the expression of each selected DEG were calculated as Pi-R (I/NI), Pi-S (I/NI), I (Pi-R/Pi-S), and NI (Pi-R/Pi-S). 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)

In detail, the transcript abundance of drug-pump gene ABC1 was strikingly increased in both Pi-R (I/NI) and I (Pi-R/Pi-S), by nearly 500- and 800-folds, respectively, while remarkably decreased in both Pi-S (I/NI) and NI (Pi-R/Pi-S); the similar (but not so strikingly) changing pattern was observed for the rest drug-pump genes including MFS1 (Fig. 6a). When comparing Pi-R (I/NI) with Pi-S (I/NI) or comparing I (Pi-R/Pi-S) with NI (Pi-R/Pi-S), the obviously higher increasing-fold of transcript abundance was also validated for the other prochloraz-responsive genes, including typical drug-target genes (i.e., CYP51A/B/C) (Fig. 6b), ergosterol biosynthesis-related genes ERG2 and ERG6 (Fig. 6b), MAPK signaling-related genes (Fig. 6c), and Ca2+ signal transducer-related genes CaMK1, CaMK2 and EF-hand2 (Fig. 6d). In addition, to functionally verify particular prochloraz-responsive gene, an mfs1-knockout mutant (Δmfs1) was constructed from its parental strain Pi-R, exhibiting obviously lower prochloraz-resistance (i.e., lower prochloraz EC50 value) as compared to the Pi-R wild-type (Additional file 12: Figure S4). This was a sort of preliminary observation from the present RNA-seq analysis, and further biological support is in process.

Discussion

In the past decades, conventional synthetic DMI-fungicides, 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 DMI-resistance in P. italicum species is still not clear, might due to rare opportunity to find highly DMI-resistant P. italicum strain(s). The EC50 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 DMI-fungicides 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 over-expression of specific drug-efflux pumps such as ABC and MFS transporters [8, 42,43,44,45,46,47,48,49,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., ABC12 and MFS14) 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 drug-efflux pumps. Interestingly, regarding these drug-efflux pump-encoding genes, over-expressed in the prochloraz-induced 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 DMI-fungicide 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 drug-pump-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 (MATE13) 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 transporters are more associated with bacterial drug-resistance [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 DMI-fungicides 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), 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 prochloraz-induced 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 up-regulated (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 DMI-resistance. Previous reports have documented azole-induced mRNA increase of ERG2 and ERG6 [39, 81], correlating these ergosterol-biosynthesis genes with acute response of fungal membrane rigidity caused by azole-fungicides. Such ERGs expression profile was observed in the DMI-resistant P. italicum, but not in the prochloraz-resistant P. digitatum strain [11] (Additional file 13: Table S9). These results suggested different strategies in ERGs response of the two Penicillium species to develop DMI-resistance.

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 prochloraz-treated 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 MAPKKK- and 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 DMI-resistance. 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 Ca2+-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 prochloraz-treated 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 (EC50 = 30.2 ± 1.5 mg·L− 1), while Pi-S is prochloraz-sensitive (EC50 = 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 EC50 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 (107 spores mL− 1) as previously described [11]. Afterwards, for each fungal strain, 200 μL of conidial suspension (approximately 2 × 106 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 EC50 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 prochloraz-induced 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 double-stranded 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 Ca2+ 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]. Double-joint 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.

Availability of data and materials

RNA-seq data discussed in this article are publicly available at the NCBI database under accession number PRJNA421419.

Abbreviations

ABC:

ATP-binding cassette

BP:

Biological process

CaM:

Ca2+/calmodulin

CaMK:

Ca2+/calmodulin-dependent kinase

CC:

Cellular component

CWI:

Cell wall integrity

DEG:

Differentially expressed gene

DMI:

Demethylation inhibitor

ERG:

Ergosterol

FPKM:

Fragments Per Kilobase per Million mapped reads

GAPD:

Glyceraldehyde-3-phosphate dehydrogenase

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MAP:

Mitogen-activated protein

MAPK:

Mitogen-activated protein kinase

MAPKK:

Mitogen-activated protein kinase (MAPK) kinase

MAPKKK:

Mitogen-activated protein kinase (MAPK) kinase kinase

MATE:

Multidrug and toxic compound extrusion

MF:

Molecular function

MFS:

Major facilitator superfamily

PCR:

Polymerase chain reaction

PDA:

Potato dextrose agar

PDB:

Potato dextrose broth

qPCR:

Quantitative real-time PCR

RNA-seq:

RNA sequencing

References

  1. 1.

    Eckert JW, Eaks IL. Postharvest disorders and diseases of citrus fruits. Citrus Industry. 1989;5:179–260.

  2. 2.

    Harding PRJ. R23979, a new imidazole derivative effective against postharvest decay of citrus by molds resistant to thiabendazole, benomyl and 2-aminobutane. Plant Dis Rep. 1976;36:37–41.

  3. 3.

    McKay AH, Förster H, Adaskaveg JE. Toxicity and resistance potential of selected fungicides to Galactomyces and Penicillium spp. causing postharvest fruit decays of citrus and other crops. Plant Dis. 2012;96:87–96.

  4. 4.

    Thapa N, Dewdney MM. Efficacy assessment of demethylation inhibitor (DMI) fungicides against Guignardia citricarpa, the causal agent of Citrus black spot (CBS). Phytopathology. 2014;104:117.

  5. 5.

    Erasmus A, Lennox CL, Korsten L, Lesar K, Fourie PH. Imazalil resistance in Penicillium digitatum and P. italicum causing citrus postharvest green and blue mould: impact and options. Postharvest Biol Technol. 2015;107:66–76.

  6. 6.

    Ruan R, Wang M, Liu X, Sun X, Chung KR, Li H. Functional analysis of two sterol regulatory element binding proteins in Penicillium digitatum. PLoS One. 2017;12:e0176485.

  7. 7.

    Kinay P, Mansour MF, Gabler FM, Margosan DA, Smilanick JL. Characterization of fungicide-resistant isolates of Penicillium digitatum, collected in California. Crop Prot. 2007;26:647–56.

  8. 8.

    Nakaune R, Adachi K, Nawata O, Tomiyama M, Akutsu K, Hibi T. A novel ATP-binding cassette transporter involved in multidrug resistance in the phytopathogenic fungus Penicillium digitatum. Appl Environ Microbiol. 1998;64:3983–8.

  9. 9.

    Hamamoto H, Hasegawa K, Nakaune R, Lee YJ, Makizumi Y, Akutsu K, et al. Tandem repeat of a transcriptional enhancer upstream of the sterol 14α-demethylase gene (CYP51) in Penicillium digitatum. Appl Environ Microbiol. 2000;66:3421–6.

  10. 10.

    Hamamoto H, Hasegawa K, Nakaune R, Lee YJ, Akutsu K, Hibi T. PCR-based detection of sterol demethylation inhibitor-resistant strains of Penicillium digitatum. Pest Manag Sci. 2001;57:839–43.

  11. 11.

    Liu J, Wang S, Qin T, Li N, Niu Y, Li D, et al. Whole transcriptome analysis of Penicillium digitatum strains treatmented with prochloraz reveals their drug-resistant mechanisms. BMC Genomics. 2015;16:855.

  12. 12.

    Ma H, Sun X, Wang M, Gai Y, Chung KR, Li H. The citrus postharvest pathogen Penicillium digitatum depends on the PdMpkB kinase for developmental and virulence functions. Int J Food Microbiol. 2016;236:167–76.

  13. 13.

    Ramón-Carbonell MD, Sánchez-Torres P. The transcription factor PdSte12, contributes to Penicillium digitatum, virulence during citrus fruit infection. Postharvest Biol Technol. 2017;125:129–39.

  14. 14.

    Bossche HV, Marichal P, Gorrens J, Bellens D, Moereels H, Janssen PAJ. Mutation in cytochrome P-450-dependent 14α-demethylase results in decreased affinity for azole antifungals. Biochem Soc Trans. 1990;18:56–9.

  15. 15.

    Tyndall JD, Sabherwal M, Sagatova AA, Keniya MV, Negroni J, Wilson RK, et al. Structural and functional elucidation of yeast lanosterol 14α-demethylase in complex with agrochemical antifungals. PLoS One. 2016;11:e0167485.

  16. 16.

    Sagatova AA, Keniya MV, Tyndall JD, Monk BC. The impact of homologous resistance mutations from pathogenic yeast on Saccharomyces cerevisiae lanosterol 14α-demethylase. Antimicrob Agents Chemother. 2017;62(3):e02242–17.

  17. 17.

    Sanglard D, Ischer F, Koymans L, Bille J. Amino acid substitutions in the cytochrome P-450 lanosterol 14α-demethylase (CYP51A1) from azole-resistant Candida albicans clinical isolates contribute to resistance to azole antifungal agents. Antimicrob Agents Chemother. 1998;42:241–53.

  18. 18.

    Marichal P, Koymans L, Willemsens S, Bellens D, Verhasselt P, Luyten W, et al. Contribution of mutations in the cytochrome P450 14α-demethylase (Erg11p, Cyp51p) to azole resistance in Candida albicans. Microbiology. 1999;145:2701–13.

  19. 19.

    Alvarez-Rueda N, Fleury A, Logé C, Pagniez F, Robert E, Morio F, Le Pape P. The amino acid substitution N136Y in Candida albicans sterol 14alpha-demethylase is involved in fluconazole resistance. Med Mycol. 2016;54:764–75.

  20. 20.

    Wu Y, Gao N, Li C, Gao J, Ying C. A newly identified amino acid substitution T123I in the 14α-demethylase (Erg11p) of Candida albicans confers azole resistance. FEMS Yeast Res. 2017;17.

  21. 21.

    Edlind TD, Henry KW, Metera KA, Katiyar SK. Aspergillus fumigatus cyp51 sequence: potential basis for fluconazole resistance. Med Mycol. 2001;39:299–302.

  22. 22.

    Camps SM, Rijs AJ, Klaassen CH, Meis JF, O'Gorman CM, Dyer PS, et al. Molecular epidemiology of Aspergillus fumigatus isolates harboring the TR34/L98H azole resistance mechanism. J Clin Microbiol. 2012;50:2674–80.

  23. 23.

    Paul RA, Rudramurthy SM, Meis JF, Mouton JW, Chakrabarti A. A novel Y319H substitution in CYP51C associated with azole resistance in Aspergillus flavus. Antimicrob Agents Chemother. 2015;59:6615–9.

  24. 24.

    Cools HJ, Hawkins NJ, Fraaije BA. Constraints on the evolution of azole resistance in plant pathogenic fungi. Plant Pathol. 2013;62:36–42.

  25. 25.

    Price CL, Parker JE, Warrilow AG, Kelly DE, Kelly SL. Azole fungicides - understanding resistance mechanisms in agricultural fungal pathogens. Pest Manag Sci. 2015;71:1054–8.

  26. 26.

    Lichtemberg PSF, Luo Y, Morales RG, Jmm F, Michailides TJ, Ll MDM. The point mutation G461S in the MfCYP51 gene is associated with tebuconazole resistance in Monilinia fructicola populations in Brazil. Phytopathology. 2017;107:1507–14.

  27. 27.

    Wang J, Yu J, Liu J, Yuan Y, Li N, He M, et al. Novel mutations in CYP51B from Penicillium digitatum involved in prochloraz resistance. J Microbiol. 2014;52:762–70.

  28. 28.

    Agarwal AK, Rogers PD, Baerson SR, Jacob MR, Barker KS, Cleary JD, et al. Genome-wide expression profiling of the response to polyene, pyrimidine, azole, and echinocandin antifungal agents in Saccharomyces cerevisiae. J Biol Chem. 2003;278:34998–5015.

  29. 29.

    Da SFM, Malavazi I, Savoldi M, Brakhage AA, Goldman MH, Kim HS, et al. Transcriptome analysis of Aspergillus fumigatus exposed to voriconazole. Curr Genet. 2006;50:32–44.

  30. 30.

    Liu X, Jiang J, Shao J, Yin Y, Ma Z. Gene transcription profiling of Fusarium graminearum treated with an azole fungicide tebuconazole. Appl Microbiol Biotechnol. 2010;85:1105–14.

  31. 31.

    Sun N, Fonzi W, Chen H, She X, Zhang L, Zhang L, Calderone R. Azole susceptibility and transcriptome profiling in Candida albicans mitochondrial electron transport chain complex I mutants. Antimicrob Agents Chemother. 2013;57:532–42.

  32. 32.

    Teymuri M, Mamishi S, Pourakbari B, Mahmoudi S, Ashtiani MT, Sadeghi RH, et al. Investigation of ERG11 gene expression among fluconazole–resistant Candida albicans: first report from an Iranian referral paediatric hospital. Br J Biomed Sci. 2015;72:28–31.

  33. 33.

    Ghosoph JM, Schmidt LS, Margosan DA, Smilanick JL. Imazalil resistance linked to a unique insertion sequence in the PdCYP51 promoter region of Penicillium digitatum. Posthavest Biol Technol. 2007;44:9–18.

  34. 34.

    Chen F, Lin D, Wang J, Li B, Duan H, Liu J, Liu X. Heterologous expression of the Monilinia fructicola CYP51 (MfCYP51) gene in Pichia pastoris confirms the mode of action of the novel fungicide, SYP-Z048. Front Microbiol. 2015;6:457.

  35. 35.

    Hernández A, Serrano-Bueno G, Perez-Castiñeira JR, Serrano A. 8-Dehydrosterols induce membrane traffic and autophagy defects through V-ATPase dysfunction in Saccharomyces cerevisae. Biochim Biophys Acta. 1853;2015:2945–56.

  36. 36.

    Luna-Tapia A, Peters BM, Eberle KE, Kerns ME, Foster TP, Marrero L, et al. Erg2 and erg24 are required for normal vacuolar physiology as well as Candida albicans pathogenicity in a murine model of disseminated but not vaginal candidiasis. Eukaryot Cell. 2015;14:1006–16.

  37. 37.

    Deng GM, Yang QS, He WD, Li CY, Yang J, Zuo CW, et al. Proteomic analysis of conidia germination in Fusarium oxysporum, f. sp. cubense, tropical race 4 reveals new targets in ergosterol biosynthesis pathway for controlling Fusarium wilt of banana. Appl Microbiol Biotechnol. 2015;99:7189–207.

  38. 38.

    Lendenmann MH, Croll D, Mcdonald BA. QTL mapping of fungicide sensitivity reveals novel genes and pleiotropy with melanization in the pathogen Zymoseptoria tritici. Fungal Genet Biol. 2015;80:53–67.

  39. 39.

    Gu X, Xue W, Yin Y, Liu H, Li S, Sun X. The Hsp90 co-chaperones Sti1, Aha1, and P23 regulate adaptive responses to antifungal azoles. Front Microbiol. 2016;7:1571.

  40. 40.

    Ouyang Q, Tao N, Jing G. Transcriptional profiling analysis of Penicillium digitatum, the causal agent of citrus green mold, unravels an inhibited ergosterol biosynthesis pathway in response to citral. BMC Genomics. 2016;17:599.

  41. 41.

    Abe F, Hiraki T. Mechanistic role of ergosterol in membrane rigidity and cycloheximide resistance in Saccharomyces cerevisiae. Biochim Biophys Acta. 1788;2009:743–52.

  42. 42.

    Saidane S, Weber S, Deken XD, St-Germain G, Raymond M. PDR16-mediated azole resistance in Candida albicans. Mol Microbiol. 2006;60:1546–62.

  43. 43.

    Tsao S, Rahkhoodaee F, Raymond M. Relative contributions of the Candida albicans ABC transporters Cdr1p and Cdr2p to clinical azole resistance. Antimicrob Agents Chemother. 2009;53:1344–52.

  44. 44.

    Jr BL, Gast CE, Mao Y, Wong B. Fluconazole transport into Candida albicans secretory vesicles by the membrane proteins Cdr1p, Cdr2p, and Mdr1p. Eukaryot Cell. 2010;9:960–70.

  45. 45.

    Sanglard D, Coste AT. Activity of isavuconazole and other azoles against Candida clinical isolates and yeast model systems with known azole resistance mechanisms. Antimicrob Agents Chemother. 2015;60:229–38.

  46. 46.

    Liu Z, Myers LC. Mediator tail module is required for Tac1 activated CDR1 expression and azole resistance in Candida albicans. Antimicrob Agents Chemother. 2017;61:e01342–17.

  47. 47.

    Pourakbari B, Teymuri M, Mahmoudi S, Valian SK, Movahedi Z, Eshaghi H, Mamishi S. Expression of major efflux pumps in fluconazole-resistant Candida albicans. Infect Disord Drug Targets. 2017;17:178–84.

  48. 48.

    Fraczek MG, Bromley M, Buied A, Moore CB, Rajendran R, Rautemaa R, et al. The cdr1B efflux transporter is associated with non-cyp51a-mediated itraconazole resistance in Aspergillus fumigatus. J Antimicrob Chemother. 2013;68:1486–96.

  49. 49.

    Paul S, Diekema D, Moye-Rowley WS. Contributions of Aspergillus fumigatus ATP-binding cassette transporter proteins to drug resistance and virulence. Eukaryot Cell. 2013;12:1619–28.

  50. 50.

    Hagiwara D, Miura D, Shimizu K, Paul S, Ohba A, Gonoi T, et al. A novel Zn2-Cys6 transcription factor AtrR plays a key role in an azole resistance mechanism of Aspergillus fumigatus by co-regulating cyp51A and cdr1B expressions. PLoS Pathog. 2017;13:e1006096.

  51. 51.

    Sorbo GD, Andrade AC, Nistelrooy JGMV, Kan JALV, Balzi E, Waard MAD. Multidrug resistance in Aspergillus nidulans involves novel ATP-binding cassette transporters. Mol Gen Genet. 1997;254:417–26.

  52. 52.

    Semighini CP, Marins M, Goldman MH, Goldman GH. Quantitative analysis of the relative transcript levels of ABC transporter Atr genes in Aspergillus nidulans by real-time reverse transcription-PCR assay. Appl Environ Microbiol. 2002;68:1351–7.

  53. 53.

    Nakaune R, Hamamoto H, Imada J, Akutsu K, Hibi T. A novel ABC transporter gene, PMR5, is involved in multidrug resistance in the phytopathogenic fungus Penicillium digitatum. Mol Gen Genomics. 2002;267:179–85.

  54. 54.

    Sánchez-Torres P, Tuset JJ. Molecular insights into fungicide resistance in sensitive and resistant Penicillium digitatum, strains infecting citrus. Postharvest Biol Technol. 2011;59:159–65.

  55. 55.

    Prasad R, Kapoor K. Multidrug resistance in yeast Candida. Int Rev Cytol. 2005;242:215–48.

  56. 56.

    Roohparvar R, De Waard MA, Kema GH, Zwiers LH. MgMfsl, a major facilitator superfamily transporter from the fungal wheat pathogen Mycosphaerella graminicola, is a strong protectant against natural toxic compounds and fungicides. Fungal Genet Biol. 2007;44:378–88.

  57. 57.

    Wang JY, Sun XP, Lin LY, Zhang TY, Ma ZH, Li HY. PdMfs1, a major facilitator superfamily transporter from Penicillium digitatum, is partially involved in the imazalil-resistance and pathogenicity. Afr J Microbiol Res. 2012;6:95–105.

  58. 58.

    Wu Z, Wang S, Yuan Y, Zhang T, Liu J, Liu D. A novel major facilitator superfamily transporter in Penicillium digitatum (PdMFS2) is required for prochloraz resistance, conidiation and full virulence. Biotechnol Lett. 2016;38:1349–57.

  59. 59.

    Kuroda T, Tsuchiya T. Multidrug efflux transporters in the MATE family. Biochim Biophys Acta. 1794;2009:763–8.

  60. 60.

    Schindler BD, Kaatz GW. Multidrug efflux pumps of gram-positive bacteria. Drug Resist Updat. 2016;27:1–13.

  61. 61.

    Lai SJ, Tu IF, Wu WL, Yang JT, Luk LY, Lai MC, et al. Site-specific his/asp phosphoproteomic analysis of prokaryotes reveals putative targets for drug resistance. BMC Microbiol. 2017;17:123.

  62. 62.

    Schlunk I, Krause K, Wirth S, Kothe E. A transporter for abiotic stress and plant metabolite resistance in the ectomycorrhizal fungus Tricholoma vaccinum. Environ Sci Pollut R. 2015;22:19384–93.

  63. 63.

    Hayes BM, Anderson MA, Traven A, van der Weerden NL, Bleackley MR. Activation of stress signalling pathways enhances tolerance of fungi to chemical fungicides and antifungal proteins. Cell Mol Life Sci. 2014;71:2651–66.

  64. 64.

    Sanz AB, García R, Rodríguez-Peña JM, Arroyo J. The CWI pathway: regulation of the transcriptional adaptive response to cell wall stress in yeast. J Fungi. 2017;41:1.

  65. 65.

    Levin DE. Cell wall integrity signaling in Saccharomyces cerevisiae. Microbiol Mol Biol Rev. 2005;69:262–91.

  66. 66.

    Saito H. Regulation of cross-talk in yeast MAPK signaling pathways. Curr Opin Microbiol. 2010;13:677–83.

  67. 67.

    Nishida N, Jing D, Kuroda K, Ueda M. Activation of signaling pathways related to cell wall integrity and multidrug resistance by organic solvent in Saccharomyces cerevisiae. Curr Genet. 2014;60:149–62.

  68. 68.

    Chen LH, Yang SL, Chung KR. Resistance to oxidative stress via regulating siderophore-mediated iron acquisition by the citrus fungal pathogen Alternaria alternata. Microbiology. 2014;160:970–9.

  69. 69.

    Yang SL, Yu PL, Chung KR. The glutathione peroxidase-mediated reactive oxygen species resistance, fungicide sensitivity and cell wall construction in the citrus fungal pathogen Alternaria alternata. Environ Microbiol. 2016;18:923–35.

  70. 70.

    Yu PL, Chen LH, Chung KR. How the pathogenic fungus Alternaria alternata copes with stress via the response regulators SSK1 and SHO1. PLoS One. 2016;11:e0149153.

  71. 71.

    Chen LH, Tsai HC, Yu PL, Chung KR. A major facilitator superfamily transporter-mediated resistance to oxidative stress and fungicides requires Yap1, Skn7, and MAP kinases in the citrus fungal pathogen Alternaria alternata. PLoS One. 2017;12:e0169103.

  72. 72.

    Wang M, Chen C, Zhu C, Sun X, Ruan R, Li H. Os2 MAP kinase-mediated osmostress tolerance in Penicillium digitatum is associated with its positive regulation on glycerol synthesis and negative regulation on ergosterol synthesis. Microbiol Res. 2014;169:511–21.

  73. 73.

    Ramón-Carbonell M, Sánchez-Torres P. PdSlt2 Penicillium digitatum mitogen-activated-protein kinase controls sporulation and virulence during citrus fruit infection. Fungal Biol. 2017;121:1063–74.

  74. 74.

    Segmueller N, Ellendorf U, Tudzynski B, Tudzynski P. BcSAK1, a stress-activated mitogen-activated protein kinase, is involved in vegetative differentiation and pathogenicity in Botrytis cinerea. Eukaryot Cell. 2007;6:211–21.

  75. 75.

    Liu WW, Soulie MC, Perrino C, Fillinger S. The osmosensing signal transduction pathway from Botrytis cinerea regulates cell wall integrity and MAP kinase pathways control melanin biosynthesis with influence of light. Fungal Genet Biol. 2011;48:377–87.

  76. 76.

    Nguyen TV, Schafer W, Bormann J. The stress-activated protein kinase FgOS-2 is a key regulator in the life cycle of the cereal pathogen Fusarium graminearum. Mol Plant Microbe In. 2012;25:1142–56.

  77. 77.

    Dudgeon DD, Zhang NN, Ositelu OO, Kim H, Cunningham KW. Nonapoptotic death of Saccharomyces cerevisiae cells that is stimulated by Hsp90 and inhibited by calcineurin and Cmk2 in response to endoplasmic reticulum stresses. Eukaryot Cell. 2008;7(12):2037–51.

  78. 78.

    Rodriguez-Caban J, Gonzalez-Velazquez W, Perez-Sanchez L, Gonzalez-Mendez R, Rodriguez-del-Valle N. Calcium/calmodulin kinase1 and its relation to thermotolerance and HSP90 in Sporothrix schenckii: an RNAi and yeast two-hybrid study. BMC Microbiol. 2011;11:162.

  79. 79.

    Ding XH, Yu QL, Zhang B, Xu N, Jia C, Dong YJ, et al. The type II Ca2+/calmodulin-dependent protein kinases are involved in the regulation of cell wall integrity and oxidative stress response in Candida albicans. Biochem Biophys Res Commun. 2014;446:1073–8.

  80. 80.

    Kumar R, Tamuli R. Calcium/calmodulin-dependent kinases are involved in growth, thermotolerance, oxidative stress survival, and fertility in Neurospora crassa. Arch Microbiol. 2014;196:295–305.

  81. 81.

    da Silva Ferreira ME, Malavazi I, Savoldi M, Brakhage A, Goldman MH, Kim HS, et al. Transcriptome analysis of Aspergillus fumigatus exposed to voriconazole. Curr Genet. 2006;50:32–44.

  82. 82.

    Bolton MD, Ebert MK, Faino L, Rivera-Varas V, de Jonge R, de Peer YV, et al. RNA-sequencing of Cercospora beticola DMI-sensitive and -resistant isolates after treatment with tetraconazole identifies common and contrasting pathway induction. Fungal Genet Biol. 2016;92:1–13.

  83. 83.

    Hellin P, King R, Urban M, Hammond-Kosack KE, Legrève A. The adaptation of Fusarium culmorum to DMI fungicides is mediated by major transcriptome modifications in response to azole fungicide, including the overexpression of a PDR transporter (FcABC1). Front Microbiol. 2018;9:1385.

  84. 84.

    Whaley SG, Caudle KE, Simonicova L, Zhang Q, Moye-Rowley WS, Rogers PD. Jjj1 is a negative regulator of Pdr1-mediated fluconazole resistance in Candida glabrata. mSphere. 2018;3:e00466–17.

  85. 85.

    Holmes GJ, Eckert JW. Sensitivity of Penicillium digitatum and P. italicum to postharvest citrus fungicides in California. Phytopathology. 1999;89(9):716–21.

  86. 86.

    Bus VG, Bongers AJ, Risse LA. Occurrence of Penicillium digitatum and P. italicum resistant to benomyl, thiabendazole, and imazalil on citrus fruit from different geographic origins. Plant Dis. 1991;75:1098–100.

  87. 87.

    Bus VG. ED50 levels of Penicillium digitatum and P. italicum with reduced sensitivity to thiabendazole, benomyl and imazalil. Postharvest Biol Technol. 1992;1:305–15.

  88. 88.

    Del Sorbo G, Schoonbeek H, De Waard MA. Fungal transporters involved in efflux of natural toxic compounds and fungicides. Fungal Genet Biol. 2000;30:1–15.

  89. 89.

    Campoy S, Adrio JL. Antifungals. Biochem Pharmacol. 2017;133:86–96.

  90. 90.

    Irie K, Takase M, Lee KS, Levin DE, Araki H, Matsumoto K, Oshima Y. MKK1 and MKK2, which encode Saccharomyces cerevisiae mitogen-activated protein kinase-kinase homologs, function in the pathway mediated by protein kinase C. Mol Cell Biol. 1993;13:3076–83.

  91. 91.

    Reinoso-Martín C, Schüller C, Schuetzer-Muehlbauer M, Kuchler K. The yeast protein kinase C cell integrity pathway mediates tolerance to the antifungal drug caspofungin through activation of Slt2p mitogen-activated protein kinase signaling. Eukaryot Cell. 2003;6(2):1200–10.

  92. 92.

    Dirr F, Echtenacher B, Heesemann J, Hoffmann P, Ebel F, Wagener J. AfMkk2 is required for cell wall integrity signaling, adhesion, and full virulence of the human pathogen Aspergillus fumigatus. Int J Med Microbiol. 2010;300:496–502.

  93. 93.

    Mehrabi R, van der Lee T, Waalwijk C, Kema G. MgSlt2, a cellular integrity MAP kinase gene of the fungal wheat pathogen Mycosphaerella graminicola, is dispensable for penetration but essential for invasive growth. Mol Plant Microbe In. 2006;19(4):389–98.

  94. 94.

    Cyert MS. Genetic analysis of calmodulin and its targets in Saccharomyces cerevisiae. Annu Rev Genet. 2001;35:647–72.

  95. 95.

    Ballester A, Marcet-Houben M, Levin E, Sela N, Selma-Lázaro C, Carmona L, et al. Genome, transcriptome, and functional analyses of Penicillium expansum provide new insights into secondary metabolism and pathogenicity. Mol Plant Microbe In. 2015;28(3):232–48.

  96. 96.

    Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.

  97. 97.

    Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008;18(9):1509–17.

  98. 98.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

  99. 99.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple hypothesis testing. J R Stat Soc B. 1995;57:289–300.

  100. 100.

    de Hoon MJL, Imoto S, Nolan J, Miyano S. Open source clustering software. Bioinformatics. 2004;20:1453–4.

  101. 101.

    Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

  102. 102.

    Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21:3787–93.

  103. 103.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2–ΔΔCt method. Methods. 2001;25:402–8.

  104. 104.

    Zhao S, Yan Y, He Q, Yang L, Yin X, Li C, et al. Comparative genomic, transcriptomic and secretomic profiling of Penicillium oxalicum HP7-1 and its cellulase and xylanase hyper-producing mutant EU2106, and identification of two novel regulatory genes of cellulase and xylanase gene expression. Biotechnol Biofuels. 2016;9:203. https://doi.org/10.1186/s13068-016-0616-9.

Download references

Acknowledgements

We sincerely thank Dr. Huazhong Shi (Texas Tech University, USA) for his help in English language revision.

Funding

This work was supported by the National Natural Science Foundations of China (No. 31371893 and 31071653), the Natural Science Fund of Hubei Province (No. 2018CFB676), and the Fundamental Research Funds for the Central Universities (No. CCNU19TS079). The funding body didn’t play any roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

YY and DL designed this study, drafted the manuscript, acquired funds, and supervised all research activities. DL, TZ and NL isolated and characterized the P. italicum strains. NL and TZ made EC50 assays. TZ and QC prepared RNA-seq samples. QC and YY performed bioinformatics analysis and figure production. TZ, DL and YY designed and conducted qRT-PCR experiments. YY substantially revised the manuscript. All authors read and approved this final manuscript.

Correspondence to Deli Liu or Yongze Yuan.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1. PDA-based EC50 determination against prochloraz for two Penicillium italicum strains (Pi-R and Pi-S) with different response to the DMI prochloraz.

Additional file 2: Table S1. EC50 values of Pi-R and Pi-S against the fungicides belonging to DMI, benzimidazole and phenylpyrrole classes.

Additional file 3: Table S2. Primers used in the present qRT-PCR validation.

Additional file 4: Figure S2. Genome mapping analysis for the clean reads from the present 4 RNA-seq libraries, i.e., Pi-R-I (a), Pi-R-NI (b), Pi-S-I (c), and Pi-S-NI (d).

Additional file 5: Table S3. Summary of DEGs in the comparison between prochloraz-induced and no-induced Pi-R strains.

Additional file 6: Table S4. Summary of DEGs in the comparison between prochloraz-induced and no-induced Pi-S strains.

Additional file 7: Table S5. Summary of DEGs in the comparison between Pi-R and Pi-S strains under prochloraz-induced conditions.

Additional file 8: Table S6. Summary of DEGs in the comparison between Pi-R and Pi-S strains without prochloraz induction.

Additional file 9: Table S7. Summary of distribution (hit records) of 19 selected DEGs towards GO enrichments in the present 4 comparisons.

Additional file 10: Table S8. FPKM values of the additional 4 prochloraz-responsive unigenes in the present comparisons.

Additional file 11: Figure S3. GAPD-based qPCR validation of 19 prochloraz-responsive DEGs including drug transporter genes (a), ergosterol biosynthesis-related genes (b), MAPK signaling pathway genes (c), and Ca2+ signal transduction genes (d).

Additional file 12: Figure S4. Comparison of prochloraz EC50 values between Pi-R (wild type) and its mfs1-knockout mutant (Δmfs1).

Additional file 13: Table S9. Comparative analysis of prochloraz-responsive gene expression profiles between the P. italicum and P. digitatum transcriptomes.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zhang, T., Cao, Q., Li, N. et al. Transcriptome analysis of fungicide-responsive gene expression profiles in two Penicillium italicum strains with different response to the sterol demethylation inhibitor (DMI) fungicide prochloraz. BMC Genomics 21, 156 (2020). https://doi.org/10.1186/s12864-020-6564-6

Download citation

Keywords

  • Transcriptome
  • Penicillium italicum
  • Demethylation inhibitor (DMI)-resistance
  • Prochloraz-responsive genes