Deciphering the modes of action of Golubevia sp., an antagonist against the causal agent of powdery mildew in wheat, using an mRNA-based systems approach

Biocontrol agents are living organisms with the potential to suppress populations of plant pathogens or pests in a cropping system. The complex interplay between the different players and the changing environment, results in a combination of different modes of action. Here, we applied an mRNA-based systems approach to gain insight into the antagonist-pathogen-host interaction of Golubevia sp. isolates BC0812 and BC0850 with the causal agent of wheat powdery mildew, Blumeria graminis f. sp. tritici, in planta over time. Bioassays were performed on potted wheat plants (water-treated control, antagonist, pathogen, antagonist+pathogen) under controlled conditions. A significantly higher percentage of mildew conidia were parasitized after treatment with Golubevia sp. BC0812 with 26% and BC0850 with 16% compared to the water control with 1%. Differential gene expression analysis of antagonists, pathogen and host 5, 6, 7, and 11 days after inoculation (dai) with the antagonist pointed to a combination of different modes of action: An interplay of modulating plant defense responses, impairing conidiogenesis of the pathogen by scavenging H2O2, facultative hyperparasitism and nitrogen competition. Microscopic observations supported the suggested hyperparasitism as thin mycelium could be observed on Bgt conidia at 6 dai and later. Taken together the results allowed the formulation of new hypothesis regarding modes of action and the interplay between antagonist, pathogen and host. It showed that a solid molecular understanding of the antagonist-pathogen relationship over time is essential for less biased mode of action studies. Understanding this complex interplay is the basis for targeted optimization strategies and allows discovery of new potential targets and markers for future biocontrol development.


Introduction
Powdery mildews are obligate biotrophic fungal pathogens that reproduce only in interaction with living host cells, thereby indispensably linking their lifecycle to that of their host. This dependency results in an arms race between mildew pathogens and their hosts (Kemen and Jones, 2012). The genomes of Blumeria graminis f. sp. tritici (hereafter Bgt) causing powdery mildew in wheat, (Wicker et al., 2013;Müller et al., 2019), and its close relative the barley-infecting Blumeria graminis f. sp. hordei (Bgh) (Spanu et al., 2010) reflect this lifestyle: They show a very high number of repetitive and transposable elements and many species-specific candidate (secreted) effectors. The dynamic nature of their genomes also involves a high risk to develop resistance to fungicides (Fungicide resistance action committee, www.frac.info). Therefore, powdery mildew diseases have been assigned a high research priority within Europe for the development of biological control products (Lamichhane et al., 2017).
A recent screening for potential antagonists against Bgt identified fungal isolates that significantly reduced the number of powdery mildew pustules per flag leaf (Köhl et al., 2019a). The isolates belong to the new genus Golubevia spp. and are closely related to G. heteromorpha and G. pallescens (previously known as Tilletiopsis pallescens) (Wang et al., 2015;Richter et al., 2019). Representatives of this genus are inhabitants of the phyllosphere and their antagonistic potential against rose and cucumber powdery mildews has been shown previously (Urquhart et al., 1994;Ng et al., 1997). Studies on a G. pallescens-like isolate and its potential mode of action so far focused on in vitro measurements of lytic enzymes and bioassays with cell extracts that resulted in plasmolyzed conidia and collapsed hyphae and conidiophores suggesting a mycoparasitic lifestyle (Urquhart and Punja, 2002).
Biocontrol agents of plant pathogens are defined as living organisms that facilitate suppression of populations of plant pathogens in a biological system (Heimpel and Mills, 2017). It is often impossible to pinpoint the mechanism of disease suppression to a single mode of action. The complex interplay between antagonist, pathogen, host and environment over time often results in a combination of different modes of action (Köhl et al., 2019b). This can be the indirect or direct interactions through induction of plant resistance genes, competition for resources, modulation of signaling cascades, competitive exclusion, hyperparasitism or antibiosis. Until now, many studies still limit the selection of potential biocontrol agents and the investigation of their mode of action to assaying known mechanisms in vitro with single timepoint measurements (Shakeel et al., 2018;Bunbury-Blanchette and Walker, 2019;Soenens and Imperial, 2019). This approach neglects the complexity of the biological system in which a biocontrol agent is supposed to fulfill its role, often leading to disappointing and inconsistent results under environmentally relevant conditions. Therefore, selection of biocontrol agents should be based on the direct screening of strains for their antagonistic potential in planta. Similarly, mode of action studies on selected strains should be performed under relevant conditions resembling targeted future application strategies.
Here, we applied an mRNA-based systems approach on potted wheat plants to gain insight into the antagonist-pathogen-host interaction of Golubevia sp. isolates BC0812 and BC0850 with Bgt in planta over time (Köhl et al., 2019a). Transcriptomic data was supported by microscopy and spore-counts to evaluate efficacy of the assay.

