Genome-wide transcriptional analyses of Clarireedia jacksonii isolates associated with multi-drug resistance

Emerging multidrug resistance (MDR) in Clarireedia spp. is a huge challenge to the management of dollar spot (DS) disease on turfgrass. Insight into the molecular basis of resistance mechanisms may help identify key molecular targets for developing novel effective chemicals. Previously, a MDR isolate (LT586) of C. jacksonii with significantly reduced sensitivities to propiconazole, boscalid, and iprodione, and a fungicide-sensitive isolate (LT15) of the same species were isolated from creeping bentgrass (Agrostis stolonifera L.). The present study aimed to further explore the molecular mechanisms of resistance by using genome-wide transcriptional analyses of the two isolates. A total of 619 and 475 differentially expressed genes (DEGs) were significantly down and upregulated in the MDR isolate LT586, compared with the sensitive isolate LT15 without fungicide treatment. Three hundreds and six and 153 DEGs showed significantly lower and higher expression in the MDR isolate LT586 than those in the sensitive isolate LT15, which were commonly induced by the three fungicides. Most of the 153 upregulated DEGs were xenobiotic detoxification-related genes and genes with transcriptional functions. Fifty and 17 upregulated DEGs were also commonly observed in HRI11 (a MDR isolate of the C. jacksonii) compared with the HRS10 (a fungicide-sensitive isolate of same species) from a previous study without and with the treatment of propiconazole, respectively. The reliability of RNA-seq data was further verified by qRT-PCR method using a few select potentially MDR-related genes. Results of this study indicated that there were multiple uncharacterized genes, possibly responsible for MDR phenotypes in Clarireedia spp., which may have important implications in understanding the molecular mechanisms underlying MDR resistance.


