Transcriptome and Metabolite Profiling of the Infection Cycle of Zymoseptoria tritici on Wheat Reveals a Biphasic Interaction with Plant Immunity Involving Differential Pathogen Chromosomal Contributions and a Variation on the Hemibiotrophic Lifestyle Definition  

Abstract

The hemibiotrophic fungus Zymoseptoria tritici causes Septoria tritici blotch disease of wheat (Triticum aestivum). Pathogen reproduction on wheat occurs without cell penetration, suggesting that dynamic and intimate intercellular communication occurs between fungus and plant throughout the disease cycle. We used deep RNA sequencing and metabolomics to investigate the physiology of plant and pathogen throughout an asexual reproductive cycle of Z. tritici on wheat leaves. Over 3,000 pathogen genes, more than 7,000 wheat genes, and more than 300 metabolites were differentially regulated. Intriguingly, individual fungal chromosomes contributed unequally to the overall gene expression changes. Early transcriptional down-regulation of putative host defense genes was detected in inoculated leaves. There was little evidence for fungal nutrient acquisition from the plant throughout symptomless colonization by Z. tritici, which may instead be utilizing lipid and fatty acid stores for growth. However, the fungus then subsequently manipulated specific plant carbohydrates, including fructan metabolites, during the switch to necrotrophic growth and reproduction. This switch coincided with increased expression of jasmonic acid biosynthesis genes and large-scale activation of other plant defense responses. Fungal genes encoding putative secondary metabolite clusters and secreted effector proteins were identified with distinct infection phase-specific expression patterns, although functional analysis suggested that many have overlapping/redundant functions in virulence. The pathogenic lifestyle of Z. tritici on wheat revealed through this study, involving initial defense suppression by a slow-growing extracellular and nutritionally limited pathogen followed by defense (hyper) activation during reproduction, reveals a subtle modification of the conceptual definition of hemibiotrophic plant infection.

Fungi are responsible for many of the most economically damaging diseases of plants that threaten future food security (Chakraborty and Newton, 2011; Fisher et al., 2012). The largest class of plant pathogenic fungi accounting for many of these diseases are the Dothideomycetes, members of which are predicted to utilize diverse and unique mechanisms for host infection leading to pathogen reproduction (Ohm et al., 2012).

For plant pathogens to reproduce effectively in host tissues, each invading species must overcome several layers of plant immunity. After surpassing any physical barriers, pathogens must then overcome the host plant’s molecular surveillance system, which detects microbe- or pathogen-associated molecular patterns (PAMPs), which trigger PAMP-triggered immunity. Successful pathogens suppress PAMP-triggered immunity by secreting protein effectors that interfere with perception or defense signaling. However, plants also can detect the function of these effectors using resistance proteins, which monitor the function of effector-targeted proteins or recognize the effector directly. This gives rise to effector-triggered immunity, which is multifaceted and strong, frequently involving localized plant cell death (Jones and Dangl, 2006). These inducible defenses also are associated with wide-ranging transcriptional and hormonal reprogramming in the plant (Hammond-Kosack and Jones, 1996). This mechanism is particularly effective at preventing the proliferation and reproduction of biotrophic pathogens (which live and feed in association with living cells) in plant tissues.

Fungi from the genus Mycosphaerella, the largest genus of plant pathogenic fungi, belong to the Dothideomycetes class. These fungi typically have long periods of symptomless intercellular (extracellular with respect to host cells) growth in infected plants (frequently described as biotrophy) followed by a rapid switch to necrotrophy associated with host tissue death. Almost all Mycosphaerella spp. fungi reproduce asexually (sporulate) in plant tissues undergoing this cell death to initiate polycyclic diseases. Plant pathogenic Mycosphaerella spp. are thus often referred to as hemibiotrophic due to having long periods of initial colonization without stimulating host cell death, where some form of biotrophic interaction is proposed (Luttrell, 1974). This initial asymptomatic phase is subsequently followed by host cell death associated with rapid fungal growth and subsequent reproduction (sporulation), a feature more commonly associated with successful colonization by necrotrophic pathogens.

Currently, one of the best studied Mycosphaerella spp. fungi is Zymoseptoria tritici (also known as Mycosphaerella graminicola or Septoria tritici), which is exclusively pathogenic toward leaves of wheat (Triticum aestivum), causing Septoria tritici blotch disease. Typical of most Mycosphaerella spp. plant pathogens, Z. tritici enters leaves via their natural openings, stomata, and completes its full asexual reproductive cycle without physically penetrating host cells (Kema et al., 1996; Pnini-Cohen et al., 2000; Deller et al., 2011; Dean et al., 2012). This strictly intercellular (or apoplastic) mode of colonization and reproduction differentiates Z. tritici from almost all other well-studied plant pathogenic fungi with available genomic resources, including all species of Magnaporthe, Fusarium, Ustilago, Colletotrichum, Blumeria, and Puccinia, all of which have at least one period of intracellular (invasive) growth within plant cells during infection (Brown and Hammond-Kosack, 2014). This lifestyle suggests that communication between the Z. tritici hyphae and the neighboring wheat cells must rely heavily on the exchange of secreted or cell surface-localized molecules but raises questions regarding the exact nature of the dialogue throughout infection.

A single asexual reproductive cycle of Z. tritici on wheat leaves lasts at least 2 to 3 weeks from spores alighting on the leaf surface. During the early phase of colonization (lasting 7 d or more depending on the particular isolate-cultivar combination), the fungus grows extremely slowly between the mesophyll cell layers of the leaf. This growth is so slow that it has often proven difficult/impossible to detect significant increases in fungal biomass in infected tissue during this period (Kema et al., 1996; Pnini-Cohen et al., 2000; Keon et al., 2007). It was demonstrated recently that the function of a secreted LysM domain-containing Z. tritici effector protein, Mg3LysM, plays an important role in establishing this initial symptomless colonization through the suppression of chitin (PAMP) recognition by the wheat chitin elicitor binding protein (TaCEBiP) and chitin elicitor receptor kinase1 (TaCERK1; Marshall et al., 2011; Lee et al., 2014). However, at some point after 7 to 10 d post infection (dpi), and for reasons that are currently unclear, wheat cells then suddenly begin to die. This response is strictly localized to the inoculated leaf area and shows many biochemical features of programmed cell death (PCD; Keon et al., 2007; Rudd et al., 2008). The onset of host cell death coincides with a sharp increase in fungal growth rate and biomass in infected leaves. This may be supported by the increased apoplastic nutrient availability as a consequence of the loss of host cell membrane integrity (Keon et al., 2007). Z. tritici then ultimately enters its reproductive mode in the dying leaf tissue, forming masses of asexual spores, produced by pycnidia formed in substomatal cavities. During high humidity and rainfall, these spores are extruded through the stomatal aperture of dead leaf areas and are propagated via rain splash to establish further cycles of infection.

The fully sequenced 38-Mb genome of the reference isolate of Z. tritici, IPO323, comprises approximately 11,000 predicted genes spanning 21 chromosomes (Wittenberg et al., 2009; Goodwin et al., 2011). The eight smallest chromosomes harboring 646 genes have been shown to be conditionally dispensable for plant infection and are now referred to as either dispensable or accessory chromosomes. It is currently unknown what role these chromosomes have in the lifestyle of this fungus, although it has been suggested that they may function to facilitate more rapid evolutionary processes (Croll and McDonald, 2012; Croll et al., 2013). A notable feature of the overall gene composition of Z. tritici is the unusually low numbers of genes encoding putative plant cell wall-attacking proteins, perhaps as a consequence of its strictly extracellular mode of colonization. The Z. tritici genome, however, does encode a large number of predicted secreted proteases and other classes of secreted proteins, including many putative Cys-rich effector candidates (Goodwin et al., 2011; Morais do Amaral et al., 2012). In addition, despite the fungus triggering infection phase-specific host cell death, no toxin (protein or metabolite) has been described thus far.

In contrast to Z. tritici, the extremely large and repeat-rich genome sequence of hexaploid wheat is still being finalized (Eversole et al., 2014). However, large numbers of ESTs, in addition to various individually sequenced chromosomes, are available to support molecular analyses (Wise et al., 2007; Eversole et al., 2014). In addition, a limited number of studies have begun to address the metabolomic responses of wheat toward pathogens and/or pathogen effectors (Du Fall and Solomon, 2013). Many conserved plant responses to pathogens, including the production of the phytohormones jasmonic acid (JA) and salicylic acid (SA), also can be measured in wheat. While various exceptions now exist, SA is often considered to trigger defenses effective against biotrophic pathogens involving localized hypersensitive cell death and the production of PATHOGENESIS-RELATED (PR) proteins, while JA triggers defenses usually more effective against necrotrophic pathogens (Hammond-Kosack and Parker, 2003; Glazebrook, 2005). Other conserved characteristic transcriptional defense responses of plants include the activation of secondary metabolite pathways, such as the phenylpropanoid pathway, which can serve to strengthen plant cell walls via lignification, thus preventing extensive cell wall hydrolysis induced by invading pathogens (Sticher et al., 1997). Monocotyledonous cereal plants also have various unique metabolite pathways that are transcriptionally activated during abiotic stresses. For wheat, one particular example is the biosynthesis of fructans, which can serve as mobile carbon sources and have been implicated in plant tolerance of cold, drought, and oxidative stresses (Hare et al., 1998; Gaudet et al., 1999; Valluru and Van den Ende, 2008; Van den Ende and Valluru, 2009; Keunen et al., 2013). However, to date, no metabolomics analyses have been reported on the interaction of wheat with Z. tritici.

RNA sequencing (RNAseq) technology has become more frequently used to profile fungal pathogen gene expression at particular phases of plant infection (de Jonge et al., 2012; Kleemann et al., 2012; O’Connell et al., 2012; Soanes et al., 2012; Cantu et al., 2013; Garnica et al., 2013; Hacquard et al., 2013; Meinhardt et al., 2014; Zhang et al., 2014). However, comparatively few studies (Kawahara et al., 2012; Yang et al., 2013; Yazawa et al., 2013) have simultaneously analyzed the responses of both pathogen and host, and very few have undertaken such an analysis throughout a complete pathogen reproductive cycle. One of the current most comprehensive studies of this type was performed, using a composite microarray, on the hemibiotrophic fungus-like filamentous oomycete Phytophthora capsici infecting tomato (Solanum lycopersicum; a dicotyledonous plant), where many infection phase-specific transcriptional programs were identified for both pathogen and plant (Jupe et al., 2013). However, unlike the true fungus Z. tritici, P. capsici has a relatively short reproductive cycle on tomato, and the initial biotrophic phase involves the generation of specific haustoria-like structures inside living plant cells. As representative members of two distinct taxonomic kingdoms, it is probable that the suggested hemibiotrophic infection strategies differ for Z. tritici and P. capsici.

