Comparative de novo transcriptome profiles in Asparagus officinalis and A. kiusianus during the early stage of Phomopsis asparagi infection

Asparagus kiusianus, an important wild relative of cultivated asparagus (A. officinalis), exhibits resistance to stem blight disease caused by Phomopsis asparagi. However, the mechanisms underlying this resistance are not understood and no transcriptomic or genetic resources are available for this species. De novo transcriptome sequencing of A. officinalis and A. kiusianus stems was performed 24 h after inoculation with P. asparagi. In total, 35,259 and 36,321 transcripts were annotated in A. officinalis and A. kiusianus, respectively. 1,027 up-regulated and 752 down-regulated transcripts were differentially expressed in the two Asparagus species. RNA sequencing data were validated using quantitative real-time reverse transcription PCR. Several defense-related genes including peroxidase 4, cationic peroxidase SPC4-like, pathogenesis-related protein-1-like, and jasmonic acid biosynthesis and signaling-related genes including phospholipase D alpha 1, 12-oxophytodienoate reductase and jasmonate-induced protein 23 KD were up-regulated in A. kiusianus relative to A. officinalis. In addition, infected A. kiusianuns exhibited a substantial increase in jasmonic acid and methyl jasmonate relative to A. officinalis. Peroxidase activity was significantly elevated in infected A. kiusianus compared with infected A. officinalis. Our transcriptomic database provides a resource for identifying novel genes and molecular markers-associated with Phomopsis disease resistance and will facilitate breeding and improvement of cultivated asparagus varieties.


Results
Stem blight disease occurrence in wild and cultivated Asparagus species. Differences in disease occurrence between susceptible cultivated A. officinalis 'Mary Washington 500W' and resistant wild A. kiusianus (AK0501 strain) were observed after inoculation with P. asparagi spores under greenhouse conditions. Primary signs of infection, seen as dark brown lesions were first visually detected on susceptible A. officinalis at 7 days post-inoculation. The fungus then spread through the stem, resulting in a fully diseased stem by 14 days post-inoculation (Fig. 1a). By contrast, typical disease symptoms were not seen on resistant A. kiusianus and the fungus was unable to spread (Fig. 1a), demonstrating that A. kiusianus exhibited resistance to P. asparagi. Disease Values represent the mean ± standard error (SE) of three independent replicates (n = 3). Asterisks indicate significant difference between AOI and AKI plants as determined by Student's t-test (***P < 0.001).
Transcriptome sequencing and de novo assembly. To identify genes involved in the response to stem blight disease, stems of susceptible cultivated A. officinalis 'Mary Washington 500W' and resistant wild A. kiusianus (AK0501 strain) were inoculated with P. asparagi. Samples were collected from SDW-treated control stems and inoculated stems 24 h after infection. Four RNA-Seq libraries were generated: A. officinalis SDW-treated control (AOC), A. officinalis-inoculated (AOI), A. kiusianus SDW-treated control (AKC), and A. kiusianus-inoculated (AKI). Libraries were sequenced using an Illumina HiSeq 2500 platform. In total, 44,65 and 41,55 Gbp nucleotide sequence were generated from A. officinalis and A. kiusianus, respectively. After removing adapter sequences and low-quality sequences, 94,440,525; 101,235,105; 90,926,718; and 94,362,439 high-quality cleaned reads with 100 bp were obtained from AOC, AOI, AKC, and AKI, respectively. Details of the Asparagus transcriptome data are shown in Supplementary Table S1. The high-quality reads were de novo assembled into unique transcripts using Trinity software. In total, 206,164 and 213,950 assembled transcripts (average length 973 bp) were obtained from A. officinalis and A. kiusianus, respectively. Many transcripts were short; nevertheless, 64,094 and 66,330 transcripts were assembled from A. officinalis and A. kiusianus, respectively, that were >1,000 bp in length. The quality of transcriptome assemblies was assessed, and the length distribution of the assembled transcripts in both Asparagus species is shown in Fig. 2a. The GC-content of the assembled Asparagus transcripts was similar in both species, with a distribution peak at ~35% GC (Fig. 2b).
Differentially Expressed Genes (DEGs). Gene expression profiles were derived from the RNA-Seq data, and normalized DEGs were determined. To find genes that were induced by stem blight disease, pairwise comparisons were made between AOI and AOC and between AKI and AKC. A total of 1,779 differentially expressed transcripts were identified, of which 1,027 were up-regulated and 752 were down-regulated, in the two Asparagus species after inoculation with P. asparagi (Fig. 4a). Of the 1,027 up-regulated transcripts, 515 were up-regulated only in A. kiusianus, 352 were up-regulated only in A. officinalis, and 160 were up-regulated in both Asparagus species (Fig. 4b,c). These results demonstrated that the gene expression responses to P. asparagi infection differed between the two Asparagus species. GO term enrichment analysis indicated that genes involved in a range of biological processes such as transcriptional regulation, stress and defence response, protein kinase activity, phenylpropanoid and hormone biosynthesis and signaling, and cell wall assembly and organization were significantly enriched in the set of genes that was up-regulated in wild A. kiusianus but not in cultivated A. officinalis (Fig. 5).