Introduction
Dollar spot (DS) is one of the most economically important diseases of turfgrass worldwide (Salgado-Salazar et al., 2018). Multiple species in the genus Clarireedia (formerly Sclerotinia homoeocarpa) have been reported to cause DS, including C. homoeocarpa, C. bennettii, C. jacksonii, C. monteithian, C. paspali, and C. hainanense (Salgado-Salazar et al., 2018;Hu et al., 2019a,b;Zhang et al., 2022). Management of DS requires repeated inputs of fungicides with OPEN ACCESS EDITED BY Vijay Soni, NewYork-Presbyterian, United States REVIEWED BY different chemical groups. Benzimidazole carbamates (MBCs), dicarboximides (DCFs), 14α-demethylase inhibitors (DMIs), and succinate dehydrogenase inhibitors (SDHIs) are among the most frequently used fungicide classes for DS control in the last few decades Stephens and Kaminski, 2019). Due to widespread, repeated uses of fungicidal compounds as well as many other factors such as management location, environmental conditions, genetic regulation, and turfgrass species (Koch et al., 2009;Latin, 2011;Stephens and Kaminski, 2019), resistance to these chemical groups in Clarireedia spp. have been reported in many regions around the world (Miller et al., 2002;Jo et al., 2006;Koch et al., 2009;Lee et al., 2017). Moreover, multiple fungicide resistant (MFR) or multi-drug resistant (MDR) populations have emerged in the United States of America, Japan, as well as China, recently (Sang et al., 2015Zhang et al., 2021), which makes DS disease increasingly difficult to control.
In the genus Clarireedia, resistance mechanisms to different chemical groups of fungicides have been systemically revealed in C. jacksonii. Hulvey et al. (2012) found that overexpression of two genes, ShCYP51B and ShatrD in the DS pathogen C. jacksonii, was responsible for reduced DMI fungicide sensitivity. Sang et al. (2016) found that nonsynonymous polymorphisms in codon 366 (isoleucine to asparagine) in histidine kinase gene (Shos1) and overexpression of ShPDR1 led to DCFs resistance in C. jacksonii isolates. Hu et al. (2018) found that nonsynonymous polymorphisms at codon 198 and 200 in β-tubulin were responsible for high and medium resistance to MBC fungicide thiophanate-methyl, respectively. The H267Y mutation in SdhB was determined to be responsible for resistance to boscalid and penthiopyrad in C. jacksonii isolates from Japan, and the G91R mutation in SdhC was confirmed to be a direct factor conferring resistance to boscalid, fluxapyroxad, isofetamid, and penthiopyrad in C. jacksonii isolates from the United States (Popko et al., 2018).
There are two different types of resistance to multiple fungicides: MFR and MDR. Generally, MFR is often caused by accumulation of different mutations at the corresponding fungicide target genes (e.g., CYP51B,os1 and Sdh genes in Clarireedia spp.), while the mechanisms of MDR in plant pathogenic fungi were often associated with overexpression of single or multiple drug efflux transporters, including ATP-binding cassette (ABC) and major facilitator superfamily (MFS) transporters (Nakaune et al., 2002;Omrane et al., 2015;Sang et al., 2015). In C. jacksonii, Sang et al. (2015) demonstrate that MDR is partially attributable to the pleiotropic drug resistance (PDR) transporter gene, ShPDR1. Furthermore, they found that a gain-of-function mutation (ATG to ACG at codon 853, M853T) in the fungus-specific transcription factor ShXDR1 led to constitutive and induced overexpression of cytochrome P450s (CYP450s) and ABC transporters, which are responsible for MDR . In the previous study, we found that multiple isolates of C. jacksonii showed decreased sensitivities to fungicides with different mode of actions . Preliminary studies revealed that no consistent mutations were found on the fungicide target genes in the resistant isolates, which suggested that they were MDR rather than MFR isolates. To further characterize if similar mechanisms of MDR responsible for the above-mentioned resistant isolates of C. jacksonii found in China, we sequenced ShXDR1 gene in the sensitive and resistant isolates, but no M853T mutation on ShXDR1 was found. The mechanisms of MDR in C. jacksonii or other species in Clarireedia may be complex and need to be thoroughly investigated.
Comparative transcriptomics is a powerful tool to elucidate the molecular mechanisms of fungicide resistance. RNA-seq analyses have been used to uncover the molecular mechanisms of resistance to DMI, SDHI, or MDR in many pathogenic fungi, including Clarireedia spp., Aspergillus fumigates (Fan et al., 2021), Fusarium virguliforme (Sang et al., 2021), and Cercospora beticola (Bolton et al., 2016). In this study, the goal is to identify the genes that are differentially expressed between fungicide-sensitive and presumably MDR isolates of C. jacksonii under exposures to three different fungicides, propiconazole, boscalid, and iprodione. Comparative transcriptomics was applied to glean insight into the potential variations at transcriptomic level. The results provide additional molecular clues that will aid in elucidating the mechanisms and metabolic pathways related to MDR in Clarireedia spp.

Fungal isolates
Fungal isolates LT15 and LT586 were recovered from symptomatic leaf blades of creeping bentgrass (Agrostis stolonifera L.). Sensitivities of the two isolates to iprodione, propiconazole and boscalid were determined in vitro using the respective discriminatory concentrations described in the previous study . Agar plugs (5 mm in diameter) were inoculated on PDA plates amended with iprodione (1.0 μg/mL), propiconazole (0.1 μg/mL), and boscalid (10 μg/mL). The non-amended PDA plates served as controls. After 3 days' incubation, the diameter of colonies was perpendicularly measured twice. The percentage (%) of relative mycelia growth was calculated as the average diameter on fungicide-amended medium divided by the average diameter on non-amended medium and multiplied by 100.
There are several species in the genus of Clarireedia. The two isolates tested here were identified to the species level, C. jacksonii by sequencing the nuclear ribosomal internal transcribed spacer (ITS) with the universal primer of ITS4 and ITS5 (Supplementary Table S1). ITS has proved useful to differentiate different species in the genus Clarireedia in previous studies, and the phylogenetic tree constructed from ITS produced a topology like the combined multi-locus dataset (Salgado-Salazar et al., 2018;Aynardi et al., 2019). Maximum likelihood method was used for constructing phylogenetic tree of select known and published isolates using the Tamura 3-parameter model in MEGA7 and tested with 1,000 bootstrap replicates (Kumar et al., 2016;Hu et al., 2019a,b).