The unique temporal features of the Z. tritici-wheat interaction, and the available genomic and EST resources, make this an attractive model to investigate infection phase-specific transcriptional and physiological changes, where clearly defined periods can be relatively easily separated (Keon et al., 2007; Shetty et al., 2007; Kema et al., 2008). Three previous reports have used next-generation RNAseq to study aspects of this interaction (Brunner et al., 2013; Yang et al., 2013; Kellner et al., 2014). Brunner et al. (2013) studied the expression of putative plant cell wall-attacking enzymes at three phases of infection by the isolate ST99CH3D7 of Z. tritici (not the publicly available reference isolate IPO323) and demonstrated that the majority were strongly expressed during the necrotrophic phase. Yang et al. (2013) investigated both pathogen and host at two time points during the initial symptomless interaction (4 and 10 dpi) and at the onset of first symptoms (the transitional phase), which occurred at 13 dpi in these experiments. This study confirmed many previously published fungal transcriptional responses deriving from EST analysis, microarray, and quantitative reverse transcription (qRT)-PCR (Keon et al., 2007; Kema et al., 2008; Rudd et al., 2010) and indicated that host defense responses may become activated at the disease transition. However, this analysis was somewhat limited by a number of features, including (1) the low overall depth of sequencing (less than 13 million reads were generated per sample), leading to (2) a very low coverage of fungal transcripts (less than 10% of the total predicted fungal genes had any RNAseq support), (3) the lack of analysis of fungus growing away from the plant (i.e. in culture), precluding the identification of fungal genes and processes up-regulated specifically during infection, (4) the absence of samples representing very early phases of the interaction, and (5) the lack of analysis of the later phases of advanced necrotrophic colonization and asexual pathogen reproduction. Finally, Kellner et al. (2014) described an analysis of the Z. tritici transcriptome in axenic culture and also at a single time point of plant infection (during symptomless growth at 4 dpi) of wheat and a nonhost plant. A notable observation from this work was that, in general, genes residing on the eight smallest accessory chromosomes were expressed at a relatively low level (although differential expression was reported for some) compared with those on the larger core chromosomes 1 to 13. However, the sampling of only a single in planta time point made it impossible to determine whether this was a feature conserved throughout infection and also prohibited the study of the dynamics of overall pathogen gene expression during the infection time course. There was no comparable analysis of the host transcriptome in that study.

The work presented here aimed to address the question of how the transcriptome and metabolome of Z. tritici and wheat leaves respond throughout clearly defined and distinct phases of infection, culminating in successful pathogen reproduction. The depth of RNAseq we achieved enabled us to report expression changes for more than 80% of all currently predicted fungal genes throughout infection. Analysis of two in vitro fungal culture conditions and corresponding mock-inoculated wheat leaves formed effective baselines for simultaneous transcriptional and metabolite profiling of the dynamics of the host-pathogen interaction. The large and comprehensive data sets we describe indicate that the early biotrophic phase of pathogen colonization does involve plant defense suppression but does not involve effective nutrient acquisition. Rapid fungal proliferation and asexual reproduction subsequently occur in an environment with both the transcriptional and hormonal features of an intense plant defense response. The data overall suggest a model for initial pathogen subterfuge followed by extreme activation of plant immune responses occurring to facilitate reproduction in infected tissues. Moreover, different Z. tritici chromosomes, and chromosomal regions, appear to contribute differentially to this mechanism, which offers an alternative and intriguing view on hemibiotrophy for plant pathogens that never physically penetrate host cells.

RESULTS AND DISCUSSION

Identification and Classification of Differentially Abundant Z. tritici and Wheat Leaf Metabolites

For our combined RNAseq and metabolite experimental pipeline (Supplemental Fig. S1), we used two fungal culture media and five time points of plant infection. We chose the culture medium Czapek-Dox broth (CDB) as a defined medium that presents carbon as Suc and nitrogen as nitrate only, while a second medium, potato dextrose broth (PDB), presents a rich and complex variety of carbon and nitrogen sources. Both media support fungal growth via a form of yeast-like budding, which occurs approximately three to four times faster in PDB than in CDB (Keon et al., 2005). Five time points of leaf infection were selected to span key phases of the host-pathogen interaction (Fig. 1). These included 1 dpi, where fungal spores had just germinated on the leaf surface; 4 dpi, which represents early poststomatal penetration and slow symptomless intercellular growth inside leaves; 9 dpi, representing the transition phase associated with the onset of more rapid fungal intercellular growth and the first appearance of macroscopic disease symptoms on leaves; and 14 dpi, representing accelerating plant cell death with rapid fungal necrotrophic growth and the onset of sporulation. Finally, we sampled tissue at 21 dpi, when leaf material is fully defeated (having a dry necrotic appearance) and mature sporulation structures (pycnidia) are present containing abundant asexual spores.

Figure 1.

Time course of infection of Z. tritici on wheat. Leaf materials were collected and subsequently analyzed by RNAseq, GC-MS, and liquid chromatography-mass spectrometry (LC-MS).

A nontargeted metabolite profile analysis was performed on fungal growth in two culture media and mock-treated and infected wheat leaves at every time point post inoculation, conducted utilizing a combination of three independent platforms: ultrahigh performance liquid chromatography (UHLC)-tandem mass spectrometry (MS/MS) optimized for basic species, UHLC-MS/MS optimized for acidic species, and gas chromatography-mass spectrometry (GC-MS). This led to the detection and quantification of 309 differentially abundant metabolites (Supplemental Table S1). Thirty-eight compounds were identified only as present in fungal cultures and were not detected at all in either mock-treated or infected wheat leaves. Ten of these (over one-quarter) were implicated in fatty acid and/or sterol metabolism (Supplemental Table S2), highlighting that these are prevalent processes occurring in Z. tritici spores. We identified 56 metabolites that were found only in leaves (mock treated or inoculated) and not detected in cultured fungus (Supplemental Table S2). These included many known plant-specific metabolites, including Suc, shikimate, salicylic acid, and fructans, the latter of which function in abiotic stress responses and also as dynamic carbon stores in wheat (Hare et al., 1998; Gaudet et al., 1999; Valluru and Van den Ende, 2008; Van den Ende and Valluru, 2009). All the remaining metabolites were detected in both fungal culture and plant (mock-treated and inoculated) samples. However, there were significant variations in the abundance of members of this group of metabolites between fungal culture and mock-inoculated leaves. For example, six metabolites were detected in at least 100-fold higher levels in fungal spores than in untreated leaves (for mean abundance data, see Supplemental Table S1). This group of metabolites included various sugar alcohols (arabitol and mannitol) and also stachydrine (Pro betaine) and trehalose (Table I).

Metabolites detected in fungal culture at levels more than 100-fold mock-inoculated wheat leaves

Table I.

Metabolites detected in fungal culture at levels more than 100-fold mock-inoculated wheat leaves

Biochemical	Platform	Mean Values, Mock Leaves	Mean Values, Fungal Culture	Ratio Fungal Culture:Mock Leaves
Stachydrine	LC-MS positive	0.027	61.733	2,294.9
Mannitol	GC-MS	0.162	109.166	672.3
Arabitol	GC-MS	0.170	87.219	510.4
Trehalose	GC-MS	0.089	42.160	469.7
Threitol	GC-MS	0.626	185.669	296.5
Valerate	LC-MS negative	0.604	71.809	118.7

Accumulation of Abundant Fungal Metabolites in Infected Leaves and RNAseq Read Mapping to the Fungal Genome Chart the Dynamics of Disease Progression

RNAseq reads from each sample of fungal culture and the time course of leaf infection were then aligned to the fully sequenced Z. tritici genome template and the filtered (see “Materials and Methods”) wheat National Center for Biotechnology Information Unigene Build 60 set. A minimum of 80% of reads originating from fungal cultured material aligned to the fungal genome template (Table II). Conversely, very small percentages of the total number of reads aligned to the fungal genome at both 1 and 4 dpi (less than 4%), again highlighting the very low levels of quantifiable biomass during this period of leaf colonization. By 9 dpi, at the onset of first symptoms, an increased proportion of reads mapped to the fungal genome (less than 7%). Subsequently, by 14 dpi, the number of reads mapping increased significantly (less than 60%), and by 21 dpi, the mapping obtained was comparable to the fungal cultured samples (approximately 80%), indicating that the host RNA had been almost completely destroyed by this point (Keon et al., 2007). Particular subsets of metabolites showed a similar pattern of accumulation throughout infection, closely matching the overall RNAseq read mapping to the fungal genome. These metabolites included those that had been detected only in fungal cultures as well as those that we had determined to be more than 100-fold more abundant in fungal culture than in untreated leaves, including arabitol, mannitol, trehalose, and stachydrine (Fig. 2). These data suggest that quantities of these metabolites might accurately reflect the disease progression of Z. tritici-infected leaves and confirmed that they accumulate to high levels in the asexual spores of this species. These probably represent important metabolites in Z. tritici, as has been demonstrated for other fungal pathogens of wheat (Solomon et al., 2006; Lowe et al., 2009).

RNAseq read counts and percentage mapping statistics to the Z. tritici genome

Table II.

RNAseq read counts and percentage mapping statistics to the Z. tritici genome

Sample	No. of Reads	No. of Reads Mapped to the Z. tritici Genome	Percentage of Reads Mapped to the Z. tritici Genome
CDB-1	32,770,793	28,316,903	86.41
CDB-2	67,034,959	58,476,386	87.23
CDB-3	23,162,180	20,887,165	90.18
PDB-1	38,371,061	31,325,941	81.64
PDB-2	38,516,569	33,946,087	88.13
PDB-3	20,672,941	19,301,735	93.37
Day1-1	35,526,015	453,509	1.28
Day1-2	41,492,156	727,796	1.75
Day1-3	27,254,131	181,633	0.67
Day4-1	54,948,633	269,050	0.49
Day4-2	46,752,435	334,890	0.72
Day4-3	24,856,079	924,376	3.72
Day9-1	47,320,630	1,335,339	2.82
Day9-2	41,382,074	2,122,459	5.13
Day9-3	24,052,623	1,517,352	6.31
Day14-1	49,844,627	8,818,097	17.69
Day14-2	39,071,063	16,476,967	42.17
Day14-3	21,605,374	13,004,622	60.19
Day21-1	33,043,199	26,201,287	79.29
Day21-2	43,409,156	36,939,935	85.10
Day21-3	19,798,698	15,420,296	77.89

Figure 2.

Abundant fungal metabolite accumulation and RNAseq read mapping to the fungal genome chart the progress of disease. Red lines represent Z. tritici-infected leaves. Blue lines represent mock-treated leaves. Asterisks indicate significant increases at adjusted P < 0.05.

Specific Whole Z. tritici Chromosomes, and Regions of Individual Chromosomes, Contribute Differentially to Overall Gene Expression throughout the Course of Infection and Asexual Reproduction

In contrast to a previous study (Yang et al., 2013), the increased depth of our sequencing allowed us to detect significant expression of fungal genes in all samples and time points of infection, despite differing overall read mapping percentages. In summary, of the 11,017 predicted genes in the Z. tritici genome sequence, the following numbers were detected as expressed in each sample: CDB, 10,541; PDB, 10,437; day 1, 9,740; day 4, 9,236; day 9, 10,041; day 14, 10,464; and day 21, 10,485 (Supplemental Table S3). Therefore, we were able to detect the expression of at least 80% of all predicted fungal genes in all samples.

Prior to performing differential gene expression analyses between samples, we first calculated a mean value for relative expression (fragments per kilobase of exon per million fragments mapped [FPKM]) of all genes whose expression was detected and displayed them graphically across the 21 Z. tritici chromosomes using the OmniMapFree software (Antoniw et al., 2011). Figure 3A, top, displays the physical distribution of all 11,017 currently predicted genes across the 21 chromosomes. We then sorted all genes by mean FPKM values (Supplemental Table S3) and subsequently displayed those with very low (mean FPKM < 1.0) or no (FPKM = 0) expression throughout the replicated experiments (Fig. 3A, bottom). In total, 260 currently predicted genes gave FPKM values of 0 in all samples, indicating that their expression level was below the limit of detection throughout this experiment. The majority of these genes lay in close physical proximity on chromosome 7, between nucleotides 1,676,643 and 2,552,085. PCR on genomic DNA spanning three independent intervals confirmed that the region was present in the experimental isolate (Fig. 3A, inset). This region currently contains 139 predicted genes and also was seen to have low overall expression at 4 dpi in an independent study (Kellner et al., 2014). The lack of expression support from this region across all replicated plant infections suggests that high-level expression of these genes is not essential for any phase of colonization or asexual reproduction by Z. tritici.

Figure 3.

