The effect on the transcriptome of Anemone coronaria following infection with rust (Tranzschelia discolor).

In order to understand plant/pathogen interaction, the transcriptome of uninfected (1S) and infected (2I) plant was sequenced at 3'end by the GS FLX 454 platform. De novo assembly of high-quality reads generated 27,231 contigs leaving 37,191 singletons in the 1S and 38,393 in the 2I libraries. ESTcalc tool suggested that 71% of the transcriptome had been captured, with 99% of the genes present being represented by at least one read. Unigene annotation showed that 50.5% of the predicted translation products shared significant homology with protein sequences in GenBank. In all 253 differential transcript abundance (DTAs) were in higher abundance and 52 in lower abundance in the 2I library. 128 higher abundance DTA genes were of fungal origin and 49 were clearly plant sequences. A tBLASTn-based search of the sequences using as query the full length predicted polypeptide product of 50 R genes identified 16 R gene products. Only one R gene (PGIP) was up-regulated. The response of the plant to fungal invasion included the up-regulation of several pathogenesis related protein (PR) genes involved in JA signaling and other genes associated with defense response and down regulation of cell wall associated genes, non-race-specific disease resistance1 (NDR1) and other genes like myb, presqualene diphosphate phosphatase (PSDPase), a UDP-glycosyltransferase 74E2-like (UGT). The DTA genes identified here should provide a basis for understanding the A. coronaria/T. discolor interaction and leads for biotechnology-based disease resistance breeding.

Introduction

The genus Anemone (Ranunculaceae) harbor over 120 species, distributed over the temperate zones of both hemispheres; many of these species are cultivated as ornamentals. The poppy anemone (Anemone coronaria), of Mediterranean origin, is the progenitor of most of the cut flower, pot and garden plant varieties currently cultivated [1]. In nature, seeds produced in late spring usually germinate in autumn and offspring starts flowering the following year. Under cultivation practice, rhizomes, derived from seed by one growth cycle in the nursery, are planted by growers, after vernalization, in order to shorten the time from planting to harvest.

A major biotic constraint for Anemone producers is the rust disease, caused by the basidiomycete Tranzschelia discolor [2], which has become an aggressive pathogen in recent years, following the widespread exploitation of tetraploid cultivars of A. coronaria. The pathogen infects Prunus spp. as its primary host and some members of Ranuncolaceae as its alternate host. During the production of A. coronaria rhizomes, seedlings frequently are infected by inoculum which has developed on Prunus spp. foliage, although the infected plants remain asymptomatic until the following vegetative cycle. The disease has a major impact on flower yield and quality and finally plants became rusted and dies.

The breeding of resistant varieties of A. coronaria has been hampered by poor state of knowledge regarding the host/pathogen interaction.

Rust pathogen fungi are obligate biotrophic parasites [3]. A successful infection requires that effectors, coded by avirulence (Avr) genes, are secreted into infected tissues to repress and manipulate host defense [4]. In turn, plants possess several hundreds of resistance (R) genes that trigger strong defense responses [5]. The ability of pathogen effectors to manipulate host functions and escape R protein recognition is thought to be the key of compatibility [6]. Specific recognition is thought to be mediated by ligand receptor binding [7]. In order to survive, plants have engaged a co-evolutionary battle engendering a wide range of constitutive and inducible defenses [8]. Constitutive defenses include many preformed barriers such as cell walls, waxy epidermal cuticle, tricomes and bark. In addition, plants have developed two innate immune systems for defense [8,9]. The primary innate immunity, is driven by pattern recognition receptor (PRRs), that recognize microbe-associated molecular patterns (MAMPs) and triggers primary defense responses, such as cell wall alterations, deposition of callose and the accumulation of defense-related proteins including chitinases, glucanases and proteases [5,9]. Virulent pathogen are able to suppress basal defense activated in the primary innate immune system, developing mechanisms to escape recognition of MAMPs [8,10]. Therefore plants have developed a secondary defense response through resistance proteins (RPs) that monitor the effectors or their perturbations of host targets and often culminate in a hypersensitive response (HR). The hypersensitive response is characterized by localized cell and tissue death at the site of infection [11]. This strong defense reaction is characterized by the accumulation of reactive oxygen species (ROS), antimicrobial proteins and phytoalexins that lead to a local cellular suicide, which stops biotrophic pathogens from further growth. Plants are also protected by a mechanism called systemic acquired resistance (SAR), which is induced simultaneously with local primary and secondary immune response [12], providing durable protection against challenge infection by a broad range of pathogens [13-15]and is dependent on different plant hormones, as salicylic acid (SA), jasmonic acid (JA), ethylene (ET), abscisic acid (ABA). Next-generation sequencing (NGS) technology has revolutionized the acquisition of nucleic acid sequence and made major contributions to our understanding of genome structure, gene expression and regulation [16-18]. The RNA-Seq platform provides a direct count of the number of specific transcripts present in an mRNA sample and thus gives an informative means not only of acquiring transcriptomic sequence, but also of identifying differential transcription [19]. The accuracy of its measurement of transcript abundance is as high, if not higher than is possible using microarray technology [20-22]. An RNA-Seq variant, that consider sequencing of 3’ ends only, permit detection of rare transcript even in case of a low number of reads [23]. As a result, the approach has been widely employed to study transcription in fungal, plant and animal genomes [24,25]. The NGS FLX 454 pyrosequencing technology (Roche, Brandford, CT, USA) has been widely used for de novo sequencing and analysis of trascriptomes in non-model organisms, such as olive [26], chestnut [27], Artemisia annua [28], ginseng [29], blueberry [30] bracken fern [31] and in switchgrass [32].

Here, we report the use of FLX 454 technology to analyze the transcriptome of A. coronaria and in particular to determine what change in transcription are induced when plant is infected by T. discolor.

Materials and Methods

RNA extraction