Metabolic pathways by KEGG analysis of differentially up-regulated genes in
qRT-PCR. Fifteen up-regulated transcripts related to stress and defence, hormone biosynthesis and signaling, and ribosome were selected for validation of RNA-Seq data using qRT-PCR. Prior to qRT-PCR analysis, RT-PCR amplification was performed to confirm the suitability of the primer pairs. Single amplicons of the expected size were produced for each primer pair ( Supplementary Fig. S1). The qRT-PCR results confirmed that the selected transcripts were all strongly up-regulated in AKI plants compared with AKC plants and down-regulated in AOI plants compared with AOC plants (Fig. 6a). In addition, the qRT-PCR results positively correlated (r = 0.78) with the RNA-Seq data (Fig. 6b), validating the sequencing results.

Phytohormone accumulation in wild and cultivated Asparagus species 24 h post-inoculated
with P. asparagi. Because JA biosynthesis and signaling-related genes were markedly elevated in AKI plants relative to AOC plants, we examined the JA, MeJA, SA and ABA levels in the AOC, AOI, AKC and AKI plants 24 h post-inoculation using LC-MS/MS (Fig. 7). JA and MeJA exhibited a significant (P < 0.001) increase at 65.92and 14746.97-fold, respectively in AKI plants in comparison with AKC plants. However, there was no significant differences between AOI and AOC plants (Fig. 7a,b). SA levels exhibited non-significant increase in the two infected Asparagus species relative to control plants (Fig. 7c). In addition, ABA levels displayed a significant increase (P < 0.05) at 1.98-fold in AOI plants relative to AOC plants (Fig. 7d). However, non-significant differences were observed between AKI and AKC plants (Fig. 7d).
Effects of P. asparagi on total protein content and stress-related enzyme activities in wild and cultivated Asparagus species. In A. kiusianus, total protein content was significantly higher (1.52-fold increase; P < 0.008) in plants inoculated with P. asparagi (AKI) than in control plants (AKC). Conversely, in A. officinalis, although protein content was slightly higher in inoculated plants, no significant difference in total protein content was seen between inoculated (AOI) and control (AOC) plants (Fig. 8a). POX and CAT enzymatic assays were carried out to investigate the nature of the defence response in the two Asparagus species. No significant difference in POX activity was detected between inoculated (AOI) and control (AOC) A. officinalis plants. However, in A. kiusianus, POX activity was significantly higher (1.32-fold increase; P < 0.002) in inoculated (AKI) plants than in control (AKC) plants (Fig. 8b). Inoculation stimulated significant increases in CAT  Values represent the mean on three independent biological replicates (n = 3) ± standard deviation (SD). Significance levels are given as: *(P < 0.05) and **(P < 0.01) and according to analysis of variance (ANOVA). activity in both A. officinalis and A. kiusianus. CAT activity was 1.64-fold higher (P < 0.001) in AOI plants than in AOC plants and 1.55-fold higher (P < 0.001) in AKI plants than in AKO plants (Fig. 8c).

