Pathogen recognition in compatible plant-microbe interactions

Microbial infections in plant leaves remain a major challenge in agriculture. Hence an understanding of disease mechanisms at the molecular level is of paramount importance for identifying possible intervention points for their control. Whole-transcriptome changes during early disease stages in susceptible plant species are less well-documented than those of resistant ones. This study focuses on the differential transcriptional changes at 24 hours post inoculation (hpi) in tomato leaflets affected by three pathogens: (1) Phytophthora infestans, (2) Botrytis cinerea, and (3) Oidium neolycopersici. Grey mould (B. cinerea) was the disease that had progressed the most by 24 hpi, both in terms of visible symptoms as well as differential gene expression. By means of RNA-seq, we identified 50 differentially expressed tomato genes specifically induced by B. cinerea infection and 18 specifically induced by P. infestans infection at 24 hpi. Additionally, a set of 63 genes were differentially expressed during all three diseases when compared by a Bayesian approach to their respective mock infections. And Gene expression patterns were found to also depend on the inoculation technique. These findings suggest a specific and distinct transcriptional response in plant leaf tissue in reaction to B. cinerea and P. infestans invasion at 24 hpi, indicating that plants may recognize the attacking pathogen.


Results
Infection with B. cinerea is evident at 24 hpi. Nine pathogen-infected samples (PI: P. infestans, BC: B. cinerea and ON: O. neolycopersici) and their respective mock-inoculated variants (PIm, BCm and ONm) were subjected to Illumina HiSeq sequencing of the transcribed RNA at 24 hpi, which resulted in more than 1.13 billion total read pairs with a Q30 of 96% (project available at http://www.ebi.ac.uk/ena/data/view/PRJEB21223). Approximately 63 million read pairs were generated per sample. The reads were mapped against the available genomes of S. lycopersicum Heinz 1706, B. cinerea T4 and P. infestans. For O. neolycopersici, there was no genome sequence available. On average 80% of the reads per sample mapped to the host S. lycopersicum. 28 million read pairs (13.7%) of the B. cinerea and 1.4 (0.7%) of the P. infestans inoculated samples mapped to their respective genomes. One mock-inoculated sample of P. infestans infections yielded only 35% reads for S. lycopersicum and a Penicillium spp contamination was identified. Nevertheless the S. lycopersicon-derived reads clustered with the other mock-infections and therefore the sample was considered healthy. Pairwise comparison using edgeR 22 on all the samples suggests similar patterns for all the mock-inoculated samples as well as for those inoculated with P. infestans and O. neolycopersici. In contrast, samples inoculated with B. cinerea clustered separately (Fig. 1), indicating a more advanced progression of the disease within the first 24 hpi.
Eleven genes that were differentially expressed in pairwise comparison of the pathogen and its corresponding mock inoculation were used to validate the RNA-seq expression. Trends in up-and down-regulated RNA-seq-logFC (log2 fold change) correspond to the mock sample-normalized expression measured by means of qPCR (quantitative real-time PCR; Supplementary Table S1). Three genes with no indicated differential expression in RNA-seq and a FDR > 0.69 were up-regulated in a qPCR evaluation. The R 2 of the linear correlation in the RNA-seq and qPCR data was 0.92 (Supplementary Figure S1). Time-course experiments with sampling in three-hours intervals from 0-24 hpi and then at 48 hpi lead to a correlation at 24 hpi with the RNA-seq expressions of R 2 = 0.8 (Supplementary Figure S2 and S3).

Inoculation techniques alter gene expression. Nine patterns of cluster type [PI&PIm] vs.
[BC&BCm&ON&ONm] comprising 668 genes were identified by summarizing genes that were differentially regulated between abaxial (PI and PIm) and adaxial (BC, BCm, ON and ONm) drop inoculation (all patterns can be found in Supplementary Data S1 with an explanatory file in Supplementary Data S2). Since B. cinerea was inoculated in half-strength grape juice, we tested if genes were clustering together in the pattern types [BC&BCm] vs. [PI&PIm&ON&ONm] and were able identify four patterns comprising 609 genes. Each sequenced gene can only be present in one pattern, hence there is no overlap between the genes expressed by different inoculation positions, inoculation media, and the ones identified as disease specific candidate genes that are described below. This suggests that both the location of the droplet deposition and the medium used for inoculation may impact gene expression.
Genes are regulated by disease in general. Using the R package EBSeq 23 we identified one gene expression pattern (Fig. 2) that separates 63 genes by their differential regulation when diseased or healthy (Supplementary Table S2).
The most DE genes were found in B. cinerea inoculated tissues. Gene expression following infection with each disease and their corresponding mock inoculations were analysed by pairwise differential expression analysis using edgeR 22 Table S3). Comparing the differential expression between the diseases, we found 511 up-regulated genes in both the B. cinerea and the P. infestans infected tissues, and 24 and 3 down-regulated genes in B. cinerea-P. infestans and B. cinerea-O. neolycopersici infected leaf tissue, respectively (Supplementary Figure S4). Grey mould and late blight-specific regulated genes were identified. In an attempt to identify disease-specific genes, the data were analysed with EBseq, and genes of informative patterns (i.e. sample combinations clustering for their gene expression against other samples which indicate whether they are "healthy" (BCm&PIm&ONm) or "diseased" (BC&PI&ON)) were extracted and compared to a pairwise comparison build  with edgeR using a false discovery rate (FDR) < 1% and a logFC > 0.9. Genes belonging to a disease-specific pattern (only expressed for one disease and not for any other tested pathogen or mock inoculated samples) were extracted and analysed for their putative function using (1) Gene Ontology (GO) enrichment analysis, (2) an ortholog gene search of A. thaliana and Solanum tuberosum derived genes, and (3) protein database comparison.
EBseq resulted in 46 patterns for B. cinerea, 3 for P. infestans and 3 for O. neolycopersici, include gene candidates that are potentially disease-specifically regulated ( Table 1). The 46 B. cinerea-specific patterns comprised 10009 genes, while the three patterns specific for P. infestans and O. neolycopersici included 178 and 80 genes, respectively (Table 1). These genes selected with EBseq were compared with the results of the DE of the pairwise comparison of the following pairs: BC vs. BCm, PI vs. PIm, ON vs. ONm, BC vs. PI, BC vs. ON and PI vs. ON. If the pattern-identified disease-specific genes were 'DE' for the specific disease-mock comparison and 'not DE' for the other pairs, they were selected for further characterization. This allowed 28 up-and 22 down-regulated genes to be identified that were specific for B. cinerea infections at 24 hpi, 18 up-regulated and no down-regulated genes with late blight specificity, and no DE genes for powdery mildew ( Table 2). The B. cinerea infection specific up-regulated genes had a maximum logFC of 8.2 (Solyc04g028460.1.1, unknown protein) and the down-regulated genes a maximum of 5.3 (Solyc12g096490.1.1, GDU1). The strongest up-regulated gene during P. infestans infection was measured with a logFC of 4.4 (Solyc02g068670.1.1, Ankyrin repeat-containing protein).
18 of B. cinerea-up-regulated genes, 18 of the B. cinerea-down-regulated genes and 16 of the P. infestans-up-regulated genes were characterized by an ITAG functional description. Among these 52 genes, we identified one carrier protein, one glycosylase, one ion binding protein, three kinases, three ligases, one lipase, one oxygenase, one proteinase, one receptor protein, one reductase, one RNA polymerase, two transferases, one nitrate transporter, one leucine zipper and three transcriptions factors ( Table 2).
According to KEGG pathway information generated using STRING 24 1 (all coding for unknown proteins) were only found in the Solanaceae family. Since the orthologues have, to our knowledge, not yet been well characterized, the orthologue and co-occurrence analyses only provided information on the genes to be compared in future work.

