Global transcriptome changes in perennial ryegrass during early infection by pink snow mould

Lack of resistance to pink snow mould (Microdochium nivale) is a major constraint for adaptation of perennial ryegrass (Lolium perenne L.) to continental regions with long-lasting snow cover at higher latitudes. Almost all investigations of genetic variation in resistance have been performed using cold acclimated plants. However, there may be variation in resistance mechanisms that are functioning independently of cold acclimation. In this study our aim was to identify candidate genes involved in such resistance mechanisms. We first characterized variation in resistance to M. nivale among non-acclimated genotypes from the Norwegian cultivar ‘Fagerlin’ based on relative regrowth and fungal quantification by real-time qPCR. One resistant and one susceptible genotype were selected for transcriptome analysis using paired-end sequencing by Illumina Hiseq 2000. Transcriptome profiles, GO enrichment and KEGG pathway analysis indicate that defense response related genes are differentially expressed between the resistant and the susceptible genotype. A significant up-regulation of defense related genes, as well as genes involved in cell wall cellulose metabolic processes and aryl-alcohol dehydrogenase (NADP+) activity, was observed in the resistant genotype. The candidate genes identified in this study might be potential molecular marker resources for breeding perennial ryegrass cultivars with improved resistance to pink snow mould.

blotch. It may also be increasingly important with the predicited climate change 12 . The projected climatic conditions 12 during fall may not allow plants to go through this process, and plants may therefore be covered with snow before they are cold hardened. Our overall aim was to identify genotypes that are able to resist/ tolerate a snow mold infection without cold hardening before exposure to the winter stress factors.
The use of molecular techniques for precise quantification of the fungal biomass in infected plants can facilitate the selection of resistant genotypes 13 . Usage of real-time PCR in quantification of plant pathogen infestation has increased in the last two decades. It is quicker, more specific and sensitive compared to traditional methods based on symptom assessment or plant dry weight 14 . Therefore, elongation factor 1a (EF-1a) gene was used in our study for accurate quantification of M. nivale DNA during snow mould infestation. The EF-1a gene has formerly been used to recognize M. nivale and M. majus as separate species 15 . The EF-1a has also been used to study the genetic variation among isolates 16 . Additionally, competitive PCR methods have also been developed for M. nivale and M. majus quantification in infected tissues 17 .
The use of whole transcriptome sequencing (RNA-seq) is an important tool in the analysis of complex plant responses 18 , and provides a more comprehensive understanding of transcription initiation sites, improved detection of alternative splicing events and the detection of gene fusion transcripts 19 . The technology is able to handle de novo sequencing of large genomes, revealing individual genome differences within the same species and quantify gene expression 20 . In particular, it now enables global transcriptome studies to be performed in non-model species that have lacked many of the array based assays that are successfully used to study gene expression in the model species.
In the present work, we have taken advantage of high throughput RNA-seq to study the global transcriptome changes in perennial ryegrass leaves during early infection by pink snow mould, more spesifically after four days incubation under artifical snow cover ( Supplementary Fig. S1). The aim of the present study was to identify molecular responses involved in snow mould resistance independent of cold acclimation in perennial ryegrass.