Discussion
Asparagus stem blight, which is caused by P. asparagi, is a serious disease that affects asparagus production worldwide, and there is an urgent need to produce asparagus cultivars with strong resistance to this disease. After artificial inoculation with P. asparagi, A. kiusianus, a wild relative of cultivated A. officinalis, exhibited significantly (P < 0.001) reduced disease severity compared with A. officinalis (Fig. 1), consistent with a previous study examining disease in these Asparagus species 4 . These results suggest that wild A. kiusianus is a potential genetic reservoir for Phomopsis disease resistance improvement in cultivated A. officinalis. Indeed, interspecific hybridization between A. kiusianus and A. officinalis resulted in F 1 hybrids with strong Phomopsis disease resistance 4 . However, the mechanisms underlying Phomopsis disease resistance in A. kiusianus remain unknown. Candidate A. kiusianus genes involved in Phomopsis disease resistance were identified by comparing transcription in susceptible A. officinalis and resistant A. kiusianus 24 h after infection with P. asparagi. Samples from inoculated and control A. officinalis and A. kiusianus plants were sequenced using the Illumina HiSeq 2500 platform, and de novo transcriptome assemblies were generated and compared. In total, 1,779 differentially expressed transcripts, 1,027 of which were up-regulated and 752 of which were down-regulated, were detected in the two Asparagus species in response to P. asparagi infection (Fig. 4). GO term enrichment and KEGG metabolic pathway analyses revealed that transcripts that were highly expressed in A. kiusianus and exhibited low expression in A. officinalis were enriched for metabolic pathways, biosynthesis of secondary metabolites, plant-pathogen interaction, protein kinase signaling and plant hormone biosynthesis and signaling transduction (Fig. 5 and Supplementary Table S3).
Expression of stress and defence response-related genes was strongly up-regulated in resistant wild A. kiusianus in comparison with susceptible A. officinalis (Figs 5 and 6) in response to P. aspargi infection. Up-regulation of PR genes, including PR-1, is a prerequisite for activation of systemic acquired resistance (SAR), Values represent the maximum, third quartile, median, first quartile and minimum of three independent replicates (n = 3). Different letters indicate statistically significant difference at P < 0.05 according to Tukey's honest significant difference (HSD) post hoc test.
which plays a central role in plant basal defence against pathogen infection 15,16 . Plant chitinases and POXs, which are also classified as PR proteins, play direct roles in plant resistance by hydrolyzing fungal cell walls and oxidating phenolic residues in the infected tissues, respectively 17-19 . A similar pattern of up-regulation of stress and defence response-related genes was reported in leaves of Withania somnifera after stimulation of salicylic acid-induced defence mechanisms 19 . Up-regulation of POX, chitinase, and PR transcripts was also reported in Fusarium-resistant wheat (Triticum aestivum) compared with a susceptible genotype 20 . Similarly, rice (Oryza sativa) variety 'Selenio' , which was resistant to Fusarium fujikuroi, exhibited up-regulation of several PR genes compared with the susceptible variety 'Dorella' 21 . Our data and results from previous studies clearly suggest that POX, chitinase, and PR-1 play important roles in plant basal defenses and SAR induction during early plant responses to infection in resistant genotypes 22 . SAR is responsible for triggering host defence mechanisms to increase de novo expression of defense-related genes, leading to enhanced expression and de novo synthesis and accumulation of chitinases, POXs, and other pathogenesis-related proteins in uninfected tissues, thereby defending them against any further pathogen attack 23 . It is therefore likely that similar defence pathways are involved in responses to P. asparagi infection in wild A. kiusianus.
Our results also indicate a possible role for genes encoding ribosomal proteins (RPLs) in Asparagus defence mechanisms. Several genes encoding RPLs were strongly up-regulated in inoculated A. kiusianus (AKI) plants and suppressed in inoculated A. officinalis (AOI) plants (Fig. 5d). RPLs have established roles in facilitating protein synthesis and preserving the stability of the ribosomal complex 24 . The involvement of RPLs in plant defence was recently reported in non-host resistance against bacterial pathogenic attack in Nicotiana benthamiana 25 . The study suggested that RPL12 and RPL19 played substantial roles in the development of non-host disease resistance through the induction of the hypersensitive response (HR) 25 . Similarly, strong up-regulation of 80% of rice RPL-related genes was observed in response to exposure to Xanthomonas oryzae 24 . The putative promoter regions of these RPL-related genes carry cis-elements that respond specifically to stress, suggesting that genes in the RPL family might be useful targets in strategies to develop stress tolerance in rice and other crops 24 . Large-scale Figure 8. Accumulation of (a) total protein content, (b) peroxidase, and (c) catalase activities in Asparagus officinalis inoculated with Phomopsis asparagi (AOI) and A. kiusianus inoculated with P. asparagi (AKI). A. officinalis and A. kiusianus treated with sterile distilled water (AOC and AKC, respectively) served as controls. Values represent the mean ± standard error (SE) of three independent replicates (n = 3). Different letters indicate statistically significant difference at P < 0.05 according to Tukey's honest significant difference (HSD) post hoc test.
analysis of differential gene expression in coffee (Coffea arabica) genotypes that were resistant or susceptible to leaf miner revealed up-regulation of RPL-encoding genes in the resistant genotype compared with the susceptible genotype 26 . Our results indicated that genes encoding RPLs may be involved in disease resistance in A. kiusianus in response to P. asparagi exposure, and that this response may be species-specific. Further molecular and biochemical analyses of RPLs identified in this study will be needed to understand the functional role of RPLs in the Asparagus defence mechanism.
In the present study the expression level of several JA biosynthesis-related genes including PLDα1, LOX1, OPR and JIP-23 KD was elevated in the resistant wild A. kiusianus in comparison with susceptible A. officinalis in response to P. aspargi infection (Figs 5 and 6). PLD and LOX genes are necessary for the initial steps of JA biosynthesis in Arabidopsis 27 . PLD release linolenic and α-linolenic acids from chloroplast membranes and these substrates are subsequently catalyzed by LOX, leading to the formation of hydroperoxy octadecadienoic acids 27 . During the early stage of plant-pathogen interaction the activation of PLD and LOX is essential for the production of important defense signaling molecules, such as oxylipins and JA which has been shown to modulate the activity of variety of proteins involved in defense signaling 28 . RNA-Seq analysis of F. fujikuroi-resistant rice variety 'Selenio' showed up-regulation in PLD and LOX genes in comparison with susceptible rice 'Dorella' , indicating a crucial role of JA pathway in the resistance of rice to F. fujikuroi infection 21 . JA signaling has broadly been associated to the defense against necrotrophic pathogens, inducing the accumulation of secondary metabolites and PR proteins 29 . However, recent studies have revealed that certain biotrophic fungal species can also trigger JA-mediated responses 29 . In this context, Phomopsis pathogens are nectrotrophic at least for the latent phase of infection and are therefore named hemi-biotrophs 30 . In A. kisuainus a substantial increase in OPR expression was observed (Figs 5 and 6). OPR catalyzes the final step in JA biosynthesis by reducing 12-oxophytodienoic acid to 3-oxo-2(2′-pentenyl)-cyclopentane-1-octanoic acid 28 . OPR was previously proposed as a disease resistance marker in tomato (Solanum lycopersicum) 31 . A cleaved amplified polymorphic sequences (CAPS) marker derived from tomato OPR showed putative co-segregation with a potato (S. tuberosum) quantitative trait locus (QTL) for late blight disease resistance. Likewise, gene expression profile of JA biosynthesis-related genes (LOX, AOC and OPR) in Plasmopara viticola disease resistance and susceptible grapevine (Vitis vinifera) cultivars, revealed a strong up-regulation of these genes in the resistant cultivar in comparison with susceptible cultivar 32 . Additionally, expression of OPR and LOX was substantially higher in a Fusarium head blight-resistant wheat variety 'Wangshuibai' than in a susceptible variety 'NAUH117' 24-48 h after infection 33 . These results suggest that OPR and LOX have potentially important roles in the early defence response in wheat against Fusarium head blight 33 . Another JA biosynthesis-related gene, JIP-23 KD, was expressed at substantially higher levels in inoculated A. kiusianus (AKI) plants than in AKC plants. In A. officinalis, JIP-23 KD exhibited reduced expression in inoculated plants (AOI) compared with control plants (AOC) (Fig. 5d). An early study examining the effect of JA on the interactions between barley (Hordeum vulgare) and powdery mildew suggested a possible role for JIP-23 KD in cell wall modification or in pathogen defence 34 . However, the specific roles of JIP-23 KD in pathogen defence remain to be determined 34 . Recently, JIP-23 KD was shown to be up-regulated in cadmium-tolerant barley genotype compared with a cadmium-susceptible genotype 35 . Collectively, these results suggest that JIP-23 KD may have a role in plant stress responses. The up-regulation of PLDα1, LOX1, OPR and JIP-23 KD genes in A. kiusianus inoculated with P. asparagi suggests that JA signal transduction may play a crucial role in Phomopsis disease resistance in A. kiusianus. Phytohormone analysis revealed a significant increase in JA and MeJA contents in wild-resistant A. kiusianus relative to cultivated-susceptible A. officinalis in response to P. asparagi infection (Fig. 7), providing an additional support for our hypothesis regarding the important role of JA-dependent signaling pathway in the P. asparagus disease resistance.
The fluctuation in POX and CAT activities between resistance and susceptible genotypes has been previously reported 36 . Therefore, stress-response enzyme activities were assessed in the present study. A significant (P < 0.002) increase in POX activity was seen in A. kiusianus plants after P. asparagi infection (AKI) compared with control plants (AKC). No significant difference was seen in POX activity between inoculated (AOI) and control (AOC) A. officinalis plants (Fig. 7). This suggested that POX played a substantial role in A. kiusianus defence via suppression of P. asparagi spread. The accumulation pattern of POX in inoculated A. kiusianus plants correlated with the up-regulation of several POX-related genes in the RNA-Seq data (Fig. 5c). Similar findings were reported in cabbage (Brassica oleracea var. capitata), where strong POX activities and up-regulation of POX-related genes was seen in cabbage resistant to black rot disease 37 . CAT activity was significantly (P < 0.001) induced by P. asparagi inoculation in both A. kiusianus (AKI) and A. officinalis (AOI) (Fig. 7c). Similar findings were previously reported in A. officinalis inoculated with P. asparagi 2 : CAT activities increased at all time points after infection whereas POX activities increased only initially, and declined thereafter 2 . These data indicate that CAT is involved in basal defence mechanisms in both resistant and susceptible Asparagus species. This is the first study to examine the molecular responses of Asparagus species to P. asparagi infection. RNA sequencing was used to identify DEGs in resistant wild A. kiusianus and susceptible cultivated A. officinalis 24 h after infection with P. asparagi. Functional annotation and KEGG pathway analysis showed that the group of up-regulated genes in A. kiusianus was enriched for metabolic pathways, biosynthesis of secondary metabolites, plant-pathogen interaction, transcriptional regulation, protein kinase signaling, phenylpropanoid and hormone biosynthesis and signaling. Results from RNA-Seq data and qRT-PCR were correlated, confirming the reliability of the transcriptome data. Activity of the stress-related enzyme POX was elevated in A. kiusianus compared with A. officinalis. Overall, comparative transcription profiling provided valuable insights into the mechanisms underlying Phomopsis disease resistance in Asparagus. These findings will be valuable in the future development of disease-resistant asparagus varieties. The RNA-Seq datasets generated in this study will be mined for sequence variations associated with gene structure and function, which will facilitate genetic trait mapping and marker-assisted selection in asparagus breeding programs.