Discussion
According to our knowledge, this paper represents the first differential gene expression study of S. lycopersicum affected by different diseases. Pairwise comparison of pathogen-and mock-treated samples at 24 hpi revealed that the greatest number of DE genes in the tomato specimens occurred when the plants were infected with B. cinerea. This number was less for P. infestans and even lower for O. neolycopersici. 75% of the DE genes induced by P. infestans were also found to be DE in the tomato-B. cinerea interaction. This is in contrast to a study into the compatible interaction of tomato with Cladosporium fulvum and Verticillium dahlia, which shared 454 DE genes out of a total of 4693 seven days after inoculation 27 . Furthermore, we identified genes that are specifically expressed  Table 2. Disease-specific genes that are differentially regulated in tomato leaves at 24 hpi. The genes have been selected using R package EBseq and pairwise comparison (edgeR). Full description of the selected genes is provided in Supplementary Table S4. at 24 hpi for both B. cinerea and P. infestans infected leaves. No such specific DE genes could be found for the infection with O. neolycopersici. Both the pairwise comparison and the disease-specific differential expression underline the limitations of the selected experimental setup, which used 24 hpi as a fixed time point for comparison, since each of the diseases has differing invasion and colonization strategies, as well as different development strategies over time and space. Due to its necrotrophic lifestyle B. cinerea is the fastest of the three pathogens in terms of invasion and destruction of host tissue, leading to visible symptoms within 24 hours 4 . An intermediate development speed was reported for P. infestans 7 , and this was also visible in our experiments. Powdery mildew was the last of the three diseases to display evident symptoms 12 . This is also true for other Erysiphaceae, whose symptoms develop in a range from 5.1 to 8 days, as has been shown for powdery mildew species screening in grapevines 28 . Additionally, the water-based inoculation method used may not provide optimal infection conditions for O. neolicopersici, since powdery mildew does not need complete wetness during infection 29 . Spores were often transferred by touching healthy tissue with sporulating lesions 12 , but the water-based method allows control and standardisation of the experimental setup in terms of experiment replication comparability with O. neolycopersici and other pathogen inoculations. And liquid conidia suspensions were used successfully in experiments with O. neolycopersici 30,31 . We identified 68 specific differentially expressed genes for grey mould and late blight at a relatively early disease progression state at 24 hpi. These genes were either up-or down-regulated and never lead to resistance, since the Heinz 1706 cultivar is highly susceptible to both diseases. Nevertheless, many of the up-and down-regulated candidate genes are related to environmental stress. All 68 identified gene accession numbers were checked using google scholar for citings in other publications. Referenced publications could only be found for eight genes. The disease-specific genes include six genes that are related to plant-pathogen interactions and two that are related to hormone signalling. One of these, the lipoxygenase (Solyc01g099210.2.1, up-regulated during B. cinerea infection) was reported to be up-regulated during root nematode infection 32 . SlWRKY80 (Solyc03g095770.2.1), which is up-regulated in P. infestans infected tissue, was also reported to be up-regulated during Pseudomonas syringae infiltration at 12 hpi 33 and 6 days after infection with Xanthomonas perforans race T3 34 . Solyc03g113580.1.1 (up-regulated during B. cinerea infection) coding for a Germin-like protein, was found to be expressed in tomato radicles grown under enhanced aluminium conditions 35 . Solyc06g053840.2.1, which is down-regulated during B. cinerea infection and coding for an auxin responsive protein (SlIAA4), is involved in hormone signal transduction and was up-regulated in young and old leaves and cotyledons compared to root 36 . Solyc06g075090.2.1 was down-regulated in B. cinerea infected leaf tissue and was reported to be involved in cytokinin related synthesis and signalling 37 . Solyc08g062490.2.1 (up-regulated in P. infestans infected leaf tissue) was annotated in the tomato genome (release ITAG2.4) as WRKY transcription factor 16, but is reported in literature as SlWRKY50 38 , a well-studied protein that mediates signalling of JA-and SA-pathways when the JA-pathway is repressed 39 .
Additionally, the B. cinerea inoculation down-regulated the Argonaute 1 gene (AGO1, Solyc01g010970.2.1), which may also be involved in biotic plant-environment interaction since an interaction of AGO1 with AGO2 in response to virus infection was demonstrated in A. thaliana 43 . The argonaute Piwi subfamily to which AGO1 belongs supports the silencing of mobile genetic elements 44 and antiviral RNA 45 . The down-regulated Solyc11g010340.1.1 gene is one of 152 bHLH transcription factors identified in tomato. One of these TFs (SlybHLH131, Solyc10g008270.2.1) was reported to be involved in the tomato reaction to tomato yellow leaf curl virus infection 46 . Solyc02g068670.1.1 was characterized to code for an "Ankyrin repeat-containing protein". These proteins were reported to be related to resistance as a potential negative regulator of pathogen-induced protein PR1 and antioxidation metabolism 47 . Furthermore, a number of DE genes may be involved in abiotic interaction: Solyc06g060830.2.1 down-regulated during B. cinerea infection codes for a putative homeobox-leucine zipper protein which is known to be involved in response to abiotic stresses 48 . The two genes coding for the homeobox-leucine zipper proteins ATHB12 and ATHB7 were up-regulated during drought stress 49 . The homologue gene AtRNP1 of Solyc07g048030.2.1 (Heterogeneous nuclear ribonucleoprotein A3) was reported to be involved in the abiotic stress response in A. thaliana 50 . In our experimental setup, the only remarkable abiotic difference of the disease infection setup of the B. cinerea infections was the media that was used for inoculation (half strength grape juice), which might have exerted some stress due to its higher sugar level and osmotic potential but Solyc07g048030.2.1 did not match the media-specific expression pattern. Solyc06g068960.1.1, Solyc03g115930.1.1 and Solyc06g069740.1.1 potentially code for calmodulin and calmodulin-like proteins which are Ca 2+ sensor proteins known to be involve in environmental stress responses 51,52 .
In summary, the genes mentioned in the previous sections are related to biotic and abiotic stresses. The Heinz 1706 tomato cultivar we used for our experiments is obviously highly susceptible to all three diseases, hence no single related stress pathway which might be induced by these genes leads to a resistant phenotype. Nevertheless, some of the 68 genes we identified as late blight and grey mould-specific could be of use as indicators for pre-symptomatic disease identification if they are expressed systemically.
The presented multi-disease comparison at 24 hpi revealed several major findings: First, we identified genes that are differentially regulated in tomato leaves both during B. cinerea and P. infestans infections by comparing grey mould, late blight and powdery mildew leaf diseases with their mock infections. The identified candidate genes may be of use in identifying one of these two diseases before symptoms development. Therefore, the regulation of candidates will be evaluated in future works for their temporal and spatial expression patterns. During O. neolycopersici infection, no disease-specific differentially expressed genes were identified when compared to late blight and grey mould disease. Nevertheless, a pairwise comparison of O. neolycopersici inoculated leaf tissue with corresponding mock-inoculated tissues identified some DE genes. Second, most of the late blight and grey mould specific DE tomato genes are apparently not directly related to plant-pathogen interactions. Third, the results suggest that inoculation location (abaxial and adaxial) and inoculum solution solvent (water and half-strength grape juice) both have an impact on gene expression.
In future studies we will analyse the expression of identified disease-specific candidate genes over time and space within the whole plant and assess the potential use of one or a combination of these genes for early pre-symptomatic disease detection.

