Transcriptome Analyses Shed New Insights into Primary Metabolism and Regulation of Blumeria graminis f. sp. tritici during Conidiation

Conidia of the obligate biotrophic fungal pathogen Blumeria graminis f. sp. tritici (Bgt) play a vital role in its survival and rapid dispersal. However, little is known about the genetic basis for its asexual reproduction. To uncover the primary metabolic and regulatory events during conidiation, we sequenced the transcriptome of Bgt epiphytic structures at 3 (vegetative hyphae growth), 4 (foot cells initiation), and 5 (conidiophore erection) days post-inoculation (dpi). RNA-seq analyses identified 556 and 404 (combined 685) differentially expressed genes (DEGs) at 4 and 5 dpi compared with their expression levels at 3 dpi, respectively. We found that several genes involved in the conversion from a variety of sugars to glucose, glycolysis, the tricarboxylic acid cycle (TAC), the electron transport chain (ETC), and unsaturated fatty acid oxidation were activated during conidiation, suggesting that more energy supply is required during this process. Moreover, we found that glucose was converted into glycogen, which was accumulated in developing conidiophores, indicating that it could be the primary energy storage molecule in Bgt conidia. Clustering for the expression profiles of 91 regulatory genes showed that calcium (Ca2+), H2O2, and phosphoinositide (PIP) signaling were involved in Bgt conidiation. Furthermore, a strong accumulation of H2O2 in developing conidiophores was detected. Application of EGTA, a Ca2+ chelator, and trifluoperazine dihydrochloride (TFP), a calmodulin (CaM) antagonist, markedly suppressed the generation of H2O2, affected foot cell and conidiophore development and reduced conidia production significantly. These results suggest that Ca2+ and H2O2 signaling play important roles in conidiogenesis and a crosslink between them is present. In addition to some conidiation-related orthologs known in other fungi, such as the velvet complex components, we identified several other novel B. graminis-specific genes that have not been previously found to be implicated in fungal conidiation, reflecting a unique molecular mechanism underlying asexual development of cereal powdery mildews.


INTRODUCTION
Bgt and Blumeria graminis hordei (Bgh) are economically important pathogens of cereals that can cause devastating damage to wheat and barley (Dean et al., 2012). Diseased plants display a white powdery mold on the leaves and stems and suffer from deprivation of nutrients by these fungi. B. graminis can reproduce not only by asexual conidia, but also by sexual crossing, which take places between two isolates with opposite mating types to generate ascospores. These two types of spores serve as primary inoculums and are spread by wind. Fungicides and resistant cultivars can be utilized in the management of these diseases. Both Bgt and Bgh are obligate biotrophic parasites that can complete their life cycles only on living hosts, and this characteristic limits the advancement of molecular functional analyses by genetic transformation. Nevertheless, significant progress in these analyses has been accomplished via global large-scale ("-omics") approaches such as genomics, transcriptomics, proteomics, and metabolomics (Bindschedler et al., 2016). Insights from these studies provide useful information for future genetic resistance improvement and fungicide development used in protection from these diseases.
Conidiation is one of the most important modes of reproduction in powdery mildews (Glawe, 2008). This process usually occurs several times during a growing season of plants, yet sexual recombination often takes place only once a year. Haplogroup analysis of the Bgt genome revealed that clonal reproduction occurs more frequently than sexual recombination (Wicker et al., 2013). Large numbers of the conidia of Bgh can be produced within several hours (Jørgensen and Torp, 1978), and can be transmitted by wind from continental Europe to Britain (Brown and Hovmoller, 2002), contributing to population migration and rapid spread of mildew disease. Therefore, understanding their molecular mechanism of conidiogenesis of B. graminis may be helpful in developing better strategies for control of these diseases.
The conidiogenesis process of B. graminis consists of the formation of foot cells, the development of conidiophores and then the release of conidia. Similar to other fungi, like Aspergillus nidulans and Magnaporthe oryzae, cereal powdery mildews must go through a period of vegetative hyphae growth after successful infection. Subsequently, conidiation typically begins when foot cells arise from vegetative hyphae and then conidiophores are erected by repeated elongation and septation of foot cells to produce massive conidia (Moriura et al., 2006). Moreover, new conidia successively form on the previous conidia generated from foot cells, thereby forming conidia chains (Glawe, 2008). Thus, foot cells control the shift from vegetative growth to asexual development and the speed of new conidial cell production of B. graminis. Uncovering the changes in metabolism and regulation during foot cell and conidiophore differentiation will contribute to better understanding the conidiation mechanism of powdery mildews.
Despite the agricultural importance of B. graminis, little work has been done to explore the genetic basis of asexual sporulation in powdery mildews. Bgh undergoes a series of dynamic changes in carbohydrate and lipid metabolism during asexual reproduction (Both et al., 2005). cDNA amplified fragment length polymorphism analysis identified 620 differentially expressed fragments related to metabolism and signaling during conidiation of Erysiphe necator (Wakefield et al., 2011). Conversely, significant advances in studies on asexual reproduction with an emphasis on regulatory factors have been made in several other fungi, indicating that there is a central regulatory pathway, the BrlA pathway, involved in conidiogenesis regulation of A. nidulans and A. fumigatus. Moreover, some upstream development activators (e.g., FluG, FlbC, FlbD, and FlbE) and downstream development regulators (e.g., velvet complex consisting of VosA, VelB, VelC, and LaeA) also act in regulating conidiation (Park and Yu, 2012;Alkhayyat et al., 2015). Additionally, at least two G protein mediated signal transduction pathways are involved in A. nidulans conidiation by negatively regulating the BrlA pathway (Yu, 2006). Conidiationrelated heterotrimeric G protein signaling components include G protein coupled receptors (GPCRs), heterotrimeric G proteins and regulators of G protein signaling proteins (Park and Yu, 2012;Wang et al., 2013). Recently, a few small monomeric GTPases have also been shown to be implicated in various cellular processes including conidiation and tolerance to multistressors such as H 2 O 2 (Guan et al., 2015). External signals from upstream G protein signaling can be transduced to some downstream central regulatory pathways that govern conidiation. It is known that the Ca 2+ mediated signaling pathway (Nguyen Q.B. et al., 2008;Li et al., 2015), mitogen-activated protein kinase (MAPK) signaling (Chen et al., 2016) and reactive oxygen species (ROS) signaling (Chung, 2012;Viefhues et al., 2014) are also correlated with conidiogenesis. Furthermore, the crosslink between MAPKs cascade and Ca 2+ signaling as well as that between MAPKs cascade and ROS signaling have been found in conidiation regulatory pathways (Huang et al., 2015;Wei et al., 2016).
The main purpose of this study was to identify the central metabolic and regulatory events during sporulation of Bgt. The characteristics of obligate biotrophy that underpin B. graminis development have historically made it difficult to design experiments that expand our understanding of these processes. Here, we performed a RNA-seq analysis of the epiphytic structures of Bgt during conidiation and conducted additional histological and pharmacological investigations. The resulting data provide new insights into the molecular mechanism underpinning conidiogenesis of B. graminis.

