Powdery mildew is one of the major fungal diseases of strawberry worldwide and causes significant economic loss because of reduced yield and berry quality. The application of sulphur and systemic fungicides to combat powdery mildew is expensive and hazardous to the environment. For this reason, the development of cultivars with enhanced resistance is the favoured strategy. However, traditional breeding of cultivated strawberry, an octoploid species, is inefficient and lengthy due to difficulties in introgression of resistance sources [1,2]. Significantly, the resistance levels of the diploid Fragaria species have not been studied in detail. Therefore, understanding the molecular mechanisms of host–pathogen interactions is of primary importance in devising strategies to control diseases. Powdery mildew of strawberry is caused by the obligate biotrophic fungus Podosphaera aphanis (syn: Sphaerotheca macularis f. sp. fragariae) [3,4]. The pathogen can infect all above ground parts including leaves, petioles, runners, flowers, fruits and their supporting stems [5]. The conidia produced on infected plants are wind dispersed over long distances. Under favourable conditions conidia germinate 4-6 h after infection to produce a germ tube that develops an appressorium after 12 h [6] which forms a penetration peg to breach the host cell wall by means of both enzymatic and mechanical force [7]. After successful penetration, a haustorium is developed 1 d after infection in the host epidermal cells by invading the host plasma membrane, whilst keeping the host cell still intact. However, the host is able to activate cell-wall associated defence to prevent fungal entry by suppressing the haustorial intake of nutrients. If the pathogen can avoid these defences and grow successfully, conidiophores are formed after 3-5 d [8,9] and conidiation starts 6 d after infection. For crop improvement programs, a better understanding of the molecular mechanisms of disease resistance is a prerequisite. Resistance to pathogens is classified into basal resistance and R-gene mediated resistance. Plant defence responses are elicited by the recognition of Pathogen Associated Molecular Patterns (PAMPs) which can be any part of the pathogen containing cell-wall proteins or flagella, toxins, etc. Plants recognize PAMPs by using extracellular receptor-like kinases (RLKs) and rapidly activate a signalling cascade through MAP kinases that result in basal immunity. R genes primarily encode the NBS-LRR class of proteins [10] and these NBS-LRRs directly recognize the Avr gene product (effectors) and initiate various downstream responses that include induction of PR genes, accumulation of inhibitory secondary metabolites and production of reactive oxygen species which lead to the hypersensitive response (HR), a form of programmed cell death, and induce disease resistance. These resistance responses still need to be characterised at the transcriptional, protein and cellular structure levels. Genetic analysis of immune response in model plants (eg Arabidopsis thaliana), plant transcriptome profiling, and whole genome sequencing have significantly increased the understanding of molecular mechanisms of the host-pathogen interaction [11]. Recent advances in high-throughput sequencing technologies aids in comprehensive understanding of host-pathogen interactions. The direct high throughput sequencing of cDNA (RNA-seq) using next-generation sequencing (NGS) technologies is progressively being used for global gene expression profiling. Highly specific, sensitive, cost effective, greater depth of sequencing and accurate quantitative measurements allow the NGS technologies to overcome the limitations of traditional hybridization or sequence-based approaches. Recently, RNA-seq was used to study patterns of gene expression during plant-pathogen interactions in cotton [12], soybean [13], rice [14] and Arabidopsis [15]. In the present study, two diploid strawberry cultivars were used, i.e. Fragaria vesca ssp. vesca accession Hawaii 4 (HW) and a highly inbred line, “Yellow Wonder 5AF7 (YW5AF7)” F. vesca f. semperflorens (YW). Despite recent availability of the genome sequence of HW [16], its resistance level to powdery mildew is not known, whereas, YW is reported to be susceptible to powdery mildew [17]. In this work, RNAseq was used to generate and compare large transcriptome datasets among two F. vesca cultivars upon P. aphanis infection. To the best of our knowledge RNA-seq has never been used to study plant-pathogen interactions in either cultivated or diploid strawberry species. Our study therefore provides the first broad overview of the time-course of transcriptional changes associated with P. aphanis infection in F. vesca.

Plant material and powdery mildew infection

Seed derived plants of two diploid strawberry cultivars, F. vesca ssp. vesca accession Hawaii 4 (originally provided by East Malling Research, UK) and F. vesca f. semperflorens line “Yellow Wonder 5AF7” (seed provided by Dr Janet Slovin, Beltsville, USA) were grown in hydroponic medium as previously described by Bindschedler et al. [18] and maintained in a Fisons growth chamber at a temperature of 20°C with 65% humidity, and a 10 h photoperiod with a photosynthetic photon flux density of 250-260 μmol/m²/s. Leaf material from an octoploid strawberry (cv. Eluica) infected with powdery mildew was provided by East Malling Research, UK. The four week old F. vesca plants were infected with P. aphanis by gently tapping conidia from infected leaves above the experimental plants. For microscopy, leaves were collected 1 h, 1 d (1 DAI) and 8 d (8 DAI) after infection, stained with 8-Anilino-Naphthalene-Sulfonic Acid (ANSA), and observed under an epifluorescence microscope (Leitz Dialux 20). For Illumina sequencing, leaves (100 mg in total) from four randomly selected plants were collected for each sample, prior to infection and then 1 DAI and 8 DAI. Leaf material was immediately frozen in liquid nitrogen and stored at -80°C until use.