Methods
Plant material. Tomato plants (S. lycopersicum, Heinz 1706 cultivar) were grown in standard soil in a semi-regulated greenhouse with open windows. The temperature was set to 20-26 °C with maxima during sunny summer days of up to 40 °C. On cloudy days, artificial light was used to achieve minimal constant lighting of 80 kW per square meter for 16 h per day. Once a week, cuttings were produced from the tomato plants, which were treated once with sulphur (Stulln WG, Andermatt Biocontrol, Grossdietwil, Switzerland) and then placed in approximately 100% rel. humidity for one week. Afterwards the cuttings were acclimatised to the same greenhouse conditions mentioned above. Young and fully unfolded leaflets were harvested from two-week old cuttings for inoculation trials. 15 leaflets were placed in miniature grow boxes (30 × 60 cm and 14 cm in height with a clear plastic cover) on paper towels wetted with distilled sterile water. A separate box was used for each pathogen and mock inoculation. All inoculations were repeated three times.
Inoculum preparation, inoculation and sampling for transcriptome analysis. A P. infestans strain K5276 that was isolated in Switzerland and kindly provided by Syngenta (Basel, Switzerland) was grown for 3-8 weeks on V8 medium (200 ml V8 Jus de Legume, Globus, Switzerland; 30 mM CaCO 3 , 1.5% (w/v) Agar, Sigma-Aldrich, Buchs, Switzerland). Sporangia were collected with 10 ml tap water, diluted to 4 × 10 5 sporangia per ml and stored in darkness for 2 h at 5 °C before inoculation. Slight shaking of the inoculum solution hindered sporangia and zoospore sedimentation. 10 μl of the suspension was applied to the abaxial leaf surface for inoculation. The inoculated leaves were stored in darkness at 16 °C for 24 h, followed by a 16/8 hour day/night regime. Mock inoculations were performed under the same conditions using tap water for inoculation.
A B. cinerea strain T4 that was kindly provided by Philippe Nicot, INRA Centre de Recherche PACA, Montfavet, France, was grown on 15 g/l malt agar (Fluka, Sigma-Aldrich, Buchs, Switzerland) plates for 3-8 weeks. Spores were harvested with 20 ml half-strength grape juice (Farmer, Landi, Dotzingen, Switzerland) and diluted to 1.2 × 10 6 spores per ml. The spore suspension was used directly for inoculation, with 10 μl drops placed on the adaxial leaf surface. The inoculated leaves were stored without light at 18 °C. Mock inoculations were performed under the same conditions using half-strength grape juice for inoculation.
O. neolycopersici that was isolated on tomato plants in the greenhouse at our institute was maintained on S. lycopersicum cv. Heinz 1706. Spores were harvested with a wet paint brush and diluted in water to a concentration of 4 × 10 4 spores per ml, which was used to directly inoculate the adaxial leaf surface with 10 μl drops. The inoculated leaves were stored at 22 °C with a 16/8 hours day/night regime. Mock inoculations were performed under the same conditions using tap water for inoculation.
Inoculation and Sampling: Each inoculum was applied as approximately eight 10-μl-drops to the abaxial or adaxial surface of 15 leaflets (folioles 53 ). Leaf disks (LD) of the inoculated sites were cut at 24 hours past inoculation (hpi) with a 5 mm-diameter cork borer for RNA-seq. LDs of eight leaflets were pooled, the remaining inoculum removed with a paper towel and frozen in liquid nitrogen. Samples were stored at −80 °C until further processing.
RNA extraction, transcriptome sequencing and quantitative real time PCR. Total RNA was extracted using NucleoSpin ® RNA Plant (Macherey-Nagel, Düren, Germany) following the manufacturer's instructions, including DNase treatment. The RNA quality and quantity was estimated using the Standard Sensitivity RNA Analysis Kit in a Fragment Analyzer (Advanced Analytical Technologies, Ames, USA) and analysed with the associated PROSize ® 2.0 v.1.3 software. Total RNA with an RNA quality number (RQN) >= 6 was sequenced with Illumina Highseq2500 v4 chemistry by GATC (Constance, Germany) using a strand-specific cDNA library from purified poly-A containing mRNA molecules.
Validating the RNA-seq-derived differential expressed genes primers for 11 DE genes (Supplementary Table S1) were designed using Primer3 54 . For qPCR total RNA was transcribed into first strand cDNA using an iScript cDNA synthesis kit (Bio-Rad, Hercules, CA, USA). Primers LSM7 and SlCBL1 that amplify genes that code for U6 snRNA-associated Sm-like protein LSm7 and Calcineurin B-like protein, respectively, were used as reference genes 55,56 . The qPCR was conducted on a LifeCycler 480 (Roche, Basel, Switzerland) with a Fast EvaGreen ® qPCR Master Mix (Biotium, Haywardm, CA, USA). For all primer pairs used in this study (Supplementary Table S1) the fast amplification protocol suggested by the master-mix provider was used, consisting of initial denaturation at 95 °C for 2' followed by 40 cycles at 95 °C for 5" and 60 °C for 30". A melting curve analysis was performed from 60-100 °C in 0.1 °C steps at the end of the run. All samples were analysed in duplicates. Data processing and statistics. To estimate the percentage of RNA-seq reads coming from the infection, the reads from all samples were aligned to the ENSEMBL 57 (release 30) genome of S. lycopersicum and to the pathogen genomes. The B. cinerea T4 58 and P. infestans genomes 59 were downloaded from the Broad Institute's website. Since no reference genome was available for O. neolycopersici, reads from the infected and corresponding mock samples were assembled using Oases 60 with a k-mer of 31. Resulting contigs were blasted against GenBank nt database using only contigs that matched fungal species for further analysis. The reads from the O. neolycopersici infected and mock-inoculated samples were mapped to the fungal contigs. For each sample, the percentage of reads that mapped the pathogen genome/contigs was calculated.
Differential gene expression analysis. Expression levels for each sample were estimated using the RSEM package (version 1.2.25) 61 in paired-end and strand-specific mode with bowtie2 62 to the ENSEMBL (release 30) annotation of the S. lycopersicum genome. The RSEM "rsem-run-ebseq" tool was used to calculate the expression pattern based on the estimated expression levels for each gene based on the EBSeq method 23 . Significant patterns were assigned with a PPDE (posterior probability of being DE) >= 99%, which corresponds to a False Discovery Rate (FDR) of 1%. Additionally a pairwise comparison of all samples was performed using the R package edgeR 22 .
The resulting genes of interest were linked to their GO terms with Blast2GO (v4.0.7) and further analysed using g:Profiler 63, 64 , STRING 24, 65 GOstat 25 for GO enrichment analysis.