Systemic Suppression of the Shoot Metabolism upon Rice Root Nematode Infection

Hirschmanniella oryzae is the most common plant-parasitic nematode in flooded rice cultivation systems. These migratory animals penetrate the plant roots and feed on the root cells, creating large cavities, extensive root necrosis and rotting. The objective of this study was to investigate the systemic response of the rice plant upon root infection by this nematode. RNA sequencing was applied on the above-ground parts of the rice plants at 3 and 7 days post inoculation. The data revealed significant modifications in the primary metabolism of the plant shoot, with a general suppression of for instance chlorophyll biosynthesis, the brassinosteroid pathway, and amino acid production. In the secondary metabolism, we detected a repression of the isoprenoid and shikimate pathways. These molecular changes can have dramatic consequences for the growth and yield of the rice plants, and could potentially change their susceptibility to above-ground pathogens and pests.


Introduction
Rice (Oryza sativa L.) is one of the most important staple foods in the world, feeding more than 50% of the human population. One of the most damaging pathogens, with major impact on rice yield, is the migratory endoparasitic nematode Hirschmanniella oryzae. Hirschmanniella oryzae is present in the soil and irrigation water of the majority of rice-growing areas, mainly where rice is grown under submerged conditions [1], [2]. These nematodes penetrate plant roots and move through the cortex of the root, producing cavities and channels and eventually necrosis. Roots invaded by H. oryzae exhibit discoloration, deterioration and rotting [1], while heavily infested plants also show above-ground symptoms like stunting and up to 60% reduced tillering [1], [2].
In comparison with the existing knowledge on the infection process of dicots by sedentary endoparasitic nematodes [3], far less is known about the interaction between monocot plants and nematodes or plant interactions with migratory endoparasitic nematode species. Plants generally respond to nematode infection by differential expression of genes involved in stress and defense responses, cell wall alteration, metabolism and nutrient allocation, and signal transduction [3]. Our previous transcriptome analysis on the local response of rice roots upon infection with the H. oryzae revealed that at 3 days after inoculation the infected roots accumulated mRNA of many biotic stress-related genes, oxidative stress and cell death pathways [4]. Concurrently the normal primary metabolism of the rice roots appeared impaired. Particularly the jasmonate and ethylene pathways, which are known to be involved in wound responses [5], were strongly induced upon migratory endoparasitic nematode infection in rice roots, although some suppression of specific defense responses was observed at later time points [4]. Next to the accumulating evidence on gene expression and metabolite changes in locally infected tissues, research has also provided evidence on the nematode-induced systemic changes in pathogenesis-related protein levels [6], [7], [8], primary metabolites [9], and hormonerelated defense pathways [10]. Below-ground feeding organisms such as insects, nematodes, root pathogens, and ectomycorrhizal fungi are known to influence the concentration of above-ground plant defense compounds including terpenoids, glucosinolates, and phenolics [11], [12], [13], [14], leading to an impact on the susceptibility of the above-ground tissues to subsequent pathogen attack (reviewed in [15]). However, the molecular mechanisms and associated genes driving such a systemic response remain barely investigated.
Therefore, to gain deeper insight into the systemic transcriptional changes in rice after H. oryzae infection we have performed mRNA-Seq on the shoots of rice root nematode infected plants.
The observations were independently validated using qRT-PCR and biochemical analyses. This research reveals significant modifications in the metabolism of the plant, with a general suppression of chlorophyll biosynthesis and primary metabolic processes involved in plant growth.