Mycelial preparation and total RNA extraction
The mycelia of the two isolates were prepared for RNA extraction with the method described previously with some modifications (Hulvey et al., 2012). In brief, isolates were inoculated into 25 mL of potato dextrose broth (PDB) with eight agar plugs (5 mm in diameter) and grown for 3 days at 25°C and 150 rpm. Propiconazole (EC, Banner Maxx, Syngenta Crop Protection), iprodione (EC, Rovral, FMC), and boscalid (WG, Cantus, BASF) were added to the tubes with the final concentrations at 0.1, 1, and 10 μg/mL, respectively Zhang et al., 2021). The untreated samples were added with same Frontiers in Microbiology 03 frontiersin.org volume of sterile water . Tubes were lightly shaken on a benchtop shaker for 1 h. Approximately 100 mg of mycelia for each sample was harvested by vacuum filtration, transferred into a 2-mL screw cap tube, and immediately dropped into liquid nitrogen. Frozen samples were stored at −80°C freezer until further processing. Three biological replicates of each treatment were generated. Total RNA was extracted from each sample (a total of 24 samples) using RNAsimple Total RNA Kit (TIANGEN Biotech, Beijing).

Transcriptome sequencing, assembly, and annotation
Total RNA of each sample was sent to Biomarker Technologies (Beijing, China) for cDNA library preparation and sequencing. RNA concentration and purity were measured using NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, United States). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). 1 μg RNA per sample was used for constructing cDNA library. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (New England Biolabs, United States) following manufacturer's instructions. RNA-seq was performed on an Illumina HiSeq4000 platform (Illumina, CA, United States). After sequencing, reads containing adapter, poly-N or low quality reads were removed through in-house perl scripts 1 to produce clean data sets. Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated. High quality clean data sets were mapped to the C. jacksonii reference genome (Zhang et al., 2023) with Hisat2 software (Kim et al., 2015). Only reads with a perfect match or one mismatch were further analyzed and annotated.