Sample Collection and RNA Preparation
The epiphytic structures of isolate 21-2, growing on leaf segments of Chancellor at 3, 4, 5 dpi were embedded in the acetate cellulose and dissected with a tweezer. These samples were chilled in liquid nitrogen immediately and stored at −80 • C. Total RNA of samples with three biological replicates at each time point was extracted using the RNAiso Plus Purification Kit (TAKARA, Japan) following the manufacturer's instruction. After RNA purity was checked using the NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Wilmington, DE, United States), the Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer (Life Technologies, Foster City, CA, United States) was employed to determine concentration of these RNA samples. Then RNA integrity was confirmed by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). Two sets of qualified RNA samples were obtained under the same conditions. One set was used for RNA sequencing, while the other set was used for quantitative reverse-transcriptase PCR (qRT-PCR).

Libraries Construction and Sequencing
Before mRNA sequencing, poly (A) mRNA was isolated from total RNA using the NEB-E7490 Poly (A) mRNA Magnetic Isolation Kit (NEB, United States), and then fragmented by the Fragmentation Kit (Ambion, Austin, TX, United States). cDNA libraries were constructed using the NEBNext R Ultra TM RNA Library Prep Kit for Illumina R (NEB, United States) according to the manufacturer's recommendations. Briefly, the first-strand and Development stage name: a time course of Bgt infection to Chancellor. Library name: RNA samples of the epiphytic structures growing on leaf segments of Chancellor at 3, 4, 5 dpi with three biological replicates for library construction; Total reads: number of total pair-end reads in each clean data; Total bases: total base number in each clean data; GC content: percentage of G and C in clean data; Q20(%): percentage of bases with quality value more than 20 in each clean data; Q30(%): percentage of bases with quality value more than 30 in each clean data.
the second of cDNA were synthesized using M-MuLV reverse transcriptase and DNA polymerase I (Invitrogen), respectively. After NEBNext adaptors were ligated to the adenylated 3 ends of cDNA, fragments of approximately 500 bp in length were purified and selected with the AMPure XP system (Beckman Coulter, Beverly, MA, United States). Next, these fragments were enriched by PCR amplifications and library quality was assessed by the Agilent Bioanalyzer 2100 System. RNA sequencing was performed using a paired-end 100 (each end with 100 bases) strategy on the Illumina Hiseq 2500 platform at Biomarker Technology, Co., Ltd (Beijing, China) to generate 2 billion bases per sample. We deposited the raw data in the NCBI sequence read archive under accession number SRP092295.

Reference Genome-Based Reads Mapping and Quantitative Analysis of Gene Expression
All reads of nine libraries were mapped to the reference genome of Bgt isolate 96224 (Wicker et al., 2013) using the Tophat2 software (Kim et al., 2013), allowing up to two base mismatches. To quantify the transcript abundance of each expressed gene, fragments per kilobase per million reads (FPKM) values were calculated by the Cufflinks software (Trapnell et al., 2010). Transcript abundances among different libraries were compared using the DESeq software (Anders and Huber, 2010) and the multiple hypothesis test was conducted to check the difference significance of gene expression with a p-value. The p-values were adjusted using the Benjamini and Hochberg's approach by controlling the false discovery rate (FDR) < 0.01. Genes with an adjusted p-value < 0.05 and | expression fold change| ≥ 2 were deemed with defined DEGs. In addition, the Pearson correlation coefficients (r 2 ) of FPKM between replicates were calculated using the method of Simple Error Ratio Estimate (Schulze et al., 2012).

DEGs Annotation
To characterize the detailed functions of DEGs, nucleotide sequences were aligned against various protein databases including the SwissProt database and the NCBI non-redundant (Nr) database via BLASTX (threshold set to E-value < 1e-5).
The BLASTX results were analyzed using the Blast2GO 1 software (Gotz et al., 2008) to conduct a Gene Ontology (GO) functional analysis. Annotations were enriched and refined using the TopGo package 2 . The DEGs were also mapped to the Clusters of Orthologous Groups (COG) database (Tatusov et al., 1997) to predict and classify functions, and to the KEGG database 3 to assign metabolic pathways. For enrichment analysis of GO terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, two tailed p-values were obtained by Fisher's exact test. Terms and pathways with FDR corrected p-values < 0.05 were considered to be significantly enriched.