RNA isolation and illumina sequencing

Total RNA was isolated from the collected frozen leaf samples using an RNAqueous® Kit (Ambion, UK) according to the manufacturer’s instruction. Six RNA samples, including three samples isolated from HW (control, 1 DAI and 8 DAI) and three samples from YW (control, 1 DAI and 8 DAI) were sent to Source BioScience (Nottingham, UK) for library preparation and sequencing. Briefly, one microgram of total RNA was subjected to library preparation using Illumina TruSeq RNA sample preparation kit guide v2 (September, 2012) following the manufacturer’s instruction. Libraries were amplified by 15 cycles of PCR and then purified and size selected for an average size of 350 bp. RNA quality and libraries were validated using the Agilent Bioanalyzer 2100 to confirm the concentrations and size distribution. Libraries were then sequenced on an Illumina HiSeq2000 flow cell to give 50 bp paired-end reads using three lanes i.e. two samples per lane.

Mapping of reads to the reference genome and transcript assembly

Mapping of Illumina reads and transcript assembly are the two main steps of the RNA-seq workflow [19] and were conducted by Source BioScience (Nottingham, UK). In brief, the raw reads were filtered with the FASTQ_Quality_Filter tool in which the FASTQ format provides per-base quality scores to each read. The Q30 threshold (indicating 99.9% accuracy) was set as a quality score and read quality was further assessed by using fast QC application. The spliced read mapper TopHat version 2.0.6 that uses Bowtie as an alignment ‘engine’ was used to map the filtered reads to the F. vesca genome v1.1 [16]. The F. vesca reference genome (HW) was downloaded from the Genome Database for Rosaceae (GDR) (www.rosaceae.org/species/fragaria/fragaria_ vesca/genome_v1.1). Alignment output files were provided to Cufflinks version 2.0.2 to generate transcript assembly and to estimate the relative transcript abundance for all six samples.

Differential gene expression analysis