Methods
Plant materials. Cultivated A. officinalis 'Mary Washington 500W' (AO0060 strain) and wild A. kiusianus (AK0501 strain) for RNA-Seq analysis were grown under standard greenhouse conditions (25 ± 2 °C and 14 h light/10 h dark) at Kagawa Prefectural Agricultural Experiment Station (Kagawa, Japan). Male A. kiusianus and female A. officinalis were kept to grow for 4-year-old in a commercial soil consisted of loam, well-drained soil, leaf mold and garden soil with fertilizers (2:1:1:2 and 2:1:1:4 for A. kiusianus and A. officinalis, respectively v:v) containing 50 mg N kg −1 , 500 mg P kg −1 , and 100 mg K kg −1 [pH 6-6.5, water-holding capacity ~70%, Nippi Engei Baido No.1 (Nihon Hiryo Co. Ltd)] as a root-supporting medium ( Supplementary Fig. S2). A. officinalis and A. kiusianus with P. asparagi for RNA-Seq analysis. Plants were artificially inoculated with spores of P. asparagi at a final concentration of 10 7 CFU ml −1 according to the vinyl cotton (VC) method 38 (Supplementary Fig. 2). Minimal methodological modifications were made: silicone tubes were used instead of vinyl tubes and water-retentive polyester fiber sheeting was used instead of cotton. Stems from three independent biological replicates (n = 3) were harvested 24 h after infection, immediately frozen in liquid nitrogen, and stored at −80 °C until used. Plants treated with SDW under the same conditions served as a control.

