Transcriptome and Resistance-Related Genes Analysis of Botrytis cinerea B05.10 Strain to Different Selective Pressures of Cyprodinil and Fenhexamid

The pathogen Botrytis cinerea is a very dangerous pathogen that infects many economically important crops such as grape, strawberry, tomato, and eggplant. Cyprodinil, a pyrimidine amine fungicide, and fenhexamid, an amide fungicide, are new reagents for controlling gray mold with special efficacy. It is necessary to understand the change trends in the toxicological and physiological characteristics of B. cinerea with successive selective pressures of cyprodinil and fenhexamid to elongate the serving life of these fungicides for effective disease control. The toxicities of cyprodinil and fenhexamid at successive concentrations of EC25, EC50 and EC75 on B. cinerea strain BO5.10 were assayed along with mycelial growth-inhibition capacity. The results showed that the EC50 value of the cyprodinil-treated F27 strain increased approximately 18-fold, whereas of which in the fenhexamid-treated F27 strain increased only 3-fold compared with that of the F0 strain. The conductivities and glycerinum contents of the strains resistant to cyprodinil and fenhexamid were obviously enhanced; in contrast, the oxalic acid contents were decreased compared with those in the F0 strain. The transcriptomes of the F27 control (T01), cyprodinil-treated (T02) and fenhexamid- treated (T03) strains were analyzed, and the expression levels of functional genes in the T02 and T03 strains were significantly increased compared with those in the T01 strain; these results were further validated using qRT-PCR. The results indicated that the relative expression of two genes encoding mixed-functional oxidases (MFOs) BC1G_16062 and BC1G_16084, two genes encoding transmembrane proteins BC1G_12366 and BC1G_13768, two genes encoding Zinc finger proteins BC1G_13764 and BC1G_10483,one gene encoding citrate synthase enzyme BC1G_09151, one gene encoding gluconolactonase BC1G_15612 in the T02 and T03 strains and one gene encoding lysophospholipids enzyme BC1G_04893 in the T3 strain increased substantially compared with that in the T1 strain (P < 0.01). Functional prediction analysis of upregulated gene expression and structural verification was also performed, and the results showed that BC1G_10483 was a ZnF_C2HC transcriptional regulator interacting with the Sp1 element of these genes to respond to the pressures from cyprodinil and fenhexamid. Our results could contribute to a better understanding of the resistance mechanism of B. cinerea against cyprodinil and fenhexamid.


INTRODUCTION
Botrytis cinerea is an aggressive plant disease and more than 200 plant species, such as tomato, grape, strawberry and so on, are infected resulting in large output loss (Shao et al., 2015). At present, the control of B. cinerea still mainly relies on chemical agents, with some other auxiliary methods, such as breeding disease-resistant varieties, rational fertilization and improving cultivation facilities. B. cinerea causes a disease with a high resistance risk (Jiang et al., 2009), and extensive use of the same types of fungicides for a long time has contributed to serious resistance, such as benzimidazole (carbendazol, benomyl) (Leroux et al., 2002), thiocarbamates (thiophanate) (Zhang et al., 2009), and dicarboxyl imide (Sansone et al., 2005;Sun et al., 2010) (iprodione, procymidone).
Due to the interactive resistance of antifungal mutant strains to fungicides, developing new fungicides has become increasingly difficult. The effective control methods for graymold developed in recent years include pyrimidine amine (cyprodinil), pyrrole (fludioxonil), and amides (fenhexamid). It is necessary to assay the biological activity and to evaluate the resistance risk before widely using new fungicides in order to delay the development of B. cinerea resistance against those newer fungicides to elongate their service time. A few cases have demonstrated that there was some cross-resistance risk even though their structural classes are significant difference. For example, the resistant-dicarboximides strains of B. cinerea also produced resistance to some fungicides, including dichloran, quintozene, biphenyl and so on, which maybe act on a same possible point, histidine kinase (Leroux et al., 2002).
Cyprodinil, a pyrimidine amine fungicide acting as an inhibitor of the biosynthesis of methionine, was developed by Novartis (Basel, Switzerland). The water-dispersible granule has been known since 2000 to provide protection, treatment, internal suction and transmission in the root and leaf (Bernardo and René, 2012). Cyprodinil has low toxicity, inhibits the secretion of the extracellular protease hydrolysis enzyme of B. cinerea and pathogen penetration, interferes with the fungal life cycle and destroys mycelium growth and development in plants (Mofeed and Mosleh, 2013). There is no cross-resistance among the fungicides benzimidazoles, thiocarbamate, carbamate, imidazoles and so on because of their different structures and target points.
Fenhexamid, an amide fungicide acting on sterol synthesis, was developed by Bayer Crop Science Co., Ltd., (Leverkusen, Germany) and displayed good efficacy against grape gray mold (Zocco et al., 2008). Fenhexamid can inhibit bacterial growth, such as by decreasing the growth of mycelium and inhibiting the extension of the bud tube and the conidia germination (Debieu et al., 2001). Fenhexamid has characteristics of biodegradability, low toxicity, environmental friendliness and no cross-resistance with the production of pyrimidine, carbamate, pyrrole, pyridine or other fungicides widely applied in agriculture.
Currently, there are only a few reports on the field efficacy of cyprodinil and fenhexamid, and the physiological and biochemical effects on resistant strains have rarely been studied. In this study, we assay the toxicities, physiological and biochemical characteristics and transcription of the F 27 strains of B. cinerea B05.10 continuously treated with cyprodinil and fenhexamid at concentrations of EC 25 , EC 50 and EC 75 . The results are expected to promote the comprehension of the mechanisms of resistance development of B. cinerea against cyprodinil and fenhexamid and provide a foundation for the design of more effective resistance management strategies.

