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 (Zhang et al., 2021). 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).

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 (Sang et al., 2018; Zhang et al., 2021). The untreated samples were added with same volume of sterile water (Zhang et al., 2021). 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).

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¹ 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.

Gene functions were annotated through the following databases: NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Protein family (Pfam), Clusters of Orthologous Groups of proteins (KOG/COG), Swiss-Prot, KEGG Ortholog database (KO), and Gene Ontology (GO).

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₂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.² 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., 2021). We used KOBAS software to test the statistical enrichment of the DEGs in KEGG pathways (Bu et al., 2021).

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.

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).

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). The raw RNA sequences were deposited in NCBI SRA Database under BioProject accession number PRJNA1007393.

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).

Clusters of Orthologous Groups of proteins function analysis revealed that all these DEGs were mainly involved in general function prediction (23.29–99%), secondary metabolites biosynthesis, transport, catabolism (12–51%), and energy production and conversion (9.88–42%; Supplementary Figure S2). KEGG pathway of all the DEGs were classified into four groups, which included cellular process, environmental information processing, genetic information processing, and metabolism (Supplementary Figure S3). KEGG pathway analysis revealed that these DEGs were mainly involved in protein processing in endoplasmic reticulum (n = 22), terpenoid backbone biosynthesis (n = 10), tyrosine metabolism (n = 13), diterpenoid biosynthesis (n = 3), and steroid biosynthesis (n = 10; Figure 2C). GO function analysis showed that the molecular functions of all these DEGs were mainly involved in dioxygenase activity (n = 10), monooxygenase activity (n = 19), iron ion binding (n = 26), oxidoreductase activity (n = 24), and heme binding (n = 29); the biological processes were mainly involved in steroid biosynthetic process (n = 4), DNA integration (n = 9), L-arabinose metabolic process (n = 4), tryptophan catabolic process to kynurenine (n = 5), and mycotoxin biosynthetic process (n = 11); the cellular components were mainly involved in integral component of membrane (n = 257; Figure 2D).

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 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 (Sang et al., 2018). 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 (Sang et al., 2018). Seven genes were involved in transmembrane transporter activity, with five genes encoding MFS proteins (GE09014_g, GE09100_g, GE09161_g, GE00853_g, and GE09173_g), one encoding sugar and other transporter protein (GE02436_g), and one encoding POT family protein (GE10772_g). Three of the MFS genes (GE09014_g, GE09100_g, and GE09173_g) were highly homologous with YLL028W (Tpo1), YLL028W (Tpo1), and YNL065W (Aqr1) in yeast (Supplementary 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 the functions of zinc finger domain and fungal specific transcription factor domain.

The transcriptome data from a previously published paper (Sang et al., 2018) were re-analyzed and mapped to the reference genome used in this study (Supplementary Table S2). Three hundred and sixty-nine genes were significantly upregulated in the MDR isolate HRI11, as compared with the sensitive isolate HRS10 under no fungicide-treatment (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 (Sang et al., 2018), 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 (Sang et al., 2018; 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).

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₂-transformed relative expression values determined by qRT-PCR and RNA-seq were positively correlated (R² = 0.7776, p < 0.001; Figure 5B). The result of qRT-PCR analysis was consistent with the RNA-seq data.