Inoculation of
Determination of stem blight disease severity in cultivated and wild Asparagus species. In a separate experiment, A. officinalis and A. kiusianus plants cultivated and inoculated as described above were kept for 3 weeks after inoculation and disease severity on the infected stems was recorded. Disease severity was scored on a 1-5 scale, where 1 = healthy and 5 = heavily infected. Disease severity percentage was calculated based on the following formula: RNA-Seq data processing and de novo transcriptome assembly. ]ABA (OlchemIm) were added to the extracts to serve as internal standards. After overnight extraction at 4 °C, solids were separated by centrifugation and re-extracted for 30 min in 5 ml of the same extraction solution. The extracts were combined, and passed through BOND ELUT C18 column (500 mg, Agilent Technologies), equilibrated with 1% acetic acid. The eluate was evaporated and dissolved in water:methanol:acetic acid (89.9:10:0.1, v-v:v) and analyzed by LC-MS/MS system. The LC-MS/MS system consisted of a Prominence 20A Series HPLC (Shimadzu) equipped with a 3200 QTrap LC/MS/MS System (AB Sciex), using an electrospray interface. For quantification of JA, SA and ABA, the purified samples were injected onto a Cadenza CD-C18 column (3 µm, 150 × 3.0 mm; Imtakt) at 45 °C and eluted at a flow rate of 0.2 ml min −1 . Hormones were separated with a gradient of mobile phase A (water:methanol:acetic acid, 89.9:10:0.1) and B (methanol). The initial conditions were 100% A, maintaining for 2 min, changing linearly to 40% A and 60% B in 5 min, changing to 100% B in 10 min, and finally maintained at 100% B for 6 min. The column was equilibrated with the starting composition of the mobile phase for 12 min. For quantification of MeJA, the purified samples were injected onto a Shim-pack XR-ODS column (2.2 μm, 75 × 2.0 mm; Shimadzu). The mobile phases and gradient conditions were described above. Quantification was obtained by multiple reaction monitoring (MRM) of the selected precursor ions and a specific product ions as described in Supplemental Table S4.
Preparation of crude enzyme extract. Stem tissue (200 mg) was ground to a fine powder in liquid nitrogen using a pre-cooled mortar and pestle, and 2 ml of extraction buffer [0.2 M phosphate buffer (pH 7.2), 0.1 mM EDTA, 1 mM DTT, and 2 U protease inhibitor cocktail] was added. The macerated suspension was centrifuged at 10,000 rpm for 5 min at 4 °C. The supernatant was collected and used as the source of enzyme. Total protein content of the crude enzyme extract was determined using a Bradford assay 48 .

Determination of stress-related enzyme activities. POX activity was estimated using a UV/Vis
Beckman DU 730 spectrophotometer (Beckman Coulter Inc. California, USA) as described previously 49 . The reaction mixture, which consisted of 0.8 ml of 0.2 M phosphate buffer (pH 7.2), 1 ml of 15 mM guaiacol, 1 ml of 3 mM hydrogen peroxide, and 0.2 ml of crude enzyme extract, was incubated at room temperature for 3 min. Absorbance of the colored product was monitored at 470 nm. POX activity expressed as ∆470 g −1 fresh weight (FW) min −1 was calculated using the following formula: (2) 1 CAT activity was determined spectrophotometrically at 240 nm as described previously 50 . CAT activity expressed as ∆240 g −1 FW min −1 was calculated using the following formula, modified with hydrogen peroxide coefficient ε240: 1