Fungicides and Strain
The technical fungicides cyprodinil and fenhexamid (a.i.98%) were purchased from Shan Dong Yi Jia Agricultural Chemical Co., Ltd., and the BO5.10 strain was kindly provided by Prof. Li Guoqing of the Plant Pathology Laboratory, Huazhong Agricultural University in January 2013.

Toxicity Changes in the BO5.10 Strain to Different Pressures of Cyprodinil and Fenhexamid
The toxicities of cyprodinil and fenhexamid on the BO5.10 strain (F 0 ) were assayed with mycelial growth inhibition as described by Kuang et al. (2011). Five concentrations were set for each fungicide. One milliliter of prepared solution was added to a sterile petri dish, followed by the addition of 9 ml of PDA and rapid shaking. After complete cooling, one 5-mm hyphal disk cut from a 3-day-old PDA culture was placed in the center of each plate. Each subsequent concentration increased twofold, and 1 ml of distilled water was used as a control. Then, incubated those colonies about 72 h at 23 • C in a biochemical incubator in the dark and measured their diameters to assess the inhibition of mycelial growth. The toxicity regression equation of the inhibition rate against the log 10 concentrations of cyprodinil and fenhexamid was calculated though SAS PROC REG (version 9.1, SAS institute, Cary, NC). Inhibition rate (%) = (colony diameter of control treatment -colony diameter of fungicide treatment)/(colony diameter of control treatment -5-mm hyphal disk diameter) (Zhou et al., 2016).
According to the toxicity regression of cyprodinil and fenhexamid on the BO5.10 strain (F 0 ), the concentrations of EC 25 , EC 50 , and EC 75 were calculated, and the BO5.10 strain (F 0 ) was successively cultured in three replicates under the three concentrations thus generating the F 3 strains. The effects of the EC 25 , EC 50 , and EC 75 concentrations of cyprodinil and fenhexamid on the F 3 strains were assayed with the same methodology, and additional screenings were continually preformed procedure until we obtained the F 27 strains. Those strains with clearly enhanced toxicities were stored in a mixture of glycerinum-saline (1:3, v/v) at −20 • C to investigate the resistance-related physiological and biochemical characteristics.

Membrane Permeability and Glycerinum and Oxalic Acid Contents of the Strain With Significantly Enhanced Toxicity
Those strains with clearly enhanced toxicities were incubated in 100 ml of PD nutrient solution in 250 ml laboratory flasks (approximately five to seven 5-mm hyphal disks per flask), and then cultivation was performed in a constant-temperature shaker (23 • C, 120 rpm). The BO5.10 strain was used as the blank control. After the tested strains had been cultivated for approximately 7 days, the hyphae were repeatedly flushed with deionized water and vacuum filtered for approximately 15 min. A 0.5-g sample for each treatment was added to 20 ml of double distilled water; the conductivities were detected at 0,5,10,20,40,60,80,100,120,140,160, and 180 min after treatment; and the final conductivities were measured via boiling the mycelium in water for 5 min. Three repetitions were performed for each treatment, and the relative conductivity was calculated as following: the conductivity of different times/the final conductivity (Li et al., 2015).
Glycerine copper colorimetry (Yan and Qiu, 2004;Duan et al., 2013) was used to assay the glycerinum contents of the hyphae under different treatments. Three kinds of solutions, including 1 ml of CuSO 4 (0.05 g/l), 3.5 ml of NaOH (0.05 g/ml) and 10-ml of glycerinum in double distilled water (0, 0.0025, 0.003, 0.004, 0.005, 0.006, 0.008, 0.01, 0.015, and 0.02 g/ml) were transferred into different 50-ml flasks, respectively. Then quickly shook the mixtures, transferred into a 50-ml centrifuge tube, and centrifuged at 100 rpm for 12 min. The absorbance of the supernatant was recorded at 630 nm to establishing the standard curve. A 0.5-g sample from each treatment was ground with 20 ml of sterile water and some quartz sand. The supernatant was transferred into a 50-ml centrifuge tube, heated in an 80 • C water bath for 15 min, and centrifuged at 12,000 × g for 10 min. The glycerinum contents in the supernatant were detected using distilled water as a control. All treatments were performed in triplicate.
The standard curve of oxalic acid was established as follows: 2 ml of FeCl 3 solution (0.5 g/ml), 20 ml of HCl-KCl buffer solution (KCl 50 mM, pH = 2) and 1.2 ml of sulphosalicylic acid solution (0.5 mg/ml) were added to a 25-ml volumetric flask in order; then, aliquots of 0, 0.1, 0.2, 0.4, and 0.8 ml of sodium oxalate solution (2 mg/ml) were transferred into different 25-ml volumetric flasks, brought to volume with doubledistilled water and shaken. The oxalic acid contents were determined based on changes in absorbance at 510 nm in a UV 2000-Spectrophotometer [Unic (Shang Hai) Instruments Incorporated] with double-distilled water as a control (Duan et al., 2014). The pretreatment of each treatment was the same as the approach when assaying for the glycerinum contents. The supernatant was transferred into a 50 ml centrifuge tube and centrifuged at 380 × g for 10 min. The absorbance values of the supernatant were recorded at 510 nm, and the oxalic acid contents were calculated with the standard curve.
Transcriptome Analysis of the BO5.10 and Fungicide-Treated Strains