The small accessory chromosomes and a region on chromosome 7 have low transcriptional activity throughout infection. A, The top displays the positions of all current predicted Z. tritici genes across the 21 chromosomes of isolate IPO323, and the bottom shows all genes with very low or no (FPKM < 1.0) detected expression throughout infection. The inset displays PCR on genomic DNA, confirming that the low/nontranscribed region on chromosome 7 was present in the experimental isolate IPO323. B, Mean value of expression per gene from each chromosome averaged from the entire data set.

Genome-wide analysis also indicated that many genes present on the eight smallest dispensable or accessory chromosomes, 14 to 21, have low expression relative to those present on the first 13 core chromosomes. The calculated mean FPKM value per gene for each chromosome reinforced this observation, highlighting that a gene present on any of the dispensable chromosomes was expressed, on average, 4- to 5-fold lower than an equivalent gene from the core chromosomes (Fig. 3B). This observation was true for all samples at every time point (Supplemental Table S3), thus confirming that the low expression previously reported for growth in one axenic culture and one time point of infection (Kellner et al., 2014) is in fact seen throughout the entire asexual infection cycle (as sampled here) and in at least two additional culture media. Overall, this suggests that the expression of genes residing on the accessory chromosomes is significantly weaker at all phases of the pathogen’s asexual life cycle in host tissue.

The Most Significant Transcriptional Reprogramming Detected in Z. tritici Occurs in Response to the Plant and Involves Genes Primarily Residing on Core Chromosomes

We next calculated the number of differentially expressed fungal genes between samples. The raw CuffDiff2 differential expression output was subjected to further filtering to increase the stringency of the analysis (for details, see “Materials and Methods “). We analyzed all differential gene expression changes relative to the levels found in the defined culture medium, CDB (Table III). This demonstrated that the least number of differentially expressed genes was between CDB and the other culture medium, PDB (607). In contrast, the highest number was between CDB and 9 dpi (1,389). Of note were over 1,000 fungal genes differentially expressed within 1 d of alighting on the leaf surface, the majority being up-regulated. This indicated that major transcriptional reprogramming takes place soon after contact with the host plant, as was suggested for the biotrophic fungal pathogen of corn (Zea mays), Ustilago maydis (Lanver et al., 2014). This is reinforced by the fact we were able to identify more than 200 fungal genes that demonstrated a peak in expression levels as early as 1 dpi (Table III).

Number of Z. tritici differentially expressed genes relative to levels in CDB culture

Table III.

Number of Z. tritici differentially expressed genes relative to levels in CDB culture

The number of genes showing their peak up-regulation/time point is also indicated.

Sample	PDB	1 dpi	4 dpi	9 dpi	14 dpi	21 dpi
Total	607	1,083	1,183	1,389	1,114	1,331
Up	395	627	771	819	640	719
Peak Up	202	224	342	516	110	487
Down	212	456	412	570	474	612

In total, 3,046 genes exhibited significant differential expression across the experiment, representing 28% of the total genes predicted in the sequenced genome (Supplemental Table S4). The majority of these genes located on the core chromosomes 1 to 13, with 2,955 genes showing differential expression (28% of all genes on chromosomes 1–13). In contrast, only 79 differentially expressed genes were present on the accessory chromosomes, 14 to 21, which represents 12% of the genes residing on these chromosomes. Taken together with the low overall relative expression of genes on chromosomes 14 to 21 (Fig. 3), this would indicate that the fungus invests fewer transcribing sequences from the accessory chromosomes than from the core chromosomes at all phases of colonization and asexual reproduction.

Combined RNAseq and Metabolite Analysis of Z. tritici Growth in Axenic Cultures and in Planta Suggests Limited Hexose and Nitrate Assimilation during Early Leaf Infection

In order to test correlations between changes in metabolite levels and changes in the expression of genes, we performed an initial study focusing solely on in vitro growth by the fungus in CDB and PDB. To reiterate, growth in both broths is morphologically indistinguishable; however, growth occurs approximately four to five times faster in PDB. CDB presents only Suc as the sole carbon source and nitrate as the sole nitrogen source. Materials were collected during log-phase growth from each broth.

The levels of almost all detected amino acids were significantly lower in cells grown in CDB than in PDB (Fig. 4A), with the exception of Glu and Gln. This probably reflects slower fungal growth as a consequence of a single nitrogen source (NO₃  ²⁻). Coupled with this observation, two predicted plasma membrane transporters of the putative nitrate-transporting nrt2 type (Galvan and Fernández, 2001) were strongly transcriptionally up-regulated in CDB relative to PDB. This was paralleled by an up-regulation of a putative Z. tritici nitrate reductase homolog (Fig. 4B). Fungal growth on Suc requires the activity of secreted invertase enzymes to generate extracellular Glc and Fru for uptake. One secreted invertase homolog was strongly up-regulated in CDB versus PDB, as were three putative transmembrane hexose transporters (Fig. 4B). The uptake of large amounts of hexoses derived from Suc gave rise to large accumulations of both Glc-6-P and its isomer Fru-6-P in CDB-grown cells (Fig. 4B). In addition, growth in CDB culture gave rise to increased levels of trehalose and mannitol in spores. A pathway to trehalose synthesis from Glc-6-P also was supported by increased UDP-Glc, suggesting that this sugar may act as a storage compound for carbon in situations where nitrogen is limiting. The accumulation of high levels of Glc-6-P and Fru-6-P in CDB-grown cells potentially occurs due to a limitation on nitrogen to fuel downstream biosynthetic processes, including amino acid biosynthesis. Perhaps most significantly, the putative transporters of nitrate and hexose, which were highly expressed in CDB, were strongly down-regulated during most phases of plant infection (Fig. 4C). This probably indicates that more diverse or alternative nutrient sources are being utilized during infection, suggesting that simple hexose and nitrate uptake are not major extracellular nutrient sources, particularly not during the early symptomless phase of colonization (Fig. 4C).

Figure 4.

Metabolite and transcriptome analysis of fungal culture illuminates metabolism during early plant infection. A, Fungal growth in CDB reduces levels of amino acids. Blue signifies decreased levels in CDB, and yellow indicates increased levels (adjusted P < 0.05). F6P, Fru-6-P; G6P, Glc-6-P; PEP, phosphoenolpyruvate. DHA, Dihydroxyacetone. B, Growth in CDB increases levels of hexose 6-phosphates, trehalose, and mannitol, accompanied by the up-regulation of genes encoding hexose and nitrate transporters. C, Low-level expression of genes involved in nitrate and hexose uptake and assimilation during phases of plant infection. Numbers represent mean FPKM values at each time point of infection.

Transcriptome and Metabolome Analyses of Symptomless Colonization from 1 to 4 dpi

For Z. tritici

In order to begin to understand the dynamics of fungal transcriptional changes during infection, we performed some global hierarchical clustering (K-means clustering) on expression. We separated all the genes with statistically significant differential expression into 12 clusters (clusters 1–12 displayed in Table IV). We then analyzed the composition of each cluster for the total numbers of genes, numbers of secreted proteins, and numbers of different functional classes (represented by InterPro [IPR] annotations).

Analysis of expression clusters of Z. tritici genes (K-clusters) throughout the course of infection

Table IV.

Analysis of expression clusters of Z. tritici genes (K-clusters) throughout the course of infection

K-Cluster Identifier	CDB	PDB		Total No. of Genes in Cluster	Percent Secreted	Percent Glycosylphosphatidylinositol Anchored	Percent with Annotation	Top Three IPR Descriptions (Percent of Genes in the Cluster)
1				308	14	1.3	67	Glycoside hydrolase (3.2)
								Glc/ribitol dehydrogenase (2.9)
								Zinc-containing alcohol dehydrogenase (2.6)
2				187	21	0	72	Glycoside hydrolase (7.5)
								Major facilitator superfamily (6.4)
								General substrate transporter (5.9)
3				367	10	1.4	69	AAA ATPase (2.7)
								Major facilitator superfamily (1.9)
								Glc/ribitol dehydrogenase (1.9)
								Cytochrome P450 (1.9)
4				383	14	2.1	73	Glc/ribitol dehydrogenase (3.6)
								Major facilitator superfamily (2.3)
								Null (1.8)
5				345	9	1.2	74	Glc/ribitol dehydrogenase (4.1)
								Major facilitator superfamily (2.3)
								Cytochrome P450 (1.7)
								Enoyl-CoA hydratase/isomerase (1.7)
								Acyl-CoA dehydrogenase (1.7)
6				126	17	0.8	56	Peptidase (3.2)
								Glycoside hydrolase (3.2)
								Major facilitator superfamily (2.4)
								Glc/ribitol dehydrogenase (2.4)
7				180	19	0.5	76	Chloroperoxidase (4.4)
								Cytochrome P450 (4.4)
								Glc/ribitol dehydrogenase (3.3)
								Major facilitator superfamily (3.3)
8				66	20	0	32	Zinc finger (3)
								Remainder all appear once on the list
9				346	13	2.0	61	Major facilitator superfamily (1.7)
								Zinc finger (1.4)
								Basic-Leu zipper transcription factor (1.4); fungal transcriptional regulatory protein (1.4)
10				221	37	0	57	Peptidase (8.6)
								General substrate transporter (3.6)
								Cytochrome P450 (3.2)
11				321	27	1.9	61	General substrate transporter (4)
								Glycoside hydrolase (3.4)
								Peptidase (3.1)
12				197	14	0.5	40	Zinc finger (2)
								Cytochrome P450 (2)
								General substrate transporter (1.5)
								Glc/ribitol dehydrogenase (1.5)

K-clusters 2, 5, 7, 10, and 11 (Table IV) all contained genes significantly up-regulated by 1 dpi relative to levels in CDB culture. The genes contained in each of these clusters are listed in Supplemental Table S4. With respect to physiological processes up-regulated at 1 dpi, the gene compositions of clusters 5 and 7 were notable. The IPR annotations Acyl-CoA dehydrogenase, AMP-dep synthetase/ligase, and Enoyl CoA hydratase were numerically more prevalent among the total and peak up-regulated genes at 1 dpi. These functional groups are predominantly involved in the β-oxidation of fatty acids and lipids, which we also identified to be abundant in Z. tritici spores (Supplemental Table S2). Two genes encoding proteins with the IPR annotation Secretory lipase also were represented in clusters 5 and 7 (Mycgr3G52069 and Mycgr3G43288). All genes in these clusters subsequently decrease in expression, particularly at 9 and 14 dpi. Inspection of the entire composition of genes in clusters 5 and 7 revealed that fatty acid generation from lipids, and their subsequent transport and β-oxidation, were prevalent processes occurring during early infection at 1 and 4 dpi (Supplemental Table S4). Several predicted secreted cutinases also displayed strong up-regulation during this early period, suggesting some degree of attack on host lipids. Moreover, two key enzymes of the glyoxylate cycle, malate synthase (Mycgr3G70677) and isocitrate lyase (Mycgr3G83726), were present in expression cluster 5 with peak expression at 1 and 4 dpi, indicating that the use of acetate, derived from fatty acid β-oxidation, is providing carbon skeletons for anabolic processes in Z. tritici during this early phase of infection.

An overall scheme, based on expression and K-cluster analysis, for the use of lipids as an initial energy source feeding the glyoxylate cycle and supporting the slow initial colonization of wheat leaves by Z. tritici is displayed in Figure 5. The importance of fatty acid β-oxidation and the glyoxylate pathway has been functionally demonstrated via gene deletion analyses in a number of other plant- and animal-infecting fungi, including Leptosphaeria maculans (Idnurm and Howlett, 2002), Magnaporthe oryzae (Wang et al., 2003, 2007), Phaeosphaeria nodorum (Solomon et al., 2004), and Candida albicans (Lorenz and Fink, 2001). In fact, it has been suggested that leaf surfaces represent a nutritionally poor environment for fungal pathogens (Tucker and Talbot, 2001), indicating that the use of stored energy sources is important to facilitate early differentiation, which ultimately supports colonization. The facts that these transcriptional responses persist for at least 4 d post inoculation, and that during this period fungal growth is slow and biomass is low, suggest that Z. tritici is not able to significantly reprogram host physiology to gain additional nutrients during this period, suggesting that, in nutritional terms at least, it is a relatively poor biotroph.