qRT-PCR Validation
To confirm the reliability of transcript levels that were quantified based on the ratio of FPKM values, the relative expression levels of 22 DEGs were checked using qRT-PCR evaluation. After having been digested by DNase I (RNase free, TAKARA, Japan), first-strand cDNAs were synthesized from 1 µg total RNA using the PrimeScript TM First-strand cDNA Synthesis Kit (TAKARA, Japan) following the manufacturer's protocol. Specific primers for each gene used in qRT-PCR were listed in Supplementary  Table S1. All qRT-PCR reactions were conducted using the SYBR Green I Nucleic Kit (Life Technologies, United States) and the ABI 7500 real-time PCR system (Applied Biosystems, United States) with the following programs: 95 • C for 10 min, followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s. The housekeeping gene encoding beta-tubulin was used as the internal control (Hacquard et al., 2013). All PCR amplifications were run in triplicates. The 2 − CT method (Livak and Schmittgen, 2001) was used to calculate the relative expression levels.

Identification of Genes Involved in Primary Metabolic Processes and Regulatory Pathways and B. graminis-Specific Proteins
To identify genes acting in primary metabolism among DEGs, we extracted genes that are affiliated with annotations in the COG database corresponding to carbohydrate transport and metabolism (E), energy production and conversion (C) and lipid transport and metabolism (I) at first. Additional genes were identified based on associated functional annotations in the GO, the KEGG, SwissProt, and Nr databases. All these genes were combined and checked by BLASTP search in the GenBank manually.
To identify putative signaling proteins, we analyzed major fungal signaling pathways including heterotrimeric G proteins/GPCRs, small G proteins, Ca 2+ , PIP, and H 2 O 2 mediated signaling. Since almost 92% of predicted Bgt genes have homologs in Bgh genome (Wicker et al., 2013), protein sequence similarity alignments against putative signaling proteins in Bgh genome (Kusch et al., 2014) were performed by BLASTP. All hits with a cutoff of E-value < 1e-10 were considered significant. Subsequently, additional candidates (such as members involved in PIP and H 2 O 2 mediated signaling) were identified by searching associated annotations in the GO, KEGG, SwissProt, Nr, and InterPro databases. GPCR candidates were further corroborated by the presence of six to eight trans-membrane domains 4 . For transcription factor (TF) mining, query sequences for BLASTP were obtained from the Fungal Transcription Factor Database 5 (Park et al., 2008). Additional TFs were also predicted based on their annotations in the database as described above. Finally, a manual inspection using BLASTP search in the NCBI website was conducted to check the combined results gene by gene.
To identify the proteins specific to Bgt or Bgh, we conducted BLASTP searches against the Nr database using standard parameters. The proteins without significant hits (E-value ≥ 1e-5 to proteins in other fungi, except Bgh) were considered as B. graminis-specific proteins. To obtain more annotations for these B. graminis-specific proteins, BLASTP alignments against the proteins of Bgh isolate DH14 (Spanu et al., 2010) were performed (E-value cutoff 1e-10).

EGTA and TFP Treatment
EGTA (GEN-VIEW company, United States), was dissolved in distilled water (pH8.0). TFP (TCI company, Shanghai, China) was added to distilled water to make a 100 mM stock solution. The surfactant Tween 20 was added to make the final concentration 0.025%. The same volume of 0.025% Tween 20 was used as a negative control. Detached leaf segments were inoculated and incubated under conditions described above. A 2.4 mL solution was sprayed on the infected leaf segments using the Potter spraying tower at pressure of 0.25 MPa. After air-drying for 1 h at 25 • C, incubation of the treated leaf segments resumed under the condition used before treatment.

Staining and Histological Observation
Leaf segments of Chancellor infected by isolate 21-2 were stained with alcoholic lactophenol trypan blue (Yang et al., 2008). The Lugol reagent (KI and I 2 ) staining method (Both et al., 2005) was employed to detect glycogen generation. After decoloration in 95% ethyl alcohol, inoculated leaf segments were stained with the Lugol reagent for 10 min at room temperature and observed immediately. To observe the generation of H 2 O 2 in this pathogen, Bgt-infected leaf segments were stained with 3, 3diaminobenzidine (DAB, Sigma-Aldrich, St. Louis, MO, United States) (Thordal-Christensen et al., 1997). Leaf segments were briefly soaked in 1mg/ml DAB solution (pH3.8), at 25 • C for 8 h and decolored with clearing solution (ethanol: acetic acid: ddH 2 O = 94:4:2, v/v) for 1 h.
To count the number of foot cells at 5 dpi and conidiophores at 6 dpi in single colonies, 30 colonies on each of three leaf segments were randomly selected and observed through a light microscope at 400× magnification. Mycelium expansion width of 30 colonies on each of three leaf segments at 5 dpi was measured using a micrometer at 100× magnification. To evaluate the impact of the treatments with EGTA or TFP on H 2 O 2 production, the percentage of colonies with stained conidiophore (s) at 6 dpi were investigated by scoring the number in 30 randomly chosen colonies on each of three leaf segments. To measure conidia production, each of three 2.5 cm long leaf segments collected at 10 dpi was put into 2.5 mL 0.025% Tween 20 to produce a suspension of spores. Concentration of spores was estimated using a hemocytometer at 100× magnification (Gadoury et al., 2012). These experiments were conducted two times independently.