We employed Cufflinks package version 2.0.2 to identify Differentially Expressed Genes (DEGs). Cufflinks (http://cufflinks.cbcb. umd.edu) was downloaded to a Linux system. The transcript GTF files produced after assembling were merged together using the Cuffmerge utility. The reference GTF (fvesca_v1.1_genemark_hybrid.gff3.gz) and FASTA files (fvesca_v1.1_pseudo.fna.gz) downloaded from GDR were provided along with the transcript assemblies of all six samples to create a single merged transcriptome file. The Cuffdiff program was performed to identify DEGs and was run using the merged transcriptome assembly file along with the accepted_hits.bam files from TopHat for all samples [19]. The expression level of each gene was normalized by the Fragments per Kilobase of transcript per Million mapped fragments (FPKM) value. The FPKM value of each gene was used to calculate fold change and log2 ratio ≥2 was used as threshold to identify DEGs.

Data access

The raw RNA-seq data have been deposited in the Sequence Read Archive at the European Nucleotide Archive under the study accession number PRJEB4896 (http://www.ebi.ac.uk/ena/data/view/PRJEB4896).

Powdery mildew infection and symptom development

The P. aphanis fungus was used to infect plants of HW and YW. Successful infection was shown by the appearance of visual symptoms on infected plants (Figure 1A). The observation of powdery mildew infected leaves under an epifluorescence microscope confirmed the presence of conidia and their development stages 1 DAI and 8 DAI (Figure 1B). After 1 h of infection the leaves showed the presence of conidia. The conidia started germinating and primary haustorium was observed 1 DAI and conidiation at 8 DAI.

Figure 1: (A) Symptom development on F. vesca ssp. vesca accession Hawaii 4 and F. vesca f. semperflorens line “Yellow Wonder 5AF7”, (a, b, c) HW control, 1 d and 8 d after infection, respectively, (d, e, f) YW control, 1 d and 8 d after infection, respectively. (B) Epifluorescence microscopic images of powdery mildew conidia and their development stages. Powdery mildew conidia stained with 8-Anilino-Naphthalene-Sulfonic Acid (ANSA), (a) Conidia on the leaves after 1 h of inoculation; (b) Primary haustoria developed at 1 d after infection; (c) Conidiation stage at 8 d after infection.

Analysis of RNA-seq data

Six cDNA libraries, HW (control, 1 DAI, 8 DAI) and YW (control, 1 DAI, 8 DAI) were generated and sequenced using the Illumina HiSeq2000 platform and the raw reads obtained successfully passed the quality filtering process with 95% accuracy. In total, more than 1,000 million paired-end reads of 50 bp in length were generated (Table 1). The sequence reads were mapped using Bowtie/TopHat to the reference genome sequence of F. vesca ssp. vesca accession Hawaii 4. Of the total reads, about 999 million reads (92%) mapped to the F. vesca genome, which is sufficient for quantitative analysis of gene expression and 120 million (11%) were unmapped. The percentage of mapped reads of all six libraries was in range of 90 to 94%.

Table 1: Summary of read numbers based on RNA-seq data.

Global analysis of transcriptome changes in F. vesca cultivars during powdery mildew infection

To investigate the underlying molecular changes in HW and YW during P. aphanis infection, we used RNA-seq to generate transcriptome profiles for each of the three stages (control, 1 DAI and 8 DAI). A total of 21,707 (HW control), 21,774 (HW 1 DAI), 22,146 (HW 8 DAI), 21,555 (YW control), 21,699 (YW 1 DAI) and 22,192 (YW 8 DAI) genes were detected in the libraries (Figure 2A, Data S1). These genes from all six libraries were combined separately, and a Venn diagram was generated to reveal unique or commonly expressed genes (Figure 2B). Among these libraries, 20,344 and 20,277 genes were commonly expressed in all three stages of HW and YW, respectively. The commonly expressed genes of HW and YW libraries were further combined and in total, 19,572 genes were expressed among all three stages of HW and YW. To identify DEGs in response to infection in HW and YW leaf tissue, FPKM values of each gene was used to calculate fold change between control and 1 DAI (1DAI/Control), control and 8 DAI (8DAI/1DAI) and 1 DAI and 8 DAI (8DAI/1DAI). Out of 34,809 predicted F. vesca genes [16], a total of 23,470 and 23, 464 genes were identified in HW and YW, respectively at all three stages and analysis based on log2 fold change ≥2 set as a threshold between stages detected the number of differentially expressed genes. As shown in (Figure 3A and 3B; Data S2 and S3), a greater number of up-regulated genes, 1,567 (1DAI/Control), 1,846 (8DAI/Control) and 1,145 (8DAI/1DAI) were identified in HW compared to YW, 1,336 (1DAI/Control), 1,619 (8DAI/Control) and 968 (8DAI/1DAI). The number of down-regulated genes, 646 (1DAI/ Control), 1,098 (8DAI/Control) and 624 (8DAI/1DAI) were also higher in HW than in YW, 571 (1DAI/Control), 754 (8DAI/Control) and 627 (8DAI/1DAI). To highlight major transcriptomic changes, a log2 fold change of greater than 5 was further examined to analyze DEGs between control, 1DAI and 8DAI (Figure 4A). In comparison, between 1DAI/Control, more up-regulated genes were identified in YW (139) compared to HW (106), whereas the number of up-regulated genes was higher in HW between 8DAI/Control (400) and 8DAI/1DAI (110) than in YW, 361 and 77, respectively. A higher number of down-regulated genes, 100 (1DAI/Control), 100 (8DAI/Control) and 70 (8DAI/1DAI) was identified in HW compared to YW 75 (1DAI/Control), 75 (8DAI/ Control) and 66 (8DAI/1DAI).

Figure 2: (A) Graph showing total number of genes expressed at all three stages (Control, 1DAI and 8DAI) of HW and YW. (B) Venn diagram showing the number of common and uniquely expressed genes in each of the three stages (Control, 1DAI and 8DAI) of HW and YW. Among these genes, 20,344 and 20,277 are expressed commonly at all three stages of HW and YW, respectively and of these, 19,572 are co-expressed in HW and YW.

Figure 3: Global views of number of differentially expressed genes (log2 fold change >2) at different times after powdery mildew infection in HW and YW. (A) Up-regulated genes. (B) Down-regulated genes.

Figure 4: (A) Global views of number of differentially expressed genes (log2 fold change > 5) at different times after powdery mildew infection in HW and YW. (a) Up-regulated genes. (b) Down-regulated genes. (B) Venn diagram showing the combination of differentially expressed genes (log2 fold change >5) between 1DAI/Control in HW and YW.

Analysis of transcriptomic changes between control and 1 DAI (1DAI/Control)

We investigated the genes that showed major transcriptional changes (log2 fold change >5) in response to P. aphanis infection in HW and YW between 1DAI/Control as major events related to defence or susceptible responses occur during this time period. Gene ontology (GO) annotation was performed by using Plaza2.5 [20] to study the biological process, molecular function, and cellular component in which DEGs (log2 fold change >5) between 1DAI/Control were involved. Analysis of categories showed that the DEGs were grouped into 27 different biological processes, of which ‘metabolic process’, ‘cellular process’, and ‘primary metabolic process’ were identified in all libraries (Figures S1 and S2). About five categories of biological processes, including ‘response to stress’, ‘macromolecule biosynthetic process’, ‘defense response’, ‘cellular nucleobase and nucleic acid metabolic process’ and ‘cellular macromolecule biosynthetic process’ were found to be up-regulated in HW, whereas, in YW, ‘oxidationreduction process’, ‘protein metabolic process’, ‘cellular protein metabolic process’, ‘phosphate metabolic process’ and ‘phosphorous metabolic process’ were up-regulated. Four categories, namely ‘transmembrane transport’, ‘establishment of localization’, ‘localization’ and ‘transport’ were detected in HW but were only down-regulated. Additionally, two categories, ‘carbohydrate metabolic process’ and ‘response to abiotic stimulus’ were down-regulated in YW. Furthermore, two biological processes included ‘response to stimulus’, and ‘nitrogen compound metabolic process’ were expressed in the opposite manner when compared between HW and YW; that is they were up-regulated in HW but down-regulated in YW. It showed that normal developmental processes were strongly affected after pathogen infection. Further analysis, involving the functional classification of genes expressed in HW and YW, showed 21 different categories of molecular functions of which two major categories demonstrated upand down-regulation in HW and YW; these were ‘binding’ and ‘catalytic activity’ while ‘hydrolase activity’ was also found in all libraries but was not dominant (Figures S3 and S4). For up-regulated genes, three molecular functional categories were represented in both HW and YW including, ‘transition metal ion binding’, ‘oxidoreductase activity’, and ‘transferase activity, transferring phosphorus-containing groups’. The molecular functions related to ‘protein binding’, ‘adenyl ribonucleotide binding’, ‘purine ribonucleotide binding’, and ‘ribonucleotide binding’ were found to be down-regulated in HW, whereas ‘hydrolase activity, hydrolysing o-glycosyl compounds’ and ‘hydrolase activity, acting on glycosyl bonds’ were down-regulated in YW. As with biological process, several molecular functional categories were up-regulated in HW and down-regulated in YW; these included ‘cation binding’, ‘ion binding’, ‘metal ion binding’, and ‘transferase activity’. Additionally, three categories namely ‘nucleotide binding’, ‘purine nucleoside binding’, and ‘nucleoside binding’ were up- and down-regulated in HW but, only upregulated in YW. Furthermore, the ‘purine nucleotide binding’ category was found to be expressed in both the libraries of HW whereas; ‘nucleic acid binding’ was up-regulated in HW but down-regulated in YW. However, the cellular component was also investigated and showed 15 different categories of which ‘cell part’, ‘cell’, and ‘intracellular’ as the major categories in all libraries (Figures S5 and S6). The ‘nucleus’ category was found to be up-regulated in HW while ‘membrane part’ and ‘integral to membrane’ were up-regulated in YW; the ‘cytoplasmic part’ category was down-regulated in both HW and YW.

Analysis of the combination of DEGs common to HW and YW between 1 DAI/Control

The genes identified as being differentially expressed (log2 fold change >5) in HW and YW between 1DAI/Control were combined to distinguish the number of unique or common genes that were shared among up- and down-regulated subgroups (Figure 4B). Of 106 and 139 up-regulated genes in HW and YW, respectively, eighteen genes were found to be shared. In contrast, four genes were identified to be similar among 100 and 75 down-regulated genes of HW and YW, respectively. Furthermore, there were some genes that were shared in an opposite manner when comparing between the up- and down-regulated genes in HW and YW. Four genes were up-regulated in HW but downregulated in YW. Similarly, nine genes were down-regulated in HW but up-regulated in YW.

Genes commonly up-regulated in HW and YW

The eighteen up-regulated genes that were common to HW and YW were analyzed (Table 2). Figure 5 displays the comparison of gene expression using log2 fold change value between 1DAI/Control in HW and YW. By using absolute and log2 fold change value, it was noticed that three genes were highly expressed in HW compared to YW; these included:- gene3227 (most highly up-regulated gene in HW but 3rd most up-regulated gene in YW) encoding a flavonoid 3´-monooxygenase, involved in secondary metabolite biosynthesis, transport, catabolism and oxidation reduction of biological process; gene16607 (2nd most upregulated gene in HW) encoding a PHD finger protein, a transcription factor that plays major role in regulation of transcription; and gene17283 (5th most up-regulated gene in HW) encoding the serine protease HTRA2, involved in heat-shock response, chaperone function and apoptosis. Gene05674 involved in molecular function (ATP binding) and encoding polyadenylate-binding protein2 was highly expressed with an absolute fold change about 6837-12745 and a log2 fold change of 13 and 14 in both HW and YW, respectively. Moreover, gene10449 (2nd most up-regulated gene in YW) involved in molecular function (transferase activity), encoding a BAHD acyltransferase, was highly expressed in YW compared to HW. Gene06825, which encodes a flagella-related protein G related to cellular component process involved in cell motility and secretion, was highly expressed in both HW and YW. Furthermore, twelve genes which were up-regulated with a log2 fold change value of less than 10, both in HW and YW were identified. These included: gene 03919 (putative ribonuclease H protein) involved in catalytic activity, gene13298 (probable molybdenum cofactor synthesis protein 2A) involved in biosynthesis of molybdenum cofactor that is essential for a diverse group of redox enzymes; gene09059 (probable a patatin-T5, precursor) with phospholipase activity involved in biological process and having a role in plant defense mechanism; gene 26471 (probable inhibitor of apoptosis protein) encodes for an apoptotic suppressor protein; gene17605 (probable ADP-ribosylation factor GTPase-activating protein AGD11) functions in ARF GTPase activator activity, zinc ion binding and involved in regulation of ARF GTPase activity; gene24405 (similar to peroxidase 15) involved in biological process of response to oxidative stress; gene18228 (probable alkyldihydroxyacetonephosphate synthase, peroxisomal (alkyl- DHAP synthase)) involved in ether type lipid metabolism; gene07051 (probable methionyl-tRNA synthetase) that plays an important role in protein biosynthesis; gene01644 (similar to a germin-like protein subfamily 1 member 18) cupin superfamily protein, which functions in manganese ion binding, nutrient reservoir activity and is involved in biological process of plant defence response; gene10952 (similar to a putative molybdate metabolism regulator); gene31927 (probable ketol-acid reductoisomerase) involved in amino acid biosynthetic process; and gene07076 (probable major allergen Pru ar 1) response to biotic stimulus of biological process category. Comparison of FPKM values of eighteen genes in control leaves of both HW and YW, showed gene31927 found to be highly expressed.

Figure 5: Graph showing the expression level of genes up-regulated (log2 fold change >5) between 1DAI/Control in both HW and YW.

Table 2: Genes (log2 fold change >5) expressed in HW and YW and showing similar up-regulation between 1DAI/Control.

Genes commonly down-regulated in HW and YW

Further analysis of four down-regulated genes shown to be common to both HW and YW was performed and it was observed that these genes were down-regulated with log2 fold change values of less than -10 and there was little difference in the level of expression of these genes between HW and YW (Table 3 and Figure 6). These included: gene06512 (probable mucin-2, precursor similar to a human intestinal protein); gene13823 (probable formate-tetrahydrofolate ligase (FHS)), which functions in formate-tetrahydrofolate ligase activity, copper ion binding, ATP binding and is involved in response to cadmium ion, folic acid and derivative biosynthetic process; gene16447 (similar to an acyltransferase-like protein Atg54570, choloroplastic, precursor) predicted to be involved in hydrolases or acyltransferases; and gene18993 (similar to the Os01g0513500 protein fragment) with unknown function. Gene16447 was found to be expressed highly in control leaves of both HW and YW compared among these four downregulated genes.

Figure 6: Graph showing the expression level of genes down-regulated (log2 fold change >5) between 1DAI/Control in both HW and YW.

Table 3: Genes (log2 fold change >5) expressed in HW and YW and showing similar down-regulation between 1DAI/Control

Common genes up-regulated in HW and down-regulated in YW

Four genes which were expressed in an opposite manner, that is upregulated in HW and down-regulated in YW, were investigated (Table 4 and Figure 7). Gene17646 showed a high log2 fold change of 8 in HW whereas it was down-regulated in YW with a log2 fold change of -7. This gene is likely to encode a Tau-tubulin kinase 1, which is similar to a human brain-specific protein kinase involved in tau phosphorylation. One gene, gene31970 (probable potassium-transporting ATPase C chain) involved in transporting potassium ion in plasma membrane was down-regulated in YW with a log2 fold change of -10 but was upregulated in HW with a log2 fold change of 7.3. The other two genes included in this category are gene20901 (probable glucan endo-1,3- beta-glucosidase 2) and gene05145 (putative lichenase); both have same functions in hydrolase activity, hydrolysing O-glycosyl compounds, and are involved in carbohydrate metabolic process. Among these four genes, gene05145 showed high expression in YW control leaves.

Figure 7: Graph showing the expression level of genes up-regulated in HW and down-regulated in YW (log2 fold change >5) between 1DAI/Control.

Table 4: Common genes (log2 fold change >5) expressed in an opposite manner in HW and YW between 1DAI/Control.

Common genes down-regulated in HW and up-regulated in YW

Furthermore, nine genes expressed in an opposite manner, that is down-regulated in HW but up-regulated in YW, were also studied (Table 4 and Figure 8). These included: gene06603 (4th most downregulated gene in HW) encoding a mitogen-activated protein kinases ANP1 involved in the molecular function of protein serine/threonine kinase activity; gene03770 (highly expressed in YW with an absolute fold change value about 5074) likely to be an aluminium activated malate transporter) expressed in response to chemical stimulus; gene06811 which is expressed in response to stress and known as a universal stress protein found to be highly down-regulated in HW; gene33420 encoding a biotin synthase involved in biotin/thiamine synthesis-associated protein; gene23225 (a serine/arginine repetitive matrix protein 2) involved in RNA binding; gene25958 (with unknown function) was found to be highly expressed in YW, gene34745 (LysM domain-containing GPI-anchored protein 2) involved in cell wall macromolecule catabolic process; gene00527 (similar to TC10/CDC42 GTPase-activating protein) involved in enhancing the GTPase activity; and gene22892 (ABC transporter B family member 11) has a molecular function of ATP binding and is involved in transport of various molecules across extra- and intra-cellular membranes

Figure 8: Graph showing the expression level of genes down-regulated in HW and up-regulated in YW (log2 fold change >5) between 1DAI/Control.

Gene regulation in response to P. Aphanis

Several subsets of genes, including transcription factors, receptorlike kinases, F-box proteins and putative disease resistance genes were identified from the DEGs (log2 fold change >5) in HW and YW between 1DAI/Control.

Transcription factors (TFs)

The F. vesca TFs induced upon P. aphanis infection were identified by manually checking against the recently reported 1,493 TFs in the F. vesca genome [21] and the NCBI database. In total, 34 differentially expressed TFs were detected in both HW and YW that belong to several families, including bZIP, C2H2, NAC, ERF, MYB, MADS, HSF, bHLH, B3, and other families of TFs (Table S1). Out of these 34 genes, one TF, PHD finger protein (gene16607) was found to be up-regulated in both HW and YW. Nineteen TFs were differentially expressed in HW with seven transcriptionally down-regulated. Conversely, YW exhibited sixteen differentially expressed TFs with eight down-regulated.

Receptor-like kinases (RLKs)

A total of 24 differentially expressed RLKs were identified in HW and YW, including LRR-RLKs, MAP kinases, WAK kinases, cysteinerich RLKS, serine/threonine-protein kinases, lectin-RLKs, HSK and CBL-kinases (Table S2). Out of 24 RLKs, Mitogen-activated protein kinase ANP1 (gene06603) was found to be down-regulated in HW but up-regulated in YW. In HW, four RLKS were up-regulated and four down-regulated; while in YW the majority of RLKs (13) were upregulated and four down-regulated.

F-box proteins

Seven F-box proteins were identified upon powdery mildew infection (Table S3) out of which, only one was up-regulated in HW, while five were up-regulated in YW. No F-box proteins were downregulated in HW, whereas one F-box protein was found to be downregulated in YW.

Putative disease resistance genes

Genes involved in plant defence mechanisms upon pathogen infection were also identified. In HW, four genes, including one NBS-LRR, two GLPs and one calmodulin gene were found to be upregulated and interestingly one gene, an NB-ARC domain containing disease resistance protein (gene25859) was down-regulated (Table S4). Conversely, YW exhibited up-regulation of five genes, including two GLPs, one LRR and NB-ARC domain containing disease resistance protein, one Pathogenesis-related protein (PR-1A) and a calmodulin (CaM) gene. Gene01644 (germin-like protein subfamily 1 member 18) was detected to be up-regulated in both HW and YW.

Despite several studies made in strawberry disease research, progress has been limited in providing molecular information to engineer strawberry plants for durable and broad-spectrum disease resistance. Gene expression profiling provides important data in understanding the molecular interaction between host and different classes of microbial pathogens. Transcriptome analyses of detrimental plant-pathogen interactions have been performed mostly using microarrays [22-24]. Recently, one study used RNA-seq in rice to identify defence related genes as well as genes involved in signalling and pathogen-related gene expression regulation in response to the rice blast pathogen M. oryzae [14]. In the present study, RNA-seq was employed to generate large transcriptome data sets of two diploid strawberry varieties, F. vesca ssp. vesca accession Hawaii 4 (HW) and F. vesca f. semperflorens line “Yellow Wonder 5AF7” (YW) upon powdery mildew infection at different time points, between control and 1 DAI (1 DAI/Control), control and 8 DAI (8 DAI/Control) and 1 DAI and 8 DAI (8 DAI/1 DAI). Out of 34,809 predicted F. vesca genes [16], 23,470 and 23,464 genes were detected in HW and YW, respectively. Comparing the number of expressed genes with a previous study, Kang et al. [22] identified 22,611 genes in the fruit tissue of Yellow Wonder 5AF7. Analysis of DEGs based on log2 fold change ≥2 showed that a greater number of up- and down-regulated genes were identified in HW compared to YW between all stages. Further detailed analysis of DEGs (log2 fold change >5) was performed in HW and YW between 1 DAI/ Control as crucial events such as conidial germination, penetration of the epidermal cells, host perception and defence, transduction and exchange of signals usually occur within 24 h after the initiation of infection [7,25]. Several common and uniquely expressed genes were identified and the analysis showed that the DEGs were involved in secondary metabolism, signal transduction, transcription regulation, and defence response. Cytochrome p450 in the phenylpropanoid pathway controls lignin synthesis, flower pigments, signalling molecules and several compounds involved in plant defence against microbes and UV light [26]. In this study, a gene encoding a member of cytochrome p450, flavonoid 3´-monooxygenase involved in flavonoid biosynthesis was strongly up-regulated in HW and YW upon powdery mildew infection. This finding was in agreement with the previously reported induction of cytochrome p450 monoxygenase gene induction in response to C. acutatum in strawberry [22] indicating that flavonoid metabolites play an important role in defence response to pathogen attack. The pathogenesis-related (PR) proteins have antimicrobial properties through hydrolytic activities on cell walls, contact toxicity, are involved in defence signalling and are induced upon stress, pathogen attack, and application of signalling compounds such as salicylic acid, jasmonic acid, and ethylene [27]. Four PR proteins were regulated in response to P. aphanis infection in HW and YW; these included peroxidase 15 (PR-9 protein family) and major allergen Pru ar 1 belong to group of PR proteins (PR-10 protein family), which were upregulated in both HW and YW. Two other PR proteins were lichenase and glucan endo-1, 3-beta-glucosidase 2 (PR-2 protein family), which can hydrolyze (1-3, 1-6) branched ß-glucans of fungal cell walls [28,29] were up-regulated in HW but interestingly, down-regulated in YW. Previous microarray analysis of barley-powdery mildew interaction [30] and grapevine-downy mildew interactions [31] showed that these PR protein families were found to be up-regulated and are involved in antifungal activity. It has been reported that more rapid induction of these PR protein transcripts could slow down fungal penetration by weakening or damaging the penetration peg formed within 15 h after infection [32]. Transcription factors (TFs) are reported to be important factors during plant development, signal transduction, response to biotic and abiotic stresses [33]. Induction of stress response genes occurs mainly at the transcription level and regulates the temporal and spatial expression patterns of specific stress genes. These TFs belong to large gene families. Several transcription family proteins including C2H2, NAC, bZIP, ERF, MYB, MADS, bHLH, HSF, B3, HSF, FAR1, WD, Trihelix, CAMTA and HB-other were up- and down-regulated in both HW and YW between 1 DAI/Control (Data analyzed for DEGs of log2 fold change >5). As recent studies have shown the up-regulation of several Zn finger, AP2/ERF and NAC family TFs in response to powdery mildew in barley suggested that these are involved in positive regulation of powdery mildew penetration resistance [30,34]. Conversely, the majority of TFs including bZIP/HD-ZIP, MYB, C2H2 Zn finger, bHLH, NAC and MYC were down-regulated based on the RNA-seq analysis after V. dahliae infection in cotton [12]. Detailed analysis of the TFs identified in this study will provide an insight on the resistance or susceptible response to powdery mildew infection in strawberry. Plant receptor-like kinases (RLKs) are transmembrane proteins containing cytoplasmic serine/threonine kinase domains and different extracellular domains which can potentially bind to a range of molecules, including carbohydrates, polypeptides, lipids, microbial cell-wall constituents and steroids, based on their extracellular domain [35]. Plant innate immunity can be achieved by pattern recognition receptors (PRRs), comprising well characterized LRR-RLKs [36]. The groups of differentially expressed LRR-RLKs identified in the present study indicate putative roles for LRR-RLKS in mediating early events in the response of strawberry to P. aphanis. Two LRR-RLKs were detected in HW, one up-regulated and one down-regulated. Three LRR-RLKs identified in YW were up-regulated. The recent studies involving interaction of B. graminis in barley and V. dahliae in cotton also induced several LRR-RLKs [12,30]. Wall-associated kinases (WAKs) represent a novel sub-family of plant RLKs [37]. WAKs in Arabidopsis are involved in cell expansion, resistance to pathogens and tolerance to heavy metal stress [38]. Three WAKs were detected in this study and were upregulated in YW. Mitogen-activated protein kinase (MAPK) cascades are associated with a wide range of developmental, hormonal and signalling responses. The MAPK cascade involving the triple kinases, MAPKKK was expressed commonly in both HW and YW, but strongly down-regulated in HW and up-regulated in YW. It has been recently reported that several WAKs and MAPKs were differentially expressed in response to blast fungus in resistant and susceptible rice cultivars based on RNA-seq analysis [14]. Detailed functional investigations would allow the identification of RLKs mediated strawberry immunitybox proteins contain the F-box domain and are involved in regulation of several plant developmental processes including, floral meristem and floral organ development, photomorphogenesis, circadian clock regulation and self-incompatibility [39,40]. The C terminus of F-box proteins comprise of one or several protein-protein interaction domains, such as LRR, Kelch repeat, WD40 repeat and tetratricopeptide repeat which interact with specific targets [41] . Almost 700 F-box proteins have been identified in Arabidopsis [42] and 687 in rice [43]. Previously, Woo et al. [44] reported that ORE9, an F-box protein functions in early senescence of leaf in A. thaliana through ubiquitindependent proteolysis. In the present study, gene00915 (similar to F-box protein ORE9) was up-regulated in HW and another gene, gene09652 (probable F-box protein ORE9) was also up-regulated in YW suggesting that these proteins play a major role in the early senescence of leaf in strawberry in response to powdery mildew infection. In a recent study, Gou et al. [45] showed that an F-box protein CPR30, functions as a negative regulator of the defence response in A. thaliana and the loss of CPR30 function exhibited a dwarf phenotype and increased resistance to both virulent and avirulant bacterial pathogens. One gene, gene17272 encoding the F-box protein CPR30 was found to be up-regulated in YW in our study, indicating that establishment of the strawberry resistance to P. aphanis might require the loss-of-function of CPR30 to modulate the plant immunity. Several other F-box proteins containing LRR and kelch repeat were differentially expressed in YW. A comprehensive expression analysis of strawberry F-box protein encoding genes will help in understanding their roles in plant defence response. In plants, the majority of the disease resistance (R) genes encode the ‘nucleotide-binding site plus leucine-rich repeat’ (NBSLRR) class of proteins characterized by NBS and LRR domains as well as variable amino- and carboxy-terminal domains. The NBS domain is also called an NB-ARC domain since it is composed of two sub domains, NB and ARC (nucleotide binding adaptor shared by NODLRR proteins, APAF-1, R proteins and CED4) [46]. These abundant proteins are involved in detection of diverse pathogens and disease resistance [47]. Several NBS-LRR genes have been identified in many crops. Recently, twelve R genes belong to the NBS-LRR gene family of disease resistance genes have been identified during rice-M. oryzae interaction [48]. In grapevine, Dry et al. [49] identified seven RGAs encoding NBS-LRR resistance proteins in response to powdery mildew. In our study, one NBS-LRR protein was up-regulated in HW and interestingly, one NB-ARC domain containing protein found to be down-regulated in HW. One LRR and NB-ARC domain containing protein was up-regulated in YW. These resistance genes serve as potential candidates in mediating powdery mildew resistance in strawberry. Germins and germin-like proteins (GLPs) are characterized as glycoproteins are capable of oxidoreductase activity which results in the production of H2O2, and are believed to play significant role in plant signalling mediated defence against pathogens [50]. GLP encoding genes can be induced by fungal pathogen infection [51], and several studies have demonstrated that GLPs are involved in plant defence responses against fungal pathogens [52]. Previously, Knecht et al. [53] showed the expression of BvGLP-1, encoding a GLP from sugar beet in transgenic A. thaliana conferred significant resistance to two pathogenic fungi, R. solani and V. longisporum. A recent study based on a comparative proteome analysis of cultivated strawberry in response to F. oxysporum showed the down-regulation of a gene09299 encoding a GLP [54] and in our study, the same gene was up-regulated in YW in response to P. aphanis. Two other genes encoding GLPs were upregulated and one was expressed in both HW and YW. No GLP genes were found to be down-regulated in this study (Data analyzed for DEGs log2 fold change >5). It appears that GLPs contribute to disease resistance in this pathosystem. In summary, the large transcriptome data set generated for two diploid strawberry varieties in response to powdery mildew infection provides a wealth of genomic information to better understand the molecular mechanism of the strawberry defence response to P. aphanis and enabled new insights into the identification of genes expressed in terms of both up- and down-regulation during powdery mildew interaction. Detailed studies on the expression patterns of relevant genes (including genes log2 fold change >2) would provide potential candidates for genetic improvement of strawberry. Furthermore, our study provides molecular basis for the selection of certain strawberry disease related genes, such as peroxidase 15, GLPs, glucan endo-1, 3-beta-glucosidase 2, MAP kinases, RLKs, NBS-LRR and TFs. Further characterization of these genes may provide potential candidates for genetic manipulation to generate transgenic strawberry plants that are resistant to powdery mildew and other fungal diseases. Moreover, the availability of our raw RNA-seq data (via Sequence Read Archive at the European Nucleotide Archive under the study accession number PRJEB4896) paves the way for future functional analysis of genes and networks that regulate disease resistance in strawberry. In addition the unmapped reads (>120 million) produced in this transcriptomic study can also be used in the future to identify transcripts of fungal origin. This will be dependent on the production of a pathogen genome, the aim of present studies in other laboratories.

We would like to thank Felix trust for financial support, and Dr Janet Slovin, Beltsville, USA and Dr Richard Harrison, East Malling Research, UK for providing materials utilised for the study.