Figure 5.

Early leaf infection (1–4 dpi) induces the expression of Z. tritici genes involved in the β-oxidation of lipids and fatty acids and up-regulation of the glyoxylate pathway. Red coloring indicates up-regulation, and no shading represents no significant change relative to the expression levels in CDB culture (adjusted P < 0.05). FA, Fatty acids; PM, plasma membrane.

We have previously characterized, in silico, the predicted protein secretome of Z. tritici (Morais do Amaral et al., 2012). From a total of 970 genes encoding predicted secreted proteins, 513 (more than 50%) showed differential expression in at least one time point or condition within this experiment (Supplemental Table S5). Three hundred sixty-six of these showed increased expression during at least one time point of leaf infection. The largest number of up-regulated genes of this type was detected at 9 dpi (235 genes), which fall within K-clusters 10 and 11 (Table IV), with the least detected at 21 dpi (122). However, 155 genes encoding putative secreted proteins were up-regulated as early as 1 dpi (Supplemental Table S5), again suggesting an early up-regulated protein secretion response on the host. The number of genes encoding predicted secreted proteins showing increased expression across multiple time points of infection is shown in Figure 6. Only 40 of the 366 were up-regulated at every time point of infection (from 1 to 21 dpi; Supplemental Table S5), indicating that the majority are produced at distinct phases of the interaction. The largest overlap in the identity of up-regulated genes occurred between 9 and 14 dpi (152 of the 366), followed by 4 and 9 dpi (140) and then 1 and 4 dpi (124). Also of note are 49 genes up-regulated at 1, 4, and 9 dpi, suggesting an early and sustained up-regulation of gene expression throughout early infection for this group (Fig. 6).

Figure 6.

Expression profiles of 366 Z. tritici genes encoding putative secreted proteins up-regulated during leaf infection.

Recent genome and secretome comparative analysis involving more than 80 fungal species identified a family of secreted chloroperoxidases (possessing PFAM PF01328), which is expanded in number, in all analyzed Mycosphaerella spp. fungi (Morais do Amaral et al., 2012). The predicted secretome of Z. tritici contains 19 members of this family. These enzymes have been suggested to be involved in the synthesis of halogenated compounds, including, for example, particular antibiotics (Hofrichter and Ullrich, 2006). In addition, chloroperoxidases can also potentially act upon harmful reactive oxygen species (hydrogen peroxide) that may be generated during plant resistance responses to pathogens (Wojtaszek, 1997; Apel and Hirt, 2004). Nine of the putative secreted fungal chloroperoxidases exhibited differential expression in this experiment. Eight of these fell into K-cluster 7 (Table IV; Supplemental Table S5) and displayed up-regulation by 1 dpi with peak up-regulation at 4 dpi, during the postpenetration symptomless phase of leaf colonization. This indicates that this expanded family of genes may have functions important for the initial symptomless intercellular growth phase of Z. tritici.

For Wheat

The most notable aspect of the wheat response to early intercellular colonization by Z. tritici was the down-regulation of specific gene classes characteristically associated with plant defense responses. At 1 dpi, while 101 unigenes out of 334 that were significantly down-regulated had no functional annotation, of the remaining 233, more than 60% had predicted roles in plant defense responses (Table V; Supplemental Table S6). It is likely that this early down-regulation of defense-associated processes is a consequence of the initial fungus-plant recognition and communication. For example, at 1 dpi, down-regulation of a number of wheat genes encoding different pathogenesis-related (PR) protein classes, including PR1, glucan endo-1,3-β-glucosidases (PR2), basic chitinases (PR3), chitin-binding wheat win proteins (PR4), and thaumatin/osmotin-like proteins (PR5), was observed (Supplemental Table S6).

Numbers of differentially expressed wheat genes in Z. tritici-infected leaves versus mock-inoculated leaves at 1 to 14 dpi

Table V.

Numbers of differentially expressed wheat genes in Z. tritici-infected leaves versus mock-inoculated leaves at 1 to 14 dpi

Data for 21 dpi were not analyzed due to the lack of wheat RNA in the heavily diseased necrotic leaves.

Sample	1 dpi	4 dpi	9 dpi	14 dpi
Up-regulated	352	36	2,401	3,056
Down-regulated	334	51	181	3,552
Total	686	87	2,582	6,608

In addition to the typical PR genes, we also observed down-regulation of other gene classes that have been linked with roles in activating plant defenses toward pathogens. This includes genes encoding WRKY transcription factors, subtilases (predicted plant equivalents of animal PCD-inducing caspases), and signaling receptor kinases of different classes, such as Leu-rich repeat receptor-like kinases, wall-associated Ser/Thr kinases, Cys-rich receptor-like kinases, and lectin receptor kinases (Supplemental Table S6).

Oxylipins (e.g. 12-oxo-phytodienoic acid [OPDA] and JA) are signaling molecules derived from α-linolenic acid that play multiple functions in plants, including defense against pests and pathogens (Dong, 1998; Kunkel and Brooks, 2002; Bari and Jones, 2009). This includes C18 members (e.g. OPDA) as well as the C12 molecules JA and methyl jasmonate. Peroxisomal OPDA reductase is a critical checkpoint in the octadecanoid biosynthesis pathway, as it controls metabolite flow from the C18 to the C12 group of compounds (Fig. 7A). It is represented by a small gene family (n = 3–10) in every plant species analyzed to date. Five wheat unigenes annotated as OPDA reductase were found to be significantly down-regulated at 1 dpi (Fig. 7B; Supplemental Table S6).

Figure 7.

Fungus-induced changes in expression of the wheat jasmonate biosynthesis pathway. A, Schematic pathway of JA biosynthesis in plants. B, Relative expression changes of JA pathway components in wheat leaves from 1 to 14 dpi with Z. tritici. Green shading indicates significant down-regulation (adjusted P < 0.05), while red shading indicates up-regulation (adjusted P < 0.05). DAD1, DEFECTIVE IN ANTHER DEHISCENCE1; LOX, lipoxygenase; 13-HPOT, 13(S)-hydroperoxylinolenic acid; AOS, allene oxide synthase; 12,13-EOT, 12,13(S)-epoxylinolenic acid; AOC, allene oxide cyclase; OPDA, 12-oxo-cis-10,15-phytodienoic acid; OPR3, 12-oxophytodienoate reductase; OPC:8, 3-oxo-2-(cis-2′-pentenyl)-cyclopentane-1-octanoic acid; OPCL, 3-oxo-2-(cis-2′-pentenyl)-cyclopentane-1-octanoic acid CoA Ligase; ACX, acyl-CoA oxidase; MFP, multifunctional protein; KAT, 3-ketoacyl-CoA thiolase; CTS, COMATOSE; ACH, acyl-thioesterase.

Lignin synthesis is an important plant defense mechanism against pathogens often induced at the site of pathogen attack. This response is most commonly effective during interactions with pathogens that seek to breach the plant cell wall. Cinnamyl alcohol dehydrogenase (CAD) is a key enzyme in lignin biosynthesis, catalyzing the final step in the synthesis of monolignols (Fig. 8A). CAD deficiency in grasses is known to result in alterations in lignin structure and its overall content (Saathoff et al., 2011; Fornalé et al., 2012). Interestingly, three wheat unigenes, Ta.24122, Ta.62207, and Ta.70656, annotated as CAD, were down-regulated at 1 dpi in addition to one unigene (Ta.9712) corresponding to caffeic acid 3-O-methyltransferase, another important lignin biosynthetic enzyme (Fig. 8B; Supplemental Table S7).

Figure 8.

Fungus-induced changes in expression of the wheat lignin biosynthesis pathway. A, Schematic of pathways leading to lignin formation in plants. Metabolite analysis of Z. tritici-infected wheat leaves identified increased levels of Phe and ferulate during the late stages of infection. B, Relative expression changes of lignin biosynthetic pathway components in wheat leaves from 1 to 14 dpi with Z. tritici. Green shading indicates significant down-regulation (adjusted P < 0.05), while red shading indicates up-regulation (adjusted P < 0.05). PAL, Phenylalanine ammonia lyase; C4H, cinnamate 4-hydroxylase; C3H, p-coumarate 3-hydroxylase; COMT, caffeic acid O-methyl transferase; F5H, ferulate 5-hydroxylase; 4CL, 4-coumarate:CoA ligase; CCR, cinnamoyl-CoA reductase; CCoA-3H, cinnamyl alcohol dehydrogenase; CCoA-OMT, caffeoyl-CoA O-methyltransferase; CAD, cinnamyl-alcohol dehydrogenase.

Fluctuations in the intracellular calcium ion concentration represent a key second messenger pathway during the plant immune response to pathogens (Nürnberger and Scheel, 2001; Rudd and Franklin-Tong, 2001), and calmodulin- and calcium-transporting ATPase proteins are important sensors and mediators of Ca²⁺-dependent signals. In total, we also detected 11 unigenes belonging to this class that were significantly down-regulated at 1 dpi (Supplemental Table S6).

Finally, at 1 dpi, we also found coordinated down-regulation of genes encoding a multitude of proteins involved in the potential detoxification of xenobiotics. This included several members of the glutathione transferase family (n = 19), cytochrome P450s (CYP450; n = 13), ATP-binding cassette (ABC) transporter/multidrug resistance-associated protein family (n = 9), and UDP glucosyl/glucoronyl transferases (n = 5; Supplemental Table S6).

Notably, all the genes discussed above were subsequently up-regulated at 9 dpi (or 9 and 14 dpi), when Z. tritici transitions from asymptomatic to necrotic growth (described below). By contrast, at 4 dpi, only a very small number of wheat genes (n = 87) showed significant change in expression compared with the corresponding mock-inoculated control plants (Table V), suggesting that the plant cells have only limited perception of, and response to, the fungus. One possible explanation for this may be that, by 4 dpi, Z. tritici has effectively deployed the range of effector proteins that were transcriptionally activated as early as 1 dpi (Supplemental Table S5). It is conceivable that these may now be suppressing further recognition by the plant surveillance machinery and/or blocking defense signaling. In support of this, the functionally characterized chitin-binding effector protein Mg3LysM (Marshall et al., 2011) was already significantly up-regulated at both 1 and 4 dpi (Supplemental Tables S3 and S4).

Overall Summary for the Symptomless Phase

Our experimental approach using fungal culture as an expression baseline identified over 1,000 fungal genes that were already differentially expressed within 24 h of spores alighting on the plant surface, suggesting rapid adaptation to this environment. The most notable biological processes up-regulated at 1 dpi and also still at 4 dpi in Z. tritici involve the metabolism of lipids and fatty acids, suggesting that this, as opposed to the uptake of external carbon sources, may support the early phase of colonization. For Z. tritici to still require these pathways as long as 4 dpi suggests that it remains in a nutrient-limiting environment within the wheat leaf, which may be (at least in part) responsible for the limited increase in biomass seen during symptomless colonization. Fungus-secreted proteins are also prevalent among the up-regulated genes during this period, including previously characterized LysM effectors that function to suppress PAMP-triggered immunity (Marshall et al., 2011; Lee et al., 2014). Consequently, at 1 dpi, plant transcriptome analysis identified many defense-associated processes to be suppressed, possibly as a result of the functions of putative fungal-secreted protein effectors. This suppression involved PR proteins, putative JA and Ca²⁺ signaling responses, and cell wall strengthening. The enlarged gene family of secreted Z. tritici chloroperoxidases was rapidly up-regulated by the fungus during this period, suggesting that it had also minimized the risk of exposure to oxidative stress in the plant environment. The overall scheme suggests a form of pathogen subterfuge with minimal utilization of host-derived nutrients during this initial phase of colonization.