Statistical Analysis
The analysis of variance was performed to determine whether the error variances of treatment and untreated control were homogeneous. The significance of the differences was evaluated using Student's t-tests.

Asexual Development of Bgt
The asexual life cycle of isolate 21-2 on Chancellor was observed under a light microscope (Figure 1). At the early infection stage, conidia with appressorium generated primary haustoria underneath the infected sites by 1 dpi and haustoria matured by 2 dpi. Then secondary haustoria began to form in the epidermis and vegetative hyphae proliferated on leaf surface by 3 dpi. Subsequently, some swollen cells, called foot cells, appeared on vegetative hyphae by 4 dpi and many foot cells could be observed by 5 dpi. The upper portion of foot cells then began to elongate and repeated septation separated a few conidial cells, which compose conidophores, the specialized developmental structures for conidia generation. By repeating this process, several conidophores started to erect on foot cells in single colonies by 5 dpi and the upper conidial cells enlarged during maturation. By 6 dpi, clustered conidiophores were seen in single colonies. Finally, mature spores were gradually released from the conidiophores.

High-Quality RNA-seq Data Obtained from Epiphytic Tissues of Bgt
To analyze changes in the transcription profiles of Bgt in the transition from vegetative growth to asexual reproduction, the epiphytic materials at 3, 4, and 5 dpi were sampled for transcriptome sequencing. Three biological replicate samples were collected for each phase. For mRNA sequencing, nine cDNA libraries were constructed and then a total of 23.05 Gb of clean data composed of 91,495,902 pair-end reads was generated. Q30% of these clean data from each library were more than 87.8% (Table 1), indicating that these Illumina reads were of a  Table S2). Among the matched reads more than 65% of the reads from the combined data were mapped to exonic regions and almost all of the rest were mapped to intergenic regions (Supplementary Figure S1). Saturation analysis indicated these libraries were sequenced to saturation and represented the transcripts in the conditions tested (approximate 1G clean data, Supplementary Figure S2). Correlation clustering analysis based on FPKM values suggested that our RNA-seq data were highly reproducible (Supplementary Figure S3). To further confirm the reliability of our RNA-seq expression data measured by FPKM fold changes, we selected 22 DEGs and monitored their expression levels using qRT-PCR. These proteins included nine functional proteins, 10 regulatory factors, and three hypothetical proteins with unknown functions. The relative transcript levels determined by qRT-PCR and RNA-seq were positively correlated with a Pearson coefficient R 2 = 0.777 (Supplementary Figure S4). These results support the validity of our RNA-seq data and can be used for downstream analysis.

Differential Expression Analysis of Genes Regulated during Conidiation
To analyze the changes in gene expression across different developmental stages, we performed differential expression analysis at three different time points (3, 4, and 5 dpi) using the DESeq software ( Table 2). 556 and 404 genes were differentially expressed in the comparisons of C4 dpi/C3 dpi, and C5 dpi/C3 dpi, respectively (Supplementary Table S3). Thus we found a total of 685 DEGs, of which 317 genes were up-regulated and 369 were down-regulated (Figure 2). Within these genes, 275 (including one gene coding a protein, EPQ64082.1, down-regulated at 4 dpi but up-regulated at 5 dpi) were differentially expressed in both comparison pairs. Moreover, 281 and 129 specific in the library pairs of C4 dpi/C3 dpi and C5 dpi/C3 dpi, were related to the foot cell and conidiophore development, respectively.
To reveal the potential functions of the 685 DEGs, the COG, GO, and KEGG databases were used to analyze functional classification and pathways. Among the 25 COG categories, except for the cluster for "general function prediction only, " the largest group was the cluster for "carbohydrate transport and metabolism" in both pair-comparisons (26 genes in C4 dpi/C3 dpi and 21 genes in C5 dpi/C3 dpi), followed by the clusters for "post-translational modification, protein turnover, chaperones, " "lipid transport and metabolism, " "amino acid transport and metabolism, " and "energy production and conversion" (Supplementary Figure S5).