A. coronaria plants (cv ‘Tetraelite’ blue) were grown under shade netting. Thirty healthy and thirty T. discolor infected plants were monitored throughout their life cycle for disease symptoms. Early infected plants were easily identified by plethoric vegetation, robust, erect leaf stems and thick, slightly curled leaf lamina. Leaves of infected plants were harvested as soon as plant showed disease symptoms. This time point covers leaf invasion by hyphae from the plant rhizomes under real field condition. Healthy leaves of the same age were harvested from uninfected plants (Fig. 1). Leaf tissues was snap-frozen in liquid nitrogen and stored at -80°C until required. Total RNA was isolated from 100 mg of frozen leaf using an RNeasy Plant Mini kit (QIAGEN GmbH, Hilden, Germany) and treated with recombinant DNase I (QIAGEN) within the column, following the manufacturer’s protocol. The concentration of recovered RNA was estimated using a Nanodrop 2000 device (Thermo Fisher Scientific Inc., Wilmington, DE, U.S.A.) and its integrity assessed using a Total RNA Stdsens chip (Experion system, Biorad, Hercules, CA, USA). High quality RNAs from five uninfected plants were combined to form the “1S” pool and similarly from five infected ones to form the “2I” pool.

Fig 1

Healthy and Tranzschelia discolor infected plant.

In comparison to healthy plant (A), during biotropich relationship, infected plant displays plethoric vegetation with robust leaf stems and curly leaf lamina. Flowering is strongly repressed.

454 titanium sequencing