Infection
Oryza sativa cv. 'Nipponbare' (GSOR-100, USDA) was germinated on wet filter paper for 6 days at 30uC and then transferred to SAP-substrate (Sand and Absorbent Polymer; [16]) and further grown at 26uC under a 16 h/8 h light-regime. Hirschmanniella oryzae infected rice roots were collected in rice fields in Myanmar and sent to Ghent University by Zin Thu Zar Maung (Plant Protection Division, Yangon, Myanmar). This sampling was done at Padan Village, Hlaingtharyar Region in Yangon Division, Myanmar. As this did not involve endangered or protected species, no specific permissions were required. The nematode-infected rice roots were imported in Belgium, issued by permit BE/LOA/2012/033. All further experiments were done at Ghent University, authorized by biosafety permit T38-0428, activity number 13.
The nematodes were extracted using the Baermann funnel method [17]. Eighteen-day-old rice plants were inoculated with a mixed population of H. oryzae at the rate of 400 per plant, or mock-inoculated with water. This was done by inserting a 1 ml pipette tip just adjacent to the plant root system and releasing the nematode suspension (dissolved in water). One day after inoculation the plants were transferred to a hydroponic culturing system with Hoagland solution [16] to synchronize the infection process. Shoots of inoculated and non-inoculated plants were collected 3 and 7 days post inoculation (dpi). At each time point two independent biological replicates, each containing a pool of 6 plants, were collected for RNA-Sequencing. A third independent biological replicate of both inoculated and non-inoculated plants, again each containing a pool of 6 different plants, was used for qRT-PCR validation.

RNA extraction, library preparation and sequencing
RNA was extracted using the Qiagen RNeasy Plant Mini Kit (Qiagen), with an additional sonication step after addition of buffer RLT. RNA integrity was checked using the Agilent BioAnalyzer 2100 (Agilent). Approximately 2 mg of total RNA was used for mRNA-Seq library construction according to the manufacturer's recommendations (Illumina).
We used the multiplexing sequencing adapters provided in the Multiplexing Sample Preparation Oligo Kit (Illumina). Size selection of the library was performed on a 2% agarose gel (Low Range Ultra Agarose, Biorad 161-3107). The denatured library was diluted to a final concentration of 6 pM and loaded on a paired-end read flow cell (TruSeq v5 kit, Illumina). To minimize lane effects the samples were multiplexed. Each sample was sequenced in duplicate in 2 different lanes (4 lanes total with 8 MID tags per lane). After cluster generation, the multiplexed library was sequenced on an Illumina Genome Analyzer IIx (36 cycles, paired end).

Mapping reads to genome data and annotated transcripts
Reads were mapped to the Oryza sativa subsp. japonica reference genome (build MSU7.0) in two phases using TopHat version 1.3.1 [18] and Cufflinks, version 1.0.3 [19]. The workflow and settings used in the data analysis have previously been described [4].

Calculation, normalization and profiling of gene expression
Expression was quantified per sample and per annotated or unannotated transcript as the sum of all reads mapped to the respective gene exons with a 16 base pair tolerance on either side to compensate for potential errors in the gene annotation. Expression profiles were assessed using the R-package ''baySeq'', version 1.5.1. [20]. To compensate for artificial differences in read distributions, the original library sizes were multiplied by additional normalisation factors calculated using the Trimmed Median of M-values method described in [21] with standard settings as implemented in the edgeR package (version 2.0.3). For statistical significance an FDR cut-off of 0.05 was considered for the 'general response' analysis, where the time point after inoculation was not taken into account. For the analysis per time point (3 and 7 dpi), considering that only 2 biological replicates (although pools) per treatment were analyzed for these tissues, the statistical analysis was performed with less stringent conditions (FDR,0.1).
For all further analyses the expression level of each transcript for each condition was estimated as the fold change (FC) of mapped reads relative to the controls. The FC was calculated as follows: reads were normalized as described earlier and averaged over the biological replicates. Before calculating the base 2 log of the ratio of these averages, the number of reads was increased by 1 in each group (to avoid 0-values).

Gene ontology and enrichment analyses
Gene ontology (GO) analysis and GO enrichment were performed using agriGO [22]. Parametric Analysis of Gene Set Enrichment (PAGE) [23], based on differential gene expression levels (log 2 FC), was executed. Benjamin and Hochberg false discovery rate analysis (FDR) was performed to adjust the PAGE P-values.
In addition, MapMan [24] was used to visualize the expression of genes onto metabolic pathways and the WSR-test (with Benjamin and Hochberg correction) was used to test the statistical significance of differential expression of these pathways.

