|Home | About | Journals | Submit | Contact Us | Français|
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.
Perennial ryegrass (Lolium perenne L.,) belongs to the Poaceae family. It is a diploid species (2n=2x=14) native to Europe, Asia and Northern Africa1. It is an important forage grass in the temperate regions of the world because of its high forage quality and yield. Out of 52 million ha of grasslands in Europe, 23% is cultivated with Lolium species, with perennial ryegrass being the most widespread. Perennial ryegrass has low resistance against pink snow mould, however, tetraploid cultivars have better resistance than diploid and turf cultivars1.
Winter injury is regarded as a serious constraint for the production of winter cereals and grasses at northern latitudes2,3. The fungus Microdochium nivale (Fr.) Samuels & Hallet is considered to be the most widespread cause of biotic winter injury in these crops4. It is an opportunistic species causing pink snow mould on winter cereals, turf and forage grasses at low temperatures, with or without a snow cover. High humidity and constant low temperatures under snow cover are highly favourable for the development of pink snow mould5,6.
Resistance to pink snow mould is enhanced by cold acclimation2,3,4. During this process the plant undergoes numerous physiological and bio-chemical changes which are essential for winter survival7. Some of these changes are also thought to increase resistance to diseases, e.g. cellular dehydration and accumulation of defense-related proteins and fructans8. Genetic variation in cold-induced resistance to pink snow mould in triticale has been shown to be associated with changes in physical and chemical properties of the leaf surface and cell walls9, and with photosynthetic acclimation and peroxidase activity10.
Previous studies on pink snow mould resistance has almost exclusively been performed on cold acclimated plants. However, some genetic variation in resistance is also present in non-acclimated winter wheat11, and this resistance may be masked when testing cold acclimated plants. Inherent resistance that is independent of cold acclimation may be more specific to M. nivale than cold-induced disease resistance, and is likely to be important for preventing diseases caused by M. nivale during the growing season, such as microdochium patch and leaf blotch. It may also be increasingly important with the predicited climate change12. The projected climatic conditions12 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 genotypes13. 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 weight14. 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 species15. The EF-1a has also been used to study the genetic variation among isolates16. Additionally, competitive PCR methods have also been developed for M. nivale and M. majus quantification in infected tissues17.
The use of whole transcriptome sequencing (RNA-seq) is an important tool in the analysis of complex plant responses18, and provides a more comprehensive understanding of transcription initiation sites, improved detection of alternative splicing events and the detection of gene fusion transcripts19. The technology is able to handle de novo sequencing of large genomes, revealing individual genome differences within the same species and quantify gene expression20. 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.
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).
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).
A total of 178 million reads and 165 million reads of 100bp 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,784bp, and 188,355 contigs for the R genotype with N50 of 1,672bp (Table 1). The longest assembled contigs in the S and R genotype were 17,632 and 12,882bp, 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 pipeline21 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.
The expression levels of each assembled transcript were estimated at three different treatments, i.e. (1: Non-inoculated and non-incubated plants (NI-NI) (control plants kept in ambient temperature), 2: Non-inoculated plants 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) (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 transcriptome22. 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).
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 Blast2GO23. 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 WEGO24. 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 Blast2GO23 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 from I-NI vs. I-I conditions of the two genotypes (Fig. 6, Supplementary 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).
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 (r2 value in range of 0.82 to 0.97), while two genes (CHI2, PR5) were not well correlated (r2 values are 0.52 and 0.39) (Supplementary Fig. S3).
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 genotypes25,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 experiments27. 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 18hours 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 studies28.
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 PR1429. The higher expression levels of these PR-proteins are often considered as markers for activation of the salicylic acid (SA) signalling pathway25,29,30,31. Pociecha et al.31 also demonstrated that resistant genotypes of Festulolium are characterized by high SA concentrations during snow mould infection. PR proteins seem to be part of a larger set of SA and jasmonic acid (JA)-dependent defense responses in which each PR protein may contribute differently to the snow mould infection. For instance the A. thaliana mutant ein2, which is defective in ethylene/JA signalling, showed low expression level of PR12, PR3 and PR4 and high susceptibility to B. cinerea (necrotrophic pathogens). Conversely, the salicylic acid induction–deficient mutants of Arabidopsis expressed PR2 and PR5 and accumulated high levels of camalexin after pathogen inoculation32.
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 Arabidopsis33.
The plant immunity system consists of two main levels34. 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 pathogens35. The transcription level of WRKY genes are up-regulated by several stress factors, in particular pathogen infection36. In Arabidopsis, 49 out of 72 WRKY genes tested responded to bacterial infection or salicylic acid37, and 8 Arabidopsis WRKY genes (WRKY 18, WRKY 38, WRKY 53, WRKY 54, WRKY 58, WRKY 59, WRKY 66, and WRKY 70) were characterized as direct targets of NPR1, a key regulator of SA signalling38. 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 pathway37. 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 pathway36,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 activated9,10.
The cross talk between cell morphogenesis and plant-pathogen interactions plays a crucial role in disease development40. 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 disease40. 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 formation41, and genetic experiments have indicated a role for this complex in the early stages of low temperature signalling42 and in the response to salt stress43. 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 nucleus44, where they initiate defense signalling45. 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 mechanisms46.
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 responses47. High K+ status decreases the occurrence of many diseases48. Furthermore, K+ affects plant hormonal pathways, i.e. the salicylic acid (SA) and jasmonic acid (JA) pathways48, which are involved in hypersensitive responses or acquired systemic resistance to pathogens. Recent studies49 showed that overexpression of GmAKT2 encoding a K+ transporter significantly increased K+ concentrations and consequently resistance to soybean mosaic virus in transgenic soybean49.
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 higginsianum50. Infection of A. thaliana with C. higginsianum showed an increase in G3P levels and a concomitant enhanced resistance in the host50. Prefoldin proteins are required for the cytoplasmic folding of actin and tubulin monomers during cytoskeleton assembly51. Recent studies52 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 Arabidopsis53. 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.
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 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.
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 18hours 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 80g/L KRISTALON fertilizer 9-5-25 (N-P-K) and 60g/L of YARALIVA CALCINIT 15.5-0-0 (Yara International ASA, Oslo, Norway), diluted to a conductivity of 2mS/cm. For the transcriptome studies, plants from the eight genotypes were randomly selected and placed on four trolleys (100×60cm) 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.
Isolate 200231 of M. nivale (isolated from L. perenne at Ås, Norway (59°N) was obtained from the fungal culture collection at the Norwegian Institute of Bioeconomy Research (NIBIO), Ås, Norway. The inoculum was prepared according to Tronsmo56 and Hofgaard et al.57. Briefly, the fungus was incubated at 9°C, in darkness, on potato dextrose agar (PDA) for two weeks. An Erlenmeyer flask containing 100ml 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 430nm. Plants were inoculated by spraying (1ml 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 18hours 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 (5cm above soil surface) and kept at −20°C for DNA extraction 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 N2 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), 300nM of each primers (Invitrogen Ltd, UK) and 2μl of 10x diluted template DNA. Specific real time PCR primers for the plant housekeeping gene LpGAPDH58 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 2min, 95°C for 10min, 40 cycles of 95°C for 15s and 60°C for 1min 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 5ng for fungal DNA and 100ng 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.
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 (Nanodrop Technologies, Wilmington, DE, USA) and BIOANALYZER (Agilent Technologies, Palo Alto, CA, USA) equipment.
Twelve RNA samples with RIN (RNA Integrity Number) values above seven were used to construct separate cDNA libraries with fragment lengths of 200bp (±25bp). Then, paired-end sequencing was performed by GATC Biotech Ltd., Germany (http://www.gatc-biotech.com/en/index.html) using the Illumina sequencing platform (HISEQ 2000). Real time analysis (RTA) output was analysed using the CASAVA software (version 1.6, Illumina), generating pass filtering FastQ files with Qphred +64 quality values. Paired-end reads with a length of 100bp were generated. The reads were deposited in the EMBL-EBI ArrayExpress Archive, under accession number E-MTAB-4459. The quality of the reads was analysed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).
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.1160. 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 samples22, and estimated transcript abundance as described above.
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 assemblies61,62. The software was run with default parameters with the included reference dataset of 248 ultra-conserved Core Eukaryotic Genes (CEGs).
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 program64 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 edgeR64 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 procedure65, 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 edgeR64 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 program66 was used to generate Venn diagrams showing up-, down- and contra-regulated transripts.
The DEGs were annotated using Blast2GO23. An E-value threshold of 10e-06 was used for the BLASTx search, and 10e-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 WEGO24. Pathway analysis was performed using the KEGG function implemented in the Blast2GO23 tool.
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 20sec, 40 cycles of 95°C for three sec (denaturation) and 60°C for 35sec (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 22.214.171.124) 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 point67.
Accession codes: The raw data generated in this study were deposited in the EMBL-EBI ArrayExpress Archive, under accession number (E-MTAB-4459). Snow mould resistant (R) and susceptible (S) de novo transcriptome assemblies generated by trinity program are deposited in DRYAD Digital Repository along with the sequences (ESTs) of candidate genes detected for snow mould resistance and sequences for the genes under GO terms significantly over represented in resistant genotype in this study (http://datadryad.org/resource/doi:10.5061/dryad.mc7c1).
How to cite this article: Kovi, M. R. et al. Global transcriptome changes in perennial ryegrass during early infection by pink snow mould. Sci. Rep. 6, 28702; doi: 10.1038/srep28702 (2016).
This work was funded by the Research Council of Norway (NFR) project “VARCLIM: Understanding the genetic and physiological basis for adaptation of Norwegian perennial forage crops to future climates” (project number 199664). Mohamed Abdelhalim was funded by a PhD scholarship from the Norwegian University of Life Sciences. The authors kindly thank Øyvind Jørgensen for excellent technical support in handling the plant material.
Author Contributions Å.E., M.B.B., I.S.H., A.M.T. and O.A.R. conceived and designed the experiments. M.A. and M.B.B. performed fungal quantification and validation of transcripts by qPCR. A.K. and M.A. performed snow mould resistance test. A.K. performed RNA extractions. K.M.R. was responsible for RNA sequencing and expression analysis with inputs from A.K., T.A. and M.A. K.M.R. drafted the manuscript with inputs from Å.E., A.M.T., M.A. and O.A.R. All authors read and approved the final manuscript.