Ten μg of RNA of each pool were sent to Eurofins MWG-Operon (Ebersberg, Germany – http://www.eurofinsdna.com/home.html) for preparation of two 3’ c-DNA libraries (“1S” and “2I”) and sequencing using GS FLX 454 Titanium system (454 Life Sciences, a Roche company, Branford, CT, USA).

Assembling, annotation and functional analysis

The Roche 454 high quality (HQ) reads generated in this study were deposited in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession SRX447797. After trimming of adapter sequences and removing reads shorter than 40 nt, libraries were mass assembled into a set of transcript contigs using CLC Genomics Workbench 5.0 operating with its default minimum identity setting of 0.8. Unigenes (contigs and remaining unique singletons) were annotated using a BLASTx search of the NCBI non-redundant protein database (nr), with the help of Blast2GO v2.5 (http//:www.blast2go.org) applying an e-value threshold of 1e-3. Blast2GO was also used to obtain gene ontology (GO) information. Sequences were annotated with respect to GO term applying the default e-value of 1e-6 and the “Augment Annotation by ANNEX” function was used to refine annotation. An InterPro search was performed for the sequences which were unsuccessfully annotated by BLASTx analysis [33]. The Plant GO-Slim algoritm was used to assign GO terms. Pathway assignment was processed with KEGG database. To identify GO categories represented differentially between the two libraries, an enrichment analysis was performed using two tail Fisher’s exact test as implemented in Blast2GO, applying a False Discovery Rate (FDR) of 0.05 and the Benjamini and Hochberg [34] Multiple Testing Correction; for this purpose, the annotated sequences of the library generated from 2I pool (contigs and singletons) was used as the test set and those of the library from the 1S pool as references set [35]. To identify genes involved in the disease response, the predicted unigene products were queried using BLASTx (as implemented in the CLC Genomics Workbench 5.0 with the default parameters) using the peptide sequences of 50 known disease resistance (R) gene [36], covering each the five major R gene classes.

Assessment of transcriptome coverage

Coverage of the transcriptome was estimated using the web-based ESTcalc tool [37]. The number of reads (610,561), the average read length of 330 nt and one run for 454 GS FLX technology determined the input parameters.

The number of eukaryotic ultraconserved orthologs (UCOs) represented in the dataset was obtained from a tBLASTx query based on the 357 Arabidopsis thaliana UCOs available at http://compgenomics.ucdavis.edu/compositae_reference.php; the chosen e-value threshold was 1e-10.

Analysis of differential transcript abundance (DTAs)

DESeq package [38] was chosen to identify gene with DTA. It integrates several statistical methods, that can estimate a theoretical replicate when an experimental one is not provided and has been routinely used [39,40].

The number of reads contributing to each contig was compared for each gene of the 1S and 2I libraries. The FDR threshold was set at 0.05.

Phylogenetic analysis and alignment of DTA genes

The selected unigenes for the alignment and phylogenetic analysis were blasted against the NCBI nr protein database, using a BLASTx search (http://blast.ncbi.nlm.nih.gov/blast/). The full length amino acid sequences with higher "Max score" and "Identity" percentage were selected for analysis.

Clustal Omega program (http://www.ebi.ac.uk/Tools/msa/clustalo/) was used to obtain sequences alignment (Gonnet matrix). The maximum likelihood (ML) method from the MEGA program (version 6.06) [41] was used to create phylogenetic trees of selected DTA genes. The reliability of each branch was tested by bootstrap analysis with 100 replications.

Validation of DTAs using qPCR

Ten plant DTAs genes (five up- and five down-regulated) putatively involved in the response to T. discolor infection, along with three fungal genes, were subjected to RT-PCR and qPCR analysis as described by Laura et al. [42]. The templates compared were the pools of uninfected and infected plants analyzed by pyrosequencing (sample A). Primers were designed using Primer 3 plus software (http://primer3plus.com/cgi-bin/dev/primer3plus.cgi) and the A. coronaria 18S rRNA gene was employed as the reference (Table 1). For each gene, three biological replicates and three technical replicates were performed. Two additional sample (B and C) each composed of five of uninfected and infected plants, not included in the libraries, were analyzed further.

Table 1

Primer sequences of differential transcript abundance (DTA) genes and Tranzschelia discolor genes used for qPCR analisys.

Contig	Sequence	Origin	Fwd primer (5’-3’)	Rev primer (5’-3’)
2476	18S Rrna (reference gene)	Plant	gagagagggggaagaaaagg	gtggaacagtgagggacaag
9288	NDR1-like	Plant	ccctaccaccgttgtttcac	cactctacccatcagcacca
4448	inner membrane protein PPF-1	Plant	atctgacgttgctgaagcct	gcaccaaaaatgagctaagcc
6057	late embryogenesis abundant	Plant	catccgttcttccagcgac	tgtgccgacttccataccc
9205	myb family transcription factor	Plant	tggcaatcaaggtcaacaac	cctcggtcttttcttcatcg
17626	XTH	Plant	tgatggctgtgaggagtctg	attctggtggcactgtaggg
3404	acidic chitinase	Plant	tggtagcactgtaggcctcg	tcttgatagccggcaaaact
12518	COMT	Plant	acttgtgaacctgctctgcc	actcttcctctgtgcgatgc
18538	metallothionein protein	Plant	cgtcaggctgtaactgtggt	tgctccaacacccatctctg
7408	PDF1.1	Plant	tgggacatggtcaggagtct	ggcggggaatttgtagttgc
21204	PGIP	Plant	cggtgatgcttcggttctgt	catttggtggtcaaggggg
21597	bax Inhibitor family protein	Fungi	ccctcgttttgcccgtca	cggcttcatcctaccttgttgc
21934	nucleoside diphosphate kinase	Fungi	cctttgcggtttgtggattg	atgtgctttgaaggggtgag
20893	inorganic phosphate transporter 1–7	Fungi	gcagcctctggaaaacttgg	gtcagcagtcccgtaaacatca

After normalization, transcript abundances were compared using the 2−ΔΔCt method [43]. The Wilcoxon-Mann-Whitney [44] test was applied as implemented in the GraphPad InStat v3.10 package (www.graphpad.com). Transcript abundance data were expressed in the form mean ± standard error (SE).

Results and Discussion

EST sequencing and assembly

The GS FLX 454 output yielded 304,487 (1S library) and 306,074 (2I library) raw reads of average length, 330 and 322 nt, respectively. The size distribution of reads is given in S1 Fig. and a summary of sequencing and assembly outcomes presented in Table 2. Respectively, 12,052 and 12,488 sequences were discarded on the basis of shortness of length (< 40bp) or low quality score (Mira version 4.0.2), resulting in the acquisition of, respectively, 292,435 and 293,586 HQ reads. These HQ reads were assembled into 27,231 contigs leaving 37,191 singletons in the 1S library and 38,393 in the 2I library. The number of contigs specific to one library was 5,802 (1S library) and 9,085 (2I library). Among the contigs, their length varied from 40 to 1367 bp (mean 377 nt) and 507 were longer than 800 nt (S2 Fig.); the length of the singletons, ranged from 40–696 nt (mean 247.5 nt) in the 1S library and from 40–767 nt (mean 248 nt) in the 2I library. Sequencing coverage, as estimated from the mean number of reads per contig [45], was 19.52.

Table 2

Overview of the sequencing and assembly.

454 pyrosequencing terms	Sequences (n)		Bases (bp)
	1S	2I	1S	2I
Sequencing
High-quality (HQ) reads	304,487	306,074	100,541,607	98,494,613
Average sequence length of HQ reads	330.2 bp	321.8 bp	-	-
N° of reads used in assembling	292,435	293,586	92,327,408	90,081,755
Contigs
Reads mass assembled in contigs	255,244	255,193	83,123,220	80,564,794
N° of contigs	27,231		10,262,491
Average length	377 bp		-	-
Range of length	40–1367 bp		-	-
N50	443 bp		-	-
Singleton
N° of singletons	37,191	38,393	9,204,188	9,516,961
Average length	247.5 bp	248 bp	-	-
Range of lenght	40–696 bp	40–767 bp	-	-
Unigenes	64,422	65,624	19,466,679	19,779,452

Transcriptome coverage

Since the genomic sequence of the A. coronaria is unavailable, the true size and composition of its transcriptome is unknown. Thus the simulation-based ESTcalc tool [37], was used to estimate the coverage of the transcriptome produced by the RNA-Seq data set. This exercise suggested that 71% of the transcriptome had been captured, with 99% of the genes present being represented by at least one read (Table 3). With respect to the UCOs, 348 of the 357 tested sequences were represented in the A. coronaria contigs.

Table 3

EST calc-based transcriptome coverage.

Input parameters	EST calc
Number of technologies	1
Technology	454 GSFLX
Library type	non-normalized
Number of Plates	0.5
Reads/Plate	610,561
Mean read lenght (bp)	330
Predicted assembly	EST calc
Total Sequence Amount (MB)	100.7
Total Assembled Sequence (MB)	21.1
Unigene count	37,822
Mean unigene length (bp)	558
Mean unigene length (longest unigene per gene, bp)	812
Singleton yield (%)	33
Percent transcriptome (%)	71
Percent of genes tagged (%)	99
Percent of genes with 90% coverage (%)	41.6
Percent of genes with 90% coverage by largest unigene (%)	29.7
Percent of genes with 100% coverage (%)	9.7
Percent of genes with 100% coverage by largest unigene (%)	9

Sequence annotation

Unigenes annotation, through use of the Blast2GO tool, showed that 50.5% of the predicted translation products shared significant homology with known protein sequences deposited in GenBank and 1.7% with hypothetical proteins, leaving 47.8% of the sequences unannotated. The proportion of sequences lacking any BLASTx alignment and shorter than 250 nt was 42.1% (S3 Fig.). Short sequences are thought to derive predominantly from the highly divergent 3’ untranslated regions (UTRs), so may account for the high proportion of the low homology sequences. An InterPro search of 13,028 unannotated contigs identified 5,891 as harbouring known protein domains.

The BLASTx positive contigs identified Vitis vinifera (grape) as the most frequently occurring species, followed by Populus tricocarpa (black cottonwood), Ricinus communis (the castor oil plant), Glycine max (soybean) and Puccinia graminis (cereal stem black rust) (S4 Fig.). The number of fungal contigs identified was 1203 (S2 Table) of which 1194 were not represented in the 1S library. The presence of eight fungal contigs (11 reads) in the 1S library is thought to reflect field-based aeciospores contamination. None of the 1S library singletons matched sequences in the P. graminis proteome.

Gene Ontology (GO) annotation

Sequences showing significant similarity to previously annotated proteins were assigned GO terms based on their associated biological processes (P), molecular functions (F) and cellular components (C). Plant specific GO slim terms were associated with 10,362 (38%) of the contigs, of which 8,433 were given an F, 6,726 a C and 6,451 a P assignation. The GO categories represented showed no significant bias and were distributed similarly to what has been described in other plant species [31,46,47].

The Predominant P categories were biosynthetic process, catabolic process, carbohydrate metabolic process and protein modification process (Fig. 2A). The C assignation of most of the contigs was to the plastids or mitochondrion, but many were associated with the ribosome (Fig. 2B). Nucleotide binding, protein binding, kinase activity and transporter activity were the major F categories present. Genes involved in the responses to stress (354), abiotic (14) and biotic stimulus (77) and signal transduction (413) are also well represented (Fig. 2C).

Fig 2

GO multilevel-Pie.

Pie chart representation of Gene Ontology classification of (A) biological process (B) cellular component (C) molecular function, using a sequences cutoff of 5.0.

Fisher’s Exact test confirmed that the distribution of GO categories differed between the two libraries. Specifically, 447 GO categories were differentially represented (S1 Table), of which 304 involved P, 67 F and 76 C. In all, 267 GO categories were under-represented in the 2I library and 180 over-represented. Among the latter, were genes encoding transferase activity (GO:0016740), hydrolase activity (GO:0016787), RNA binding (GO:0003723), kinase activity (GO:0016301), lysozyme activity (GO:0003796), chitinase activity (GO:0004568), peroxidase activity (GO:0004601) and hydro-lyase activity (GO:0016836). The 3,4% of the over-represented categories showed no unigene sequences in the reference group.

R genes homologs in A. coronaria

A tBLASTn-based search of the sequences using as query the full length predicted polypeptide product of 50 R genes [36] identified 84 contigs, along with six (1S library) and one (2I library) singletons. The 91 unigenes recognized, related to 16 R gene product (those derived from Vf1, Fls2, Pbs1, Xa21, Xa26, Rps5, Ssi4, Rpg1, Mlo, Hm1, Hs1, Cf-2, Cf-5, Cf-9, Pto and Vrgl1), distributed between the R gene classes NBS-LRR (3), LRR (2), LRR-TM (10), LRR-PK (19), PK (54) and TM (2) with one showing a high level of similarity to R genes carrying a Toxin reductase domain (Table 4).

Table 4

Anemone coronaria unigenes with homology to known Resistance genes.

Unigenes ID	R genes	Structure	E-value	Protein ID	Unigenes ID	R genes	Structure	E-value	Protein ID
19281	Cf-2	LRR-TM	7,29E-12	AAC15779	8128	Pbs1	PK	5,95E-19	ABR46085
26249	Cf-2	LRR-TM	1,21E-14	AAC15779	18170	Pbs1	PK	2,19E-16	ABR46085
5613	Cf-2	LRR-TM	1,57E-14	AAC15779	9610	Pbs1	PK	8,74E-10	ABR46085
21204	Cf-2	LRR-TM	2,19E-14	AAC15779	15776	Pbs1	PK	9,34E-10	ABR46085
11915	Cf-5	LRR-TM	3,84E-21	AAC78591	15403	Pbs1	PK	1,25E-18	ABR46085
3264	Cf-5	LRR-TM	4,46E-16	AAC78591	24986	Pbs1	PK	2,67E-18	ABR46085
4230	Cf-5	LRR-TM	3,26E-26	AAC78591	5452	Pbs1	PK	3,67E-45	ABR46085
4420	Cf9	LRR-TM	1,70E-17	CAA05274	15392	Pbs1	PK	5,29E-23	ABR46085
21205	Fls2	LRR-PK	1,08E-12	BAB11088	20766	Pbs1	PK	6,85E-15	ABR46085
2470	Fls2	LRR-PK	7,78E-12	BAB11088	15954	Pbs1	PK	3,22E-10	ABR46085
16108	Fls2	LRR-PK	1,35E-10	BAB11088	15959	Pbs1	PK	3,26E-14	ABR46085
191	Fls2	LRR-PK	8,25E-18	BAB11088	917	Pbs1	PK	8,35E-16	ABR46085
9548	Fls2	LRR-PK	1,67E-19	BAB11088	9947	Pbs1	PK	7,38E-26	ABR46085
19550	Fls2	LRR-PK	3,72E-15	BAB11088	GG9DABN02FRDGE	Pbs1	PK	1,79E-14	ABR46085
26189	Hm1	tox reduc	1,03E-10	AAC04333	GG9DABN02HL1AG	Pbs1	PK	2,89E-20	ABR46085
22173	Hs1	LRR-TM	2,84E-19	AAW03319	19538	Pto	PK	1,76E-10	AAB47421
GG9DABN02F9DDG	Hs1	LRR-TM	2,92E-13	AAW03319	7481	Pto	PK	2,38E-13	AAB47421
17299	Mlo	TM	1,55E-13	CAB06083	2969	Pto	PK	6,53E-23	AAB47421
GG9DABN02G8QRM	Mlo	TM	5,53E-12	CAB06083	13334	Pto	PK	4,86E-10	AAB47421
14622	Pbs1	PK	5,85E-14	ABR46085	2466	Pto	PK	4,21E-14	AAB47421
14845	Pbs1	PK	1,44E-23	ABR46085	9388	Pto	PK	1,01E-11	AAB47421
3503	Pbs1	PK	4,16E-10	ABR46085	7553	Pto	PK	4,19E-12	AAB47421
12107	Pbs1	PK	6,30E-13	ABR46085	14902	Pto	PK	8,90E-13	AAB47421
4272	Pbs1	PK	1,37E-21	ABR46085	8930	Pto	PK	1,17E-15	AAB47421
25490	Pbs1	PK	3,40E-15	ABR46085	26150	Pto	PK	3,55E-11	AAB47421
GG9DABN02F3PO5	Pbs1	PK	3,98E-10	ABR46085	18600	Pto	PK	1,21E-12	AAB47421
8545	Pbs1	PK	3,46E-12	ABR46085	GG9DABN02HEM5G	Pto	PK	8,53E-16	AAB47421
19451	Pbs1	PK	1,03E-15	ABR46085	14229	Rps5	NBS-LRR	1,11E-14	AAC26126
16933	Pbs1	PK	1,48E-19	ABR46085	2985	Ssi	NBS-LRR	8,07E-12	AAN86124
17238	Pbs1	PK	1,17E-12	ABR46085	18516	Vf1	LRR	1,89E-14	CAC40825
22527	Pbs1	PK	1,23E-16	ABR46085	19604	Vf1	LRR	6,83E-17	CAC40825
8763	Pbs1	PK	2,31E-25	ABR46085	GG9DABN02HVY07	VRGL1	NBS-LRR	8,78E-12	AAF19148
13820	Pbs1	PK	2,04E-24	ABR46085	919	Xa21	LRR-PK	6,04E-11	AAC80225
10598	Pbs1	PK	1,55E-10	ABR46085	16298	Xa21	LRR-PK	4,41E-10	AAC80225
6752	Pbs1	PK	7,86E-13	ABR46085	4299	Xa21	LRR-PK	8,13E-10	AAC80225
13774	Pbs1	PK	1,21E-12	ABR46085	1845	Xa21	LRR-PK	9,17E-11	AAC80225
14120	Pbs1	PK	5,07E-20	ABR46085	11988	Xa21	LRR-PK	1,40E-21	AAC80225
12619	Pbs1	PK	1,15E-13	ABR46085	18834	Xa21	LRR-PK	1,87E-13	AAC80225
12545	Pbs1	PK	5,78E-12	ABR46085	9452	Xa21	LRR-PK	5,85E-13	AAC80225
13256	Pbs1	PK	5,97E-11	ABR46085	9691	Xa21	LRR-PK	3,10E-10	AAC80225
17036	Pbs1	PK	5,79E-15	ABR46085	8221	Xa26	LRR-PK	6,73E-13	ABK51312
17518	Pbs1	PK	1,34E-19	ABR46085	16817	Xa26	LRR-PK	1,35E-18	ABD36512
22913	Pbs1	PK	1,19E-12	ABR46085	5926	Xa26	LRR-PK	2,17E-13	ABD36512
14431	Pbs1	PK	5,44E-14	ABR46085	11237	Xa26	LRR-PK	3,58E-11	ABK51312
15636	Pbs1	PK	1,84E-11	ABR46085	15836	Xa26	LRR-PK	1,51E-16	ABD36512
8046	Pbs1	PK	5,37E-14	ABR46085

Unigene products were queried using BLASTx as implemented in the CLC Genomics Workbench 5.0 with the default parameters.

Identification of differential transcript abundance (DTAs)

In all, 305 DTA genes were identified by comparing transcript abundances between the two libraries. Of these, 253 were present in higher abundance in the 2I library and 52 in lower abundance. In the former set, 234 were not detected in the 1S library and their read number per transcript in the 2I library varied from 15 to 456; similarly 25 of the down-regulated DTAs were not represented in the 2I library, whereas they were present in 15–49 copies in the 1S library (S2 Table). About a half (128) of the DTAs present in the 2I library were of fungal origin, 49 were clearly plant sequences and 7 could have encoded a fungal or a plant protein; the remaining DTAs could not be functionally assigned using BLAST.

Fungal up-regulated genes

Among the likely fungal sequences, 118 had homologs in either P. graminis f. sp. tritici or Melampsora larici-populina and 78 were associated with a likely function (S2 Table). Of the 30 genes encoding ribosomal proteins (RPs), 26 were likely to have been of fungal origin, reflecting the active protein synthesis exhibited by fungi during the early phase of infection [48]. Genes hydrolytic enzymes acting on plant biopolymers (cellulase), proteinase (subtilase-type proteinase psp3, vacuolar protease A, proteasome subunit 1) and several carbohydrate-active enzymes (glycoside hydrolase, glyceraldehyde-3-phosphate dehydrogenase, enolase, glucose-repressible protein) were well represented, as would be predicted since the invading fungus penetrates the host cells by degrading enzymes [49,50]. Apart from these, a fungal chitinase gene was recognized; this enzyme is used to remodel fungal cell wall during infection, either to promote hyphal invasion and/or to avoid recognition by the host’s defense system [51]. Other strongly represented fungal genes encoded histones, an argonaute-like protein, thiamine synthesis, a mitochondrial thiazole synthesis enzyme and the ubiquitin-conjugating enzyme E2, as was also the case during the infection by rust of both Populus sp. and wheat [49]. P. triticina genes encoding a thiamine synthesis protein and a cyclophilin are also induced in planta, as reported by Thara et al. [48]. The first protein is a cofactor controlling the activity of several enzymes involved in the central carbon metabolism [52], whereas cyclophilin is involved in a wide variety of cellular processes, including the response to abiotic stress, the control of cell cycling, the regulation of calcium signaling and control of transcriptional repression [53-55]. The transcription of a gene encoding a thaumatin-like protein (TLP) may reflect the fungus’ attempt to interfere with the host’s defense signaling apparatus [56,57]; the presence of this protein has been noted in the Melampsora secretome [58,59]. The phylogenetic tree of TLP resolved the entries into two major branches. One includes rust fungal proteins only (P. graminis and M. larici-populina and T. discolor), the other groups proteins of all other fungal taxa. T. discolor TLP is well separated from P. graminis and M. larici-populina proteins with a bootstrap of 99% (Fig. 3A). Amino acid sequence of T. discolor TLP shares with rust TPLs the 16 conserved cysteine residues that characterize the large type TLPs [60].

Fig 3

Phylogenetic tree and aligment of thaumatin-like protein.

(A) Anemone coronaria thaumatin-like protein (contig 13789) cluster with Puccinia graminis and Melampsora larici-populina proteins. TLPs of no-rust fungal taxa cluster in a separate group. Bootstrap values are indicated in relevant nodes. Arabidopsis thaliana TLP was used as out-group. (B) Amino acid sequence alignment of five rust TLPs. The conserved 16 cystein residues are highlighted in the boxes.

Plant up-regulated genes

Among the A. coronaria genes up-regulated in the 2I library, 17 had homologs in V. vinifera, five in A. thaliana and the remainder in another species. The BLAST assignment of these DTAs is given in S2 Table. In what follows, the function of plant DTAs with potential relevance for host/pathogen interaction is explored.

Ribosomal proteins (RPs) genes. Four genes encoding 40S or 60S RPs were up-regulated in the infected plants, suggesting not only an increased level of protein synthesis induced by the infection process, but also the promotion of a suite of extra ribosomal activities, such as DNA repair, apoptosis, inflammation, tumorigenesis and transcriptional regulation [61]. The ribosomal protein S3 (RPS3), which is a component of the eukaryotic 40S ribosome is known to be involved in certain host–pathogen interactions [62].

Genes involved jasmonate (JA) signaling. Among the genes up-regulated in the infected plants were six which encode various proteins involved in JA signaling; these included those encoding 12-oxophytodienoate reductase 2 [63] and acyl-CoA oxidase [64], which are both components of JA synthesis, two JA-induced proteins, one JA signaling repressor (TIFY 3B, also known as JAZ12) and strictosidine synthase 1 [65]. JAZ proteins were degraded on perception of jasmonyl-isoleucine (JA-Ile, active form of JA) allowing the JA-Ile dependent gene expression [66,67]. Strictosidine synthase 1, a key enzyme in alkaloid biosynthesis, was induced by plant defence signalling compounds, such as salicylic acid (SA), ethylene and methyl jasmonate [65]. The Kegg analysis of the jasmonate biosynthetic pathway proves that six genes, in addition to the two up regulated, were identified by the trascriptome sequencing (Fig. 4).

Fig 4

KEGG pathway of α-linoleic acid metabolism.

Anemone coronaria transcrips involved in jasmonic acid metablolic pathway are highlighted by grey tone. 12-oxophytodienoate reductase 2 (EC: 5.3.99.6) and acyl-CoA oxidase (ACX) are upregulated during Tranzschelia discolor infection.

Among R genes identified by tBLASTn analysis, only polygalacturonase inhibitor-like (PGIPs), which harbor a leucine rich repeat and a transmembrane domain (LRR-TM) was significantly up-regulated. PGIPs can inhibit fungal endopolygalacturonases (PGs) which are responsible for breaking down the host cell wall and their encoding genes are typically induced by pathogen infection [68]. The levels of PGIP were correlated with an increased resistance to fungi in raspberry fruits [69], in older bean hypocotyls [70] and in tomato transgenic plants [71].

A. coronaria PGIP clusters together with Monocotyledon. This unexpected result draws a parallel with the phylogenetic classification of the species into the early diverging Eudicotyledon clade [72]. PGIP of core Eudicotyledon clusters into two well separated subgroups (Fig. 5).

Fig 5

Phylogenetic tree of PGIP.

Anemone coronaria PGIP (contig 21204) clusters with those of Monocotyledon species. Dicotyledon PGIPs cluster in two separate groups. Bootstrap values are indicated in relevant nodes. Entamoeba invadens PGIP protein was used as out-group.

Genes encoding pathogenesis related protein (PR). Ten PR genes were up regulated in the infected plants; they encoded either a chitinase (three genes), a bacterial-induced peroxidases (two genes), a defensin-like protein 13, a major latex protein (MLP28), an S-norcoclaurine synthase-like (NCS) enzyme, a thionin and a metallothionein. Chitinases are an important group of PR proteins because chitin is the major component of many fungal cell walls [73,74]. Duplessis et al. [75] have shown that the early expression of chitinase is needed for an incompatible Populus-Melampsora interaction.

Phylogenetic analysis resolved fungal and plant family 18 chitinases [76] into two main branches. The first branch contains plant chitinases, comprehensive of A. coronaria contigs 3404, 3405 and 12644. The second branch bring together fungal chitinases and include T. discolor contigs 26906 and 20409, that is strongly over expressed during plant infection (Fig. 6).

Fig 6

Phylogenetic tree of chitinases.

Analysis resolves family 18 chitinases into two main branches: the first includes plant chitinases and the second bring together fungal chitinases. Escherichia coli chitinase was used as out-group.

The large family of peroxidase represent enzymes [77] which contribute to plant disease resistance in several ways: they act to strengthen the host cell wall via deposition of lignin, which acts as a physical barrier against pathogen ingress [78] and also produce toxic radicals such as hydrogen peroxides [79,80]. The defensin-like protein 13 (PDF1.1) belong to a family of antimicrobial peptides which are intimately involved in determining innate immunity [81,82]. MLP28 and NCS were homologous to PR10 proteins, that are thought to participate in the defense of plants against microorganisms and fungi [83]. The MLP protein family has been associated with pathogen defense, although how they act remains unknown. NCS catalyzes the first committed step in the synthesis of benzylisoquinoline alkaloids [84]. Thionin (PR 13) is a well studied compound known to be able to permeate pathogen membranes [85]; the presence of these compounds is frequently induced in the leaf and they are present at high levels in floral tissue [86]. Finally, the metallothioneins are small cystein-rich proteins involved in correcting for imbalances in metal ions and the regulation of homeostasis under various stresses. Their participation in plant defense is thought to involve the induction of reactive oxygen species (ROS) and the suppression of ROS scavenging enzymes [87-89]. Nishimura et al. [90] have recently proposed that, they could also be used by the plant to control the synthesis of pathogen toxins via inhibition of zinc absorption by the pathogen.

Other up-regulated genes putatively involved in defense response. Nine other genes associated with the defense response were up-regulated in the infected plants. These encoded caffeic acid 3-O-methyltransferase (COMT), Cytochrome P450 (CYP 450), Early responsive to dehydration (ERD), Flavonol synthase (FLS), Heat shock proteins (HSPs), Lipid binding protein (LTPs), SNARE-interacting protein KEULE (KEU) and UDP-glucose transglucosylase-like protein. Tremblay et al. [91] used the up-regulation of COMT genes (their product is an important component of phenylpropanoid synthesis) as a marker for the activation of the plant defense response. CYP 450 contributes to oxidative metabolism and the production of ROS and is reportedly involved in the hypersensitive response (HR) to pathogen infection [92]. Some CYPs participate in the synthesis of the defense-associated compounds: lignin, phytoalexins and anthocyanins [93]. The ERD gene family comprises at least 21 members in Arabidopsis and have been identified as part of the immediate response to drought stress. Altering the level of ERD15 transcript not only had an effect on the plant’s abiotic stress tolerance but also on its level of disease resistance [94]. FLS converts both flavanones and dihydroflavonols to their related flavonols; the enzyme is a bifunctional dioxygenase, with certain hydroxylation and desaturation activities [95,96]. While many studies have indicated a role for flavonoids in disease resistance, the multi-functionality of these compound complicates the interpretation of results. [97]. Two HSP genes were up-regulated: HSP90 product has a role in signal transduction during the plant defenses response [98], while HSP23.6 accumulates during the systemic infection of Actinidia chinensis with P. syringae pv. Actinidiae [99]. Silencing of HSP90 in N. benthamiana compromises not only induction of the HR, but also non host resistance [100]. HSP90, in conjunction with other proteins, is also known to modulate N gene-mediated resistance to Tobacco mosaic virus in tobacco and RPS2 and RPM1-mediated resistance to P. syringae in Arabidopsis [101,102]. LTPs are able to transfer phospholipids between membranes and to bind fatty acids in vitro and are putatively involved in cutin synthesis, surface wax formation, defenses against pathogen and adaptation to environmental changes [103]. KEU interacts with the SNARE domain present in certain genes active in plant defense [104]. Loss-of-function of gene encoding SNARE enable elevated levels of host cell entry either by non adapted fungal species and delay in the formation of localized cell wall appositions [105]. UDP-glucose transglucosylase is thought to be involved in the synthesis of cell wall polysaccharide [106] and is active in grapevine plants exposed to pathogen infection [107]. A total of 16 other genes with no known involvement in pathogen defense were also represented by enhanced transcription the 2I library (S2 Table).

Down-regulated genes

Among the down-regulated genes, 32 were of plant origin, one of viral origin and one shared homology with both bacterial and plant proteins. The remaining 18 genes gave no BLAST hit. The genes had homologs in A. thaliana (nine), V. vinifera (eight) and R. comunis (four), with the other nine related to genes from G. max, P. tricocarpa, M. truncatula (S2 Table). In what follows, the function of plant DTAs with potential relevance for host/pathogen interaction is explored.

Cell wall associated genes. Six genes encoding components of the constitutive defense response were down-regulated in the infected plants. Two involved cell wall-associated hydrolases which act to degrade and reorganize the cell wall [108,109]; one was a cellulose synthase-like protein (CesA superfamily) which synthesizes cellulose and is required for secondary cell wall formation, one was a xyloglucan endotransglucosylase (XTH), one was a white-brown-complex (WBC) ATP-binding cassette (ABC) transporters family and one was a protein containing a galactose-binding domain-like fold (lectins). Affecting the cell wall’s integrity by inhibiting cellulose synthesis induces the activation of a number of host defense mechanism designed to produce an environment enriched with respect to antimicrobial compounds [110]. XTH restructures and loosens the xyloglucan network in the cell wall, thereby enabling cell expansion [111]. In A. thaliana, Gruner et al. [112] have shown that several genes encoding XTH, arabinogalactans, expansin- and extension-like proteins and polygalacturonase are all strongly down-regulated during the development of the systemic acquired resistance (SAR) process. Certain ABC transporters are known to be important for assuring the movement of cutin monomers [113] and others in the resistance to a number of fungal pathogens in wheat [114]. Together with other defense genes work in a sequential and concerted manner to result in a hypersensitive response to Puccinia striiformis infection [115]. Lectins, which contain a galactose-binding domain-like fold, act to bind specific ligands (such as, for example, cell surface-attached carbohydrate) and represent the only plant proteins capable of recognizing and binding the glycol-conjugates present on the outer surface of bacteria and fungi [116]. While the down regulation of CesA superfamily and XTH may activate certain defense pathways, the suppression of genes encoding cell wall associated hydrolases, ABC transporters and lectins is quite conceivably one of the means whereby T. discolor overcomes the host’s constitutive defense machinery.

Two genes encoding a product with significant homology to NON-RACE-SPECIFIC DISEASE RESISTANCE1 (NDR1), a plasma membrane-localized protein were down regulated. NDR1 is involved in the maintenance of the integrity of the cell wall/plasma membrane connection and represents a key signaling component during pathogen infection [117]. In Arabidopsis a member of the CC-NBS-LRR R protein require the signaling gene NDR1 for full activity [118]. Alignment of A. coronaria NDR1 deduced protein with the NCBI nr database proteins identifies motif 2 and 3 of the three NDR1/ HIN1-like (NHL) protein superfamily [119,120]. Motif 1 was not covered by the A. coronaria sequence (Fig. 7). Phylogenetic analysis together with down regulation during T. discolor infection provide evidence that A. coronaria NDR1 genes are credible candidate of the fungus to establish biotrophic relationship.

Fig 7

Aligment of NDR1 proteins.

Anemone coronaria contigs 6057 and 9288 are aligned with selected member of NDR1 proteins. Motif 1, motif 2 and 3 of NDR1/HIN-like (NHL) protein superfamily are highlighted in the boxes.

Other down-regulated genes putatively involved in defense response. Five genes encoding defense-associated proteins were down-regulated in the infected plants: one was a myb transcription factor and the others encoded a presqualene diphosphate phosphatase (PSDPase), a UDP-glycosyltransferase 74E2-like protein (UGT), a peptide transporter with homology to PTR3 and a secologanin synthase-like protein. The large family of myb factors includes many involved in regulating the defense response [121]. PdMYB3, for example, is more strongly activated in disease susceptible than in disease resistant Prunus domestica cultivars [122]. Fukunaga et al. [123] have shown that a human PSDPase (which converts PSDP to a monophosphate form) is important for maintaining cell function in the face of disease pressure. In the A. coronaria / T. discolor interaction, the down-regulation of the gene encoding PSDPase would likely shift the PSDP pathway in the direction of oxidosqualene, the precursor of membrane sterols, brassinosteroids, saponins and other defense compounds [124]. Most pathogen-induced SA is glycosylated by UGT to form the non-toxic SA 2-O-β-D-glucoside. The combination of SA methylation, amino acid conjugation and glycosylation forms an intimate part of the plant defense response [125,126]. A UGT loss-of-function mutant has been shown to express an enhanced level of SAR [127]. PTR3 is regulated by both SA and JA. Its A. thaliana homolog AtPTR3 is induced by the presence of the P. syringae pathogen [128], while loss-of-function mutants show accentuated susceptibility to both Erwinia carotovora and P. syringae. Secologanin synthase catalyzes the oxidative cleavage of loganin into secologanin [129], a component of the terpenoid indole alkaloids proposed to be involved in plant defense [130]. The down-regulation of myb, PSDPase and UGT is suggestive of the activation of host defense against T. discolor, while the down-regulation of PTR3 and the gene encoding the secologanin synthase-like protein may reflect pathogen growth and the establishment of a compatible host/pathogen interaction.

Validation of DTA by real-time quantitative PCR (qPCR)

The qPCR analysis based on ten plant and three fungal target results (unigenes) confirmed that the 3’ sequencing of non-normalized libraries was informative with respect to recognizing DTAs. qPCR data of sample A (plants analyzed by pyrosequencing) were compared with qPCR data of samples B and C, each composed of five distinct uninfected and infected plants. Significant differences in expression was observed for the 13 genes tested (Fig. 8 and Table 1). The up regulated genes showed the same behavior, whereas down regulated genes varied significantly likely as a result of sample bias at low expression levels. In addition, divergences in level of gene expression may reflect time points of the plant/pathogen interaction or the genetic eterogenicity in A. coronaria population.

Fig 8

Gene expression in Anemone coronaria infected with Tranzschelia discolor.

Expression analysis was conducted among sample A (plants analyzed by pyrosequencing) and samples B and C, each composed of five distinct uninfected and infected plants. Ten plant DTAs genes (five up- and five down-regulated) putatively involved in the response to Tranzschelia discolor infection and three fungal genes, were tested. The data were normalized using Anemone coronaria 18s rRNA gene as the reference. Expression analysis was performed in triplicate on three biological replicates. Transcript abundance data were expressed in the form mean ± standard error (SE).

Conclusions

Until now, the amount of genomic information for A. coronaria in the public domain has been limited to one EST and 12 DNA sequences. This has now been rectified by the acquisition of 600,000 cDNA sequences, assembled into over 27,000 contigs. The estimated coverage of the gene content of the species was 71%, with almost all genes being represented by at least one read. ESTcalc and UCO analysis also estimated that almost all gene were represented by at least one read. Taken together these data demonstrate the potential of 3’ sequencing, although an half 454 plate only was sequenced.

Biotrophic fungi require the presence of living host tissue for their survival. Rusts such as P. graminis, Melampsora spp. and T. discolor are obligate biotrophs which often require two phylogenetically non-related hosts [49]. They have evolved specialized structures, haustoria, formed within host tissue to efficiently acquire nutrients and suppress host defense responses [131].

In cultivated A. coronaria susceptible plants, a compatible interaction occurs when seedling challenge T. discolor teliospores formed on Prunus leaf. In this phase the pathogen overcome constitutive defenses including many preformed barriers such as cell walls and waxy epidermal cuticle The A. coronaria transcriptome included the products of 16 of the 50 R genes described in A. thaliana [36], but in infected A. coronaria leaf tissue, only one of them was up-regulated. In the meanwhile two NDR1 genes involved in activation of CC-NBS-LRR R genes were down regulated. As previously reported for several plant / rust interaction [6], the ability of T. discolor effectors to escape A. coronaria R protein recognition and activation is likely the key of compatibility. During leaf colonization by fungal hyphae, field grown A. coronaria plants, activate their own immune systems and overexpress PR proteins as chitinases, involved in degradation of fungal cell wall chitin, peroxidases that may have a role in inhibiting the hyphal extention, and several additional protein (defensin-like, metallothionein, MLP-like protein 28-like, S-norcoclaurine synthase-like, thionin precursor).

The response of the plant to fungal invasion includes in addition the up-regulation of several genes associated with cellular defense. Some of these, encoding peroxidase, CYP 450, superoxide dismutase Cu/Zn chloroplast and metallothionein are involved in HR that lead to cell death and stops biotrophism. SAR is an important component of the defensive armoury of plants; it provides protection against infection by a broad range of pathogens [11,13-15,132]. The phenomenon is co-ordinated by various phytohormones. Several genes involved in JA signaling were induced in the infected A. coronaria plants, which suggests the activation of SAR, as does the induction of the genes ERD, PDF1.1, ABA glucosyl transferase and the down-regulation of XTH and UGT. Taken together these data show that either constitutive or R gene mediated defense are overcome in A. coronaria by T. discolor. A. coronaria activate both primary and secondary immune system that trigger HR. SAR is induced simultaneously. Despite plant reaction T. discolor strongly affect A. coronaria gene expression to support mating, sporulation and completing its life cycle. Transcriptome sequencing is a convenient choice to investigate a complex traits as plant pathogen interaction despite the wide genome size of A. coronaria (165.28 MB corresponding to 137X A. thaliana). To fulfill an exhaustive set of knowledge, transcriptome data provided in this work need to be implemented with sequencing of full coding and regulatory regions together with an analysis of RNA interference. Next-generation re sequencing of selected genomic regions of a large amount of accession represents a powerful approach to identify the complete spectrum of DNA sequence variants [133]. This technology is a powerful approach to discover resistance alleles in candidate genes selected among the DTA s during A. coronaria / T. discolor interaction. A short cut strategy to bread resistant genotype was proposed after the advent of genetic engineering [134]. Coding sequences of the differentially over expressed A. coronaria gene can be expressed at early stage of pathogen infection by constitutive, inducible or tissue specific promoters to effectively counteract the disease. On the other side, silencing of genes that are activated directly by pathogen effectors or indirectly by the guardee proteins may result in an attenuated virulence. Targeted mutagenesis is the most recent tool to disrupt Avr gene targets [135] and to confer new recognition specificities [136]. The major constrain for utilization of genetic engineering and targeted mutagenesis in A. coronaria, is the lack of a reliable transformation method and transient system for expression of nucleases respectively.

The Ranunculaceae family belongs to an ancient eudicotyledonous clade [72] which includes a number of both ornamental and medicinal species. The present study represents the first analysis of the transcriptome of such an early diverging species [137]. The identification of gene sequences of the pathogen T. discolor will enable its interaction with its primary host (Prunus spp.) to be investigated: the latter genus of trees and shrubs is much utilized both for its fruit and flowers. The DTA genes identified here should provide a basis for understanding the A. coronaria / T. discolor interaction and leads for biotechnology-based disease resistance breeding.

GO categories differentially represented between the 2I (test set) and 1S (references set) libraries.

(XLS)

Click here for additional data file.

Genes differentially expressed on the base of transcript abundance between the 2I and 1S libraries.

Sheet one: up regulated genes of Tranzschelia discolor; sheets two and three: up and down regulated genes of Anemone coronaria respectively.

(XLSX)

Click here for additional data file.

Size distribution of the 454 raw reads.

2I represents infected library and 1S represent uninfected library.

(TIF)

Click here for additional data file.

Size distribution of the contigs.

The contigs were mass assembled from the two libraries; the mean length is 377 nt.

(TIF)

Click here for additional data file.

Size distribution of sequences with or without BLASTx.

The 50.5% of the predicted translation products shared significant homology with known protein sequences deposited in GenBank and 1.7% with hypothetical proteins, leaving 47.8% of the sequences unannotated. The proportion of sequences lacking any BLASTx alignment and shorter than 250 nt was 42.1%.

(TIF)

Click here for additional data file.

Top-hit species distribution.

Vitis vinifera (grape) is the most frequently occurring species, followed by Populus tricocarpa (black cottonwood), Ricinus communis (the castor oil plant), Glycine max (soybean) and Puccinia graminis (cereal stem black rust).

(TIF)

Click here for additional data file.