Wheat plants and fungi
Winter wheat seedlings cv. Julius (moderately resistant against Bgt) and Blumeria graminis f. sp. tritici (Bgt) FAL92315 were produced according to Köhl et al. (2019a). 30 leaves with Bgt pustules were cut from the newly infected plants and used within 1 h for inoculation of wheat plants in the bioassay.
Golubevia sp. isolates were produced in T2 medium (30 g L -1 glucose, 1 g L -1 yeast extract (Oxoid), 10 g L -1 bacteriological peptone (Oxoid), pH 7.8) in 30 ml medium in 100 ml-conical flasks on a shaker with 100 rpm during 15 days at 16°C. The obtained mycelium was harvested by centrifugation at 10000 g and homogenized in tap water with 0.01% Tween 80 with a sterile blender (at high speed for 5 s). Concentrations of obtained mycelial fragments were determined using a hemocytometer and adjusted by adding 0.01% Tween 80 in tap water to 1x10 7 mycelial fragments ml -1 for Golubevia sp. isolate BC0812 and 5x10 6 mycelial fragments ml -1 for BC0850 (due to differences in growth).

Bioassay
A bioassay was conducted in a randomized block design with four treatments (Bgt, Bgt+BC0812, Bgt+BC0850, Tween control) arranged in six blocks (replicates). Each experimental unit consisted of one pot with six wheat plants. Wheat plants were spray-inoculated using an atomizer with suspensions of the two Golubevia sp. isolates or water containing 0.01% Tween 80 as control treatment until run-off. After 24 h plants were dry inoculated in a closed inoculation box with Bgt conidia produced on wheat leaves using compressed air (Cantu et al., 2013). The obtained density was 3.6 Bgt conidia mm -2 on a hemocytometer placed in the box between the pots. Subsequently the plants belonging to the same block were transferred a plastic box (710 x 440 x 380 mm L x W x H) with a layer of water and kept at 15°C with 12 h light. At 13 dai with Bgt, the first leaves of the six plants per pot were cut and placed in a 50 ml-falcon tube that contained 10ml water with 0.01% Tween, followed by 1-minute vortex. The concentration of the conidial suspension was determined using a hemocytometer to determine the total amount of Bgt conidia and for healthy Bgt conidia not showing symptoms of collapse or other damage per cm -2 leaf surface. A sub-sample of the conidial suspensions was stained using aniline blue solution (20% lactic acid, 40% glycerol, 25% demi water, 15% of a 1% solution of aniline blue in demi water) and 2 x 25 arbitrarily chosen conidia per experimental unit were assessed for signs of hyperparasitism with visible attached or penetrated thin mycelium of Golubevia sp.. The percentage hyperparasitized Bgt conidia and numbers of nonhyperparasitzed Bgt conidia per leaf and of hyperparasitized Bgt conidia per leaf were calculated.
In parallel, four sets (replicates) of four pots each with six wheat plants were treated (BC0812, BC0850, Tween-treated control), incubated for 24 h and half of them inoculated with Bgt as described above. This resulted in the following treatments for RNASeq analysis: BC0812, BC0850, Bgt+BC0812, Bgt+BC0850, Bgt, Tween-treated control. Pots that obtained the same treatment were incubated together in a plastic box (710 x 440 x 380 mm L x W x H) under the same conditions as described above (to avoid any cross contamination of Golubevia sp. or Bgt during the experiment). First leaves were sampled after five, six, seven and 11 days after inoculation with Bgt. On each sampling date, five leaves from the same pot sampled and pooled and immediately stored in liquid nitrogen for RNA isolation and subsequent gene expression studies.

Statistics
Effects of treatments with the different Golubevia sp. isolates were compared with the water-treated control for log 10 -transfromed data of numbers of Bgt conidia per leaf and for percentages of parasitized Bgt conidia using unprotected LSD-tests (P = 0.05).