Transcriptome and Metabolome Analysis of the Transition to Disease Symptoms at 9 dpi

For Z. tritici

Nine days after inoculation, the first symptoms of disease became visible (Fig. 1), and the overall RNAseq read mapping to the fungal genomic template began to increase (Fig. 2; Table II). The fungus was no longer up-regulating genes implicated in lipid metabolism and catabolism (Fig. 5; Supplemental Table S4), presumably as a consequence of the increased availability of extracellular (plant-derived) nutrients.

The genome sequencing of Z. tritici isolate IPO323 previously predicted that it contains a large number of proteases (Goodwin et al., 2011), suggesting that proteolysis may be providing some essential nutrients at some point during infection. K-clusters 10 and 11 (Table IV) contain genes with peak expression at 9 dpi. The IPR classification protease is the most prevalent annotation in cluster 10, and the third most prevalent in cluster 11. In total, 60 genes with the IPR annotation Protease were present in the entire differentially expressed data set (Supplemental Table S4). Twenty-nine of these fell into K-clusters 10 and 11. Notably, 17 of these are predicted to be secreted (Supplemental Table S5), suggesting that, at 9 dpi, at the transition to disease symptoms, the fungus was producing its maximal secreted protease response. This may be in order to either attack host cells directly or to begin the hydrolysis of host proteins being released from dying cells.

A very high level of expression of many secreted proteins was a general and notable feature of the transition phase of infection. This was evident by simply sorting genes by transcript abundance from highest to lowest FPKM (Supplemental Table S3) and then counting the number of predicted secreted proteins in the top 50 most abundant transcripts. Figure 9 highlights that, at 9 dpi, almost 50% of all transcripts in this top abundance category encoded secreted proteins, of which 35% have no functional annotation. The two functionally characterized LysM domain effector proteins Mg3LysM and Mg1LysM lie in 12th and 27th place in the top 50 genes sorted by transcript abundance at 9 dpi (Supplemental Tables S3 and S7). Of the functionally unknown secreted proteins in this list are 13 genes encoding proteins with more than 5% Cys, indicating that they are Cys rich, a hallmark of apoplastic fungal effector proteins (De Wit et al., 2009). In total, 68 genes encoding functionally unannotated Cys-rich effector proteins were differentially expressed in this study (Fig. 9B; Supplemental Table S8), with the majority falling into K-clusters 10 and 11 shown in Table IV.

Figure 9.

Expression profiles of Z. tritici-secreted protein subclasses throughout plant infection. A, Percentage of secreted and unannotated proteins present in the top 50 most abundant transcripts at each time point. B, Expression profile of 68 unannotated small Cys-rich secreted proteins. C, Expression profile of 31 differentially expressed secreted proteases.

Five genes encoding candidate Cys-rich effectors were selected for functional analysis through generating gene deletion/disruption strains for each and testing for altered virulence on susceptible wheat leaves. All five genes were strongly up-regulated at 9 dpi with high overall relative expression (highlighted by the light gray shading in Supplemental Table S8). The genes selected all encoded mature proteins of between 47 and 127 amino acid residues containing between six and 10 Cys residues. Supplemental Figure S2A displays the amino acid sequence of each protein along with the gene deletion strategy and demonstrates that many independent deletion strains were obtained for each gene targeted. However, none of the mutant strains showed any detectable defect in virulence toward the susceptible wheat cultivar (Supplemental Fig. S2B), indicating that a high level of genetic redundancy between these in planta-expressed Z. tritici secreted protein effectors may exist. A similar problem in reverse genetics-based functional studies of effectors has been encountered in other plant pathosystems; for example, in M. oryzae, 78 different secreted fungal effector gene deletion stains were generated, with only one (MC69) impacting upon virulence (Saitoh et al., 2012). That this level of redundancy may also occur in Z. tritici is perhaps emphasized by the fact that 115 genes encoding functionally unknown secreted proteins were shown to be up-regulated, often displaying peak up-regulation at 9 dpi (Supplemental Table S5).

We then further analyzed the fungal genes present within K-clusters 1 to 12 of Table IV to identify whether any were physically linked and thus potentially coregulated gene clusters (or loci). We set a minimum cutoff of three adjacent genes falling into a single K-cluster as evidence for coregulation at a locus. This identified 16 putative regions (Fig. 10; Table VI). In contrast to observations made for the plant pathogenic fungus U. maydis (Kämper et al., 2006), none of the putative physical clusters were predominantly composed of small secreted proteins. Instead, many of the clustered genes had functions (IPR) related to metabolic processes, and in particular secondary metabolism. For example, the largest putative coexpressed physical cluster of genes (locus 11 in Table VI) was surrounding a polyketide synthase, PKS7 on chromosome 7, with a total of eight physically adjacent genes all falling into K-cluster 10 (Table IV; Supplemental Table S3). This concerted up-regulation of all genes in the PKS7 cluster, and their peak expression at 9 dpi, were confirmed by independent qRT-PCR analysis (Supplemental Fig. S3). Interestingly, the PKS7 cluster includes the ABC transporter ATR3, which has been shown to be dispensable for fungal virulence on wheat leaves (Stergiopoulos et al., 2003). Despite this, the PKS7 cluster represents a strong candidate region for the production of an as yet unknown secondary metabolite generated during the period of disease symptom formation. The production of other currently unknown secondary metabolites at 9 dpi is also supported by the identification of two further, albeit smaller, physical clusters/loci of coregulated genes. This included a three-gene cluster containing a hybrid PKS-nonribosomal peptide synthetase (HPS1) on chromosome 10 and a three-gene cluster surrounding an isopentenyltransferase-like protein on chromosome 2 (physical clusters/loci 15 and 12 in Table VI and Fig. 10, respectively).

Figure 10.

Physical positions of putative coexpressed Z. tritici gene clusters. A, All currently predicted Z. tritici genes. B, Sixteen putative coexpressed physical gene clusters shown in Table VI.

Physical clusters of coexpressed Z. tritici genes

Table VI.

Physical clusters of coexpressed Z. tritici genes

N/A, Not available.

Physical Cluster/Locus	Expression Pattern	Ensembl Identifier	Chromosome Location	IPR Description	Secreted?
1		Mycgr3G103423	2:2981357–2982117	N/A	No
		Mycgr3G91101	2:2982419–2983427	N/A	No
		Mycgr3G68955	2:2984329–2985314	Glc/ribitol dehydrogenase	No
2		Mycgr3G81110	6:840277–841035	N/A	No
		Mycgr3G72903	6:842527–843565	Phospholipid/glycerol acyltransferase	No
		Mycgr3G59769	6:843950–844804	N/A	No
3		Mycgr3G97784	17:103356–104467	N/A	No
		Mycgr3G97785	17:104930–105694	N/A	No
		Mycgr3G97786	17:106684–107683	N/A	No
4		Mycgr3G37732	3:1230247–1231582	Endonuclease/exonuclease/phosphatase	No
		Mycgr3G108533	3:1231827–1233151	N/A	No
		Mycgr3G108534	3:1233680–1234574	N/A	No
5		ERG8	9:1739391–1741142	Galacto-homoserine mevalonate and phosphomevalonate kinase	No
		Mycgr3G101297	9:1741654–1745867	Pentatricopeptide repeat	No
		Mycgr3G75855	9:1746497–1747319	N/A	No
6		Mycgr3G37216	2:1767213–1769408	Alkaline phosphatase	Yes
		Mycgr3G99122	2:1769493–1771089	β-Lactamase-like	No
		Mycgr3G68371	2:1771599–1772547	Peptidoglycan-binding LysM	No
7		Mycgr3G31366	1:5296665–5298411	General substrate transporter	No
		Mycgr3G107684	1:5300985–5303217	Esterase/lipase/thioesterase	No
		Mycgr3G67335	1:5303477–5304917	Aromatic ring hydroxylase	No
8		Mycgr3G110348	8:436741–438158	S1/P1 nuclease	Yes
		Mycgr3G94962	8:438836–439859	N/A	No
		Mycgr3G110349	8:440581–442729	N/A	Yes
9		Mycgr3G91263	2:3674776–3675478	N/A	No
		Mycgr3G103509	2:3675857–3676433	N/A	No
		Mycgr3G36641	2:3677612–3678173	N/A	No
10		CYP-49	7:869030–870993	Cytochrome P450	No
		Mycgr3G28893	7:872083–872413	N/A	No
		Mycgr3G105151	7:874465–875117	N/A	No
11		Mycgr3G45489	7:2577095–2579147	Amino acid permease-associated region	No
		PKS7	7:2581919–2589469	Acyl transferase region	No
		Mycgr3G45256	7:2590016–2591228	Taurine catabolism dioxygenase TauD/TfdA	No
		Mycgr3G45561	7:2591336–2592281	Glc/ribitol dehydrogenase	No
		Mycgr3G45185	7:2592692–2594087	α/β-Hydrolase fold	No
		MGATR3	7:2597607–2601809	AAA ATPase	No
		CYP-10	7:2603447–2605040	Cytochrome P450	No
		Mycgr3G44943	7:2605772–2606813	Glc/ribitol dehydrogenase	No
12		Mycgr3G67803	2:532615–533741	Conserved hypothetical protein 730	No
		Mycgr3G107854	2:535734–536815	Isopentenyl transferase-like	No
		Mycgr3G67804	2:537679–540118	FAD-dependent pyridine nucleotide-disulfide oxidoreductase	No
13		Mycgr3G38371	3:1729334–1732135	Peptidase M3A and M3B, thimet/oligopeptidase F	Yes
		Mycgr3G38460	3:1733828–1734746	Glc/ribitol dehydrogenase	No
		Mycgr3G69836	3:1734977–1737927	Fungal transcriptional regulatory protein, N terminal	No
		Mycgr3G91716	3:1739597–1743676	N/A	No
14		Mycgr3G92614	4:1689339–1690777	N/A	No
		Mycgr3G40021	4:1691331–1691995	N/A	No
		MGSCP2	4:1692025–1693749	Esterase/lipase/thioesterase	No
		Mycgr3G22352	4:1694200–1695604	Major facilitator superfamily	No
		Mycgr3G17135	4:1696669–1697425	S-adenosyl methionine (and some other nucleotide)-binding motif	No
15		HPS1	10:448502–457362	Acyl transferase region	No
		Mycgr3G87728	10:459258–460865	General substrate transporter	No
		Mycgr3G105888	10:461405–463044	Isoflavone reductase	No
16		MGLAC5	12:157367–159218	Glycoside hydrolase, family 35	No
		Mycgr3G50352	12:159957–161713	General substrate transporter	No
		Mycgr3G50311	12:161783–162926	N/A	Yes
		Mycgr3G101694	12:163489–164961	Sulfatase	No

Overall, these data suggest that Z. tritici is producing and perhaps secreting a diverse set of secondary metabolites at the onset of disease symptoms. We did not detect or resolve any candidate molecules in our metabolite analysis, possibly due to dilution, low quantities, and/or unique structures. Thus, the role in infection of these putative metabolites is currently unknown; however, functional studies on MgATR3 may indicate that, similar to what we have shown for candidate protein effectors (Supplemental Fig. S2), they also have either redundant or perhaps cultivar-specific functions.

For Wheat