Results
Snow mould resistance test. Eight genotypes of L. perenne cv. Fagerlin were randomly selected and evaluated for resistance to M. nivale. There was a significant effect of genotype on resistance, measured as relative regrowth after inoculation and incubation under artificial snow cover for 6 and 8 weeks (Fig. 1). The differences between genotypes were most pronounced after 6 weeks. Overall, genotype M had the highest relative regrowth, while there were small differences among the other genotypes (Fig. 1). Based on the results of this resistance test we selected two genotypes for studying global changes in the transcriptome during snow mould infection; genotype F as a "susceptible" genotype (later termed S) and genotype M as a "resistant" genotype (later termed R).
Quantification of M. nivale DNA. Visual assessment of disease severity and quantification of M. nivale DNA in leaf and stem tissue of the eight genotypes showed that genotypes with severe symptoms of injury (such as genotype F) contained significantly more fungal DNA than the other genotypes, especially after 6 weeks of incubation (Fig. 2). Genotypes that contained most fungal DNA, e.g. genotype F, were also most symptomatic based on visual scoring (Fig. 2). There was a significant positive correlation (r = 0.489, P ≤ 0.05) between disease severity and the amount of M. nivale DNA, indicating the genotypes containing most fungal DNA also showed most disease symptoms (Supplementary Table S1). The correlations between relative regrowth and the amount of M. nivale DNA, and relative regrowth and visual scoring of symptoms were non-significant. Tissue samples collected from the "resistant" genotype M(R) and the "susceptible" genotype F(S) showed small differences in the amount of fungal DNA 1 day after inoculation, but 4 days after inoculation the differences were more significant in the "susceptible" genotype F(S) (Supplementary Fig. S2). De novo based susceptible (S) and resistant (R) transcriptome assemblies. A total of 178 million reads and 165 million reads of 100 bp were generated for the S and R genotypes, respectively (Table 1). Separate transcriptome assemblies were generated for each genotype using all their respective reads. The de novo assembly yielded 261,978 contigs for the S genotype, with N50 of 1,784 bp, and 188,355 contigs for the R genotype with N50 of 1,672 bp ( Table 1). The longest assembled contigs in the S and R genotype were 17,632 and 12,882 bp, respectively. To estimate the quality of the assemblies, we compared them to the Brachypodium distachyon coding sequence consisting of 31,029 entries. There were 27,135 B. distachyon sequences (87.45 percent) that had a significant hit in the S transcriptome assembly and 27,399 (88.30 percent) that had a significant hit in the R transcriptome assembly. Further, we used the CEGMA pipeline 21 to evaluate the completeness of our assemblies. The percentage of complete core eukaryotic genes (CEGs) in R and S assemblies are 82.66 and 93.95, respectively, and the percentage of partially complete CEGs ranged from 90.73 to 98.79 ( Table 2). The average number of orthologs per CEG in the R and S assemblies is 3.72 and 3.83, respectively, and the percentage of detected CEGs that had more than one ortholog was 96.1 and 97.0, respectively.    (I-NI) , and 3: Inoculated plants that were incubated under artificial snow-cover for four days (I-I) (Artificial snow cover created by covering the plants with moist cellulose tissue paper and black plastic sheets in a cold chamber at 2 °C in darkness. The reads from each sample were mapped onto their respective genotype specific (S and R) de novo assemblies and to the reference inbred L. perenne transcriptome 22 . In the case of each sample, more than 83-90 percent of the reads mapped onto the assembled transcripts. Using the genotype-specific assemblies in a series of pairwise comparisons between samples, 2,354 and 3,748 differentially expressed transcripts were identified with false discovery rates (FDR) < 0.05 between NI-NI and I-I samples; 1,602 and 3,080 between NI-NI and I-NI samples; and 83 and 275 between I-NI and I-I samples in the S and the R genotype, respectively, with several up-, down-and contra-regulated transcripts (Fig. 3A1,B1). When using the reference based assembly mapping, 880 and 1,391 differentially expressed transcripts were identified between NI-NI and I-I samples; 755 and 1,050 between NI-NI and I-NI samples; and 95 and 210 between I-NI and I-I samples in the S and the R genotype, respectively, with several up-, down-and contra-regulated transcripts (Fig. 3A2,B2). In addition, heat maps were generated for each genotype based on the differential expression data from edgeR in order to determine the sample relationships (Fig. 4). A clear separation was seen between non-incubated (NI) and incubated (I) samples in both genotypes, whereas incubated inoculated (I-I) and incubated non-inoculated (I-NI) grouped together (Fig. 4). Even the expression data generated from reference based mapping clearly differentiated between incubated and non-incubated samples. Both S and R incubated grouped together and were separated from non-incubated samples (Fig. 4).

Annotation and GO of differentially expressed transcripts.
Approximately 75 percent of the differentially expressed transcripts had blast hits to the Viridiplantae database extracted from NCBI. The top hit species are Brachypodium distachyon followed by Hordeum vulgare, which are most closely related to L. perenne. Among the transcripts with blast hits, 40-52 percent of the differentially expressed transcripts were annotated using Blast2GO 23 . Putative descriptions and functions were assigned to the transcripts predominantly based on annotations from H. vulgare and B. distachyon and Arabidopsis thaliana. Gene Ontology classifications of DEGs at I-NI vs. I-I conditions in genotypes R and S were generated using WEGO 24 . The results are summarized in three main GO categories: cellular component, molecular function and biological process (Fig. 5). Comparisons of the functional categories of genotype R with those of genotype S reveal differences in terms of the biological processes. DEGs responses to stress, biotic stimulus were highly represented in the R genotype, while DEGs response to death is only seen in S genotype Fisher's exact test from Blast2GO 23 was used for GO enrichment analysis between R and S to determine if any gene ontology (GO) terms were over-or under-represented in the various sets of differentially expressed transcripts. A total of seven GO terms were enriched when comparing the differentially expressed transcript sets    Table S4). Out of these, five were overrepresented in the R genotype, in terms related to cell wall cellulose metabolic process, cell wall pectin metabolic process, cell morphogenesis, actin nucleation and organelle epidermal cell differentiation. Transcripts assigned to aryl-alcohol dehydrogenase (NADP+ ) activity and phycobilisome were present only in R genotype.
Several genes involved in the initiation of pathogen-associated molecular pattern (PAMP) immunity, like cysteine-rich receptor-like protein kinase (CRK), cyclic nucleotide gated channel (CNGC), calcium-dependent protein kinase (CDPK), respiratory burst oxidase homolog (Rboh), calcium-binding protein CML (CaM/CML), and NADPH oxidase were detected in these studies. Several pathogen related genes like PR1, β-1,3-Glucanase (PR2), chitinase II/V (PR3), thaumatin-like (PR5), and lipid-transfer protein (PR14) are upregulated in genotype R compared with genotype S under I-I conditions ( Table 3, Supplementary Table S5). We also found several potential pathogen resistance candidate genes like chitinase V, lipid transfer protein, serine-glyoxylate aminotransferase and WRKY 75 highly upregulated in the R genotype under I-NI treatment compared with the I-I treatment ( Table 3). All potential candidate genes involved in the response of L. perenne to inoculation with M. nivale are listed in Table 3 with homologues in A. thaliana and B. distachyon. A hypothetical model for gene regulation in the plant-pathogen interaction pathway after four days of incubation with the pink snow mould pathogen, based on the pathogen related DEGs identified in this study, is presented in Fig. 7.
Furthermore, the KEGG database (http://www.genome.jp/kegg/) was used to detect different pathways in response to M. nivale in the S and R genotypes. Blast to the KEGG database showed that 5009 DEGs were involved in 135 pathways (Supplementary Table S2). Pathways with highest representation among the genes were involved in purine metabolism (5.19 percent, 260 genes), biosynthesis of antibiotics (5.09 percent, 255 genes), thiamine metabolism (4.25 percent, 213 genes), starch and sucrose metabolism (2.91 percent, 146 genes) and aminobenzoate degradation (2.61 percent, 131 genes).

Validation of transcripts by real-time PCR.
In order to validate the expression profiling by Illumina sequencing, the expression levels of six genes, including two chitinase genes, three pathogen-related genes and one WRKY family gene were further analysed by qRT-PCR. All the genes showed differential expression levels between S and R genotypes. Correlation analysis was performed between the RNA-seq and qPCR log transformed data for each gene. Among the six genes, four genes (WRKY, CHI5, PR1 and THI) were highly correlated (r 2 value in range of 0.82 to 0.97), while two genes (CHI2, PR5) were not well correlated (r 2 values are 0.52 and 0.39) ( Supplementary Fig. S3).

Discussion
A significant positive correlation was found between the amount of M. nivale DNA in leaf samples and visual assessment of disease severity after 6 weeks of incubation under artificial snow cover. After 8 weeks of incubation, disease severity was similar across genotypes, while significant differences in the quantity of M. nivale DNA were detected. The plants incubated for a longer period (8 weeks) had a higher amount of fungal DNA than the plants incubated for 6 weeks, which is in agreement with other studies reporting that longer incubation period increases snow mould infestation even in resistant genotypes 25,26 . In general, neither visual assessment of disease severity, nor the amount of M. nivale DNA in leaves were good indicators of snow mould resistance in L. perenne. In the present study, plant regrowth, after inoculation with M. nivale and incubation for several weeks, was not correlated with fungal biomass nor disease severity. Some genotypes such as M and K showed severe symptoms on their leaf tissues, but still had good regrowth, possibly because the lower stem was not infected. On the other hand, genotype C had a poor regrowth despite limited symptoms and M. nivale DNA detected in the leaves (Fig. 2, Supplementary Table S1).  A qPCR test could in principle be a useful method for breeders in the selection of snow mould resistant materials. For such a test to be useful, it should be based on infestation of the lower stem tissue of the plants. More importantly, the application of quantitative real-time PCR will facilitate the detection of latent infections and early diagnosis of disease. For such purposes the test described in this study would need thorough validation with respect to sensitivity and specificity.
High throughput sequencing capabilities have made the process of assembling a transcriptome easier, even for non-model organisms without a reference genome. But the quality of a transcriptome assembly must be good enough to capture the most comprehensive catalogue of transcripts and their variations, and to carry out further transcriptomic experiments 27 . The CEGMA analysis (Table 2) showed high coverage of ultra-conserved CEGs in the assemblies of the S and R genotypes, demonstrating their completeness in terms of gene content. However, a common question is whether reference based assembly gives better results than a de novo based assembly.
In this study we detected a larger number of differentially expressed transcripts in a pairwise comparisons between NI-NI and I-I; NI-NI and I-NI conditions, than the pairwise comparison between I-NI and I-I condition both by de novo and reference based assembly mapping (Fig. 3). It was expected that there would be a larger number of transcripts differentially expressed when plants were transferred from growth (non-incubation) conditions at 20-22 °C and 18 hours of light, to incubation conditions at 2 °C and darkness due to the significant changes in temperature and light. This was seen both for the susceptible (S) and the resistant (R) genotypes, as several differentially expressed transcripts involved in rapid responses to cold stress and light, in addition to abiotic  stress related transcripts were detected. On the other hand, few differentially expressed transcripts were observed between treatments I-NI and I-I in both genotypes. The annotation results of the detected transcripts between I-NI and I-I conditions in both de novo and reference based mapping identified similar genes involved in biotic stress, immune response and cell death. This shows the potential of de novo method in capturing the essential transcripts even in the absence of reference genome, which also has been demonstrated in raspberry studies 28 . To our knowledge, this is the first transcriptome study using RNA-seq to understand the response of L. perenne to the early infection of pink snow mould (M. nivale). Several of the differentially regulated genes such as disease related proteins, calmodulin binding proteins, lipid transfer proteins, and flavonoid biosynthesis (Table 3,  Supplementary Table S5) detected in the R and S genotypes between the I-NI and I-I conditions are involved in different defense-response mechanisms. The R genotype showed higher expression levels of several pathogenesis related genes such as PR1, PR2, PR3, PR5, PR13 and PR14. These results are similar to those from winter wheat, where snow mould resistance was associated with the accumulation of PR1a, PR2, PR5, and PR14 29 . The higher expression levels of these PR-proteins are often considered as markers for activation of the salicylic acid (SA) signalling pathway 25

. Hypothetical modules for plant-pathogen interaction after 4 days of incubation with snow mould pathogen M. nivale derived by KEGG plant-pathogen interaction pathway (http://www.genome.jp/kegg/) and network of WRKY transcription factors (Eulgem & Somssich (2007). Red color indicates down regulated genes and green color indicates up-regulated genes. Intensity of the colors indicates the fold change.
Circle represents the fold change under non-inoculated and incubated, condition (I-NI) a . The square represents the fold change under inoculated, incubated (I-I) b conditions. The recognition of pathogen-associated molecular pattern (PAMP) initiate PAMP trigger immunity via the activation of cysteine-rich receptor-like protein kinase (CRK), cyclic nucleotide gated channel (CNGC), calcium-dependent protein kinase (CDPK), respiratory burst oxidase homolog (Rboh), calcium-binding protein CML (CaM/CML), and NADPH oxidase. The activation of PAMP trigger immunity initiate the production of reactive oxygen species (ROS), which might activate the plant hypersensitive response (HR), cell wall reinforcement, as well as stomata closure. Defense responses are also instigated upon recognition of the fungal effectors in the host cell by serine/threonine-protein kinase PBS (PBS) and the activation of MAP kinase cascades such as mitogen-activated protein kinase kinase kinase (MAPKKK), mitogen-activated protein kinase kinase (MAPKK), and mitogen-activated protein kinase (MAPK). Effectors triggered immunity (ETI) initiate the production of several pathogen related proteins such as PR-1, β -1,3-glucanase (PR-2), chitinase II/V (PR-3), thaumatin-like (PR-5), and lipid-transfer protein (PR-14). Both PAMP triggered immunity and effectors triggered immunity alternate the production of salicylic acid (SA) and jasmonic acid (JA) by the action of distinct transcription factors WRKY such as WRKY 75, WRKY 70, WRKY 18, and WRKY 33. Pathogen-triggered SA signaling also by the activation of serine/threonine-protein kinase2 (SRK2), auxin receptors, and abscisic acid responsive element binding factor (ABF). Conversely, the salicylic acid induction-deficient mutants of Arabidopsis expressed PR2 and PR5 and accumulated high levels of camalexin after pathogen inoculation 32 .
WRKY proteins, another important defense related group of proteins, constitutes a superfamily of transcription factors, involved in the regulation of different physiological platforms in plants, including pathogen defense, trichome development and senescence. In this study, WRKY65, WRKY70 and WRKY75 were upregulated after inoculation with snow mould (Table 3). It is also reported that WRKY transcription factors contributed to the defense against Pseudomonas syringae in tomato and play a partially conserved role in basal defense in tomato and Arabidopsis 33 . The plant immunity system consists of two main levels 34 . The first level is based on the perception of pathogen-associated molecular patterns (PAMPs), which activates the PAMP-triggered immunity pathway (PTI). The second level is the recognition of pathogen effectors, which activates pathogen related PR genes in a process called effector-triggered immunity (ETI). In the present study, the transcriptome analysis of the snow mould resistant genotype showed that the PTI pathway was activated (Fig. 7), particularly by the up-regulation of the expression level of calcium-dependent protein kinase CDPK, respiratory burst oxidase homolog Rboh and calcium-binding protein CaM/CML. Therefore, the activation of the PTI inhibits the snow mould pathogen from colonizing the plant tissues by increasing the production of reactive oxygen species and cell wall reinforcement. These results are similar to studies in Festulolium, where the resistant genotypes are characterized by high peroxidase activity, intensive lignification, callus formation and high concentrations of reactive oxygen species during the stage of early infection (within 6 days from inoculation) 31 .
Transcription factors, such as WRKY, play important roles in defense responses towards several plant pathogens 35 . The transcription level of WRKY genes are up-regulated by several stress factors, in particular pathogen infection 36 38 . In the present study, the resistant genotype showed high transcription levels of several WRKY genes such as WRKY 70 and WRKY 75. Therefore, it is expected that the up-regulation of these genes will lead to the activation of the salicylic acid pathway 37 . Furthermore, our results also showed down-regulation of WRKY 18 and WRKY 33, which are responsible for the activation of the JA pathway and the deactivation of the SA pathway 36,39 . In a study by Gaudet et al. 25 , the expression levels of WRKY 34 and WRKY 16 were up-regulated in snow mould resistant genotypes of winter wheat, which led to the activation of the JA pathway. Other studies showed that M. nivale infection is usually influenced by the physical and the chemical conditions of the plant tissue, thus the fungus behaves as a biotroph when the plant defense system is induced and the SA pathway is activated 9,10 .
The cross talk between cell morphogenesis and plant-pathogen interactions plays a crucial role in disease development 40 . Plants have developed a system for sensing pathogens and monitoring the cell wall integrity, upon which they activate defense responses that lead to a dynamic cell wall remodelling required to prevent disease 40 . Genes responsible for actin nucleation, aryl-alcohol dehydrogenase (NADP+ ) activity and cell differentiation were significantly enriched in the R genotype based on GO enrichment analysis by Fisher's exact test (Fig. 6,  Supplementary Table S4). Under GO term actin nucleation, we found that genes encoding actin-related protein 2 (ARP2), importin-β and serine threonine-protein kinase (TOR) were over-represented in the R genotype during infection. ARP2 in complex with ARP3 plays a central role in actin cytoskeletal formation 41 , and genetic experiments have indicated a role for this complex in the early stages of low temperature signalling 42 and in the response to salt stress 43 . The importin-β subunit belongs to nuclear import receptors which play an essential roles in transferring defense proteins, such as nucleotide-binding and leucine-rich repeats (NB-LRRs), from the cytoplasm to the nucleus 44 , where they initiate defense signalling 45 . Moreover, protein kinases such as serine threonine-protein kinase (TOR) play a key role in signalling during pathogen recognition and subsequent activation of plant defense mechanisms 46 .
Under GO term aryl-alcohol dehydrogenase activity, a gene encoding a voltage-gated potassium channel beta subunit like was over-represented in the R genotype. Potassium (K + ) plays many important regulatory roles in plant development and stress responses 47 . High K + status decreases the occurrence of many diseases 48 . Furthermore, K + affects plant hormonal pathways, i.e. the salicylic acid (SA) and jasmonic acid (JA) pathways 48 , which are involved in hypersensitive responses or acquired systemic resistance to pathogens. Recent studies 49 showed that overexpression of GmAKT2 encoding a K+ transporter significantly increased K+ concentrations and consequently resistance to soybean mosaic virus in transgenic soybean 49 .
Under cell differentiation gene ontology, we detected that the genes encoding glycerol-3-phospahate-1 (G3P) transporter and prefoldin (PFD) were over-represented in the snow mold infected R genotype. The G3P transporter is an important component of carbohydrate and lipid metabolic processes. G3P levels in A. thaliana plants were previously associated with defense to the hemibiotrophic fungal pathogen Colletotrichum higginsianum 50 . Infection of A. thaliana with C. higginsianum showed an increase in G3P levels and a concomitant enhanced resistance in the host 50 . Prefoldin proteins are required for the cytoplasmic folding of actin and tubulin monomers during cytoskeleton assembly 51 . Recent studies 52 showed that prefoldin 6 interacts with two P. syringae effectors and defense regulatory protein EDS1 (enhanced disease susceptibility 1). Additionally, prefoldins 3 and 5 have been shown to play essential roles in tolerance to salt stress in Arabidopsis 53 . Significant enrichment of these GO terms in the R genotype show that these gene systems are involved in defense responses to pink snow mould infection in perennial ryegrass.

Conclusions
In this study, a susceptible (S) and a resistant (R) genotype of L. perenne cv. Fagerlin were found to be significantly different in resistance, as measured by relative regrowth, and accumulation of M. nivale DNA, as quantified by real-time qPCR. RNA sequencing of transcriptome responses of non-cold acclimated plants to early infection by Scientific RepoRts | 6:28702 | DOI: 10.1038/srep28702 an aggressive M. nivale isolate identified differentially expressed genes between the S and R genotype. Many pathogen related genes were found to be upregulated during snow mould infection in the resistant genotype. Further GO enrichment analysis confirmed that specific GO terms related to plant defense are over-represented in the resistant genotype. The list of putative candidate defense associated genes and cell morphogenesis associated genes identified in this study might provide a scientific basis for further investigations to obtain more in-depth understanding of host-pathogen interactions and development of resistant cultivars by marker-assisted breeding.

Materials and Methods
Plant materials and growth conditions. For screening of snow mould resistance and quantification of fungal DNA, mother plants of eight randomly selected genotypes (A-F, M and K) of perennial ryegrass, cultivar Fagerlin, were divided into multiple ramets and planted in a fertilized soil mixture (Gartnerjord, TJERBO) and grown in the greenhouse at 20-22 °C (day/night) and 18 hours light period with light intensity (CONSTANTCOLOR CMH lamps 400W) at about 200-250 μ mol m −2 s −1 . The plants were fertilized weekly with a mixture of 80 g/L KRISTALON fertilizer 9-5-25 (N-P-K) and 60 g/L of YARALIVA CALCINIT 15.5-0-0 (Yara International ASA, Oslo, Norway), diluted to a conductivity of 2 mS/cm. For the transcriptome studies, plants from the eight genotypes were randomly selected and placed on four trolleys (100 × 60 cm) and moved to controlled growth chamber at 18/20 °C (day/night) temperature, with light intensity of 220-240 μ mol m −2 s −1 for four weeks. Perennial ryegrass is often infected with endophyte Epichloë festucae var. lolii. However, the cultivar Fagerlin used in this study, is not infected with endophyte Epichloë festucae var. lolii (pers. communication with the forage grass breeder Petter Marum, Graminor AS, the owner of the variety). Further to be sure, we performed PCR to screen the eight genotypes of Fagerlin for endophytes with the β -tubulin (TUB2) 54 , translation elongation factor 1-a (TefA) 55 and the soft (SO) 55 genes (which are essential in screening for endophyte infection). We did not detect any traces of these genes in our genotypes, indicating that Fagerlin is free of endophyte infection.  57 . Briefly, the fungus was incubated at 9 °C, in darkness, on potato dextrose agar (PDA) for two weeks. An Erlenmeyer flask containing 100 ml potato dextrose broth (PDB) was then inoculated with four plugs (5mm diameter) of fungal mycelium and incubated at 15 °C in darkness. Fungal mycelium was harvested after 10 days by filtering through cheesecloth. The mycelium was homogenized in distilled water containing 0.01% TWEEN 20 (SIGMA) using an ULTRA TURRAX. The inoculum was diluted to an optical destiny of 0.5 at 430 nm. Plants were inoculated by spraying (1 ml inoculum per plant on average) and the control plants were sprayed with distilled water. After inoculation, the plants were incubated under artificial snow cover by covering the plants with moist cellulose tissue paper and black plastic sheets, then placed in a cold chamber at 2 °C in darkness for six and eight weeks. Each week during incubation, the plants were repositioned in the room.
After incubation, the plants were moved to a greenhouse at 20-22 °C and 18 hours of light for recovery. The plants were cut at five cm above the soil surface and allowed to regrow for two weeks. The regrown plants were harvested (all parts above soil surface) and dried at 60 °C for three days in order to measure dry weight (g DW plant −1 ). Relative regrowth was calculated for each inoculated plant as the dry weight divided by the average dry weight of non-inoculated plants within the same genotype. Relative regrowth values approaching 1 represents resistant plants. Disease severity was visually scored two days later according to the following scale: 0 = no green tillers, 1 = some green tillers visible, 2 = green tillers found in less than half of the total plant area, 3 = green tillers found in more than half of the plant area, and 4 = green tillers observed in the whole plant area. After visual assessment of the symptoms, plant leaves and stems were harvested (5 cm above soil surface) and kept at − 20 °C for DNA extraction for fungal quantification.
Real-time PCR for fungal quantification. Plant materials (leaves and stems above five cm from soil) were collected from the eight genotypes at two different time points (6 and 8 weeks after inoculation), from the genotypes F and M also at 1 and 4 days after inoculation. For the samples collected 6 and 8 week after inoculation, plant materials were stored at − 20 °C until DNA was extracted. For samples collected after 1 and 4 days, the same plant materials used for gene expression analysis (see section below) were utilized for fungal DNA quantification. For DNA extractions, samples were frozen quickly in liquid N 2 and ground using mortar and pestle. DNA was extracted from the ground plant tissue using DNeasy Plant Mini Kit (QIAGEN) according to manufacturer's protocol (Qiagen Inc., Germany). The quality of the extracts was measured using a NANODROP ND-1000 UV-Vis Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and visualized by electrophoresis through 1.5% agarose gels.
Real-time PCR primers specific for M. nivale were designed based on the EF-1a gene sequence by Glynn et al. 15 using PRIMER EXPRESS software version 2 (Applied Biosystems, Foster City, USA) based on the following parameters; amplicon Length of 50 to 150 bases (for optimum PCR efficiency), primer length of 20 bases, melting temperature (Tm) of 58 °C to 60 °C (Optimal 59 °C), G + C content being between 30 and 80%; and the last five nucleotides at the 3′ end do not contain more than two G + C residues. The following primer set was chosen based on the regions of identity within the EF-1a gene sequence between isolates, forward primer EF1-F: 5´-GGTCTTGGCTTGCACAAACA-3´ and reverse primer EF1-R: 5´-AGCACAACAGGCGTGGATAAG -3´. Quantification of plant and fungal DNA was carried out by real time PCR in a total volume of 25 μ l, using 2x SYBR Green PCR master mix (Applied Biosystems), 300 nM of each primers (Invitrogen Ltd, UK) and 2 μ l of 10x diluted template DNA. Specific real time PCR primers for the plant housekeeping gene LpGAPDH 58 were used as internal control for plant DNA. PCR was performed on Applied Biosystems 7900HT instrument with a standard 96-well block (Applied Biosystems). For all the PCR reactions the following cycling parameters were used: 50 °C for 2 min, 95 °C for 10 min, 40 cycles of 95 °C for 15 s and 60 °C for 1 min followed by dissociation curve analysis at 60 °C-95 °C. The data was analysed using SEQUENCE DETECTION SOFTWARE (SDS) Version 2.2.1 (Applied Biosystems). The amount of fungal and plant DNA in the samples were quantified by a standard curve algorithm based on cycle threshold value (Ct) using a 10-fold dilution series of known amount of DNA, starting with 5 ng for fungal DNA and 100 ng for plant DNA and three technical replicates. Samples were tested in two technical replicates. The amount of fungal DNA was calculated as pg fungal DNA per μ g plant DNA for each sample.
Tissue sampling and RNA extraction for gene expression analysis. Leaf samples were collected from genotype F, a susceptible genotype (hereafter termed S) and genotype M, a resistance genotype (actually less susceptible, here after termed R), which had been exposed to three different treatments: 1: Non-inoculated and non-incubated plants (NI-NI) (Control plants kept in ambient temperature), 2: Non-inoculated plants that were incubated under artificial snow-cover for four days (I-NI) , and 3: Inoculated plants that were incubated under artificial snow-cover for four days (I-I), with two biological replicates, a total number of 12 samples. The collected samples were immediately placed in liquid nitrogen and stored at − 80 °C until used for RNA extraction. The frozen leaf samples were crushed with a pestle and mortar and total RNA was extracted using the PURE LINK RNA MINI KIT (Life technologies, USA) and PLANT RNA ISOLATION AID (Life technologies, USA). On-Column DNAse kit was used to remove the DNA contamination. The concentration and quality was checked using the NANODROP ( De novo and reference based transcriptome analysis. After trimming adapter sequences and filtering low quality reads using the sickle program (https://github.com/najoshi/sickle/blob/master/README.md), the bioinformatics pipeline ( Supplementary Fig. S4) was followed for de novo assembly and further detection of differentially expressed genes. Briefly, the clean reads derived from the two genotypes susceptible (S) and resistant (R) were used to construct separate de novo assemblies for each genotype using the Trinity assembler (release 2013-02-25) 59 with the following settings; Trinity.pl --seqType fq --JM 20G --left F_1_Lolium.fq-QT.gz --right F_2_Lolium.fq-QT.gz --CPU 16 -min_contig_length 200 --SS_lib_type FR --full_cleanup --min_kmer_cov 2 --output Trinity_201 2 > &1 > logfile.lolium-F. The de novo assembled transcripts were then used as a reference to map back the individual reads by Bowtie. Further, we estimated transcript abundances in each genotype and treatment combination using RSEM version 1.1.11 60 . A maximum of one mismatch (-bowtie-n 1) was allowed in the seed region of the reads. In another approach, to facilitate comparison of the two genotypes, we aligned all clean reads from each genotype and treatment combination to a reference transcriptome of an inbred L. perenne genotype, generated from a combination of root, stem, leaf sheath, leaf and meristem samples 22 , and estimated transcript abundance as described above.
De novo assembly validation by CEGMA. CEGMA software (version 2.4) 21 was used to assess the completeness of the S and R transcriptome assembly datasets. This program assesses the presence and coverage of a set of 248 extremely conserved core eukaryotic genes (CEGs). It is routinely used for evaluating genomic assemblies, however, it has also been used for evaluating transcriptome assemblies 61,62 . The software was run with default parameters with the included reference dataset of 248 ultra-conserved Core Eukaryotic Genes (CEGs).

Identification of differentially expressed genes, BLAST and functional annotation.
Gene expression levels were measured as expected number of fragments per kilobase of transcript sequence per millions mapped reads (FPKM) 63 . The transcript matrix files derived from RSEM program were processed with the edgeR program 64 using perl script (DE_analysis.pl) in the Trinity pipeline for detecting differentially expressed genes determined with a False Discovery Rate (FDR) of 0.05. Briefly, pairwise comparisons were carried out between all the selected time points and edgeR 64 analysis was performed by fitting normalized count data with a generalized linear model (GLM) estimating a negative binomial distribution to the calculated mean values of the two biologically independent samples. For each gene, fold changes and P values (pval) as well as P values adjusted (padj) for multiple testing with the Benjamini-Hochberg procedure 65 , were used to control FDR. The sequence with padj of less than 0.05 was deemed to be significantly differentially expressed genes (DEGs). The variance stabilized data obtained from edgeR 64 was used as input for clustering, and for constructing multidimensional scaling plots using R integrated in the Trinity pipeline. The transcripts showing differential expression at any time point during snow mould infection were clustered using a K-means clustering algorithm. VennPlex program 66 was used to generate Venn diagrams showing up-, down-and contra-regulated transripts.
The DEGs were annotated using Blast2GO 23 . An E-value threshold of 10 e-06 was used for the BLASTx search, and 10 e-10 for the annotation, with a cut-off value of 55 and a GO weight Hsp-hit value of 20. The enrichment analysis for the differential gene ontology term distribution was performed with a p-value significance cut-off value of 0.01. Gene ontology classifications of differentially expressed genes in the resistant (R) and susceptible (S) genotypes were generated using the web histogram tool WEGO 24 . Pathway analysis was performed using the KEGG function implemented in the Blast2GO 23 tool.
Validation of RNAseq expression profiles by qRT-PCR. Expression patterns of five defense related genes (chitinase 2, chitinase 5, WRKY, PR3, PR1 and PR5) differentially expressed between I-NI and I-I samples identified in this transcriptome studies were analysed using qRT-PCR. The same RNA used for sequencing was used for validating the genes by qRT-PCR. Based on the transcriptome sequences of the five genes, primers (Supplementary Table S3) were designed using primer express software version 2 (Applied Biosystems, Foster City, USA). Efficiency test of the primers was performed on different samples for normalization of the expression level. The EXPRESS two-Step qRT-PCR kit, which includes the SuperScript VILO cDNA Synthesis kit, was used for generating the single-stranded cDNA that was later used for quantifying the amount of specific gene expression using forward and reverse primers, following the manufacturer's instructions. cDNA synthesis was done using up to 2.5 μ g of the total RNA in 20 μ l reaction. Five μ l of cDNA (5x diluted) was used in each well of Fast Optical 96 well plate along with other components making the total volume 20 μ l. The fast cycling program was then set at 95 °C for 20 sec, 40 cycles of 95 °C for three sec (denaturation) and 60 °C for 35 sec (annealing). Then each plate (with samples) for each gene with a bar code was placed in ABI7500 qRT-PCR machine. The SYBR Green dye was used to detect the amplified products. The expressions of the specific genes were normalized by using LpGAPDH (EC 1.2.1.12) as the reference gene. The expression of the target gene relative to the reference gene at 4 days after inoculation using the 2 −ΔΔCT method where the Δ Δ CT = (CT of target -CT of reference) 4 days after inoculation -(CT of target -CT of reference) before inoculation, which gives the mean relative expression of target genes at this time point 67 .