Analysis of Novel Transcriptionally Active Regions (nTARs)
The Cufflinks program generates a GTF file including all transcripts annotated in MSU7.0 and putative novel transcripts derived from the data. BLASTx searches were performed against Swiss-Prot and trEMBL and all predicted rice proteins (http:// rice.plantbiology.msu.edu/). Homologues of the nTARs in rice ESTs (downloaded from NCBI on February 1, 2014) were searched by BLASTn. For all analyses a bitscore cut-off of 50 was used.

Validation of mRNA-Seq by qRT-PCR
Based on potential functional importance, a subset of genes was selected for validation in an independent biological sample (pool of 6 plants) by qRT-PCR, at 3 dpi. Locus numbers of these transcripts and primer sequences are presented in Table 1. PCRs were performed under the following conditions: 10 min at 95uC and 45 cycles of [25 s at 95uC, 60 s at 58uC and 20 s at 72uC]. After the PCR reaction, a melting curve was generated by gradually increasing the temperature to 95uC to test the length as an indication of PCR specificity. Normalization of the expression level of the target genes was done using two reference genes (LOC_Os03g27010, [25]; LOC_Os07g02340.1, [26]). Data was analysed as previously described [4].

Chlorophyll measurement
Chlorophyll measurements were performed in an independent experiment. Nematode collection, as well as growth conditions and inoculation of the plants were as described above. Shoots of inoculated and non-inoculated plants were harvested at 3 and 7 dpi.  [27]. After checking for normality and homoscedasticity of the data, statistical significance was evaluated using student's t-test. This whole experiment was repeated in time, with similar results.

Results
More than 95 million RNA fragments were acquired and made publicly available in the GEO database (GSE57707). The short reads were aligned against the reference genome sequence of cv. Nipponbare (MSU7.0) and in total 81% of the sequenced reads could be mapped ( Table 2). The total number of sequenced bases was nearly 7 billion, representing 18-fold the rice genome size and about 50-fold coverage of the annotated transcriptome. The expression of a total of 89,246 different rice transcripts was detected in the analyzed tissues. Comparative gene expression profiling was performed by statistical analysis of differential gene expression levels, Gene Set Enrichment, and pathway mapping. In addition, a search was performed to detect novel transcriptionally active regions (nTARs) not annotated in rice genome assembly MSU7.0. In a first analysis, the expression level in all sampled systemic tissues of infected plants (at both time points) was compared with the shoot tissues of healthy plants to look for consistent trends, regardless of the time point after inoculation.
At an FDR cut-off of 0.05, 52 annotated genes were found to be significantly differentially expressed in the shoots after root infection with H. oryzae (Table 3). Twenty one of these differentially expressed genes (DEGs) were significantly downregulated, while 31 showed higher expression in infected plant shoots. Among the downregulated genes for instance, we identified one encoding a RAS-related protein (Small monomeric GTPase), and a gene encoding 3-oxo-5-alpha-steroid 4-dehydrogenase, involved in brassinosteroid biosynthesis. Among the upregulated DEGs, 2 genes encoding AP2 domain containing proteins, and a cupin domain containing protein were found. Next to the annotated transcripts, 12 nTARs were found to be differentially expressed in shoots of inoculated vs. non-inoculated plants (Table 3).
Parametric Analysis of Gene Set Enrichment of the differentially expressed genes in the infected versus uninfected plants revealed that the repressed genes were enriched for the secondary GOgroups 'developmental process', 'multicellular organismal process', 'cellular process', 'metabolic process', 'organelle', 'cell', 'cell part', 'catalytic activity' and 'binding' (Fig. 1A). None of the enriched gene sets showed induction in systemic tissue of infected plants in comparison with control shoot tissue.
When focusing on the general plant metabolism (Fig. 1B) the WSR-test of Mapman found a significant suppression of many primary metabolic processes, like 'light reaction', 'fatty acid synthesis and elongation', 'ATP synthesis', 'tetrapyrrole synthesis', 'nucleotide metabolism', 'amino acid synthesis', 'TCA', 'glycolysis', and 'cell wall precursor synthesis'. Not only the primary metabolism of the plants was impaired, we also detected a specific suppression of a few plant defense related pathways. For instance, the biosynthesis of chorismate through the shikimate pathway, was strongly suppressed in the systemic tissue. Additionally, a significant suppression of the 'isoprenoids' pathway was observed.
The impact of the observed transcriptome changes in the primary metabolic pathways was analyzed by measuring the length of roots and shoots of infected versus uninfected plants. At 3 and 7 dpi, plant root and shoot length was not different from the control plants (data not shown). Because 'light reaction' and 'tetrapyrrole synthesis' were among the gene sets with the strongest levels of suppression upon infection, the chlorophyll content was measured in the leaves of H. oryzae-infected and control plants at 3 and 7 dpi (Table 4). At 3 dpi, a decreasing trend in levels of both

Systemic transcriptome changes at 3 days after rice root nematode infection
When considering the 3 dpi samples alone, only 2 biological replicates (although pools) per treatment were analyzed, and therefore the statistical analysis was performed with less stringent conditions (FDR,0.1). This resulted in 195 DEGs, with 169 repressed and 26 induced (Table S1). Next to the annotated transcripts, 25 nTARs were found to be differentially expressed in this tissue versus uninfected plant shoots (Table S1). Pathway mapping with MapMan showed a general suppression of the plant primary metabolism in the systemic tissues at 3 dpi ( Figure S1). Among the most intriguing and strongest downregulated DEGs were genes encoding (1)   as non-host-specific defense response, wounding response, and symbiosis [28]. Two other suppressed DEGs were encoding enzymes belonging to the shikimate pathway, chorismate mutase (Log 2 FC = -1.51) and isochorismatase (Log 2 FC = -1.78).
To validate the gene expression levels as detected by RNA-seq, an independent validation was performed using qRT-PCR, comparing infected with uninfected plants at 3 dpi. In total, 14 out of the 16 tested gene expression profiles (87.5%) were confirmed (Table 1), in line with an expected FDR of 10%.

Systemic transcriptome changes at 7 days after rice root nematode infection
For the 7 dpi data, the statistical analysis was again performed with less stringent conditions (FDR,0.1), resulting in 17 DEGs, with 11 repressed and 6 induced. Next to the annotated transcripts, 10 nTARs were found to be differentially expressed in this tissue versus uninfected plant shoots (table S2). Most induced DEGs are annotated as 'expressed protein' or 'hypothetical protein', except for a gene annotated as heat shock protein DnaJ (LOC_Os01g42190; Log 2 FC = 0.528), and LO-C_Os04g18650 (Log 2 FC = 1.652), coding for a pathogenesisrelated transcriptional factor, with an ethylene-responsive AP2 domain (EREBP34, TSRF1, belonging to class IIIb; [29]). The tomato homolog of this AP2 gene, TSRF1 has previously been shown to enhance the osmotic and drought tolerance of rice [30]. This locus was also slightly (but not significantly) induced at 3 dpi in the systemic leaf tissue. Among the most strongly repressed

Analysis of novel transcripts
A total number of 33,638 nTARs were detected in the analysed tissues (Table S3). A BLAST search against all ESTs from O. sativa resulted in 30,214 nTARs giving a significant hit (bitscore. 50), while tBLASTx against all proteins of O. sativa resulted in significant hits for 19,735 of the nTARs, indicating that these nTARs are potential paralogues of known rice transcripts (Table  S3). A Singular Enrichment analysis (FDR,0.05) of these rice transcripts revealed four enriched GO-terms of the category 'Molecular function': kinase activity, transferase activity, transferase activity -transferring phosphorus-containing groups and catalytic activity. A SwissProt search was done, and this was successful for 6,271 transcripts (Table S3), revealing for instance nTARs showing resemblance to protein domains annotated as xylanase inhibitor, MATE efflux family protein, histone-lysine Nmethyltransferase, ribonuclease H protein, sugar transport protein, and ethylene-responsive transcription factor ABI4. A BLAST search against Trembl resulted in 16,134 hits, mainly to uncharacterized proteins from O. sativa.

Discussion
Gene expression analyses investigating the plant response upon nematode attack have mainly targeted the local reaction of the plant [4], [25], [31], [32], [33]. This localized approach is most appropriate to identify responses directly regulated by nematodes, while studies using systemic tissues can provide a global view of the infection response in host plants and offer a broader perspective on plant health in the context of plant-nematode interactions.
Species within the genus Hirschmanniella infest 58% of the world's rice fields, causing 25% yield losses [1]. Hirschmanniella oryzae is one of the few plant-parasitic nematodes that can survive in anoxic environments [34] and hence this nematode is the most common plant-parasitic nematode in rice grown under constantly flooded conditions [1]. In previous studies performed by our research group we investigated (1) the systemic transcriptional responses of selected marker genes involved in hormone related defense pathways using qRT-PCR [10] and (2) the local reaction of rice roots upon H. oryzae and Meloidogyne graminicola infection using mRNA-Seq [4], [25].
In the current study, a more general overview of the systemic responses of H. oryzae-infested rice plants is provided. Using mRNA-Seq we have acquired almost 100 million reads from infected and healthy rice shoot tissues (.5 M per sample), thereby obtaining count data for more than 80,000 annotated and unannotated rice transcripts. We have used a synchronized inoculation method, in which the nematodes only had a limited period (24 hours) to infect roots. This allowed a relatively uniform infection process, so that the transcript changes detected here are not compromised by continuous entry of juveniles. Systemic effects on defense responses Locally, roots infected by H. oryzae exhibit a fast (3 dpi) and strong induction of genes involved in cell death as well as defense response pathways, like for instance the phenylpropanoid and jasmonate pathways [4]. The systemic shoot tissues also showed an induction of some defense-related genes at 3 dpi, including genes involved in jasmonate and ethylene production [10] and genes involved in flavonoid and phytoalexin biosynthesis (this study). However, at 7 dpi the expression of many of these genes were significantly downregulated or had returned to non-infected tissue level [10]. In contrast, shoot tissue of Arabidopsis plants infected with the cyst nematode Heterodera schachtii contained continuously higher levels of SA and JA-marker genes than corresponding control plants, from 5 till 14 dpi [8]. However, in the case of root knot nematode infection, a strong and early suppression of SA, JA, ET-biosynthesis genes, as well as PR genes has been observed already at 3 dpi [4] and 5 dpi [8], in systemic tissue of the Meloidogyne graminicola-rice [4] and Meloidogyne incognita-Arabidopsis [8] interaction, respectively. Hence, the here-observed systemic suppression of the plant defense system upon H. oryzae infection seems to be slower and less consistent than what was seen in systemic shoot tissues upon root knot nematode infection in rice and Arabidopsis [8].
In the current mRNA-Seq data we found two notable exceptions to this rule, since we detected a repression of the isoprenoid and shikimate pathway already at 3 dpi in the systemic tissue. The suppression of genes from the shikimate pathway could potentially deplete the pool of chorismate and prephenate available for the production of aromatic amino acids and salicylate. Aromatic amino acids are precursors for several secondary metabolites with proven roles in plant defense, like lignin, phenolics, flavonoids and salicylic acid, produced through the phenylpropanoid pathway [35]. Interestingly, we have recently discovered two candidate effectors of H. oryzae with homology to chorismate mutase and isochorismatase [36]. These effectors are most likely involved in the deregulation of the shikimate and subsequent phenylpropanoid pathway in the host plant. A secreted chorismate mutase from Ustilago maydis has for instance been demonstrated to lower the total salicylic acid content in infected leaves [37]. If and how these effectors can also affect the plant defence system in the above-ground tissue is an intriguing route for further research.
The isoprenoid pathway, which is also suppressed at 3 dpi in systemic shoots, is responsible for the production of a variety of both primary and secondary metabolites with very diverse functions, like tocopherol, carotenoids, mevalonate and provides precursors for brassinosteroid, gibberellin and chlorophyll biosynthesis [38]. Consistent with this view, genes involved in both chlorophyll biosynthesis as well as brassinosteroid biosynthesis and signalling were impaired in the systemic plant tissue. Brassinosteroids (BRs) are polyhydroxysteroids that offer a wide range of physiological activities starting from seed germination till flowering and senescence [39]. The strong suppression of this pathway as seen in the infected plants can potentially have dramatic effects on plant development. Although we did not observe significant effects on plant growth at 3 and 7 dpi, effects at later time points can certainly be expected.
In addition to their well-known role in plant development, BRs have recently been demonstrated to also play a role in plant responses to a broad spectrum of environmental stresses [40], [41]. For instance, BRI1, which was here found to be strongly repressed in the shoots of H. oryzae infected plants, is an important player in PAMP triggered Immunity [42].

Systemic effects on primary metabolism
The data provided in the current study reveals a remarkable suppression of the whole basic metabolic machinery in the rice shoot. The transcriptome revealed changes in the pathways needed for production of different metabolites for plant growth and development, such as carbohydrates, lipids, proteins, and nucleic acids. Similarly, Hofmann et al. [9] showed that cyst nematode (Heterodera schachtii) infection was causing dramatic changes in metabolite profiles in systemic root and shoot tissues. For instance, amino acid levels decreased in shoots of cystnematode infected plants in comparison with control shoots, but levels of some organic acids were increased [9]. It is known that defense metabolites such as phenylpropanoid pathway products, PR-proteins, ROS and callose represent a major flow of carbon from primary metabolism into secondary metabolism [43]. Large carbon fluxes into secondary metabolism during the defense response cannot occur without influencing reactions in primary metabolism [43]. The local strong defense induction inflicted in the H. oryzae infected root system [10], might be responsible for a depletion of carbon and energy available for the primary metabolism of the plant shoot.

Potential consequences for interactions with other pathogens and herbivores
In the case of root feeding by insects, a general induction of defense responses in the systemic shoot tissue has been observed [12], [44], [45]. For instance in cotton (Gossypium herbaceum), root-chewing by insect larvae led to increased concentrations of non-volatile terpenoid compounds in the leaves, known to deter insect feeding [12]. Recently Wondafrash et al. [15] provided a review of the systemic effects that root feeding by nematodes can have on subsequent attack by above-ground herbivores. They concluded that the effects varied depending on the feeding style of the nematode, nematode population densities and the susceptibility of the host plant, and different outcomes have been reported. For instance, soybeans with moderate and high Heterodera glycines infection levels supported greater Helicoverpa zea larval populations than observed on control or low density H. glycines infested plants. Nevertheless, very high H. glycines population densities caused severe stunting and chlorosis, reducing the suitability of these soybeans as a food source for H. zea larvae [46]. However, in Arabidopsis, a simultaneous inoculation with the sedentary endoparasitic nematode H. schachtii and the aphid Brevicoryne brassicae was performed [47] and no effects of the nematode infection were observed on the aphid performance at day 3. Again different, Agrostis capillaris and Anthoxanthum odoratum plants attacked by ectoparasitic and migratory endoparasitic nematodes, caused a reduced fecundity of the aphids Rhopalosiphum padi feeding on the shoots [48]. A lower amino acid content was observed in the phloem sap of the nematodeinfected plants, and was probably one of the causes of this decreased aphid performance. Since our transcriptome data provides evidence for a reduced amino acid production in the above-ground tissues of H. oryzae infected rice plants, one could deduce that this tissue might be less attractive to phloem-feeding insects. Nevertheless, experimental evidence for this is currently lacking.
In conclusion, this transcriptome analysis of above-ground tissues upon rice root nematode infection unveils a major impairment of primary metabolite biosynthetic pathways in the shoots of infested plants. Biochemical investigation of chlorophyll levels in these tissues, corroborated the significant reduction of the plant's photosynthetic capacity upon H. oryzae infection. Although less prominent, some subtle changes in defence-related pathways were also detected, among which mainly the isoprenoid and shikimate pathway seem to be negatively affected. If and how these changes affect the abiotic and biotic stress tolerance of the surviving rice plants remains to be further investigated. Figure S1 Mapman visualization of the expression profiles of genes involved in the general metabolism of the rice plant. The visualization shows the observed differential expression patterns, based on the Log 2 fold changes of mRNA levels, in shoots of infected versus uninfected control plants at 3 dpi (A) and 7 dpi (B). Red dots indicate that the gene is upregulated in infected plants versus the corresponding healthy control plants, while blue indicates downregulation. (ZIP)   Table S3 List of putative novel transcriptionally active regions identified in the investigated rice shoot material, with their genomic location, BLASTn results against cDNA and genomic DNA of Oryza sativa cv. 'Japonica', BLASTx results against rice proteins, SwissProt/tremble. For all analyses a bitscore cut-off of 50 was used. (XLSX)