At 9 dpi, the earlier suppression of plant defenses appears to be relieved and significant up-regulation of 2,401 wheat transcripts was observed (Table V; Supplemental Table S6). Over one-third of these up-regulated genes currently have no functional annotation (Supplemental Table S6). Similar to the analysis of Z. tritici transcripts (Supplemental Table S8), the relative abundance of wheat transcripts was sorted from highest to lowest for all samples and the 100 most highly expressed wheat genes were analyzed (top 100). Four bona fide PR genes and one of the nepenthesin-like aspartic proteinase genes (genes known to play roles in defense against various pathogens) were among the list of 100 most abundant wheat genes detected at 9 dpi (Supplemental Table S9). Interestingly, a lipoxygenase2 gene, a candidate for the first component in the biosynthetic pathway of JA (Fig. 7), is also in the top 100 list at 9 dpi and is ranked 67th (Supplemental Table S9). So, it appears that at 9 dpi, the infected plants start investing significant resources into the synthesis of defense-related compounds and proteins.

At 9 dpi, we notably detected the up-regulation of all three key enzymes of fructan, Fru-based polymer (fructan) biosynthesis, namely 1-SST (for Suc:Suc 1-fructosyltransferase), 6-SFT (for Suc:fructan 6-fructosyltransferase), and 1-FFT (for fructan:fructan 1-fructosyltransferase; Fig. 11; Supplemental Table S6). Cereal fructans are produced in the vacuoles and are then known to be transported to the apoplast. One of the important functions of fructans in the extracellular space is the stabilization of cellular membranes during various abiotic stresses such as drought and cold (Hare et al., 1998; Gaudet et al., 1999; Valluru and Van den Ende, 2008; Van den Ende and Valluru, 2009; Keunen et al., 2013). Fructans also may serve as principal short-term energy stores in plant leaves, bridging the temporal gaps between resource availability and demands. The transcriptional data suggest that fructan biosynthesis is stimulated during Z. tritici infection, because the genes coding for 1-SST, 6-SFT, and 1-FFT were all found to be very highly up-regulated at 9 dpi (Fig. 11; Supplemental Table S6). In fact, these three biosynthetic genes reside in the top 5% of the most highly up-regulated genes detected at 9 dpi (Supplemental Table S6). Metabolite analysis also began to detect statistically significant increases in fructan levels in infected leaves at 9 dpi (Fig. 11; Supplemental Table S1).

Figure 11.

The fructan biosynthesis pathway and fructan metabolite accumulation are triggered in wheat leaves during the onset of disease symptoms induced by Z. tritici. A, Quantification of fructan metabolites during the course of infection. **, Significant at adjusted P < 0.05. B, Transcriptional up-regulation of the fructan biosynthetic pathway. Green shading indicates significant down-regulation (adjusted P < 0.05), while red shading indicates up-regulation (adjusted P < 0.05).

Up-regulation of other groups of transcripts whose annotation suggests a role in the response to pathogens and abiotic stresses also was observed (Supplemental Table S6). These include transcripts encoding other PR proteins, cytoplasmic nucleotide-binding Leu-rich repeat-like resistance proteins, receptor-like kinases (Leu-rich repeat-, lectin-, Cys-rich, Leu-rich, and wall-associated kinases), cytoplasmic kinases (e.g. mitogen-activated protein kinase3 [MPK3], which has previously been demonstrated to accumulate and be activated at the protein level in response to Z. tritici; Rudd et al., 2008) and phosphatases, potential cell death-inducing subtilases, other proteases (AAA-type, Ser-, Cys-, Asp-, and metallo-proteases), plant ubiquitination machinery proteins (ubiquitin, E2 and E3 ubiquitin ligases, and Skp, Cullin, F-box containing complex/F-box proteins), reactive oxygen species-generating peroxidases, and WRKY-type transcription factors. Genes encoding plastocyanin-like proteins, which are known to be associated with PCD, plant defense responses, and pathogen infection, also were found to be up-regulated at 9 dpi.

In addition, a large number of genes involved in the biosynthesis of JA and ethylene, hormones known for their role in defense signaling, in particular to necrotrophic pathogens, genes involved in the biosynthesis of lignin and other phenylpropanoid, isoprenoid, and flavonoid compounds, and numerous CYP450 genes potentially involved in the synthesis of a variety of other secondary (defensive) metabolites were significantly up-regulated during this infection time point (Figs. 7 and 8; Supplemental Table S6). Quite noticeable also was the up-regulation of many potential detoxification enzymes, including ABC transporters and UDP-glucosyl/glucoronyl transferases, and enzymes involved in cell wall remodeling, such as callose biosynthesis genes, xyloglucan endotransglucosylases, and fasciclin-like arabinogalactan proteins (Supplemental Table S6).

Interestingly, and in contrast to the earlier stage of the interaction (Table V), only a relatively small number of transcripts were found to be down-regulated at this stage (n = 181), and approximately 50% of these are unknowns (Supplemental Table S6). The most noticeable was the down-regulation of photosynthesis-related genes, genes involved in the biosynthesis of tetrapyrroles (active cores of chlorophyll), light signaling genes, protein biosynthesis genes, and genes encoding pentatricopeptide repeat proteins. Pentatricopeptide repeats in the model plant Arabidopsis (Arabidopsis thaliana) are known to play essential roles in the biogenesis of mitochondria and chloroplasts (Lurin et al., 2004).

Overall Summary for the Disease Transition Phase

The onset of macroscopic cell death in the leaves had a profound effect upon the transcriptomes of both the pathogen and the plant. The plant response comprised the up-regulation of many defense pathways, including defense-associated hormone biosynthetic processes, as well as the characteristic agents of a plant counterattack, secreted PR-type proteins. The fungus was also putting maximum effort into secretion during this phase. The lack of an altered infection phenotype for the five secreted candidate protein effectors, tested by generating single gene deletion strains, suggests that these may function in overlapping processes during the interaction. However, protease secretion in particular might be compensating for the cessation of intracellular lipid degradation to provide energy to Z. tritici. Many candidate effector proteins had maximal expression at the transition time point as well as various putative physical and transcriptional coregulated secondary metabolite biosynthetic clusters. We can envisage at least two potential roles for the large-scale production of effectors and secondary metabolites at this stage: (1) they might be required to defend the fungus against the individual components of the plant counterattack through inhibitory interactions; or (2) they might be actively involved in stimulating the activation of plant defenses culminating in host cell death to facilitate necrotrophic colonization and reproduction (i.e. they may be putative toxins). Finally, at this point, infected wheat leaves have up-regulated the entire fructan biosynthetic pathway either to attempt to limit the degree of cell death or to remobilize carbohydrate stores away from the infection site. Concurrent to all the activated processes, the plant also now diverts energy away from any further photosynthesis in the infected tissue.

Transcriptome and Metabolome Analysis of Rapid Necrotrophic Colonization and Asexual Reproduction by Z. tritici at 14 to 21 dpi

For Z. tritici

At 14 dpi, disease symptoms are clearly visible as the plants extensively undergo a form of widespread PCD with characteristics similar to hypersensitive cell death (Keon et al., 2007; Rudd et al., 2008; Dean et al., 2012). This process is complete at 21 dpi, when the fungus is now asexually sporulating in defeated plant tissue (Fig. 1).

Of the 115 predicted secreted protein effectors (functionally unannotated proteins) that were strongly up-regulated at 9 dpi, 37 were no longer significantly up-regulated at 14 dpi. Of the remaining 78, a reduced, but still significant, level of up-regulation (versus CDB) was observed for 63 of them, with the remainder (15 genes) showing higher expression at 14 dpi than at 9 dpi (Supplemental Table S5). By 14 dpi, only 13 genes encoding predicted secreted proteins were present in the top 50 abundant transcripts, with the remaining time points/materials having nine or less (Supplemental Table S3). By 21 dpi, 92 of the original 115 candidate effectors with peak expression at 9 dpi were no longer up-regulated relative to fungal culture, while 17 showed markedly decreased (but still significant) levels of up-regulation relative to fungal culture (Supplemental Table S5). Only six genes encoding small secreted proteins showed increased expression at 21 dpi, the majority building upon increased levels already detected at 14 dpi.

Overall, this suggests a progressive switching off of the production of many previously up-regulated unannotated candidate effector proteins in the lead up to the complete defeat of plant tissue. This pattern is supported by the overall expression profiles of the 68 Cys-rich effector candidates within this set (Fig. 9B), many of which had peak expression at 9 dpi that was decreased by 14 dpi before being almost completely down to culture levels at 21 dpi. This trend is also mirrored by a comparable pattern of progressive down-regulation of secreted proteases, suggesting a switch away from the digestion of host proteins for nutrition at 14 and 21 dpi (Fig. 9C).

In contrast to the situation for proteases, the contents of K-clusters 1 and 2 (Table IV) highlight that a significant number of glycoside hydrolases, transporters, and alcohol dehydrogenases have increased expression at 14 dpi, often with a subsequent peak at 21 dpi. Combined with our previous in silico protein secretome analysis, we detected 20 putative secreted plant cell wall-attacking enzymes that were differentially expressed in our experiment (Table VII). Fourteen of these genes were present in K-clusters 1 and 2 (Table IV), suggesting that they are deployed late in the infection cycle and are significantly associated with necrotrophic colonization. This agrees well with previous data describing the expression of various plant cell wall-attacking enzymes during wheat infection by a different Z. tritici isolate (Brunner et al., 2013). Overall, the genome sequence of Z. tritici contains relatively few genes encoding enzymes that may have the ability to attack plant cell walls to liberate simple sugars (Goodwin et al., 2011; Morais do Amaral et al., 2012). Nevertheless, the expression data here indicate a switch to putative complex carbohydrates as a more significant energy source very late in infection, potentially involving proteins whose substrate preferences are likely to include cellulose, hemicellulose, and pectin, which are key components of plant cell walls (Table VII).

Expression profiles of differentially expressed Z. tritici-secreted plant cell wall-attacking enzymes

Table VII.

Expression profiles of differentially expressed Z. tritici-secreted plant cell wall-attacking enzymes

Joint Genome Institute Identifier	Ensembl Identifier	Expression Pattern	Chromosome	bp	Description	Protein Name	CAZy Family	Putative Substrate
105871	MGEND1		10	310,310 to 311,338	Glycoside hydrolase family 12	β-1,4-Endoglucanase	GH12	Cellulose
103512	Mycgr3G103512		2	3,707,056 to 3,708,585	Glycoside hydrolase family 61		GH61	Cellulose
100252	Mycgr3G100252		5	1,526,256 to 1,527,825	Glycoside hydrolase family 7	Cellobiohydrolase	GH7	Cellulose
76589	MGEND2		10	1,542,406 to 1,543,537	Glycoside hydrolase family 45	β-1,4-Endoglucanase	GH45	Cellulose
110289	GHY		7	1,860,196 to 1,863,081	Glycoside hydrolase family 36	α-1,4-Galactosidase	GH36	Hemicellulose
106075	MGAXE1		11	312,607 to 313,691	Esterase poly(3-hydroxybutyrate) depolymerase		CE1	Hemicellulose
95636	Mycgr3G95636		9	705,132 to 706,300	Esterase poly(3-hydroxybutyrate) depolymerase		CE1	Hemicellulose
61141	MGXYN2		7	2,500,781 to 2,502,585	Glycoside hydrolase family 10	β-1,4-Endoxylanase	GH10	Hemicellulose
60105	MGXYL1		6	1,657,807 to 1,658,690	Glycoside hydrolase family 11	β-1,4-Endoxylanase	GH11	Hemicellulose
51439	MGFAE3		13	742,995 to 744,596	Tannase and feruloyl esterase		CE1	Hemicellulose
49510	TAN		11	1,328,072 to 1,329,987	Tannase and feruloyl esterase		CE1	Hemicellulose
46238	FAE		8	2,030,409 to 2,032,034	Tannase and feruloyl esterase		CE1	Hemicellulose
105323	MGABN1		8	140,967 to 142,117	Glycoside hydrolase family 43	β-1,4-Xylosidase	GH43	Hemicellulose or pectin
96505	MGXYL7		11	298,906 to 300,120	Glycoside hydrolase family 43	β-1,4-Xylosidase	GH43	Hemicellulose or pectin
75584	BXL		9	1,079,287 to 1,081,613	Glycoside hydrolase family 43	β-1,4-Xylosidase	GH43	Hemicellulose or pectin
68922	MGABFB3		2	2,908,495 to 2,909,491	Glycoside hydrolase family 62	Arabinoxylan arabinofuranohydrolase	GH62	Hemicellulose or pectin
51381	Mycgr3G51381		13	304,851 to 305,879	Glycoside hydrolase family 53	β-1,4-Endogalactanase	GH53	Hemicellulose or pectin
85457	MGPEL3		3	3,248,792 to 3,249,927	Pectate lyase	Pectin lyase	PL1	Pectin
69329	MGRHA1		3	471,147 to 473,377	Glycoside hydrolase family 78	α-Rhamnosidase	GH78	Pectin
66866	MGPME1		1	4,254,552 to 4,255,738	Pectinesterase		CE8	Pectin