Library Construction and Sequencing
Total RNAs were extracted from F 27 strains screened from BO5.10 (T01) with EC 50 concentrations of fenhexamid (T02) and cyprodinil (T03) and BO5.10 (F 0 ) using Trizol R Reagent (Invitrogen TM ) according to the manufacturer's instructions and all following of library construction and sequencing were performed as the description of Wang et al. (2018). After the quality of total RNA was checked out, eukaryotic mRNA was enriched by Oligo(dT) beads, then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcripted into cDNA with random primers. The second-strand cDNA were synthesized by DNA polymerase I, RNase H, dNTP and buffer. Then the cDNA fragments were purified with QiaQuik PCR extraction kit, end repaired, polys(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using IlluminaHiSeq TM 4000.

Sequence Alignment Between Transcriptome Data and Reference Genome Sequence
The high quality reads obtained in this study have been deposited in the NCBI SRA database (accession number: SUB4543717). Owing to reads obtained from the sequencing machines including dirty reads, reads containing more than 10% of unknown nucleotides, more than 50% of low quality (Q-value ≤ 10) bases and adapters were further filtered in order to get high quality clean reads. Clean reads of three samples were aligned with the reference genome 1 through TopHat2, and the reference genome positions and the unique sequences of the tested samples were identified (Kim et al., 2013). Clean reads that were able to align to the reference genome were referred to as mapped reads.

Statistics and Annotation of Single Nucleotide Polymorphisms
Based on the alignment results between the reads of the three treatments and the reference genome used by TopHat2, single erroneously paired nucleotides (SNPs) as determined by comparison with the reference genome were distinguished by The Genome Analysis Toolkit (GATK) (McKenna et al., 2010), and then, these SNP sites were analyzed by the program SnpEff (Cingolani et al., 2012). Furthermore, it was possible to analyze whether these SNP sites affect the gene expression level or encoded protein type.

Discover New Gene and Expressed Gene
Based on the selected reference genome sequence, mapped reads were spliced though Cufflinks software (Trapnell et al., 2010) and compared with the original genome annotation to search for original transcription areas that were not annotated, and the species of transcription and new genes were explored to complement and complete the existing genome annotation information. The newly discovered genes were aligned using BLAST software (Altschul et al., 1997) in NR, Swiss-Prot, GO, COG, KOG, Pfam, ad KEGG databases, and the KEGG orthology (Ashburner et al., 2000;Tatusov et al., 2000;Apweiler et al., 2004;Koonin et al., 2004;Deng et al., 2006;Finn et al., 2013;Kanehisa et al., 2004) results of new genes were generated using KOBAS2.0 (Xie et al., 2011). Information on new gene annotation was acquired by comparison between HMMER software and the Pfam database after the amino acid sequences of new genes had been predicted.
Transcriptome sequencing can be modeled as a random sampling process, that is, from an arbitrary sample of the transcriptome is generated a nucleic acid sequence of an independent random sequence fragment. The related to the mapped data, the length of the transcribed sequence and the transcription expression levels were used to determine a fragment number that truly reflected the transcript expression level, the number of mapped reads and the normalized length of transcription. Cuffquant and Cuffnorm used Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) as a measure of transcription or gene expression levels (Florea et al., 2013;Lei et al., 2014) as follows: cDNA fragements mapped fragements (Millions) × transcript length(kb) In the formula, cDNA fragments stand for the number of segments that were able to align on a transcript, that is, the number of two read pairs; mapped fragments (millions) represent the total number of segments that were able to align during transcription, in a denomination of 10 ∧ 6 bp; and transcript length (kb) represents the length of transcription, in a denomination of 10 ∧ 3 bp.

Differential Expression Analysis
The expressed genes of three samples (T01, T02, and T03) were analyzed using EBSeq. Because differential expression analysis of the transcriptome sequencing was performed independently on a large, statistically hypothetical set of gene expression values, some results could be false positives. Therefore, the significant P-values extracted from the original hypothesis testing were corrected using the recognized Benjamini-Hochberg correction method, and eventually, the fold change (the ratio of expression amount between two samples) and FDR (based on correcting the significant P-values) were used as the key indicators of screening differentially expressed genes (DEGs), with Fold Change ≥ 2 and FDR < 0.01 as screening criteria. According to the requirements of each group, a Venn diagram of DEGs was drawn and analyzed by hierarchical clustering and grouped genes with the same or similar expression patterns. The DEGs were aligned using BLAST software in the COG and KEGG databases, and annotated with KOBAS2.0 and aligned using BLAST software in the STRING database (a database of multiple species hypothetic and experimental protein-protein interactions), and then, homologous proteins were searched for. The protein-protein interaction networks were built based on the interactions of homologous proteins (Franceschini et al., 2013), and a protein-protein interaction network scheme was drawn using Cytoscape software for visualization (Shannon et al., 2003).