RNA isolation
Snap-frozen RNA samples were submerged in TRIzol™ reagent (ThermoFisher) and homogenized by bead-beating for 2x30sec at 8000rpm. Sample lysis and phase separation were performed according to the TRIzol™ protocol. The aqueous phase was mixed with 100% RNA-free EtOH (1:1 v/v) and transferred to an RNeasy Spin Column (Qiagen). Following steps were performed according to the RNeasy Plant Mini kit (Qiagen) protocol. RNA concentrations were determined with RiboGreen and analyzed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). Sequencing of the RNA was carried out by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China) on an Illumina HiSeq 2500 sequencer (Illumina) generating 12GB of 150bp PE reads per sample.

Gene expression analysis
Reads were cleaned by Novogene according to the following parameters: removal of reads containing adapters, reads containing N > 10% and reads containing low quality (Qscore < = 5) base. The resulting datasets were submitted to the Sequence Read Archive (SRA) at NCBI under bioproject number PRJNA577523 (SAMN13028261 -SAMN13028324). After cleaning, reads were mapped to the genomes of Blumeria graminis f. sp. tritici 96224 (Wicker et al., 2013), Golubevia sp. BC0812 and BC0850 (Russ et al., 2020) using the hisat2 software excluding discordant alignments (Kim et al., 2015). Gene expression counts were obtained with htseq-count version 0.11.1 and analyzed with DESeq2 version 1.14.1 (Love et al., 2014) in R comparing the effects of antagonist treatments in presence and absence of Bgt according to the manual. The selection of 75 Triticum aestivum (durum/ monococcum) genes involved in defense responses (Supplementary file 1 (Cantu et al., 2013) was mapped to raw transcripts with CLC genomics workbench 11.0.1 using 100% similarity mapping (Qiagen, Denmark).

Bioassay
On water-treated wheat plants, 2.82 x10 6 conidia cm -2 leaf had been produced 13 dai. Although the Golubevia sp. isolates BC0812 and BC0850 were genetically very similar (Russ et al. 2020), they did show differences in growth speed and pattern, leading to different inoculation densities.
Treatment with Golubevia sp. BC0850 reduced the number of Bgt conidia by 23 % (Table 1). Treatment with Golubevia sp. BC0812, applied at twice the concentration of mycelial fragments, reduced the number of Bgt conidia by 28 %. However, these treatment effects were not statistically significant. Development of thin mycelium as sign for hyperparasitism by Golubevia sp. isolates was observed on Bgt conidia at 6 dai and later ( Fig. 1). A significantly higher percentage of parasitized Bgt conidia were counted after treatment with Golubevia sp. BC0812 at 13 dai with 26% and BC0850 with 16% compared to the water control with 1%. The calculated numbers of non-parasitized Bgt conidia was significantly lower with 46 % after treatment with Golubevia sp. BC0812. Treatment with BC0850 significantly reduced the number of non-parasitized Bgt conidia by 35 %.

Read mapping to features
Mapping the QC trimmed reads against the genome of Bgt showed that at least 94 % of its coding sequences were expressed (Supplementary Figure 1), increasing to more than 99% until 11 dai with Bgt. The contribution of Bgt to all reads ranged from 0.4 (5 dai) to up to 14.5% (day 11) (Fig. 2A). The increase in gene expression of Bgt is time-as well as treatment-dependent: It did show a higher expression at 6 dai especially in combination with the antagonists and at 11 dai in all the treatments ( Fig. 2A).
A lower number of reads could be mapped to the coding regions of Golubevia sp., ranging between 0.17 and 1.24% (Fig. 2B). This is also reflected in the coverage of the transcriptome which is lower than for Bgt: Between 8.46% (day 11 Bgt+BC0812) and 29.09% (day 7 Bgt +BC0812) of the Golubevia sp. genes did not show any expression. At day 7 the relative read count and the genome coverage were significantly reduced in Bgt-treated samples for both isolates ( Figure 1). The percentage of reads mapped and the genome coverage at day 11 was higher for BC0812.
The overall reproducibility was confirmed by gene expression levels of the duplicates of Bgt and Golubevia sp.. In both cases the number of genes with a read count ≥1 were in a similar range (Supplementary Figure 1).

Differential gene expression in wheat defense response genes
Out of the 75 genes that are presumed to play a role induced resistance in wheat 74 were expressed in at least one of the treatments. Although there was no statistical significance a trend was visible in the differential expression of PR1.2 (AJ007349.1) (Fig. 3). Treatment with BC0812 and BC0850 increased the expression of PR1.2 in the presence as well as absence of Bgt. The gene was upregulated by 2-4-fold in the  L. Russ, et al. Biological Control 152 (2021) 104446 samples treated with antagonist and Bgt compared to treatment with Bgt alone (Supplementary Table 2). Antagonist-only treatment resulted in a 3-8-fold upregulation of gene expression compared to the watertreated control (Fig. 3). Also, a β-1,3-glucanase (Z22874.1) gave a similar expression pattern, with the non-Bgt-treated antagonist samples showing the highest gene expression (at least 3 times higher than the control and 2 times higher than the Bgt-treated samples with and without antagonist). A type II chitinase (AB029934.1, wrongly annotated as type I), was also upregulated in presence of BC0812 disregarding treatment with Bgt (Fig. 3). The other chitinases of type III (AB029936.1) and type IV (AF112966.1) were constitutively expressed at the same level in all treatments however. The genes coding for germin-like protein 2a (AJ237942.1) and allene oxide synthase (AY196004.1) were upregulated in the presence of Bgt (Fig. 3). Treatment with any of the antagonists alone did not result in an upregulated expression of these genes.

Bgt gene expression changes
The gene expression pattern in Bgt did not show strong treatmentdependent fluctuations, but rather seemed to cluster in a time-dependent manner according to the developmental stage of the fungus (Fig. 4). An exception was day 6: Here Bgt-only and antagonist+Bgt samples show more variance (Fig. 4, day 6). Underlying this separation in the PCA-plot was the downregulation (2-6x) of a high number of putative secreted effector proteins, hypothetical proteins, a sporulation-specific chitinase (BGT96224_1166), a catalase (BGT96224_717), an oligopeptide transporter (BGT96224_E5560) and an endoglucanase (BGT96224_848) (Fig. 5). The downregulation was clearly linked to the application of the antagonists, as duplicates of BC0812 and BC0850 showed similar results (Fig. 4). This shift in gene expression of Bgt was transient however, as following sampling timepoints did not show such a strong effect of the treatment with the antagonists anymore.
Only six other genes showed differential gene expression in a treatment-dependent manner. All of them were upregulated in combined treatments of antagonists and Bgt compared to Bgt alone (Fig. 5). An alpha carbonic anhydrase (BGT96224_3750), cytochrome P450 (BGT96224_2744) and a gene encoding an unknown protein of 544 amino acids (BGT96224_4015) were constitutively upregulated at all sampling timepoints (Fig. 5). Whereas the carbonic anhydrase was upregulated at similar levels (5-6x with BC0812), the cytochrome P450 showed highest differential expression at day 5 (~4x) and BGT96224_4015 at 11 dai (~4.5x). Another uncharacterized protein (BGT96224_1906) was only upregulated in the early interaction phase (5 dai, up to 3.4x), but did not seem does not seem to play a role at 11 dai, in contrast to BGT96224_2925, encoding a protein of unknown function which had increased expression levels at 11 dai. The same was true for an ammonium transporter (BGT96224_4139), which was upregulated at the later sampling timepoints (7 dai and 11 up to 3.4x).
Differential expression of Golubevia genes Golubevia sp. did not show a strong treatment-depended gene expression pattern, neither did gene expression changes cluster according L. Russ, et al. Biological Control 152 (2021) 104446 Fig . 3. Differential transcription of 74 Triticum sp. stress response genes at 5 dai with 6 different treatments. Heatmap shows log 2 (n+1) for normalized counts (generated by DEseq2).
to sampling time at 5-7 dai (Fig. 6A&B, Supplement 5&6). Although this is true for the overall gene expression during early interaction, both Golubevia isolates showed increased expression of a gene encoding an unspecific fungal peroxygenase (UPO) (BC0812_6932~10x and BC0850_407~3x) in the presence of Bgt during its conidiation phase on 6 and 7 dai. A clear treatment-dependent effect on in gene expression was visible at day 11, when Bgt treatment caused a shift in gene expression in both isolates (Fig. 6A&B). The effect was more evident in BC0812 (Fig. 6A). This could already be seen in the significantly higher number of reads that could be mapped to BC0812 in Bgt-treated samples at 11 dai (Fig. 2B). The genes that shared the same expression pattern at 11 dai in both isolates were all upregulated in a combined treatment (antagonist +Bgt) and included a metalloprotease (BC0850_6326 9x, BC0812_502 5x), a P-type ATPase (BC0850_6726 6x, BC0812_6550 9x) and a MFS transporter protein similar to the family of tetracycline resistance proteins (BC0850_2323 10x, BC0812_2492 5x). These genes shared 100% sequence identity in both strains. BC0812 also upregulated several genes involved in ribosome biogenesis at day 11 when Bgt was present. This included ribosomal proteins of both eukaryotic subunits as well as several histones (Fig. 7). A combined treatment also influenced nitrogen metabolism of BC0812 with a putative nitrite reductase (BC0812_4076, 5x) and a formate/nitrite transporter (BC0812_6863, 4.9x) both being downregulated. Accordingly, many other transporters of different families were also downregulated by 4-5-fold, including a transporter of the xanthine/uracil/vitamin C family (BC0812_1113), two putative purine permeases (BC0812_1331&4158), an uracil permease (BC0812_640), an amino acid transporter (BC0812_7306), a solute symporter (BC0812_1845) and two MFS family transporters (BC0812_2522&75). Interestingly, a similar number of transporters was upregulated, indicating a metabolic shift: a copper family transporter (BC0812_1333), the MFS-type transporter mentioned before (BC0812_2092), an Fe 2 /Zn 2 -regulated transporter (BC0812_546), an amino acid transporter (BC0812_585) and an ABC transporter (BC0812_8446).
Also, several lytic enzymes were upregulated in the presence of Bgt, including a class 3 lipase (BC0812_1601), two metalloproteases one of which mentioned before (BC0812_1897&502), a glycoside hydrolase family 26 (BC0812_6333) and a glycosyl hydrolase family 16 (BC0812_2475).

No clear evidence for induced resistance by Golubevia sp.
To investigate whether our selected biocontrol strains could potentially modulate induced resistance in wheat, reads from the different treatments at 5 dai were mapped to a selection of 75 Triticum aestivum genes that have been associated with stress response (Supplementary file). The overall differences in gene expression between the watertreated control and the plants treated with either BC0812 or BC0850 alone were not evident. When treated with either antagonist alone, several genes were upregulated in a similar pattern to a combined an-tagonist+Bgt treatment. Expression was even slightly higher than after treatment with Bgt alone. This included PR 1.2, which plays a still unknown role in early plant defense against biotic stress, (Molina et al., 1999;Breen et al., 2016), a class II chitinase and a β-1,3-glucanase. Class II chitinases are induced in the systemic response where they inhibit growth of fungal pathogens in an interplay with β-1,3-glucanases by causing lysis of hyphal tips (Grover, 2012). This pointed to a similar transcriptomic response of the plant to the antagonist alone as it did to the pathogen.
Although the role of Golubevia sp. in the phyllosphere and its interaction with the plant itself is not yet understood the increased expression of plant defense genes upon treatment with BC0812 could point to modulation of the plant defense responses. It is currently unknown whether the antagonist is acting on priming stress responses in the plant only or could possibly induce a real lasting stress in the plant. As we did not work under sterile conditions it is also possible that the effect of BC0812 on the plants gene expression is of a secondary nature: the antagonist might also have interacted with other community members leading to the resulting gene expression patterns.

Interaction of Golubevia sp. and Bgt during conidiation
Major treatment-dependent fluctuations in gene expression of Bgt were limited to day 6. This coincided with conidiogenesis of Bgt. The differential expression of genes during the conidial stage has been described previously hinting at significant physiological changes during asexual reproduction (Zeng et al., 2017;Hu et al., 2018;Pham et al., 2019b). The treatment of the host plant with either BC0812 or BC0850 in combination with Bgt seemed to delay the molecular onset of conidiation in Bgt on the transcriptome level. The expression of many putative secreted effectors as well as other genes that have been described to play a role in conidiophore erection were downregulated in Fig. 4. Principal component analysis plot showing sample-to-sample distances between the Bgt control (red), BC0812+Btg (green) and BC0850+Btg (blue) at 4 timepoints. The PCA plot is based on the differential expression analysis with the R package DEseq2 (Love et al., 2014).

Fig. 5.
Differential transcription of Bgt genes (with log 2 fold change ≥ ± 2 and p adj ≤0.001 in at least one sample). Heatmap shows the average log 2 fold change of combined treatments with Bgt and BC0812 or BC0850 at 4 timepoints relative to the Bgt-only control. samples taken on 6 dai ( Figure 2) when compared to the Bgt-only treated control. Accordingly, a catalase (BGT96224_717) catalyzing the decomposition of H 2 O 2 to water and oxygen was downregulated. Catalase is an important pathogenicity factor during Bgh infection and conidiation, to protect itself from the oxidative burst it elicits to suppress the plants immune response (Zhang et al., 2004).
We hypothesize that H 2 O 2 might also play a role in the mode of action of Golubevia sp., as both isolates upregulated an unspecified fungal peroxygenase (UPO) at day 6. This group of enzymes are hemecontaining extracellular monooxygenases that transfer one oxygen atom from a peroxide (H 2 O 2 ) molecule to diverse organic substrates (Hofrichter and Ullrich, 2006;Manoj and Hager, 2008;Hofrichter et al., 2020). So far, the only members of the heme-thiolate peroxidases have been studied intensively are the chloroperoxidase (CPO) of Caldariomyces fumago and the aromatic peroxidase of several agaric basidiomycetes (Morris and Hager, 1966;Ullrich et al., 2004). The biochemistry and potential value for organic chemistry is well known, but the biological and environmental significance of this group of enzymes remains to be resolved. They could play a role in the degradation of complex carbon and nitrogen compounds, such as plant exudates, lignin, toxins or xenobiotics (Hofrichter et al., 2010;Kellner et al., 2014;Karich et al., 2017). We assume that they also play an important role in the interaction of Golubevia sp. with the pathogen, as this group of enzymes has catalase activity, detoxifying the H 2 O 2 released by the pathogen during condiogenesis.
The genomes of BC0812 and BC0850 each encode at least seven different UPOs, presumably belonging to different subfamilies, as well as a truncated homolog (supplement 7) (Faiza et al., 2019). Although both isolates share the same set of UPO genes (100% similarity), they expressed different homologs in (BC0812_6932 and BC0850_407 showing only 48.1% similarity on amino acid level according to pairwise alignment with EMBOSS Needle. The biological replicates confirmed this observation, pointing to a general difference in the UPO regulatory network in these two very closely related strains. Whether the Golubevia isolates use the H 2 O 2 to drive the conversion of complex compounds with the expressed UPOs is unknown. However, actively decreasing the H 2 O 2 levels might explain the delay in the molecular onset of conidiation in Bgt. Lower H 2 O 2 levels have been associated with impaired conidiogenesis before (Zeng et al., 2017). The effect seemed to be transient, as the expression pattern between the treatments and the control seems to equilibrate again on day 7 (Fig. 3). It is interesting however to pursue whether this mechanism indeed adds to the biocontrol potential of the strains and how it may be exploited.

Bgt stress response to antagonist treatment
Several genes were constitutively upregulated when plants were treated with Bgt and the antagonist compared with Bgt alone disregarding the time of sampling. This includes an uncharacterized protein of 544 aa (BGT96224_4015) which also has a homolog in Bgh DH14 (bgh04406), but to which no known protein domains could be assigned. Although the function of this protein therefore remains elusive, it seems to be conserved and unique in Bgt making it an interesting target for future studies. Also, a putatively secreted alpha carbonic anhydrase of Bgt (BGT96224_3750) is upregulated in the presence of the antagonist. Like in other organisms, fungal carbonic anhydrases are at the basis of many biological processes such as respiration, pH homeostasis and CO 2 transport by interconverting CO 2 to HCO 3 - (Elleuche and Pöggeler, 2010). They have recently been described to also play an important role in cAMP signaling cascades that regulate gene expression under nutrient and oxidative stress as well as during reproduction and pathogenesis (Raffaele et al., 2010;Kusch et al., 2014;Martin et al., 2017). Also, a cytochrome P450 monooxygenase (BGT96224_2744) is constitutively upregulated in the presence of BC0812 and BC0850. This group of proteins plays a role in secondary metabolism, fatty acid metabolism and the degradation of xenobiotics (Nelson et al., 1996;Chen et al., 2014). We hypothesize that the carbonic anhydrase and the cytochrome P450 may be part of the Bgt stress response as their upregulation is linked to the presence of the antagonist in the phyllosphere.

Gene expression of Bgt and Golubevia sp. over time
Determining the percentage of the total trimmed reads assigned to either Golubevia sp. or Bgt and the number of genes without any hits in a time and treatment dependent manner, could give some insight into biological dynamics of both players. The coverage of the Golubevia sp. genome decreased only slightly over time in the absence of Bgt (Supplementary Figure 1). This supported the conclusion that Golubevia sp. is not an obligate hyperparasite on Bgt but can survive in the absence of Bgt. Furthermore, both Golubevia sp. isolates showed a strong decrease in coverage and overall read count at day 7 in combined treatments with Bgt ( Fig. 2 and Supplementary Figure 1), which might be related to the physiological changes triggered by the H 2 O 2 burst upon conidiation of Bgt at day 6. Although direct evidence is lacking, cell death of part the Golubevia sp. population is a plausible explanation for this decrease. The reason might by the oxidative burst elicited by Bgt during conidia formation. Especially BC0812 seemed to recover from this event until day 11, showing a similar percentage of reads mapped and read count ≥1 in the duplicates as on day 6. Samples treated with BC0850 generally showed a lower coverage of its genome and a higher mapping to Bgt, observations which were in line with its poorer performance in the bioassay (Tab. 1). The initially inoculated concentration of biomass for this isolate was also two times lower, probably too low to show efficacy after 11 days.
Gene expression changes in Golubevia sp. and Bgt at day 11 One of the factors that seemed to drive many of the observed changes in Golubevia sp. gene expression was a change in nitrogen Fig. 7. Differential transcription of genes with log 2 fold change ≥ ± 2 and p adj ≤0.001 for treatments of wheat plants with BC0812 (top) or BC0850 (bottom) in the presence and absence of Bgt at four timepoints. Heatmap shows log 2 (n+1) for normalized counts (generated by DEseq2). Column clustering is based on Euclidean distance.
availability. BC0812 downregulated a formate/nitrite transporter as well as a putative nitrite reductase in the presence of Bgt. Via this pathway nitrite is imported into the cell and converted to ammonium, probably the preferred nitrogen source of Golubevia sp.. Although some yeasts are known to possess only the nitrite reduction pathway (Linder, 2019) it is generally more common that nitrogen assimilation proceeds from nitrate via nitrite to ammonium or that ammonium, if available, is imported and used directly (Garrett and Amy, 1979;Marschner, 1995). In the list of genes that did not pass the 4-fold change threshold a nitrate transporter (BC0812_4077, 2.7x) and nitrate reductase (BC0812_4075, 2.9x), but also an ammonium transporter (BC0812_3003, 3.8x) which is highly expressed, showed a similar downregulation in the presence of Bgt at day 11. The group of transporters that were downregulated in the presence of Bgt were purine, uracil and amino acid permeases whose substrates can also be used as sole nitrogen source (Diallinas et al., 1995;Linder, 2019). In the absence of Bgt, Golubevia isolates depended on plant-derived nitrogen compounds leaching to the leaf surface (Tejera et al., 2006). Although it is assumed that nitrogen limitation is only a secondary issue on wellfertilized plants (Vorholt, 2012), it could be a limiting compound in our system, as plants did neither receive leaf-fertilization nor were they exposed to atmospheric deposition. The switch in the nitrogen assimilation pathway that becomes evident at 11 dai, is probably the result of ammonium accumulation in wheat leaves. It has been described previously that barley leaves infected with powdery mildew accumulate ammonium due to impaired assimilation by the infected plant (Sadler and Scott, 1974;Walters and Ayres, 1980). As nitrogen assimilation pathways via nitrate, nucleobases and amino acids are negatively regulated by the presence of ammonium, the increasing ammonium availability explains the downregulation of nitrogen assimilation pathways in Golubevia isolates. The reaction of Bgt to the changes in nitrogen availability is different: It upregulates its ammonium transporter (BGT96224_4139) in combined treatments from day 7 onwards. This could be the result of competition for ammonium with the antagonist. Bgt is an obligate biotroph with a reduced genome that strongly depends on what the host delivers, while Golubevia sp. is a known inhabitant of the phyllosphere with a versatile repertoire of nitrogen assimilation pathways. The response of both is therefore a logical consequence of their different lifestyles: Whereas Golubevia sp. is adapted to a nitrogen-limited environment and downregulates its uptake mechanisms, Bgt is less used to share and upregulates ammonium transporters to scavenge more. Differences in affinities of the ammonium transporters and growth speed of Golubevia sp. and Bgt respectively could play an additional role.
Next to the changes in nitrogen assimilation pathways a high number of lytic enzymes was differentially expressed 11 dai in combined treatments. Microscopic observations already hinted at a direct hyperparasitic interaction of Golubevia sp. and Bgt (Fig.1). An integral part of the parasite-host relationship is the production of enzymes that lyse host cell walls, including chitinase and β-1, 3-glucanase (Sivan and Chet, 1989;Sahai and Manocha, 1993;Geraldine et al., 2013). In a recent study on barley powdery mildew it was shown that the conidia of Bgh very actively express chitin deacetylase (Pham et al., 2019a). Chitin deacetylases catalyze the conversion of chitin to chitosan and play a role in cell wall formation, but also plant-pathogen interactions (Tsigos et al., 2000). Conversion of chitin to chitosan enables plant pathogens to avoid detection as chitosan does not elicit the plant immune response. The high expression of this enzyme led the researchers to hypothesize that the majority of chitin in the conidial cell wall is deacetylated to chitosan. In addition, the outer cell wall is decorated with galactofuranosyl residues and xylomannan side chains. These compounds are rarely present in fungal cell walls and point to the adaptation of Bgh and probably also Bgt to its obligate biotrophic lifestyle (Pham et al., 2019a). This implies that Golubevia sp. would require a more specialized set of genes to penetrate the Bgt cell wall. The genomes BC0812 and BC0850 did not encode a dedicated chitosanase or exo-β-d-glucosaminidase, which are able to hydrolyze chitosan. But, several non-specific enzymes such as common carbohydrases, proteases and lipases have been described to possess chitosanolytic activities as well (Thadathil and Velappan, 2014). Several candidate genes that could be involved in Bgt cell wall degradation were specifically upregulated in combined treatments: A lipase (BC0812_1601) and two metalloproteases (BC0812_1897&502) that could potentially act on chitosan, an endo-beta-1,3(4)-glucanase (B0812_2475) hydrolyzing linkages in beta-D-glucans and a gene encoding an enzyme of 831 aa (BC00812_6333). This enzyme contains a glycoside hydrolase family 26 domain at the N-terminal end. Most of proteins within this family are endo-1,4-beta-mannanases, that are less common in fungi (van Zyl et al., 2010). The C-terminal end does not contain any known domains. The encoded protein seems to be noncytoplasmic and so far, unique to Golubevia isolates: Closest blastp hits did show only~30% query coverage at~50% similarity and did belong to the Pucciniomycotina. BC0850 also expressed the gene (BC0850_3334) at day 11 in the presence of Bgt (Supplementary table), but the cutoff for significance (p adj ≤0.001) was not reached. The selective expression of this gene after conidiation Bgt points to its role in the parasitic relationship. It could be responsible for the degradation of the outer, mannan-containing cell wall of Bgt. As it is hardly expressed in the absence of Bgt, we assume that the enzyme does not play a role in active degradation of plant mannans under environmental conditions. The exact role of this enzyme remains speculative and should be investigated further.
Also, the expression of an ABC transporter (BC0812_8446) depended on the presence of Bgt at 11 dai. It showed some similarity to multidrug resistance proteins, responsible for the export of toxic compounds and possibly involved in the defense system of Golubevia sp..
In accordance with the increased relative gene expression and coverage of the transcriptome in the presence of Bgt, BC0812 also showed increased metabolic activity as it upregulated genes involved in ribosome biogenesis and several histones. This is the molecular onset for cell division and growth.

Conclusion
To unravel the modes of action of antagonist Golubevia sp. on Blumeria graminis f. sp. tritici we subjected material from our bioassays to RNA sequencing. In this way it was possible to measure the molecular first-line response of the host plant, Bgt and antagonists over time during interaction with each other. We found indications that Golubevia sp. modes of action could be an interplay of modulating plant defense responses, impairing conidiogenesis of Bgt by scavenging H 2 O 2 , facultative hyperparasitism and nitrogen competition. According to the results at hand we think Golubevia sp. to be more suitable for curative instead of preventive application. These hypotheses need further confirmation with targeted in situ measurements, microscopy and longterm bioassays.
Unlike the commonly performed experiments to determine mode of action (i.e. measuring lytic potential, antibiosis etc.) RNASeq is unbiased towards specific pathways and enzymes. This gives a much better picture of the complexity of a biological system and its many dependencies that are very closely linked to efficacy of a biocontrol agent (i.e. timing, nutrient availability, interaction partners etc.) (Köhl et al., 2019b). It also shows the importance of measuring the right parameters under environmentally relevant conditions (chitosanolytic activity instead if chitinolytic). The future development and optimization of biocontrol agents with proven efficacy could greatly benefit from such an approach, because it allows formulation of new hypotheses and optimization strategies.
Union's Seventh Framework Program for research, technological development and demonstration under grant agreement no 612713 and the Dutch Ministry of Agriculture, Nature and Food Quality (Program Durable Plant Production, project BO-25.10-001-002). We kindly acknowledge Dirk-Jan Valkenburg and Helen Goossen-van de Geijn for technical assistance, as well as Beatriz Andreo Jimenez for discussion and feedback.