The most prominent functional category of genes with peak up-regulation at 21 dpi includes a variety of Zn²⁺-containing alcohol dehydrogenases. These functional classes are most prevalent in the genes listed in expression cluster 1 (Table IV). The function of many of these is currently unknown; however, it is intriguing that the simultaneous analysis of metabolites in this tissue also identified very high levels of sugar alcohols, including mannitol and arabitol, at 21 dpi (Fig. 2; Supplemental Table S1). Various other reductases (including aldo/keto) and dehydrogenases (Glc/ribitol) with undefined substrate specificities also show peak up-regulation at 21 dpi, suggesting an involvement in metabolism related to the accumulation or interconversion of abundant spore compounds. In support of this, the Z. tritici homolog of mannitol dehydrogenase (Mycgr3G104277) showed peak expression at 21 dpi, having already displayed significant up-regulation at 14 dpi (Supplemental Table S4). Homologs of this gene, and mannitol production in general, have been shown to be involved in asexual sporulation of the related wheat-infecting fungal pathogen P. nodorum (Solomon et al., 2005, 2006).

For Wheat

The most significant transcriptional reprogramming in wheat detected in this study occurred in infected leaf tissue at 14 dpi, coincident with large-scale host cell death. As many as 6,608 genes (i.e. 13.5% of the analyzed transcriptome) were found to be significantly up- or down-regulated at this stage of infection (Table V).

Of 3,056 up-regulated genes in total, 1,055 were of unknown function (Supplemental Table S6). Interestingly, the vast majority of these (94%) appear to be up-regulated more than 2-fold specifically at 14 dpi, implying that these are transcriptional programs associated with the fungus colonizing in its more rapid necrotrophic growth phase. Analysis of differentially expressed genes with functional annotation revealed that those up-regulated at 14 dpi can be largely grouped into the same overall functional categories as those up-regulated at 9 dpi, with a few exceptions (Supplemental Table S7). Thus, it appears that the main difference between lists of functionally annotated genes up-regulated at 14 and 9 dpi is in the numbers of genes residing in each functional category. The most noticeable differences were the larger numbers of genes belonging to the following functional categories: amino acid, lipid, nucleotide, and secondary (isoprenoid) metabolism, photosynthesis, and tricarboxylic acid cycle. Perhaps most significantly, the average level of gene expression (fold change versus the corresponding mock treatment) at 14 dpi is noticeably higher (on average, approximately 30-fold) than that at 9 dpi (approximately 6-fold), suggesting that peak plant defense activation is occurring during this rapid fungal necrotrophic colonization phase. Conversely, transcripts encoding photosynthesis-related processes are strongly down-regulated at 14 dpi.

Comparative analysis of metabolite data at 14 dpi now highlights a very clear and large accumulation of wheat fructans as a consequence of continued transcriptional up-regulation of this biosynthetic pathway from 9 to 14 dpi (Fig. 11; Supplemental Tables S1 and S6). However, by 21 dpi, when plant tissue appears necrotic, the levels of fructans detected in this tissue were diminished, suggesting that the fructan pool either had been metabolized by the fungus or mobilized away from the infection court. Similarly, while transcriptional activation of plant secondary metabolism including lignification was up-regulated at 9 and 14 dpi, key metabolites such as ferulate were detected maximally in the residual 21-dpi tissue (Supplemental Table S1), suggesting that, in contrast, the fungus had not effectively degraded/remobilized these metabolites and that they are most probably embedded within the residual plant cell wall structures. Finally, a sharp increase in the plant defense hormone SA was detected at 14 dpi that persisted to a lesser extent in the residual plant tissue at 21 dpi (Fig. 12; Supplemental Table S1). This suggests that SA-mediated defenses may be transiently triggered during the period of extensive host cell death involving the rapid necrotrophic colonization by Z. tritici.

Figure 12.

SA accumulation coincides with the necrotrophic phase of fungal colonization. The red line indicates fungus-inoculated leaves, and the blue line indicates mock-inoculated leaves. **, Significant at adjusted P < 0.05.

Overall Summary of the Necrotrophic Phase and Asexual Reproduction

The maximal plant defense response was detected at 14 dpi, and this correlated inversely with a decreased overall expression of putative fungal effectors, which were ultimately switched off entirely by sporulating Z. tritici at 21 dpi. During this period, the fungus appeared less reliant on secreted proteases and began to invest more effort in producing its putative cell wall-attacking capability, despite this genome feature being unusually low in numbers for a necrotrophic pathogen (Goodwin et al., 2011; Morais do Amaral et al., 2012). This suggests some degree of effort to utilize carbohydrates from plant cell walls late in infection. Despite this, the plant continues to attempt to reinforce cell walls via lignification and continues to maintain the up-regulation of the JA pathway and potentially also SA-mediated defenses (perhaps transiently). Wheat fructan metabolites are now clearly detected at high levels at 14 dpi, suggesting use by the plant to either minimize cell death or remobilize intracellular carbohydrates. This transcriptional and metabolite analysis highlights a clear overlap with abiotic stress responses in wheat (Valluru and Van den Ende, 2008). By 21 dpi, upon complete defeat of plant tissues, major fungal metabolites accumulate that can be associated with asexual spore formation, and certain classes of enzymes that may be involved in generating and/or interconverting spore storage compounds show peak up-regulation. Overall, these data suggest that, at 14 dpi, the fungus thrives in plant tissue that in fact is undergoing a massive defense response involving localized cell death. The fact that the pathogen reproduces in this environment suggests that a form of hijack of plant cell death and defense activation is a feature of successful Z. tritici colonization, perhaps facilitating reproduction during this period.

CONCLUSION

This study describes an in-depth analysis of the interaction between Z. tritici and wheat leaves, via combined genome-wide RNAseq and global LC-MS- and GC-MS-based metabolite profiling, spanning a full asexual disease cycle of the pathogen on its host as well as two in vitro growth conditions. This analysis allowed several novel key findings on the Z. tritici-wheat interaction to be made: (1) relatively minor contributions of genes on dispensable/accessory chromosomes to overall gene expression at any infection phase; (2) rapid and large-scale in planta up-regulation of fungal genes encoding candidate effector proteins; (3) different infection phase-specific fungal nutritional programs; (4) novel coexpressed secondary metabolite clusters; and (5) early repression followed by specific activation of plant defense responses involving the biosynthesis of fructans, SA, JA, and cell wall reinforcement at the switch to necrotrophy and pathogen reproduction.

The collective temporal regulation of all these processes indicates that, while nutritionally, Z. tritici appears to be a relatively poor biotroph, it uses an effective biphasic interaction with plant immunity involving an initial subterfuge strategy that is followed by an apparent hijack through sudden acute activation. How this is achieved for a strictly extracellular colonizing (apoplastic) filamentous pathogen remains unclear, but it may involve any number of the strongly up-regulated secretion responses (proteins and/or metabolites) detected in this study, which are now open to functional analyses. This strategy may serve as an intriguing model for many other plant-infecting fungi, which are currently perceived to have hemibiotrophic or necrotrophic lifestyles.

MATERIALS AND METHODS

Biological Materials/Samples

The Zymoseptoria tritici isolate IPO323 and wheat (Triticum aestivum) ‘Riband’ were used in all experiments. All materials were freeze dried prior to RNA or metabolite isolation. Fungal cultures were propagated in shake flasks (220 rpm) at 18°C for 3 d (for PDB) or 5 d (for CDB) and then harvested via filtration. These incubation periods were determined to be within the logarithmic growth phase Z. tritici. Plant inoculation experiments were done as described previously (Keon et al., 2007) using a spore density of 1 × 10⁶ spores mL⁻¹ in 0.01% (v/v) Tween 20 in sterile water. Mock inoculations of plants were made using just the Tween 20 water solution.

Each biological replicate plant sample for either RNA or metabolite isolation was made up of five leaves collected from five independent plants randomly distributed in a single walk-in temperature- and humidity-controlled glasshouse. Leaves were immediately frozen in liquid nitrogen, freeze dried, and then ground to fine powder in liquid nitrogen.

Metabolomics

Metabolite analysis was performed at Metabolon (Research Triangle Park, NC). The global unbiased metabolic profiling platform consists of a combination of three independent platforms: UHLC-MS/MS optimized for basic species, UHLC-MS/MS optimized for acidic species, and GC-MS. These platforms have been described in detail in previous publications (Evans et al., 2009; Ohta et al., 2009). Essentially, 20 mg of each sample was thawed on ice and extracted using an automated MicroLab STAR system (Hamilton) in 400 mL of methanol containing recovery standards. ultraperformance liquid chromatography/MS was performed using a Waters Acquity UHPLC device coupled to an LTQ mass spectrometer (Thermo Fisher Scientific) equipped with an electrospray ionization source. The extracted samples were then split into three equal aliquots for analysis by the three methods. For the two liquid chromatography methods, chromatographic separation followed by full-scan mass spectra was carried out to record retention time, M  _r (mass-to-charge ratio), and MS/MS of all detectable ions presented in the samples (Evans et al., 2009). For GC-MS, bis-trimethyl-silyl-triflouroacetamide-derivatized samples were analyzed on a Thermo-Finnigan Trace DSQ fast-scanning single-quadrupole mass spectrometer operated at unit mass resolving power (Ohta et al., 2009). Compounds were identified by automated comparison of the ion features in the experimental samples with a reference library of chemical standard entries that included retention time, M  _r (mass-to-charge ratio), preferred adducts, and in-source fragments as well as their associated MS/MS spectra (DeHaven et al., 2010, 2012). This library allowed the rapid identification of metabolites in the experimental samples with high confidence. Comparison of experimental samples with process blanks (water only) and solvent blanks allowed the removal of signal derived from artifacts.

Data Imputation and Statistical Analysis