Functional analysis of differentially expressed genes
Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). The formula for FPKM calculation was: FPKM = (cDNA fragments)/mapped fragments (millions) × transcript length (kb). Differential expressed genes (DEGs) were analyzed by using the R package DESeq2 (Love et al., 2014). The value of p < 0.05, |log 2 Fold Change| ≥ 1.5 and FPKM > 0.05 were set as the threshold for the identifications of DEGs. Principal component analysis (PCA) and heat map analysis were performed using BMKCloud. 2 Enrichment analysis of the DEGs was implemented by the R package GOseq based on Wallenius non-central hypergeometric distribution (https://cran.r-project.org/; Young et al.,

Quantitative real time RT-PCR validation of RNA-seq data
The relative expressions of selected DEGs were assayed with quantitative real-time PCR (qRT-PCR) before and after treatments with propiconazole, iprodione, and boscalid as described above. Total RNA was extracted with RNAsimple Total RNA Kit (TIANGEN Biotech, Beijing), cDNA synthesis from total RNA was conducted with the one-step gDNA removal and cDNA synthesis super mix kit (TransGen Biotech, Beijing). cDNA was diluted to 10-fold, and 1 μL (100 ng) of cDNA which was used for qRT-PCR with HiScript II Q RT SuperMix (Vazyme Biotech, Nanjing). qRT-PCR was conducted on LightCycler 480 II (Roche Life Sciences, Swiss). The actin (Shact) gene of C. jacksonii was selected as the housekeeping gene, primers for quantifying Shact and ShPDR1 gene were used as reported by Hulvey et al. (2012), and primers for ShCYP68 gene quantification were used according to Sang et al. (2018). Other primers were designed with Primer5 (https://www.researchgate.net/journal/Biotech-Software-Internet-Report-1527-9162; Supplementary Table S1). The comparative threshold cycle (C T ) method was used for the calculation of relative gene expression. Linear regression analysis of the selected DEGs was performed with FPKM values from RNA-Seq and relative expression values from RT-qPCR.

Genetic background and fungicide sensitivities of the isolates for RNA-seq
Phylogenetic analysis based on ITS sequences revealed that the two isolates LT15 and LT586 belonged to the same species: C. jacksonii (Supplementary Figure S1). In vitro sensitivity assays confirmed that the isolate LT586 showed significantly reduced sensitivities to iprodione, propiconazole, and boscalid ( Figure 1A) compared with the isolate LT15 (p < 0.001) ( Figure 1B). No point mutations were found on the fungicide target genes (ShCYP51B, Shβ-tubulin, Shos1, and ShSdh genes) in LT586, suggesting that LT586 is a MDR isolate of C. jacksonii (Figure 1).

Transcriptomic assembly data
In this study, comparative transcriptomic analyses were conducted to identify the DEGs between the two isolates (LT15 and LT586) without or with the treatments of propiconazole, boscalid, and iprodione (Supplementary Table S2). A total of 39,386,718-59,568,566 clean reads were obtained from 24 libraries. GC content was ranging from 48.75 to 49.23% and the percentage of Q30 base was 93.27% or above (Supplementary Table S2). By iterative alignment, 93.56-95.79% of the clean reads were successfully mapped to C. jacksonii reference genome (Zhang et al., 2023). Of the mapped reads, 92.49-94.59% of them was mapped to unique loci, while 0.80-1.26% of the reads was mapped to multiple loci in the reference genome (Supplementary Table S2

Differential expression between the two isolates without fungicide treatments
Principal component analysis revealed that different compartments were clustered apparently, PC1 and PC2 explained 24.87 and 16.88% variation between the two isolates without and with fungicide treatment, respectively ( Figure 2A). The DEGs were identified, a total of 619 and 475 DEGs were significantly down-and upregulated in the MDR isolate LT586, compared with the sensitive isolate LT15 without fungicide treatment ( Figure 2B).

Differential expression between the two isolates induced by the treatment of three fungicides
RNA-seq analyses revealed that 398, 351, and 530 genes were significantly upregulated, while 736, 529, and 668 genes were significantly downregulated in MDR isolate LT586, compared with the sensitive isolate LT15 under the treatment of propiconazole, boscalid, and iprodione, respectively ( Figure 3A). A total of 306 and 153 genes showed significantly lower and higher expressions in the MDR isolate LT586 than those in the sensitive isolate LT15, which were commonly induced by the three fungicides ( Figure 3B). KEGG pathway enrichment analysis revealed that the 306 downregulated genes were mainly enriched in butanoate metabolism (n = 6), terpenoid backbone biosynthesis (n = 6), tyrosine metabolism (n = 6), and starch and sucrose metabolism (n = 7; Figure 3C); the 153 upregulated genes were mainly enriched in peroxisome (n = 4), glycerophospholipid metabolism (n = 4), protein processing in endoplasmic reticulum (n = 5), and ABC transporters (n = 4; Figure 3D).
Heatmap analysis showed that the 153 DEGs were clustered apparently in LT15 and LT586, and were significantly higher expressed in the MDR isolate LT586 under the treatments by the three fungicides ( Figure 3E). The function of 153 upregulated DEGs further Mycelial growth of the MDR isolate (LT586) and fungicide-sensitive isolate (LT15) of Clarireedia jacksonii on non-amended PDA and propiconazole, boscalid, and iprodione amended PDA at 0.1, 10, and 1 μg/mL, respectively. *** means significant difference between the two strains under the fungicide treatment with p value < 0.001. Microbiology  05 frontiersin.org characterized, which were mainly classified into seven molecular function terms, five biological processes, and three cellular components ( Figure 3F). The molecular functions of the 153 DEGs were mainly involved in FAD binding, monooxygenase activity, transmembrane transporter activity, heme binding, zinc ion binding, ATP binding, and oxidoreductase activity. The five biological processes were intracellular signal transduction, response to stimulus, cellular protein metabolic process, L-arabinose metabolic process, and carbohydrate metabolic process. The cellular components included extracellular region, cytoplasm, and integral component of membrane ( Figure 3F). The specific information of the genes with Pfam function annotations were shown in Supplementary Table S3. Most of them were involved in the functions of xenobiotic detoxification-related genes, which mainly included oxidoreductase activity, transferase activity, ATPase activity, and transmembrane movement of substance. There were 10 genes involved in oxidoreductase activity (GE10720_g, GE15552_g, GE05437_g, GE03831_g, GE03356_g, GE02832_g, GE09522_g, GE05251_g, GE10755_g, and GE01921_g), which mainly included the functions of alcohol dehydrogenase, zinc-binding dehydrogenase, and NADH: flavin oxidoreductase. Nine genes were involved in ATP binding, including three ABC transporter genes (GE04562_g, GE10243_g, and GE13106_g), one Cdc48 subfamily gene (GE15443_g), one asparagine synthase gene (GE01321_g), one E1-E2 ATPase gene (GE06702_g), two fungal protein kinase genes (GE09508_g, GE16043_g), and one Hsp70 protein encoding gene (GE08240_g). Of the three ABC transporter genes involved in ATP binding (Schoonbeek et al., 2001;Sang et al., 2015), two (GE04562_g = ShPDR1, and GE10243_g = ShatrD) had been shown to be associated with MDR in C. jacksonii . There were seven genes involved in heme binding, including six CYP450s genes (GE00961_g, GE07314_g, GE05983_g, GE10739_g, GE14553_g, and GE06676_g) and one glycosyl hydrolase family 20 gene (GE05305_g). Three of the six CYP450s genes (GE10739_g = CYP561, GE14553_g = CYP65, GE00961_g = CYP68) were also involved in monooxygenase activity and had been shown to be associated with MDR in C. jacksonii . Seven genes were involved in transmembrane transporter activity, with five genes encoding MFS proteins  Table S3), which had been shown to be associated with fungicide resistance by active efflux of xenobiotics (Kretschmer et al., 2009;Omrane et al., 2015). Three genes were involved in FAD binding (GE05437_g, GE13986_g, and GE13409_g). Eight genes were involved in zinc ion binding (GE05988_g, GE06509_g, GE09053_g, GE09312_g, GE09402_g, GE10720_g, GE12966_g, and GE15454_g), which mainly included

Comparison with the published data
The transcriptome data from a previously published paper  were re-analyzed and mapped to the reference genome used in this study (Supplementary Table S2). Three hundred and sixtynine genes were significantly upregulated in the MDR isolate HRI11, as compared with the sensitive isolate HRS10 under no fungicidetreatment (Supplementary Figure S4A). Fifty of these genes also showed significantly higher expression in LT586 than LT15 from this study (Supplementary Figure S4B). Most of the differentially expressed genes between the two studies were both from the same gene families, such as major facilitator superfamily, CYP450s, and fungal specific transcription factors (Supplementary Figure S4C).
In the previously published paper , C. jacksonii isolates were only treated by propiconazole, therefore the transcriptional expression analyses between the two studies were based on the data of propiconazole treatments. A total of 203 (PHRS10 vs. PHRI11) and 398 (PLT15 vs. PLT586) genes were significantly upregulated in the MDR isolates HRI11 and LT586 after treatment of propiconazole ( Figure 4A). Seventeen genes showed commonly upregulated overexpression in the two MDR isolates, which included the following five genes (CYP561, CYP65, CYP68, ShPDR1, and ShatrD), confirmed to be related with MDR function in the MBio paper  Figure 4B). These genes showed a generally constitutive overexpression and fungicide-induced expression trend in the two MDR isolates, except for CYP65, which was only induced expressed under the treatment of boscalid in this study ( Figure 4C).

qRT-PCR analysis of key DEGs involved in fungicide resistance
To validate the expression pattern of some key DEGs involved in fungicide resistance from RNA-Seq analysis, three zinc finger transcript factor genes (GE07854_g, GE09003_g, and GE09664_g), four MFS transporter genes (GE06394_g, GE09014_g, GE09100_g, and GE09173_g), one CYP450 gene (GE00961_g), and two ABC transporter genes (GE13106_g and GE04562_g) were further quantified with qRT-PCR ( Figure 5A). The log 2 -transformed relative expression values determined by qRT-PCR and RNA-seq were positively correlated (R 2 = 0.7776, p < 0.001; Figure 5B). The result of qRT-PCR analysis was consistent with the RNA-seq data.

Discussion
Clarireedia jacksonii is the most widely distributed DS pathogen species across the world (Viji et al., 2004;Hu et al., 2019a,b). Due to widespread, repeated uses of multiple fungicide classes for DS control, MFR and MDR have become an increasing issue in C. jacksonii FIGURE 4 Significantly upregulated genes in HRS10 plus propiconazole vs. HRI11 plus propiconazole and LT15 plus propiconazole vs. LT586 plus propiconazole. (A) Venn plot analysis. (B) The specific information of 17 commonly upregulated genes through Venn plot analysis. (C) Gene expressions of CYP561, CYP65, CYP68, ShPDR1, and ShatrD in the MDR isolate LT586 and the sensitive isolate LT15 according to RNA-seq data. *means significant difference between the fungicide treatment group and the PDA treatment group in the LT586 isolate value of p < 0.05.
Frontiers in Microbiology 08 frontiersin.org populations from different regions or countries in recent years (Sang et al., 2016;Stephens and Kaminski, 2019;Zhang et al., 2021). In this study, fungicide sensitivity tests and comparative analyses of fungicide target genes clearly indicated that the C. jacksonii isolate LT586 resistant to multiple fungicides was a MDR isolate. At present, MDR in fungal pathogens has become great concerns both in clinical and agricultural fields (Leroux et al., 2012;Sang et al., 2015). Several mechanisms have been reported to cause MDR in pathogenic fungi (Nakaune et al., 2002;Thakur et al., 2008;Kretschmer et al., 2009;Sang et al., 2018). Therefore, we profiled the response of two C. jacksonii isolates contrasting in sensitivity to different active ingredients of three fungicide classes: iprodione, propiconazole, and boscalid, to identify the genes involved in MDR. Moreover, we compared the results between this study and a previous study  to get a better understanding of MDR mechanisms. This study identified common MDR-related genes, which have also been described in study of Sang et al. (2018), as well as new genes potentially related to MDR. We observed that many genes were constitutively or induced over-expressed by the three fungicides in the MDR isolate LT586 (Figures 2, 3). In general, both constitutive and induced expressions of MDR-related genes can lead to MDR . Therefore, there were still a substantial number of DEG changes between LT15 and LT586 even in the absence of fungicide treatment. Moreover, the genetic differences also might contribute to the constitutive changes of the DEGs between the two isolates. In this study, both constitutive and induced DEGs were identified between the two isolates, which provided us a better reference to locate the potential MDR genes for further analyses. Among these DEGs identified, some have been previously characterized to be associated with MDR in C. jacksonii , some homologous genes were identified to be responsible for fungicide resistance in other fungi, which have not been studied in C. jacksonii. It is proposed that the MDR mechanisms in C. jacksonii are complicated and need to be further investigated.
When compared with the data from a similar study conducted by Sang et al. (2018), we found same and novel DEGs between the two studies, which suggested that the MDR mechanisms in C. jacksonii might have common as well as contrasting regulatory pathways, depending on the origins of the isolates. There are two main reasons that may contribute to the differences between the two studies. Firstly, the C. jacksonii isolates used in the two studies were different in genetic background and fungicide sensitivities. The key genes regulating MDR in these isolates might be selected differently by different ways of fungicide uses or different environment conditions (Koch et al., 2009;Latin, 2011). Secondly, the isolates in the current study were treated by three fungicides, while only propiconazole was used to treat C. jacksonii isolates in the study by Sang et al. (2018).
Modifications, over-expression of target genes, and over-expression of efflux transporters in plasma membranes are well known mechanisms of fungicide resistance in filamentous fungi . Transcriptome analysis in this study revealed that the target gene of DMI, propiconazole (CYP51) was significantly upregulated in LT586, while the target genes of SDHI, boscalid (SdhB, SdhC, and SdhD) and dicarboximide, iprodione (Shos1) did not show differential expressions between the two isolates (Supplementary Figure S5). Combined with no mutations found on these target genes, indicated that the observed MDR in LT586 was probably caused by the over-expressions of CYP51 and other genes involved in xenobiotics detoxification systems. A significantly higher level of resistance to propiconazole observed in LT586 further supported partial contribution of CYP51 to MDR. Further investigations should be oriented to characterize if there exist commonly regulatory factors responsible for over-expressions of CYP51 as well as the genes in xenobiotics detoxification systems simultaneously.
General xenobiotics detoxification systems work through pathways that can be divided into three phases. Phase I includes metabolizing enzymes such as CYPs, that catalyze xenobiotics, phase II involves conjugating enzymes (e.g., glutathione S-transferase) by adding polar molecules onto compounds (Xu et al., 2005), and in phase III secretion system, ABC transporters or other transmembrane transporters export parent and/or metabolized compounds (Ihunnah et al., 2011). Previous study revealed that the constitutive and induced over-expressions of three phase I metabolizing enzyme genes CYP561, CYP65, and CYP68, Frontiers in Microbiology 09 frontiersin.org two phase III transporter genes ShPDR1 and ShatrD were associated with MDR in C. jacksonii . However, in this study, CYP65 in the MDR isolate LT586 was only upregulated in response to the treatment of boscalid, but not propiconazole or iprodione.
Overexpression of different CYPs displayed various sensitivities to different chemical groups of fungicides, which might be due to the different substrate specificity of CYPs (Moktali et al., 2012). Beside the three above CYPs, three other CYPs were significantly upregulated in the MDR isolate LT586. Further studies are needed to characterize the full functional roles of these differentially expressed CYPs found in this study, which may gain a better understanding of the associations between these different CYPs and MDR. In this study, one phase II glutathione S-transferase (GE09048_g) showed constitutive and induced overexpression in the MDR isolate LT586, this gene was only constitutively overexpressed in the MDR isolate HRS11 in the previous study . The different expression patterns in different MDR isolates are expected and suggest the role of the gene in fungicide resistance needs to be further investigated.
Notably, the family members including ABC transporters, MFS transporters, zinc-cluster transcription factors, and CYP450s were significantly enriched in these potential multidrug resistant DEGs (Supplementary Table S3). Regarding MFS transporters, they are generally less understood compared with ABC transporters. MFS transporter genes such as GE09014_g, GE09100_g, and GE09173_g in this study (Supplementary Table S3) were highly homologous with YLL028W (Tpo1), YLL028W (Tpo1), and YNL065W (Aqr1) in yeast. Tpo1 was demonstrated to confer resistance to benzoic acid in Saccharomyces cerevisiae, and under the controls of the Gcn4 and Stp1 transcription factors involved in the response to amino acid availability (Godinho et al., 2017). However, the sequences of these three genes are consistent in MDR isolate LT586 and fungicide sensitive isolate LT15 without point mutations. Therefore, the roles of these MFS genes in fungicide resistance and the interactive mechanisms between MFS genes and other transcription factors need to be further studied.
In addition, it has been not well characterized how fungi transcriptionally regulate and coordinate MDR related genes. Transcription factors have been confirmed to play important roles in drug resistance in multiple pathogenic fungi and they generally interacted with downstream target genes or other transcription factors (Nakaune and Nakano, 2007;Thakur et al., 2008;Wang et al., 2018). To date, only the transcription factor ShXDR1 has been demonstrated to regulate downstream MDR related genes in C. jacksonii . Our results indicated that essential and specific roles for other fungal transcription factors in regulating MDR of Clarireedia spp. should be investigated (Supplementary Table S3). Further genetic transformation and protein interactive analyses will help better understand if other transcription factors are involved in MDR, and how they regulate MDR in Clarireedia spp.

Conclusion
In this study, we utilized comparative transcriptomics to glean insight into the potential MDR-related genes associated with the reduced sensitivities to multiple fungicides (propiconazole, boscalid, and iprodione) in C. jacksonii isolates, comparative analyses were also conducted with the datasets from the current and a previous study . Common as well as novel DEGs were found to be associated with MDR in this study. The molecular functions of these DEGs were mainly involved in xenobiotic detoxification-related systems and transcriptional regulations. The results obtained from this study will provide valuable information for further identifying the molecular mechanisms of MDR. The study also strongly indicated that the MDR mechanisms in C. jacksonii are complicated, which need to be thoroughly characterized through functional studies in the future.

Data availability statement
The datasets presented in this study can be found in online repositories. The raw RNA sequences were deposited in NCBI SRA Database under BioProject accession number PRJNA1007393.

Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. We thank National Natural Science Foundation of China (grant number: 32272597), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX23_0830), and Jiangsu Institute of Botany Talent Fund (JIBTF202212) for financial supports of this research.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.