Quantitative PCR (qRT-PCR)
Total RNA from the F 27 strains screened from BO5.10 treated with the EC 50 concentrations of fenhexamid (T02) and cyprodinil (T03) and from BO5.10 (T01) was extracted as described in section of Library Construction and Sequencing. For the reverse transcription reaction, the synthesis of First-Strand cDNA was used according to the kit. Then, the cDNA was stored at −20 • C for qRT-PCR.
Sixteen pairs of respective primers, including the cDNAs of fifteen genes (BC1G_12765, BC1G_12768, BC1G_16062, BC1G_16084, BC1G_04656, BC1G_04779, BC1G_12366, BC1G_13768, BC1G_04893, BC1G_09151, BC1G_15612, BC1G_02313, BC1G_04879, BC1G_13764, BC1G_10483) and β_tubulin using for a house-keeping gene (Yang et al., 2012) in the T1, T2 and T3 strains were displayed in Supplementary  Table 1 for amplifying by qRT-PCR, which was a total 20 µl reaction system, including 1 µl of cDNA, 0.4 µl of each primer, 10 µl of 2 × TransStart R TopTaq Green qPCR SuperMix (Transgene, Beijing, China) and 8.2 µl of double-distilled water. The standard cycle was performed in 94 • C for 30 s, 40 cycles of 94 • C for 30 s and 60 • C for 30 s. The initial temperature of the dissolution curve was 68 • C. The β_tubulin gene of B. cinerea was used as the internal control. The relative normalized expression of target genes was calculated using the 2 − Ct method (Livak and Schmittgen, 2001). Three independent replicates were held in all strains and target genes, and all data were analyzed via one-way ANOVA by PROC GLM of the SAS program. All means were compared by least squared difference (LSD) tests at Type I error = 0.05.

Function Prediction of Upregulated Genes and Structural Verification
The structural domains of the upregulated genes that were annotated using GO annotation were predicted to be P450 2 , their structural domains were verified using motif_scan and ProDom 3 , and their topological domains were analyzed 4 .
The cis-acting elements of the DEGs were predicted using AliBaba2.1 5 , and their trans-acting factors, which could combine Transcriptome Analysis of Botrytis cinereaB05.10 Strain with the cis-acting elements through JASPAR CORE Fungi 6 , were analyzed. Prediction analysis of transcription factors of DEGs was performed through a BLAST Search of FTFD 7 , and their structural domains were verified using SMART 8 and the PFAM SEARCH 9 .

Toxicity Changes in the BO5.10 Strain in Response to Different Concentrations of Cyprodinil and Fenhexamid
The original EC 50 value of cyprodinil on the BO5.10 strain (F 0 ) was 0.03 µg/ml, and those in the F 6 strains generated through successively screening with the EC 25 , EC 50 and EC 75 concentrations increased to 0.13 µg/ml, 0.17 µg/ml and 0.16 µg/ml, respectively, and the toxicities of these strains were increased almost 5-fold. The toxicities of strains F 9 to F 15 screened with the EC 25 , EC 50 and EC 75 concentrations slightly increased, whereas the toxicities of F 18 strains rapidly increased under the pressures of EC 50 (0.50 µg/ml); the toxicities plateaued after 20 screenings. Finally, the EC 50 values of cyprodinil on the F 27 strain screened with the EC 25 , EC 50 and EC 75 concentrations reached 0.53 µg/ml, 0.55 µg/ml, and 0.51 µg/ml, respectively, and increased almost 18-fold compared with that of the BO5.10 strain (F 0 ) ( Figure 1A).
The toxicity of fenhexamid on the BO5.10 strain (F 0 ) continuously increased as the EC 25 concentration pressure increased from strains F 0 to F 27 , which was enlarged approximately 3-fold compared with that of the F 0 strain (0.18 µg/ml). The toxicity of the strains screened with the EC 50 concentration from strains F 0 to F 15 rapidly increased, followed by a slim change from strains F 15 to F 18 , with EC 50 values from 0.37 to 0.49 µg/ml. The toxicity of the strains screened with EC 75 from strains F 0 to F 9 significantly increased. The trend was slow from strains F 9 to F 15 , whereas the toxicities of F 15 to F 27 strains continuously increased with EC 50 values from 0.39 µg/ml to 0.50 µg/ml ( Figure 1B).