(Supplementary
Using KEGG pathway enrichment analysis, we categorized the biological functions of the 685 DEGs. No significant enriched pathway with FDR corrected p-value < 0.05 was observed between library pairs (Supplementary Table S5). As with the GO analysis, pathways with non-corrected p-value < 0.05 (Fisher's exact test) were also taken in account. This analysis showed that the pathway of "peroxisome" (ko04146) was highly represented in both pair-comparisons. Six DEGs were observed in the "peroxisome" pathway, which can be classified into three functional groups, peroxisome biogenesis (two peroxisomal proteins, EPQ61877.1 and EPQ67652.1), fatty acid oxidation (two enzymes, EPQ67195.1 and EPQ66210.1), and antioxidant system (two scavengers, EPQ63972.1 and EPQ62743.1) (Supplementary Figure S7). The results indicate that fatty acid and ROS metabolism might be involved in conidiation.

DEGs of Bgt Involved in Carbohydrate Metabolism, Energy Production, and Unsaturated Fatty Acid Metabolism during Conidiation
Since many DEGs were found to be involved in basal metabolic processes, we used BLAST to search the 685 DEGs against the SwissProt and NCBI Nr databases to further identify enzymes involved in the carbohydrate, energy, and lipid metabolism. The results are listed in Supplementary Table S6. This analysis identified 18 carbohydrate metabolism-related genes. Within these genes, eight genes play roles in metabolizing various sugars into glucose or glucose phosphate. Most of them (7/8) were significantly up-regulated at 4 and/or 5 dpi (Figure 3). For example, two FGGY carbohydrate kinases, which can catalyze phosphorylation of several sugars and their products, connected to glycolysis and the pentose phosphate pathway (Zhang et al., 2011), displayed more than eightfold increases of mRNA levels in the library of C5 dpi. The substrates catalyzed by the seven enzymes included β-1,3-glucan, starch, sucrose, xylulose and fuculose, suggesting that conversion of diverse sugars into glucose become active during conidiation. Furthermore, the expression profiles of genes associated with glycolysis (EPQ67062.1, up-regulated), the pentose phosphate pathway (EPQ66083.1, up-regulated), and gluconeogenesis (EPQ64945.1, down-regulated), suggest that generated glucose might partially be utilized as fuel during Bgt conidiation (Figure 3). In addition, we found five genes related to glycogen biosynthesis to be up-regulated at 5 dpi (Figure 3). This indicates that glycogen synthesis may be activated in conidiophores. To further confirm this result, we performed an iodine staining assay to test whether glycogen was accumulated in conidiophores. As predicted, glycogen began to generate in the foot cells at 4 dpi and accumulated at high levels in the developing conidiophores at 5 dpi ( Figure 4A).
For energy metabolism, 14 genes related to energy production and transport showed differential expression during conidiation. Among these genes, 11 genes, which are mainly involved in the TAC, the glycoxylate cycle, energy transport and the ETC (Figure 3), are significantly up-regulated. These data suggest that energy metabolism is enhanced during conidiation of Bgt.
For lipid metabolism, a total of 29 DEGs were obtained in our transcriptomic analysis. Most DEGs (8/11) involved in triacylglycerol and phosphatidate degradation were significantly down-regulated during sporulation. However, seven genes related to the break down pathway of peroxisomal unsaturated fatty acid, such as peroxisomal 2,4-dienoyl-CoA reductase, enoyl-CoA isomerase (Wanders and Waterham, 2006) and the multifunctional enzymes of the peroxisomal fatty acid betaoxidation pathway (Nguyen S.D. et al., 2008) were up-regulated (Figure 3). Moreover, two genes acting in the glyoxylate cycle and two in the carnitine transport, where acetyl CoA is transported from peroxisome to mitochondria (van Rossum et al., 2016), were also activated (Figure 3). These results suggest that lipid degradation is repressed, while unsaturated fatty acid oxidation is enhanced. We did not find evidence of the activation of enzymes responsible for lipid synthesis. However, a significant increase of the expression level of one gene coding the caleosin protein (EPQ66680.1), which is related to storage lipid bodies, was detected.

Activation of Antioxidant Redox System and Accumulation of H 2 O 2 in Conidiophores of Bgt
Based on the results from the GO and KEGG analysis, we hypothesized that cell redox homeostasis-related processes might be correlated to Bgt conidiation. Therefore, we looked for DEGs with annotations associated with ROS clearance and production in the SwissProt, Nr, and InterPro databases. This analysis resulted in 14 genes coding for putative antioxidant proteins, including two superoxide dismutases, one catalase, two peroxiredoxins, five thioredoxins, two glutathione oxidoreductases, and two glutaredoxins. All these genes, except for the catalase coding gene, exhibited up-regulation patterns at 4 or/and 5 dpi, which indicates that the antioxidant system for some ROS, such as H 2 O 2 , is activated during sporulation. For ROS-producing enzymes, only one putative NADPH oxidase A (NoxA, EPQ61671.1), tightly linked with H 2 O 2 -generating reactions in conjunction with superoxide dismutases, was significantly down-regulated at 4 dpi.
To assess if H 2 O 2 accumulated in sporulating tissues, a DAB staining experiment was carried out. As shown in Figure 4B, H 2 O 2 began to appear in mycelium and foot cells at 4 dpi and bursted in the developing conidiophores at 5 dpi. These results indicate that H 2 O 2 metabolism is correlated with the conidiophore development. By contrast, low quantities of H 2 O 2 were detected in the infected epidermal cells during the timecourse of 3, 4 and 5 dpi, suggesting that little H 2 O 2 was generated by the host cells.

DEGs Involved in Major Regulatory Pathways during Bgt Conidiation
To identify DEGs involved in major signaling pathways and transcriptional regulation, we used a combined strategy of  sequence alignments against known regulators in other fungi and searches of genes with corresponding functional annotations in several databases (the GO, KEGG, SwissProt, Nr, and InterPro databases). The known regulators included those in Bgh (Kusch et al., 2014) and TFs in the Fungal Transcription Factor Database. In addition, prediction of transmembrane helices was performed in putative GPCRs. A total of 91 DEGs were identified, which are involved in a variety of regulatory pathways (Figure 5 and Supplementary Table S7). The gene expression profiles were grouped based on their FPKM using hierarchical clustering. This analysis revealed three major clusters exhibiting distinct expression patterns during conidiation (Figure 6). Cluster A contained only one up-regulated glutathione oxidoreductase coding gene. Cluster B consisted mostly of members in G-protein FIGURE 6 | A heatmap for the expression patterns of 91 putative regulator coding genes in Bgt at 3, 4, and 5 dpi. Regulatory pathways in which these regulators participated in include the G protein signaling, small G protein signaling, calcium signaling, H 2 O 2 signaling, PIP signaling and transcriptional regulation. A hierarchical clustering analysis revealed three major clusters, Cluster A, Cluster B, and Cluster C. The detail information of these genes is listed in Supplementary Table S7. G-protein-beta, heterotrimeric G protein beta subunit; RGS, regulator of G protein signaling; GPCR, G protein coupled receptor; SmallG, small G protein; GEF, guanine nucleotide exchange factor; SOD, superoxide dismutase; CAT, catalase; Prx, peroxiredoxin; Trx1, thioredoxin; GSHR, signaling and some TF genes, which were down-regulated at 4 or/and 5 dpi. In contrast, genes in cluster C were markedly induced and most of them code for regulators in Ca 2+ , H 2 O 2 , and PIP signaling. For the G protein signaling pathways, a total of 21 regulators were identified in silico. We found one heterotrimeric G protein, five GPCRs, and one RGS in the heterotrimeric G protein signaling pathway. For another type of G proteins, Small G proteins, we found nine putative GTPase activating proteins, one GTP binding protein, and four guanine nucleotide exchange factors. Among these 21 DEGs, 17 were down-regulated during sporulation and most down-regulation (11/17) took place at 4 dpi.
We also observed significant transcriptional changes associated with components of Ca 2+ and PIP mediated signaling pathways, which can sense signals from the G protein signaling pathways. In this analysis, we found up-regulation of one phospholipase C (EPQ61894.1) coding gene. Since phospholipase C hydrolyses phosphatidyl inositol 4,5bisphosphate to form inositol 1,4,5-trisphosphate (IP3) and diacylglycerol, two secondary messengers that involved in Ca 2+ influx from other organelles or extracelluar space to cytosol, we hypothesize the cytoplasmic Ca 2+ level might be increased. This is supported by two findings as follows: first, five CaM Ca 2+ binding proteins coding genes were upregulated at 4 and/or 5 dpi, because these proteins are the major transducers of cytoplasmic Ca 2+ in eukaryotic cells. Second, protein sequence alignments revealed that a serine/threonine protein kinase (EPQ63189.1) showed a significant similarity to the Cmk1p (CCU77137.1), a putative Ca 2+ /CaM-dependent protein kinase (CCaMK) in Bgh (E-value = 1e-24). This gene also exhibited significantly increased mRNA levels at 4 and 5 dpi. These results suggest the changes in Ca 2+ concentration are implicated in Bgt sporulation. However, no DEGs encoding for Ca 2+ channel proteins, Ca 2+ exchangers, Ca 2+ -transporting ATPases or calcipressin were found. In PIP mediated signaling, the increased mRNA levels of genes coding the phospholipase C and three PIP binding proteins suggest that PIP metabolism is enhanced on the cytosolic side of cell membranes.
To further analyze the DEGs involved in the cell response to H 2 O 2 and signal transduction, we searched for putative regulators related to H 2 O 2 signaling. In addition to the NoxA (EPQ61671.1) and 14 antioxidant proteins mentioned above, this analysis identified several other orthologous genes with putative roles in ROS signaling in Alternaria alternata (Chung, 2012), including one histidine kinase (HSK, EPQ67055.1), two TFs homologous with the redox-sensitive YAP1 (KZV08838.1) in yeast and a MAPK (EPQ63189.1) homologous to the protein kinase in high osmolarity glycerol pathway (HOG1) in Bgh (CCU77378.1). In total we identified 19 DEGs involved in ROS signaling, of which 14 were up-regulated at 4 or/and 5 dpi.

Effects of Ca 2+ Chelator and CaM Antagonist on Conidiation of Bgt and Production of H 2 O 2
To further analyze the role of Ca 2+ and CaMs in regulating conidiation, EGTA, a extracellular Ca 2+ chelator, and TFP, a CaM antagonist, were applied on wheat leaf infected with the powdery mildew isolate. A series of concentration of EGTA (1, 2, 3, 4, 5, 6, and 7 mM) and TFP (6.25,12.5,25,50,100,200 µM) were sprayed at 1 dpi and the haustorium formation rates were investigated at 3 dpi. We found no significant difference of the haustorium formation rates between the treatments (1, 2, 3, 4, 5 mM EGTA; 6.25, 12.5 µM TFP) and untreated control ( Figure 7A). However, when applied with 5 mM EGTA or 12.5 µM TFP at 4 dpi, the two inhibitors hampered foot cell and conidiophore formation significantly (p < 0.01, Figures 7B,C,  8A), thereby reduced conidia production significantly (p < 0.01, Figures 7D, 8C). In addition, only a slight decrease of mycelium growth was observed (p < 0.05, Figure 7E). These results indicate that extracellular Ca 2+ and intracellular CaMs are involved in Bgt conidiation. Meanwhile, to assess a potential influence of the two treatments on H 2 O 2 accumulation, DAB staining was carried out at 5 and 6 dpi. The rates of colonies with DAB stained conidiophore (s) were significantly lower in treated samples than the control (Figures 7F, 8B). Similar results from another independent experiment are presented in Supplementary Figure S8. These data indicate that chelation of Ca 2+ and antagonism of CaMs inhibit H 2 O 2 generation and that both H 2 O 2 and Ca 2+ levels are tightly linked with asexual sporulation of Bgt.

B. graminis-Specific Protein-Coding Genes
To identify proteins specific to B. graminis among the 685 proteins, we performed BLASTP alignments against the Nr database with an E-value cut off of 10 −5 . In total, 79 proteins without significant hits except from hits to other Bgh proteins were obtained (Supplementary Table S8). We did not find any annotations of these proteins in the COG, GO, KEGG and SwissProt databases. To gain more information, we aligned the 79 proteins against the Bgh protein database using BLASTP (Evalue ≤ 1e-10). We found 47 putative effector proteins and 27 unknown hypothetical proteins. Among the 27 hypothetical genes, several genes were dramatically up-regulated during Bgt sporulation, for instance, two genes (EPQ64343.1, 18.5-fold rise in C5 dpi and EPQ62428.1, 13.3-fold rise in C4 dpi).

Primary Metabolism in B. graminis during Conidiation
Asexual sporulation of Bgt is a complex process correlated with nutrient metabolism, energy consumption and fuel storage, which have been previously shown in other powdery mildews (Both et al., 2005;Wakefield et al., 2011). Our results indicate that besides glucose (reported previously; Both et al., 2005), other sugars, such as starch, sucrose, xylulose and fuculose, could serve as carbon sources for the foot cell and conidiophore development. This means that B. graminis may use specific sugars for conidiation like other fungi, for example, Hypocrea atroviridis, which has carbon-source predilection in conidiation (Friedl et al., 2008). It is well-known that energy is generated through glycolysis in cytosol and the TAC in mitochondria where acetyl-CoA is completely oxidized. Acetyl-CoA can be produced from the unsaturated fatty acid degradation and then enter the TAC so as to generate more energy. These results agree with the recent finding that deletion mutants of PEX6 in Coniothyrium minitans cannot utilize oleic acid and produce considerably fewer conidia than the wild type (Wei et al., 2013). These results indicate that, compared with vegetative growth, Bgt needs more energy, which could come from glucose and acetyl-CoA oxidation, during conidiation. A major change in anabolism was conversion of glucose into glycogen during the conidiophore development. In addition, up-regulation of the caleosin protein-related gene occurred, which suggests that storage lipid bodies might be associated with Bgt conidiation. Glycogen and lipids are universally used for energy storage in fungal conidia. For example, glycogen and lipids have been detected in Bgh conidiophores as fuel for the next infection (Both et al., 2005). Caleosins are strongly correlated with lipid body maintenance and can be used as storage compounds in plant seed (Partridge and Murphy, 2009). Although diverse fatty acids have been detected in conidia of B. graminis (Muchembled et al., 2005), we have no direct evidence from the RNA-seq data that fatty acid biosynthesis is activated. This will need further studies focused on lipid FIGURE 7 | Effect of 5 mM EGTA and 12.5 µM trifluoperazine dihydrochloride (TFP) on conidiation and H 2 O 2 accumulation during conidiation of Bgt. Wheat leaf segments of Chancellor were inoculated with isolate 21-2. A pilot experiment with a series concentration of two chemicals (dissolved in 0.025% Tween 20) was carried out by spraying at 1 dpi and the haustorium formation rates were investigated at 3 dpi (A). Based on these data, 5 mM EGTA and 12.5 µM TFP were applied at 4 dpi by spraying the solutions. Infected leaf segments were sampled at 5 and 6 dpi and stained with the trypan blue and DAB solution. Next, a microscopy observation was carried out. Number of foot cells per single colonies at 5 dpi (B), number of conidiophores per single colonies at 6 dpi (C), conidia concentration of 2.5 mL conidiospore suspension dislodged from 2.5 cm long diseased leaf segments at 10 dpi (D), mycelium expansion width of single colonies at 100× magnification at 5 dpi (E) and percentage of colonies with DAB stained conidiophore(s) at 5 dpi (F) were investigated with three replicates. For each replicate, 30 colonies were observed. All the experiments were performed two times with similar results. Data from another independent experiment are presented in Supplementary Figure S8. Values are means of data with three replicates of each treatment and error bars indicate standard error of the means. The significance of the differences between the treatments and untreated control was determined by Student's t-tests. Significant decreases in the rates of foot cell and conidiophore formation and conidia production after 5 mM EGTA and 12.5 µM TFP treatment occurred ( * * p < 0.01). Both treatments have a slight impact on the mycelium growth ( * p < 0.05). and fatty acid biosynthesis during conidia formation and release.

Comparison of the Genetic Control of Conidiation between Bgt and Other Fungi
It is known that several important signal pathways such as the BrlA pathway, G-protein signaling, ROS and Ca 2+ mediated signaling have essential roles in coordinating conidiation in other fungi (Nguyen Q.B. et al., 2008;Park and Yu, 2012;Medina-Castellanos et al., 2014). Some orthologous genes in these pathways are differentially expressed during conidiation of powdery mildews, including E. necator (Wakefield et al., 2011). Our transcriptome analysis also identified some orthologs of the regulatory factors, such as the velvet complex components, that participate in these pathways. Therefore, the genetic control of conidiogenesis in B. graminis shares conserved elements known from other fungi. In addition, oxylipins, products produced from unsaturated fatty acid oxidization, have been known to function as signal molecules in fungal conidiation (Herrero-Garcia et al., 2011;Scala et al., 2014). Moreover, in biosynthesis of oxylipins, PpoA, one of the conserved fatty acid oxygenases plays a crucial role in regulating conidiation (Brown et al., 2009). Activation of the oxylipins biosynthesis-related gene in Bgt might suggest that the oxylipin pathway is also implicated in conidiogenesis of powdery mildews. Nevertheless, powdery mildews may recruit their unique proteins with unknown function as well, based on the result that several B. graminisspecific hypothetical genes were significantly up-regulated during conidiation. Functional characterization of these genes may contribute to uncovering unusual conidiogenesis mechanism employed by powdery mildews.

Roles of H 2 O 2 in Conidiogenesis
H 2 O 2 has dual roles, acting as both a signal molecule and a stressor during Bgt asexual reproduction. Normal conidiation is correlated with a strong accumulation of H 2 O 2 , while impaired conidiation is associated with reduced H 2 O 2 levels. This suggests that H 2 O 2 is tightly linked with Bgt conidiogenesis. This observation is in accordance with the previous finding in C. minitans mutants involved in the disruption of the CmPEX6 gene result in a significant decrease in the amount of H 2 O 2 , as well as reductions in conidia production compared to the values of the wild type (Wei et al., 2013). H 2 O 2 also plays an important role in activating the oxidative stress response and triggering conidiation, partially via regulating the redoxresponsive regulators (YAP1 and SKN7) and the HOG1 mediated pathway in A. alternata (Chung, 2012). In this pathway, YAP1 acts in H 2 O 2 sensing and activation of antioxidants and the H 2 O 2 signal can be transmitted through HOG1, which is regulated by NoxA and HSK. Comparatively, based on the expression profiles of the NoxA, the HSK, two YAP1s and the HOG1 during Bgt conidiation, it seems that the H 2 O 2 signal is likely transmitted through HOG1 as well. However, unlike the situation in A. alternata, the NoxA and the YAP1s in Bgt appear not to sense and respond to H 2 O 2 and the NoxA and the HSK seem not to be associated with transcriptional activation of HOG1. Further functional analysis of these genes will test this hypothesis. It appears, however, that the concentration of cellular H 2 O 2 must be strictly regulated in time and space, since excessive H 2 O 2 is undoubtedly harmful to the membrane system of organelles in cells (Sies, 2014). According to our data, the activated antioxidant system provided a 'redox-buffer' to detoxify H 2 O 2 so as to modulate increased H 2 O 2 level, and thereby enhance resistance to this oxidative stress. This result suggests the importance of maintaining redox homeostasis during Bgt conidiation. This balance is also required during interaction between Botrytis cinerea and its host to ensure normal growth and virulence (Viefhues et al., 2014).

Source of H 2 O 2
We propose two possible explanations for the source of H 2 O 2 accumulated in the developing conidiophores. One is uptake from host epidermis cells by haustoria, and the other is production of H 2 O 2 by the pathogen itself. Because the H 2 O 2 burst was obviously suppressed in the infected cells of Chancellor (Figure 4B), it is unlikely that most of it was derived from the host. For producing by itself, O 2− can be generated through the metabolic reactions catalyzed by the Nox complex on the plasma membrane and then can be further dismutated to H 2 O 2 by superoxide dismutases (Toledano et al., 2010;Reddi and Culotta, 2013). Another source of O 2− is mitochondria, from which it can be generated via the ETC and released (Suraniti et al., 2014). Note that we found that the ETC was enhanced, while the NoxA gene expression was repressed at the foot cell forming stage. We speculate that more H 2 O 2 is produced by the pathogen itself through mitochondria. This is supported by the findings that enhanced activity of the ETC gives rise to generate more O 2− , while inhibition of Ca 2+ signaling, which plays a positive role in mitochondrial ATP synthesis (Feissner et al., 2009), markedly affects H 2 O 2 accumulation in our investigation.

Involvement of Ca 2+ , H 2 O 2 , and PIP Mediated Signaling Pathways in Conidiogenesis
A crosslink between Ca 2+ and H 2 O 2 mediated signaling pathways is present in Bgt conidiogenesis. One of the most important results of this study is the finding that disturbed Ca 2+ /CaM signaling have a negative influence on both H 2 O 2 production and asexual reproduction of the pathogen. An alternative explanation is that the treatments with the two drugs altered the defense state of the host, including pre-and post-invasion defenses (Lipka et al., 2005), which caused less conidiation. However, 5 mM EGTA and 12.5 µM TFP did not affect the haustorium formation of the pathogen significantly. This means that the preinvasion resistance is not induced. Furthermore, according to the results from DAB staining, very little H 2 O 2 was detected in the infected epidemic cells (Figure 8). This indicates that the post-invasion defense responses during the interaction between the pathogen and the host at 5 and 6 dpi are suppressed. Therefore, the significant reduction in conidiation does not result from the alteration of defense state of the host. These data suggest that the two pathways have crucial roles in regulating conidiation of Bgt and a communication between them is present. The crosslink between Ca 2+ and H 2 O 2 signaling has been reported in mammals (Feissner et al., 2009;Liu et al., 2016), plant (Gonzalez et al., 2012) and fungi, such as Trichoderma atroviride (Medina-Castellanos et al., 2014). The HOG1 with a putative role in ROS signaling showed a significant sequence similarity to the CCaMKs in Bgh, which suggests that it may have dual functions in the crosslinked signaling mediated by H 2 O 2 and Ca 2+ . In addition, since the antioxidant protein coding gene expressions can be induced by CCaMKs in plants (Gonzalez et al., 2012;Ma et al., 2012), those in Bgt are possibly induced by the HOG1 during conidiation. For initiation of Ca 2+ fluxes, it can be activated through Ca 2+ channels or IP3 sensitive channels (Wen et al., 2015). However, no significant differential expression of genes involved in Ca 2+ channels or exchangers was detected. Alternatively, activated PIP signaling can lead to release of IP3, which can induce Ca 2+ fluxes from endoplasmic reticulum to cytosol (Alzayady et al., 2016). Hence, PIP signaling is likely involved in the release of free Ca 2+ . More evidence from monitoring of cellular free Ca 2+ and IP3 levels is needed to confirm this hypothesis. Taken together, our data indicate that at least three signaling pathways mediated by Ca 2+ , H 2 O 2 , and PIP are active in regulating conidiogenesis of B. graminis.

CONCLUSION
Several events in primary metabolism including the activation of metabolism of diverse sugars, glycogen synthesis, energy production, and unsaturated fatty acid oxidation occurred during Bgt conidiation. A crosslink between H 2 O 2 and Ca 2+ signaling and involvement of some other regulators are associated with regulation of Bgt conidiation. Our findings broadened and strengthened our knowledge of conidiogenesis of powdery mildews.