The missing value for a given metabolite was imputed with the observed minimum detected value, based on the assumption that they were below the limits of detection/sensitivity of the instrument. Statistical analysis of the data was performed using JMP (SAS; http://www.jmp.com) and R (http://cran.r-project.org). Welch’s two-sample t tests were performed on the log-transformed data to compare the treatment and control groups at each time point. Multiple comparisons were accounted for with the false discovery rate method, and each false discovery rate was estimated using q values. For the convenience of data visualization, the raw area counts for each biochemical were rescaled by dividing the value for a specific biochemical in each sample by the median value observed for that specific biochemical (median scaled).

RNAseq

Total RNA was isolated from freeze-dried tissues using the Trizol procedure incorporating a final LiCl₂ precipitation. RNA was further purified using the Qiagen Plant RNA isolation kit columns. The mRNA was isolated from 4 μg of total RNA, then subjected to fragmentation followed by complementary DNA (cDNA) synthesis. The ends of the cDNA fragments were repaired, and A tails were added. Sequencing adapters containing individual index sequences were ligated to the A-tailed fragments. The adapter-ligated cDNA library was amplified using 12 cycles of PCR. All reactions were done using NEBNext modules from New England Biolabs. The libraries were run on an Agilent Technologies Bioanalyzer to determine the size and concentration and then normalized to 10 nm for sequencing. Sequencing was done using an Illumina HiSeq 2000 sequencer to generate 100-base single-end reads. Five libraries were run per lane, with an average yield of 45,000,000 reads per sample (range, 30,000,000–67,000,000 reads per sample).

In order to remove bias from targets such as ribosomal sequences (22 Mycosphaerella spp. ribosomal RNA [rRNAs] and 99 wheat unigenes), these were removed from the analysis before proceeding with data normalization to get the read count (expression values). All short reads (single end) were mapped with TopHat (version 2.0.6) against the Z. tritici genome (-G Mycosphaerella_graminicola.MG2.16.gtf; Trapnell et al., 2012). Cufflinks (version 2.1.1) was used to calculate FPKM values for reference annotations (-G Mycosphaerella_graminicola.MG2.16.gtf) but excluding genes annotated as rRNA (-M rRNA_genes.gtf). Differential expression analysis was done with Cuffdiff (cuffdiff -u -M rRNA_genes.gtf -b Mycosphaerella_graminicola.MG2.16.dna.toplevel.fa). We extracted all significant genes (P < 0.05) between CDB and all other samples (PDB, 1, 4, 9, 14, and 21 dpi). Analysis was performed on both the three-biological replicates experiment and a two-biological replicates experiment. This was due to biological replicate 3 disease tracking slightly faster than biological replicates 1 and 2 and also because biological replicate 3 was sequenced on a different platform (generating only 75-bp reads) 2 years previously. Only genes that were consistently differentially expressed in both analyses were considered further. The total output of the differentially expressed genes from Z. tritici was further filtered to remove any predicted differentially expressed gene where very low expression levels overall were detected. This involved the removal of genes where FPKM values never exceeded less than 5 in all samples. This minimum cutoff was selected following a dispersal analysis on the range of mean FPKM values for fungal growth in both media and at every time point of plant infection (Supplemental Fig. S4, A and B). A further cutoff was then applied at a 2-fold change in expression. These steps were taken to increase the stringency of identifying robust differential expression events and led to 3,034 Z. tritici candidate genes. This list of fungal differentially expressed genes was then subjected to K-means clustering using the csCluster method in cummeRbund (version 2.0.0). We tested different K-values ranging from 1 to 16. We selected K = 12, as this produced a set of distinct profiles with about 60% of the total original dissimilarity.

Transcriptomes of wheat leaves infected with Z. tritici (1, 4, 9, and 14 dpi) along with the mock (buffer)-inoculated wheat leaves of identical age (i.e. corresponding to 1, 4, 9, and 14 dpi) were also studied. We preprocessed all reads by removing reads that mapped to wheat rRNA sequences or Z. tritici transcripts (bowtie2-very-fast-ungz filtered.fastq.gz). The remaining reads were mapped with bowtie2 to a wheat transcriptome reference sequence to generate raw counts. This wheat transcriptome reference was constructed using wheat Unigene Build 60 (http://www.ncbi.nlm.nih.gov/UniGene/) by filtering out all the contaminating nonplant sequences (i.e. those with more than 80% coverage of unigene and more than 90% similarity to the genome of the known wheat pathogens [e.g. Z. tritici, Fusarium graminearum, Puccinia spp., and Blumeria spp.]). The final filtered wheat unigene data set comprised 48,878 sequences. We then functionally annotated these unigene sequences using the Mercator pipeline (Lohse et al., 2014), which uses the MapMan BIN ontology (Thimm et al., 2004), tailored for functional annotation of plant omics data.

The resulting raw counts data were subjected to in-house statistical analysis for the identification of differentially expressed genes as follows. Eight treatments were used forming a two-by-four factorial set. The two factors were treatment (mock-inoculated control or infected wheat leaves) and time (1, 4, 9, and 14 dpi). Three separate replicated biological experiments, forming three statistical blocks, were done. In the first two blocks, both control and treated samples were observed, but in the third block, only data from genes of the control samples were observed. Hence, there were three replicates of control samples but only two replicates of infected samples at the four time points, making 20 RNAseq samples altogether. Regression analysis was used on a per gene basis to assess the significance of differences between the treatments, times, and interaction between these two factors. A generalized linear model to the count data for each gene was generated, assuming a Poisson distribution with a log link and accounting for the difference in sequencing depth by way of an estimated offset variable in the model:

where count_ijk is the count for Block_i, i = 1, 2, 3; Treatment_j, j = 1 (mock) or 2 (infected); and Time_k, k = 1 (1 dpi), 2 (4 dpi), 3 (9 dpi), or 4 (14 dpi); and where d  _ijk is the offset variable with values estimated as the proportion of total counts per sample using a set of genes proposed not to be differentially expressed and therefore contributing only to the assessment of sequencing depth, based on a Poisson goodness-of-fit statistic. The analysis took account of variance due to the three experiments as statistical blocks, with 2 degrees of freedom (df), before testing the effects of treatment (with 1 df), time (with 3 df), and the interaction (with 3 df). Therefore, with 20 observations, F tests (rather than χ² tests) for the effects were based on 10 (residual) df, these F tests being appropriate to account for any overdispersion in the data above that expected for a Poisson distribution.

Following the analysis, the P values from the F test results for the differences between treatments, time, and the interaction were adjusted using the Benjamini and Hochberg (1995) method (B-H) to account for false discovery at the 5% level of significance given 48,878 genes. Predicted means back on the raw Illumina read counts scale for these genes were then also output and compared using the approximate lsd value at the 5% level of significance. This allowed 12,901 differentially expressed wheat genes with a significant interaction between treatment and time (P < 0.05, F test, and Benjamini and Hochberg [1995]-adjusted P < 0.05) to be output. The lsd values were approximated, as the model used was nonlinear. In particular, these values were likely to have been poorly estimated for those unigenes represented by very low Illumina read counts. Therefore, these very lowly expressed genes (with reads mapped per million total reads < 2) were removed from further analyses. A further cutoff was then applied at a 2-fold change in expression. This resulted in the identification of 7,441 statistically and biologically most important wheat differentially expressed genes.

Real-Time qRT-PCR

In order to confirm selected transcriptional changes in planta, quantitative real-time PCR was conducted on cDNA synthesized from RNA from Illumina-analyzed replicates combined. A 200-ng aliquot of total RNA was used for each sample, and cDNA was synthesized using the Applied Biosystems High-Capacity cDNA Reverse Transcription Kit according to the manufacturer’s instructions. In order to check for genomic contamination, a standard PCR assay was conducted using intron-spanning primers. Only cDNA samples for which a single band of the correct size was observed were used. The cDNA samples were diluted 1:3 and subjected to quantitative real-time PCR using SYBR Green JumpStart Taq ReadyMix (Sigma) in accordance with the manufacturer’s instructions. Thermocycler settings were as follows: 2 min at 95°C, followed by 50 cycles of 15 s at 95°C, 30 s at 60°C, and 45 s at 72°C. Melt curves were run for all primer pairs (Supplemental Table S10) in order to check for dimerization. Data were analyzed using the Applied Biosystems 7500 software version 2.0.1, and the Z. tritici β-tubulin gene was used as a reference.

Fungal Transformation

Agrobacterium tumefaciens-mediated fungal transformation was performed as described previously (Motteram et al., 2009) using the Ku70 modified form of isolate IPO323 (Bowler et al., 2010). Constructs designed to either delete or disrupt candidate secreted protein effector genes were generated in vector pCHYG (Motteram et al., 2009). Several independent transformants with each gene disrupted were validated by PCR on genomic DNA. All primer sequences used for the design of knockout constructs and targeting validation of transformants are shown in Supplemental Table S10.

Cryoscanning Electron Microscopy

Sections of inoculated wheat leaves (approximately 5 × 5 mm) were excised using a sterile blade and attached to cryo stubs using a smear of Optimal Cutting Temperature compound (Agar Scientific). Slotted stubs were used to enable the freeze fracturing of leaves. The mounted samples were plunged into prefrozen liquid nitrogen and then transferred under vacuum to the cooled (−180°C) cryo chamber.

Fractured surfaces were prepared and imaged using the JSM LVSEM 6360 scanning electron microscope (JEOL). Leaves were fractured using a cold blade mounted in the Alto 2100 (Gatan) cryochamber, followed by sublimation to remove any ice on the surfaces. This was performed by raising the temperature of the stage to −95°C for 1 min. The temperature was allowed to recover to −150°C before coating the samples with gold for 1 min (approximately 10 nm thickness). The samples were transferred to the scanning electron microscope with the Gatan stage temperature maintained at −150°C for examination and imaging in high-vacuum mode.

Imaging of leaf surfaces inoculated with the fungal pathogen was done under cryo mode in a JSM FEGSEM 6700 scanning electron microscope (JEOL). Sample preparation in the attached Alto 2500 (Gatan) cryochamber was performed in a similar way to remove any contaminating ice, and the samples were coated with gold-palladium (approximately 5-nm thickness).

Supplemental Data

The following supplemental materials are available.

• Supplemental Figure S1. Schematic overview of the experimental workflow.

• Supplemental Figure S2. Generation and testing of fungal secreted protein effector mutants.

• Supplemental Figure S3. qRT-PCR expression validation on the PKS7 metabolite cluster.

• Supplemental Figure S4. Fungal RNAseq statistical analysis.

• Supplemental Table S1. Full fungal and wheat metabolite data.

• Supplemental Table S2. Metabolites detected only in fungal or plant samples.

• Supplemental Table S3. Mean FPKM values for fungal genes.

• Supplemental Table S4. All differentially expressed fungal genes.

• Supplemental Table S5. All differentially expressed genes encoding fungal secreted proteins.

• Supplemental Table S6. All differentially expressed wheat genes.

• Supplemental Table S7. Differentially abundant wheat functional gene categories at 9 and 14 dpi.

• Supplemental Table S8. The most abundant fungal transcripts (FPKM) at 9 dpi.

• Supplemental Table S9. Wheat genes ranked by transcript abundance.

• Supplemental Table S10. Primer sequences.

ACKNOWLEDGMENTS

We thank members of the Rothamsted Research Horticultural support team.

Glossary

 
• PAMP

  pathogen-associated molecular pattern

 
• PCD

  programmed cell death

 
• JA

  jasmonic acid

 
• SA

  salicylic acid

 
• RNAseq

  RNA sequencing

 
• dpi

  days post infection

 
• qRT

  quantitative reverse transcription

 
• CDB

  Czapek-Dox broth

 
• PDB

  potato dextrose broth

 
• GC-MS

  gas chromatography-mass spectrometry

 
• LC-MS

  liquid chromatography-mass spectrometry

 
• FPKM

  fragments per kilobase of exon per million fragments mapped

 
• IPR

  InterPro

 
• OPDA

  12-oxo-phytodienoic acid

 
• CAD

  cinnamyl alcohol dehydrogenase

 
• ABC

  ATP-binding cassette

 
• cDNA

  complementary DNA

 
• rRNA

  ribosomal RNA

 
• df

  degrees of freedom

Author notes

1

This work was supported by the Biotechnology and Biological Sciences Research Council of the United Kingdom (BBSRC) through the Institute Strategic Program (grant no. BB/J/00426X/1, 20:20 Wheat, to Rothamsted Research) performed in collaboration with Syngenta Biosciences and by a Syngenta-sponsored Industrial Case (BBSRC) doctoral studentship (to M.D.).

*

Address correspondence to .

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Jason J. Rudd ().

[OPEN]

Articles can be viewed without a subscription.

www.plantphysiol.org/cgi/doi/10.1104/pp.114.255927