Membrane Permeability and Glycerinum and Oxalic Acid Contents of BO5.10 Under the Successive Pressure of Cyprodinil and Fenhexamid
The results indicated that the relative conductivities of BO5.10 from strains F 6 to F 27 without fungicidal pressure were 0.2816∼0.3014 (P > 0.05), whereas which in the strains screened with the EC 25 , EC 50 , and EC 75 concentrations of cyprodinil from strains F 6 to F 27 increased, reaching 0.3005∼0.4962, 0.3477∼0.5106 and 0.3712∼0.5165, respectively, and were significantly higher than that in the blank control (P < 0.01) ( Figure 2A1). The relative conductivities of strains F 6 -F 27 screened with different concentrations of fenhexamid were clearly enhanced, especially for strains F 6 -F 19 and F 24 -F 27 screened with the EC 50 and EC 75 concentrations, achieving 0.3174∼0.4489, 0.3326∼0.4691 and 0.4536∼0.4594, 0.4698∼0.4805, respectively, and were significantly different from those in the EC 25 treatment (0.2906∼0.3903 and 0.4213∼0.4354) (P < 0.05). The relative conductivities of strains F 6 -F 27 screened with the EC 25 , EC 50 , and EC 75 concentrations of fenhexamid were significantly stronger than those in the control (P < 0.01) ( Figure 2B1).
The glycerinum contents in those strains screened with fungicides were positively correlated with the screen dosage, especially for strains F 19 -F 27 screened with the EC 75 of cyprodinil, which were remarkably enhanced (3.87∼7.99 mg/ml) and significantly different from the other treatments (P < 0.05) ( Figure 2A2). The glycerinum contents of strains F 15 -F 27 screened with the EC 75 and EC 50 of fenhexamid were  T01, T02, and T03 stand for the BO5.10 strain with no fungicide treatment, F 27 strains successively screened with EC 50 concentration of cyprodinil or fenhexamid, respectively; BMK-ID, the sample number for analysis; Clean reads, the total of pair-end Reads in Clean Data; Clean bases, the total base Clean Data; GC content, the percent of G and C in total base for the Clean Data; ≥Q30%, the percent of quantity ≥30 in Clean Data. BMK-ID, the sample number for analysis; Total Reads, the number of Clean Reads; Mapped Reads, the number of Clean Reads that were able to align on the reference genome data, and the ratio of the mapped reads compared with the total Reads; Uniq Mapped Reads, the number of Clean Reads that were able to align on the only location of the reference genome, and the ratio of the uniq mapped reads compared with the total Reads; Multiple Map Reads, the number of Clean Reads that were able to align on the multiple location of the reference genome, and the ratio of the multiple mapped reads compared with the total Reads; Reads Map to "+", the number of Clean Reads that were able to align on the positive-sense strands of the reference genome data, and the ratio of these reads compared with the total Reads; Reads Map to "−", the number of Clean Reads that were able to align on the antisense strands of the reference genome data, and the ratio of these reads compared with the total Reads. The number of Mapped Reads in different regions of the reference genome (exon, intron, and intergenic region) was counted, and the distribution of Mapped Reads was drawn. In theory, Reads from mature mRNA should be compared to exon. Reads compared to the introns were the precursor of mRNA. Reads in the intergenic regions were genes not well commented on the genome.
conspicuously elevated (3.83∼7.64 and 3.87∼8.09 mg/ml) and were extremely higher than in the other treatments (P < 0.01), followed by strains F 19 -F 27 screened with the EC 25 of fenhexamid, reaching 3.75∼7.20 mg/ml, and were significantly different from those in the control (3.64∼3.81 mg/ml) ( Figure 2B2). The oxalic acid contents in the strains screened with fungicides were negatively correlated with the screen dosages, and those in the blank control were the highest, especially for strains F 6 -F 27 , with 90.67∼98.33 µg/ml. In contrast, the oxalic acid content of the F 27 strain screened with the EC 75 treatment of cyprodinil was the lowest (28.67 µg/ml) (P < 0.01), followed by those in the EC 50 and EC 25 treatments, reaching 84.0∼36.67 and 89.00∼59.33 µg/ml, respectively ( Figure 2A3). The oxalic acid contents in strains F 6 -F 27 screened with the EC 75 , EC 50 , and EC 25 concentrations of fenhexamid reached 89.00∼55.00 µg/ml, 84.00∼48.33 µg/ml and 78.33∼39.67 µg/ml, respectively. Furthermore, those in the treatments with EC 75 and EC 50 were clearly lower than those in the treatment with the EC 25 (P < 0.05) (Figure 2B3).

Transcriptional Data Statistics and Sequence Alignment Between Transcriptome Data and the Reference Genome Sequence
The results indicated that these samples had nearly 13.79 Gb of clean data, and the Q30 exceeded 95.72% (Table 1). According to the results, the contrast efficiencies of reads reached 80.62∼82.96% compared to the reference genome ( Table 2). Transcriptome sequencing of three samples was analyzed by Illumina HiSeqTM 4000, and more than 23,000,000 mapped reads were obtained after screening and filtration in each sample contracted by the reference genome. The ratios of mapped reads, uniquely mapped reads, and multiple mapped reads compared with the total reads reached over 80, 70, and 10%, respectively. The ratios of reads mapped to "+" and reads mapped to "−" were all over 35% (Figure 3A).

Statistics and Annotation of Single Nucleotide Polymorphisms
The SNP site can be divided into transition and transversion according to the different methods of base replacement. Based on the selected reference genome sequence, 3610 SNP sites in the T01 strain, 4057 SNP sites in the T02 strain, and 3797 SNP sites in the T03 strain were demonstrated (Table 3), of which 2439 SNP sites existed in all three strains, 680 SNP sites were unique to the T01 strain, 899 SNP sites were unique to theT02 strain, and 697 SNP sites were unique to the T03 strain ( Figure 3B).

Discover New Genes
Short peptide chains of less than fifty amino acid residues or containing only an exon were filtered out, and 123 new genes were discovered. Based on the messages from the COG database (Supplementary Figure 1A), twelve new genes were found to participate in carbohydrate transport and metabolism (1), post-translational modification (2), replication (3), energy production and conversion (4), amino acid transport and metabolism (5), translation (6), signal transduction mechanisms (7), general function prediction only(8), and lipid transport and metabolism (9). Seventy-four homologous genes of B. cinerea, five homologous genes of Sclerotinia sclerotiorum, one homologous gene of Sclerotinia borealis, one homologous   Figure 1B).

Gene Expression and Differential Expression Analyses in Treatments
The results indicated that 4076 of 10414 genes in the T01 strain, 3962 of 10528 genes in the T02 strain, and 4102 of 10388 genes in the T02 strain were not expressed, while the expression levels of genes in the T02 and T03 strains were higher than those in the T01 strain (Figure 4A), which demonstrated that some genes could be significantly enhanced when the strains were induced by the successive fungicide exposure. Gene expression has temporal and spatial specificity, and the expression levels of gene transcription significantly differing under two different conditions denote DEGs. According to the transcription database, there were 1382 DEGs between the T01 and T02 strains (G0), of which 555 were up-regulated and 827 were down-regulated; there were 1228 DEGs between the T01

Clustering and Enrichment Analysis of DEGs
The results shown in an interactive hot diagram indicated that the DEGs in the T02 and T03 strains were less related than those in the T01 and T02 strains or those in the T01 and T03 strains (Figure 5). A total of 196 of 1864 DEGs were annotated in the KEGG database (Figure 6), of which 26 participated in cellular processes; 4 genes participated in MAPK signaling pathway; 23 participated in genetic information processing, 143 participated in metabolism (some genes participated in two or more metabolisms), including 51 participated in amino acid metabolism, 2 in biosynthesis of other secondary metabolites, 59 in carbohydrate metabolism, 10 in energy metabolism, 37 in global and overview maps, 3 in glycan biosynthesis and metabolism, 25 in lipid metabolism, 11 in metabolism of cofactors and vitamins, 6 in metabolism of terpenoids and polyketides, and 3 in nucleotide metabolism. A total of 1041 of 1864 DEGs were annotated in the COG database (Supplementary Figure 2), including 19 translation, ribosomal structure and biogenesis genes; 26 transcription genes; 31 replication, recombination and repair genes; 3 chromatin structure and dynamics genes; 8 cell cycle control, cell division, chromosome partitioning genes; 14 defense mechanisms genes; 30 signal transduction mechanisms genes; 20 cell wall/membrane/envelope biogenesis genes; 1 cell motility gene; 8 cytoskeleton genes; 2 intracellular trafficking, secretion, and vesicular transport genes; 26 post-translational modification, protein turnover, and chaperones genes; 68 energy production and conversion genes; 132 carbohydrate transport and metabolism genes; 135 amino acid transport and metabolism genes; 6 nucleotide transport and metabolism genes; 26 coenzyme transport and metabolism genes; 64 lipid transport and metabolism genes; 92 inorganic ion transport and metabolism genes; 102 secondary metabolites biosynthesis, transport and catabolism genes; and 211 general function prediction genes only.
FIGURE 5 | Heat map of the DEGs. The column represents samples, and the row means genes; Color-scaled log 2 (fold change) values for resistant lines, the redder the color is, the higher the gene expression is, and on the contrary, the greener the color is.

qRT-PCR of DEGs
According to the transcriptome data, the obvious upregulated expression of four MFO genes, including BC1G_12765, BC1G_12768, BC1G_16062, and BC1G_16084; four genes encoding transmembrane proteins, BC1G_04656, BC1G_04779, BC1G_12366, and BC1G_13768; four genes encoding Zinc finger proteins acting as regulation factors, including BC1G_02313 BC1G_04879, BC1G_13764 and BC1G_10483; and three other important functional genes, including one gene encoding lysophospholipids enzyme BC1G_04893, one gene encoding citrate synthase enzyme BC1G_09151 and one gene encoding gluconolactonase BC1G_15612, could be related to the development of resistance of B. cinerea against cyprodinil and fenhexamid ( Table 5).
To validate the transcriptome data, 15 significantly upregulated genes in strains T2 and T3 were analyzed using qRT-PCR (Figure 7). The expression of BC1G_16062 and BC1G_16084 in strainsT2 and T3 substantially increased compared with that in the T1 strain (P < 0.01), with RNE of 10.7-, 4.3-and 14.1-, 6.7-fold, respectively. However, the expression of BC1G_12768 was significantly down-regulated in the T2 and T3 strains compared with that in the T1 strain (P < 0.05), with RNE of 0.72-and 0.30-fold, respectively ( Figure 7A). The RNE of BC1G_12366 and BC1G_13768 in the T2 and T3 strains was 41.6-, 29.4-and 4.6-, 12.0-fold (P < 0.01), respectively ( Figure 7B). The expression of BC1G_13764 and BC1G_10483 in the T2 and T3 strains clearly increased the RNE by 5.1-, 6.3-and 21.0-, 11.2-fold, respectively ( Figure 7C). The expression of other functional genes, including BC1G_09151 and BC1G_15612 in the T2 strain and BC1G_04893, BC1G_09151, and BC1G_15612 in the T3 strain, was also enhanced and significantly different from that in the T1 strain (P < 0.05) (Figure 7D).

Prediction Analysis of the Function of Upregulated Genes and Structural Verification
All typical characteristics of the P450 gene, including the sites of heme binding and iron ion binding, monooxygenase activity source, and oxidoreductase activity, were demonstrated in BC1G_16062, BC1G_1276 and BC1G_16084. The MFS_1 domains of BC1G_13768, BC1G_12366, BC1G_04779, and BC1G_04656 were predicted. Meanwhile the N-terminal ends of the proteins encoded by BC1G_04779, BC1G_12366 and BC1G_04656 were in the cell, and that of BC1G_13768 was outside of the cell. Each transmembrane domain (TMD) was linked though the ring protection structure. The reverse repetition structure and typical topological characteristics of the MFS transporter protein were found in each domain. Therefore, BC1G_12366, BC1G_04779 and BC1G_04656 belonged to the MFS gene (Supplementary Figure 4), which was widely related to multidrug resistance (MDR).
The cis-acting elements from 1000 bp upstream to 100 bp downstream of the initiation codons for the significant upregulated genes, including BC1G_12765, BC1G_16062, BC1G_16084, BC1G_04656, BC1G_04779, BC1G_12366, BC1G_09151 and BC1G_04893, were predicted and found that each had a GC-rich Sp1 element. The Sp1 element was a part of the Zinc-coordinating element and could be regulated by beta-beta-alpha zinc finger proteins. BC1G_02313 (Zn2Cys6), BC1G_04879 (Zn2Cys6), BC1G_10483 (zinc finger, CCHC-type) and BC1G_13764 (Zn2Cys6) were demonstrated to bezinc finger proteins (Supplementary Figure 5). The GAL4 (Zn2Cys6 (or C6 zinc) binuclear cluster of DNA-binding domains was discovered in the BC1G_02313, BC1G_04879 and BC1G_13764 genes and participated in the regulation of transcription with cis-acting elements. However, the ZnF_C2HC domain was found in the BC1G_10483 gene and participated in the regulation of transcription via cis-acting elements. A coiled coil, beta-betaalpha zinc finger was also anchored on the element. Therefore, we hypothesize that BC1G_10483 is a very important element in response to pressures from cyprodinil and fenhexamid. Based on the network analysis of protein interaction, there were no DEPs, and the expression of BC1G_10483 could be regulated by the synergism of multiple transcription factors.

DISCUSSION
Fungicide resistance has been coming the key problem in crop protection worldwide in the two decades, and has resulted in the greatly decreasing the efficacy, and the increasing concentration to ensure efficacy has produced other problems, such as environment pollution, excessive residue in food, successively strengthening the resistance and increasing the cost of Integrated Pest Management (IPM) strategies (Brent, 2012). Therefore, it is necessary to evaluate the resistance risk in the lab before newer chemicals are wildly put into practice. In this experiment, we continuously screened the B05.10 strain with different concentrations of cyprodinil and fenhexamid and found that the toxicities steadily increased. The cell membrane permeability and glycerinum contents were also significantly increased in those strains whose toxicities have been clearly weakened; in contrast, the oxalic acid (OA) contents markedly decreased. The results were consistent with the early report by Duan et al. (2013). According to the results of the pathway of enrichment analysis, BC1G_09151, a speed-limited enzyme of citrate synthase in the tricarboxylic acid (TCA) cycle, was found to be expressed in the T2 and T3 strains and to decrease the oxaloacetic acid content, being is the main precursor in the synthetic process of OA.OA is known to play a very important role in pathogenesis and fungal development (Godoy et al., 1990). Duan et al. (2013) found that when the resistant strain of S. sclerotiorum were dealt with fludioxonil, the OA content and virulence of pathogens decreased significantly, then finally resulted in the failure of infection. Meanwhile, BC1G_15612, a gluconolactonase gene, was also significantly up-regulated in the glycerol synthetic pathway in the T2 and T3 strains, and its upregulated expression increased the contents of 6phosphoric acid glucuronic acid and glyceraldehyde 3-phosphate, ultimately leading to the accumulation of glycerinum. In addition, the upregulated expression of BC1G_13768 promoted the degradation of glycerol phospholipids to glycerophosphate and increased the accumulation of glycerinum in the T2 and T3 strains. Glycerinum is an important factor that can modulate osmotic pressure in cells. In our experiment, the glycerinum contents accumulated in the T2 and T3 strains with the screened pressure of cyprodinil and fenhexamid, whereas the glycerol synthesis rate gradually decreased. This result was consistent with the development of resistance, which was mainly attributed to the enhanced MFO activities to promote the metabolism of fungicides and transmembrane proteins to decreasing fungicide transport into the cell. Some research has demonstrated that fludioxonil interferes with osmotic signal transduction and induces glycerol synthesis in Neurospora crassa (Pillonel and Meyer, 1997), Candida albicans (Ochiai et al., 2002), and S. sclerotiorum (Duan et al., 2013). Our results have suggested that the inhibition of cyprodinil and fenhexamid on B05.10 has induced those strains to produce more glycerinum. At the same time, the changes in the virulence of pathogens in our strains on the host are still worthy of future research. In recent years, transcriptome analyses has been widely applied in research on fungicidal resistance in some plant pathogens, including Mycosphaerella graminicola (Cools et al., 2007), Fusarium graminearum (Liu et al., 2010), and Cucumis sativus (Wu et al., 2013). Based on our transcriptome data, the expression of genes encoding MFO, transmembrane protein, and zinc finger protein in the T2 and T3 strains was clearly up-regulated through successive treatment with cyprodinil and fenhexamid.
P450 genes are superfamily, which widely distributed in proand eukaryotic organisms and were related to the biosynthesis of many biologically important compounds, such as hormones, fatty acids, and sterols (Guengerich, 2004). Becher et al. (2011) reported that the resistance in fungi to azole fungicides mainly relied on the mutation of P450 gene CYP51 encoding cytochrome P450-dependent sterol demethylase (P450 14αdm ), and resulted in the loss of demethylation activity. In our experiment, RNE of BC1G_16062 in the T1 strain screened with fenhexamid was found to be significantly promoted. The high homologies of BC1G_16062 with CYP37B1 (participating in the metabolism of drugs and steroids in Caenorhabditis elegans) (Laing et al., 2012), BC1G_12765 and BC1G_16084 with CYP4 (which played an important role in xenobiotic biotransformation and in modulating the concentrations of eicosanoids during inflammation as well as in the metabolism of clinically significant drugs) (Bylund et al., 2001;Kalsotra and Strobel, 2006) were also demonstrated (Supplementary Figure 6).
The transporter in fungi, removing the intracellular fungicide, is an important resistant mechanism for fungicides and includes two types: one is an ATP binding cassette (ABC), such as CDR genes, and the other is the major facilitator superfamily (MFS), such as the CaMDR genes. ABC transporters are involved in the energy-dependent efflux of sterol demethylation inhibitors (DMIs), which have been described for Aspergillus nidulans, C. albicans, M. graminicola, etc., (Kerr, 2002;Piddockv, 2006). If the charged molecules are unidirectionally pumped as a consequence of the consumption of a primary cellular energy source, electron chemical potential results, which could be used to drive the active transport of additional solutes via secondary carriers (Ren and Paulsen, 2005). BC1G_13768, BC1G_12366, BC1G_04779, and BC1G_04656 were widely related to multidrug resistance (MDR). The development of MDR based on a single mechanism of resistance, such as the overexpression of genes encoding drug efflux transporters (ABC and MFS), has been observed in isolates of B. cinerea resistant to different classes of fungicides (Vermeulen et al., 2001;Moyano et al., 2004;Myresiotis et al., 2007). The expression of the BcatrD gene, an ABC transporter gene, is a determinant of the sensitivity of B. cinerea to DMI fungicides (Hayashi et al., 2002). The biosynthesis of sterol in B. cinerea was stopped by cyprodinil restraining the activity of the 3-keto reductase gene ERG27 in the C-4 methylation process (Albertini and Leroux, 2004). Fenhexamid inhibits the activity of cystathionine β-lyase and hinders the biosynthesis of methionine (Leroux et al., 2002). According to the data analysis, the SNP sites encoding the genes of the 3-keto reductase gene ERG27 and cystathionine β-lyase were not found in the T02 and T03 strains. FIGURE 7 | qRT-PCR analysis of significant DEGs in T02 and T03 strains. The asterisk * and * * designate statistically significant differences (P < 0.05) and extremely significant differences (P < 0.01) for each DEG between three strains, respectively. T01, T02, and T03 stand for the BO5.10 strain with no fungicide treatment, F 27 strains successively screened with EC 50 concentration of cyprodinil or fenhexamid, respectively. The relative normalized expressions in the DEGs of coding the MFO (A), transmembrane protein (B), Zinc finger proteins (C), and other functional genes (D).
Therefore, we hypothesized that the B05.10 strain screened with cyprodinil and fenhexamid contained base mutations at target points.
Sp1 is a stress response element of −50 bp in cells and plays a key role in adapting to adverse ecological conditions, including the selective pressure of fungicides (Chang et al., 2017). SP1 can modulate the expression of surviving cells through the MAPK signaling pathway to regulate the drug resistance of leukemia stem cells . It is a transcription activating factor of zinc finger-rich glutamine that regulates the expression of target genes involved in transcription regulatory factors, contains the structure of the zinc finger in the carboxyl end and especially combines the GC box in DNA promoter (Pearson et al., 2008).The overexpression of EGR1 (zif268, zinc fingers) was able to activate a construct containing tandem of the MDR1 promoter at SP1 sites and increased its expression in human acute and chronic leukemias (McCoy et al., 1995). The active zone of the zinc finger exists in the middle, participates in regulating the expression in the target genes and combines the other transcription factors. The combination ability between zinc finger and DNA, the activity of transcription and modulation domain and the content of zinc finger in cells have played a key role in the response to external stress factors. The zinc finger protein of MSN2 or MSN4 acts as a transcriptional activator in Saccharomyces cerevisiae to activate transcription via response to stress (Fabrizio et al., 2004), and the zinc finger protein of seb1 also responds to high osmotic stress in Trichoderma atroviride (Peterbauer et al., 2002).

CONCLUSION
The toxicities, conductivities and glycerinum contents of the B. cinerea BO5.10 strain with the successive pressure of cyprodinil and fenhexamid increased. Some key genes encoding MFO, transmembrane proteins, zinc finger proteins, citrate synthase enzyme, gluconolactonase, and lysophospholipid enzymes have been significantly upregulated, and their functional prediction analysis indicated that BC1G_10483 was a ZnF_C2HC transcriptional regulator by means of interaction with the Sp1 element of these genes in response to pressures from cyprodinil and fenhexamid. Our conclusions were necessarily tenuous and further studies were required since the information about the resistant level of B. cinerea on traditional chemicals, and whether cross-resistance between those chemicals and the newer tested fungicides or not were not clearly. However, all results could contribute to a better understanding of the resistance mechanism of B. cinerea on cyprodinil and fenhexamid and provided a foundation for designing the IPM strategy for the effectively controlling gray mold in the field.

AUTHOR CONTRIBUTIONS
CG conceived and designed the experiments. CG, YZ, LS, and XW performed the experiments. XW and CG wrote and revised the paper. All authors approved the final version of the manuscript.