Chorismate mutase and isochorismatase, two potential effectors of the migratory nematode Hirschmanniella oryzae, increase host susceptibility by manipulating secondary metabolite content of rice.

Hirschmanniella oryzae is one of the most devastating nematodes on rice, leading to substantial yield losses. Effector proteins aid the nematode during the infection process by subduing plant defence responses. In this research we characterized two potential H. oryzae effector proteins, chorismate mutase (HoCM) and isochorismatase (HoICM), and investigated their enzymatic activity and their role in plant immunity. Both HoCM and HoICM proved to be enzymatically active in complementation tests in mutant Escherichia coli strains. Infection success by the migratory nematode H. oryzae was significantly higher in transgenic rice lines constitutively expressing HoCM or HoICM. Expression of HoCM, but not HoICM, increased rice susceptibility against the sedentary nematode Meloidogyne graminicola also. Transcriptome and metabolome analyses indicated reductions in secondary metabolites in the transgenic rice plants expressing the potential nematode effectors. The results presented here demonstrate that both HoCM and HoICM suppress the host immune system and that this may be accomplished by lowering secondary metabolite levels in the plant.

INTRODUCTION

Hirschmanniella oryzae is one of the most devastating pests in rice fields worldwide. Because H. oryzae is well adapted to flooded conditions, it is a major problem in flooded rice ecosystems (Babatola, 1981). This migratory nematode penetrates the root and migrates through the aerenchyma, completing its life cycle in about 33 days (Karakas, 2004). Plant‐parasitic nematodes are equipped with an arsenal of cell wall‐modifying proteins that are secreted through the stylet into the host tissue to hydrolyse cell walls (Davis et al., 2011). Sedentary nematodes establish a feeding site, whereby it is important to remain undetected by the plant to prevent an adequate defence response. Effectors are secreted by this type of nematode to alter plant cell metabolism to their advantage and to attenuate the defence mechanisms of the host (Postma et al., 2012; Hewezi and Baum, 2013; Jaouannet et al., 2013). Migratory nematodes, like H. oryzae, stay mobile during their entire life cycle. They keep on migrating through the root whilst feeding. Although this lifestyle might enable them to outrun local plant defences and effector proteins inhibiting plant defence might seem unnecessary, several potential effector proteins have been predicted from the transcriptome of plant‐parasitic migratory nematodes (Haegeman et al., 2011; Nicol et al., 2012; Bauters et al., 2014). This implies that migratory nematode species also invest energy in attenuating the defence responses of their host.

The plant hormone salicylic acid (SA) plays an important role as a signalling molecule during the defence reaction on pathogen infection. SA can be produced through two distinct pathways starting from chorismate, the final product of the shikimate pathway: the isochorismate synthase (ICS) pathway and the phenylalanine ammonia lyase (PAL) pathway (leading to SA and phenylpropanoids). Several chorismate‐derived compounds play a role in plant physiology and the defence system, for example, SA, auxin, lignin, and flavonoids. Changes in concentrations of these compounds can affect the plant defence against pathogen infection (Jones et al., 2007; Nahar et al., 2012). Interestingly, transcripts encoding chorismate mutase and isochorismatase were detected in the transcriptome of H. oryzae (Bauters et al., 2014). Both enzymes can interfere in the host production of salicylic acid and/or phenylpropanoids. Chorismate mutase catalyses the conversion of chorismate to prephenate. This enzyme is present in plants and bacteria but had not been reported in animals until its discovery in Meloidogyne javanica (Lambert et al., 1999). Since then it has been reported in several other plant‐parasitic nematodes (Jones et al., 2003; Huang et al., 2005; Vanholme et al., 2009). It has been shown that chorismate mutase, secreted by plant‐pathogenic organisms, is able to alter plant cell development by impairing the development of lateral roots and vascular tissue (Doyle and Lambert, 2003). More recent research provided experimental evidence that a secreted chorismate mutase originating from Ustilago maydis reduced SA levels in host plants (Djamei et al., 2011). A decrease in SA content was also observed in tobacco leaves transiently expressing a chorismate mutase from Meloidogyne incognita (Mi‐CM3) upon infection by the oomycete Phytophthora capsici. In addition, constitutive expression of Mi‐CM3 increased susceptibility of tobacco plants to M. incognita infection (Wang et al., 2018). The gene encoding isochorismatase has not been characterized in nematodes before, but it was reported to be present in the genome of Meloidogyne hapla and the transcriptome of Rotylenchulus reniformis (Opperman et al., 2008; Wubben et al., 2010). Moreover, isochorismatase was detected in the secretome of plant‐pathogenic fungi, but not in nonpathogenic species (Soanes et al., 2008) and it was reported to reduce SA content in plants on infection. Isochorismatase converts isochorismate to 2,3‐dihydroxy‐2,3‐dihydrobenzoate and pyruvate. Isochorismate is an intermediate in the biosynthesis of SA, hence isochorismatase is capable of reducing the pool of isochorismate available for SA synthesis (Liu et al., 2014).

In this research, the genetic structure and protein properties of chorismate mutase (HoCM) and isochorismatase (HoICM) isolated from H. oryzae are described. Their activity is demonstrated by complementation of chorismate mutase (CM) and isochorismatase (ICM) deficient Escherichia coli strains, and infection assays with H. oryzae on transgenic plants showed that both proteins increase susceptibility of the host plant to nematode infection. Moreover, HoCM, but not HoICM expression, had an effect on infection with the sedentary nematode Meloidogyne graminicola. Finally, RNA‐Seq performed on rice plants expressing HoICM or HoCM revealed that both proteins interfere with secondary metabolic pathways, a result supported by metabolome analysis.

RESULTS

Genetic structure of HoCM and HoICM

The genomic sequence of chorismate mutase (HoCM, KP297892) consists of 942 nucleotides from start to stop codon, and contains two introns of 72 and 123 nucleotides. The isochorismatase gene (HoICM, KP297893) consists of 822 nucleotides from start to stop codon and includes two introns of 69 and 81 base pairs (Figure 1a,b). All exon/intron transitions contain the consensus splice site GT/AG, except for one HoICM noncanonical splice site GT/CG. The first 60 nucleotides of HoCM code for a secretion signal, targeting the encoded protein to the classical secretory pathway. Expression in the secretory glands was shown by in situ hybridization in previous research (Bauters et al., 2014). No secretion signal was detected in HoICM. Several attempts were made to perform an in situ hybridization for HoICM transcripts, but all of them were negative or showed ambiguous results. HoICM was PCR‐amplified from genomic DNA extracted from a single nematode, three times, from an independent single hand‐picked nematode. HoICM could be amplified in all three experiments, supporting that this is an endogenous H. oryzae gene rather than being derived from any contaminating material. The fact that this type of isochorismatase (ICM) is specific to plant‐parasitic species (Bauters et al., 2014) and that it is secreted by filamentous pathogens through a nonclassical secreted pathway (Liu et al., 2014) points towards an effector function in nematodes as well.

FIGURE 1

Gene structure and partial protein sequence of HoCM and HoICM. (a) HoCM has a total length of 942 base pairs at the genomic level, containing two small introns. The protein consists of three different domains: an N‐terminal signal peptide (SP), a serine/histidine‐rich domain (S/HRD), and a chorismate mutase domain (CMD) (PF01817). The alignment shows the conserved catalytic residues in several plant‐parasitic nematode species (Mjav, Meloidogyne javanica; Mare, Meloidogyne arenaria; Minc, Meloidogyne incognita; GPLIN, Globodera pallida). Eight conserved putative catalytic residues characterized in chorismate mutase of Mycobacterium tuberculosis and/or Escherichia coli are marked with filled dots. (b) HoICM is 822 bp long and contains two small introns. It contains an isochorismatase domain (PF00857), but no predicted N‐terminal signal peptide. The three putative catalytic residues are marked with filled dots. Sequence names correspond to the protein sequences encoded by different nematodes, downloaded from the Wormbase Parasite database (Ho, Hirschmanniella oryzae, ALI53582; Dd, Ditylenchus destructor, Dd_10059; Hg, Heterodera glycines, Hetgly.G000001356; Gr, Globodera rostochiensis, GROS_g01640.t1; Mf, Meloidogyne floridiensis, scf7180000418840.g2997; Mi, M. incognita, Minc3s04290g35947; Mj, M. javanica, M.javanica_Scaff754g009582; Ma, M. arenaria, M.arenaria_Scaff10008g072685; Me, Meloidogyne enterolobii, scaffold11435_cov287.g14711; Mh, Meloidogyne hapla, MhA1_Contig1915.frz3.gene2). *Sequence adjusted (part of intron used as exon) due to probable mistake in gene prediction

Important catalytic residues are conserved in both HoCM and HoICM

The HoCM protein consists of three domains, including the signal peptide. The N‐terminal domain of the mature protein has no similarities with other domains in public databases, but it has two motifs, rich in serine and histidine, respectively, hence the name serine/histidine‐rich domain (S/HRD). The C‐terminal part is similar to the Pfam chorismate mutase (CM) type II domain (PF01817) (CMD). Alignment of HoCM and other CM sequences deduced from the genomes of cyst and root‐knot nematodes showed conservation of catalytic residues as initially characterized in E. coli and Mycobacterium tuberculosis (Lee et al., 1995; Ökvist et al., 2006). The eight conserved catalytic residues are highlighted in Figure 1a. All residues are conserved in H. oryzae, but far less in cyst and root‐knot nematodes.

An ICM domain (PF00857) is predicted in the C‐terminal part of the HoICM protein. Homologous sequences are present in the expressed sequence tag (EST) database of R. reniformis, Radopholus similis, and Heterodera glycines. No homologues were found in the genome of Globodera pallida, but several homologous sequences were retrieved from the genomes of different Meloidogyne spp. and Globodera  rostochiensis. Catalytic residues in ICMs have been identified in Oleispira antarctica and Pseudomonas aeruginosa, and are conserved in H. oryzae and other plant‐parasitic nematodes (Figure 1b) (Parsons et al., 2003; Goral et al., 2012).

Predicting the protein structure of HoCM and HoICM

CMs are classified into two main groups according to their structure: the AroH and AroQ classes. Proteins of the rare AroH class, represented by the CM of Bacillus subtilis, have both α‐helices and β‐sheets in their structure (Chook et al., 1993). On the other hand, proteins of the AroQ class are more abundant and are represented by the CM structure of E. coli, containing only α‐helices (Lee et al., 1995). Secondary and tertiary structure predictions showed that HoCM is composed of α‐helices and loops, without β‐sheets, indicating that HoCM is a member of the AroQ class. A 3D model was constructed using Phyre2 with the CM structures of M. tuberculosis (PDB code 2FP1) and Yersinia pestis (2GBB) as a template. The predicted tertiary structure of HoCM is shown in Figure 2a. The CM domain of HoCM has the same predicted topology as the global structure of an AroQγ protein, which consists of six α‐helices connected by loops. The S/HRD region (N‐terminal) of HoCM lacked homology and was modelled ab initio with only 72% of the residues modelled with over 90% confidence.

FIGURE 2

Predicted tertiary structure of HoCM and HoICM. (a) HoCM. The structure of the mature protein is shown (without secretion signal). It was constructed using homology modelling with chorismate mutase models of Mycobacterium tuberculosis and Yersinia pestis as templates. The C‐terminal chorismate mutase domain (grey) consists of two large α‐helices in the middle, surrounded by four smaller α‐helices. The N‐terminal part (black) is predicted to have two α‐helices. Eight catalytic residues are shown as stick models. (b) HoICM. The six‐stranded parallel β‐sheet (black) is in the centre of the molecule, surrounded by three α‐helices at one side and two at the other side. At the bottom (three) and the top left corner (one) in this view of the molecule there are four 310‐helices. The catalytic triad (Cys164‐ Asp42‐Lys131) is labelled

ICM‐like hydrolases (ILH) are a large family of enzymes divided into several structural subgroups by the Structural Classification of Proteins (SCOP release 2.03), that is, YecD, N‐carbamoylsarcosine amidohydrolase, phenazine biosynthesis protein (PhzD), pyrazinamidase/nicotinamidase, ribonuclease MAR1, and YcaC. These proteins belong to the class of “alpha‐beta‐alpha sandwich” according to the CATH (Class Architecture Topology/fold Homologous superfamily) structural classification, adopting a Rossman fold. The 3D model of HoICM was constructed using six different ICM structures (PDB code: 1NF9, 1NBA, 2FQ1, 3TB4, 3IRV, and 3KL2), creating a model for which 95% of the residues were modelled with more than 90% confidence with Phyre2. The model shows a molecule with a Rossman fold; a six‐stranded parallel β‐sheet flanked by three big α‐helices at one side and two at the other. Next to these α‐helices, there are four putative 310‐helices present. The presumed catalytic triad is shown in Figure 2b. The catalytic centre where isochorismate can bind has been described in Oleispira antarctica and is conserved in other plant‐parasitic nematodes (PPNs) (Goral et al., 2012). PyMOL (v. 1.6) structural alignment assigns HoICM to the structural subgroup of the phenazine biosynthesis proteins (root mean square deviation of 0.437).

HoCM and HoICM complement E. coli mutants

The activity of two different HoCM constructs was assessed by a complementation assay in a CM‐deficient E. coli strain (KA12/pKIMP‐UAUC) grown on dropout medium without phenylalanine. One construct contained the full CM coding sequence without the predicted signal peptide, the second only consisted of the catalytic region. The positive control was the E. coli mutant expressing a B. subtilis CM. The results confirmed CM activity of the HoCM protein (Figure 3a). The fact that mutants with both protein constructs grew equally well reveals that the S/HRD region is not necessary for CM activity.

FIGURE 3

Complementation tests of HoCM and HoICM. (a) Complementation assay of the chorismate mutase (CM)‐deficient Escherichia coli strain KA12/pKIMP‐UAUC with different constructs of HoCM. The left Petri dish shows bacterial growth on M9cY medium supplemented with phenylalanine, the right one bacterial growth on M9cY medium without additions. +, positive control (CM of Bacillus subtilis); −, negative control (empty vector pQE30‐UA); Cat, catalytic domain of CM; Full, mature CM protein. (b) Complementation assay of entB‐mutated E. coli AN192. AN192 cells were complemented with an empty vector (pDEST17) (white) or HoICM (grey). Cells were grown in liquid medium with different concentrations of bipyridyl to create iron‐limiting conditions. Bacterial growth was observed by measuring optical density (y axis) at three different time points (an average was taken from four bacterial cultures). Treatments were compared with the control (no bipyridyl) and statistically analysed by a Mann–Whitney test. Each experiment was performed twice with similar results. Different letters above the graph indicate the significant differences between the different treatments

The activity of HoICM was tested using an ICM (entB) deficient strain (AN192) incapable of producing the siderophore enterobactin, an iron chelator. The growth of AN192 and AN192 complemented with HoICM under iron‐limiting conditions, created by adding bipyridyl, was monitored by measuring optical density. Results showed that AN192 cells harbouring HoICM grew much better under iron‐limiting conditions compared to the empty vector control (Figure 3b). Fifteen hours postinoculation, there still was a significant difference (p value < .05) between bacteria growing under normal conditions and bacteria growing in iron‐limiting conditions, but after 24 hr the cells complemented with HoICM grew as well as the cells under conditions without bipyridyl. AN192 cells that carried the empty vector could not grow well under iron‐limiting conditions, the growth rate of this control remaining far below the growth rate of HoICM‐complemented cells under the two most stringent conditions (150 and 200 µM bipyridyl). This result proves that HoICM can complement the lack of ICM activity in the AN192 strain.

Rice lines expressing HoCM or HoICM are more susceptible to nematode infection

Transgenic rice lines constitutively expressing one of two versions of HoCM (HoCM without signal peptide (HoCM_FULL) or only the catalytic domain of HoCM (HoCM_CAT)) or HoICM were generated starting from O. sativa "Nipponbare." Lines transformed with a vector without insert (empty vector) were used as control. Five and two independent lines were obtained for HoCM_CAT and HoCM_FULL, respectively. For HoICM, eight independent lines were obtained, but only three of them produced enough viable seeds to perform further experiments. No clear phenotypes were observed between independent lines during the entire life cycle of the plant. The presence of the construct was verified at the genomic level by PCR with gene‐specific primers. Expression was validated by quantitative reverse transcription PCR (RT‐qPCR; see Table S1a,b). The transgenic lines were tested for their susceptibility to nematode infection (Figure 4). Constitutive expression of the nematode HoCM or HoICM enhanced susceptibility of the rice plants against H. oryzae. All three different constructs showed similar results. The transgenic lines were used in an infection assay with the sedentary root‐knot nematode M. graminicola as well. In contrast with the results obtained for H. oryzae, two lines of HoICM‐expressing rice plants showed no difference in M. graminicola infection levels with the control. On the other hand, both HoCM constructs enhanced susceptibility against M. graminicola, similar to the results obtained for H. oryzae (Figure 4).

FIGURE 4

Effect of expression of HoICM, HoCM_FULL, and HoCM_CAT in rice on susceptibility against Hirschmanniella oryzae (a) and Meloidogyne graminicola (b). Bars represent the average number of nematodes (H. oryzae) or galls (M. graminicola) of eight infected plants per line. The number of galls was expressed relative to the control, which was set to 1. Asterisks indicate values significantly different from the control (empty vector) according to a nonparametric Mann–Whitney test (p < .05). (*) indicates p = .052. Error bars represent the standard error of the mean. Data for M. graminicola infection are expressed as relative values. The two repetitions of each experiment are shown. The infection was conducted on one line per construct

Several genes involved in plant defence are down‐regulated in rice lines constitutively expressing the potential nematode effectors

RNA‐Seq was performed on one line per construct to study transcriptional changes in rice plants expressing HoCM_CAT, HoCM_FULL, or HoICM. Almost 150 million sequence fragments were generated. Almost 95% of all reads could be mapped to the rice genome (IRGSP‐1.0, see Table S2).

Expression levels of all genes in the transgenic lines were compared to those in the line transformed with the empty vector. Differentially expressed genes (DEGs) with a false discovery rate (FDR) <0.05 were compared using blastp to protein data for Viridiplantae available in the SwissProt database for annotation purposes. Compared to the empty vector line, the transgenic line expressing HoCM_FULL had 55 DEGs (13 up‐regulated, 42 down‐regulated). When HoCM_CAT was expressed, 196 DEGs were detected (67 up, 129 down), and 105 DEGs (36 up, 69 down) were found in plants expressing HoICM. Approximately half of these genes had a significant (bit‐score > 50) hit when compared with blast against the SwissProt database. These annotated genes were classified into six groups according to their putative function (File S1). Expression of 20 genes was down‐regulated in the transgenic plants expressing any of the three constructs, among which were two glutathione S‐transferases, a mono‐oxygenase, and geranylgeranyl pyrophosphate synthase. A nonprotein‐coding transcript was up‐regulated in all three constructs. Nine genes had the same trend in transgenic plants with either CM construct, among which were a down‐regulation of genes involved in steroid biosynthesis and auxin transport. It is interesting to note that many genes involved in signalling, stress response, and secondary metabolite production were differentially expressed upon in planta expression of both nematode genes. An overview of the significantly differentially regulated genes can be found in File S1.

The reliability of the RNA‐Seq results was validated by evaluating the expression of 12 randomly chosen DEGs by RT‐qPCR. Eleven of the genes confirmed RNA‐Seq results; only one gene showed a contrasting expression pattern: the gene coding for histone H2B.3 was significantly up‐regulated according to the RNA‐Seq data (log2FC = 1.6) in plants expressing HoCM_CAT, while it was down‐regulated according to the RT‐qPCR data. The expression of two genes, coding for glutathione S‐transferase and a receptor‐like protein kinase, was so strongly down‐regulated that it was undetectable by RT‐qPCR in the transgenic lines (see Figure S1). The clear correspondence with RT‐qPCR results for these 12 genes (p = .006, binomial test) supports the reliability of the generated RNA‐Seq data.

Gene ontology enrichment analysis shows a decrease in defence compound biosynthesis terms

Parametric analysis of gene set enrichment (PAGE) was conducted using the log2 fold changes of all genes as input. All levels of gene ontology (GO) in the category “Biological Process” were considered, with a significance level of 0.05 (FDR). GO terms with a calculated absolute Z‐score higher than 4 for at least one of the three constructs are visualized in Figure 5. The expression of each different construct induces a repression of the “secondary metabolic process” GO in rice. Looking at the child terms of this GO term, the repression was due to a reduction in the GO term “Phenylpropanoid metabolic process”. The GO term “Lignin catabolic process” is reduced in HoCM_FULL‐ and HoICM‐expressing plants. “Diterpene phytoalexin metabolic process” is suppressed in plants expressing HoCM. Rice plants expressing either one of the HoCM constructs are also repressed in the general “Defense response”, while all three lines are reduced in “Oxidation reduction”.

FIGURE 5

Gene ontology (GO) term enrichment analysis, showing the significantly enriched GOterms in the category “Biological process”. GO terms for which at least one of the three constructs resulted in plant transcripts with significant difference with an absolute Z‐score (y axis) value of 4 are shown in the graph. a and aʹ indicate an ancestor term (a) with its child terms (aʹ). All GO terms represented in this graph are significantly less abundant in the transgenic lines (Z‐score < 0) compared to control plants (empty vector line). Nonsignificant data are not shown in the graph

Metabolome analysis suggests that HoCM down‐regulates phenylpropanoid biosynthesis

No differences in SA content in shoots or roots were detected in transgenic lines compared to the control (Figure S2) using the protocol described by Haeck et al. (2018). To obtain more insight into the effects of either ICM or CM on phenolic metabolism in rice roots, a comparative liquid chromatography‐mass spectrometry (LC‐MS, see Experimental procedures) profiling was performed on the methanol extracts from HoICM‐expressing plants and the control line (referred to as the HoICM set) and from HoCM_FULL, HoCM_CAT, and the control line (referred to as the HoCM set). Taking both sets together, approximately 13,884 compounds (as estimated by the number of m/z feature groups) were profiled. The abundances of these compounds were analysed by applying an analysis of variance (ANOVA) model separately for each set, and correcting for the large number of tests by calculating the FDR (α = .01). The numbers of so‐obtained differential compounds for the HoICM and HoCM sets were 53 (Table S3) and 17 (Table S4), respectively. In the HoCM set, 1 and 16 compounds displayed higher and lower levels in the transgenic lines as compared to the control lines. However, all of them were of low abundance, hindering the recording of good‐quality MS/MS spectra for structural elucidation (Table S5). For the HoICM set, 7 and 46 compounds showed higher and lower levels in the transgenic lines as compared to the control lines, respectively. The most differentially accumulating compound in the HoICM line was structurally characterized as 1‐O‐ethyl‐6‐O‐glycoloyl hexose. In addition, another of the top accumulating compounds was very similar in the MS/MS spectrum and was identified as (1‐O‐ethyl)hexosyl(6‐>)ethylene glycol glycolate. Both compounds indicate that expression of HoICM leads to the accumulation of glycolic acid. In addition, many of the compounds that showed diminished levels in the HoICM lines were probably fatty acid‐derived molecules that contained a glycine moiety based on the chemical formula and the detection of glycine‐representing MS/MS product ions (Table S5 and Figures S3 and S4).

Because it proved to be difficult to characterize the unknown compounds that were statistically different in abundance between treatment and control, as an alternative approach only those metabolites that could be matched with a known compound were taken into consideration. To increase confidence in the proposed matches, only compounds with a fragmentation score higher than 0.8 were kept. This pipeline resulted in a list of only 17 compounds, 16 of which were linked to the phenylpropanoid pathway. Although not significantly different, for HoCM‐expressing lines 13 compounds had lower levels in the two independent lines per construct, compared to the control, suggesting down‐regulation of the phenylpropanoid pathway. A two‐way ANOVA using “relative abundance of each component” as variable and “control vs transgenic” and “independent line” as fixed and random factors, respectively, was conducted. A general lower abundance of compounds from the phenylpropanoid pathway in HoCM‐expressing lines (p = .008) was observed without any effect of the “independent line” factor (p = .5). Because there was an interaction effect between the fixed and the random factor for HoICM‐expressing plants, no clear conclusion could be drawn (see Table 1).

TABLE 1

Relative values for characterized secondary metabolite abundance in rice lines expressing HoCM or HoICM compared to the control line

	Compound ID	HoCM_FULL2	HoCM_CAT2	HoCM_CAT3	HoICM1	HoICM4
Phenylpropanoid	Tricin(4'‐O‐8)G‐5/7‐O‐hexoside 3	0.4023	0.2582	1.1936	4.6489	1.2003
Phenylpropanoid	Vanillic acid 1	0.3996	0.5069	0.4375	1.2130	0.4053
Phenylpropanoid	Syringic acid 4‐O‐hexoside	0.3837	0.4415	0.4575	1.2970	0.7224
Phenylpropanoid	Vanilloyl hexose	0.1163	0.1129	0.1869	0.6548	0.6579
Phenylpropanoid	4‐Hydroxybenzoic acid	0.3824	0.4140	0.6249	2.1176	1.6102
Phenylpropanoid	Feruloyl quinate 1	0.3513	0.6880	0.5078	1.2792	0.5786
Phenylpropanoid	Ferulic acid 4‐O‐hexoside 2 (‐dehydrated hexose)	1.0309	0.9846	0.5018	1.4663	0.5709
Phenylpropanoid	Sinapoyl hexose	0.4679	0.2185	0.6017	1.0430	1.0901
Phenylpropanoid	p‐Coumaric acid	0.3871	0.4913	0.6304	1.7322	1.0924
Phenylpropanoid	Caffeoyl shikimic acid 1	0.3424	0.6480	0.5111	0.5509	0.3112
Phenylpropanoid	t‐Ferulic acid	0.3413	0.3526	0.3254	0.6086	0.0763
Phenylpropanoid	m‐Hydroxycinnamic acid	0.3302	0.4654	0.5502	1.7059	1.1033
Phenylpropanoid	G(red8‐5)G + hexose	0.2210	0.8074	0.6627	0.6986	0.2126
Phenylpropanoid	G(8‐8)G hexoside 1	0.3728	0.9640	0.5749	0.5736	0.2960
Phenylpropanoid	G(red8‐5)G 8‐O‐hexoside	0.3414	0.3420	0.6841	3.3869	0.7431
Phenylpropanoid	G(8‐5)G [−H2O] 2	0.3202	0.2526	0.5464	1.1624	0.3864
		p = .008			p = .865
Jasmonate	12‐Hydroxy jasmonic acid sulphate	0.4016	0.5912	0.4852	0.8332	0.7970

Compounds that have a lower abundance compared to the control have a value < 1. p values show the result of a two‐way analysis of variance comparing relative abundances of compounds from the phenylpropanoid pathway between transgenic lines and control.

DISCUSSION

Although the presence of a gene encoding CM has been reported in transcriptomes of migratory PPNs before (Bauters et al., 2014; Haegeman et al., 2011), this is the first time it has been functionally characterized in a migratory nematode. Chorismate is the final product of the shikimate pathway, a pathway that has only been reported in plants and microorganisms, but not in animals (Herrmann, 1995). The first report of CM from a nematode was in M. javanica (Lambert et al., 1999), where it was thought to be involved in the establishment of a feeding site within the host. Indole‐3‐acetic acid (auxin, IAA) levels were reported to be lowered due to expression of this CM in plant roots, which had an influence on root development. This observation is in contradiction with the findings that IAA is important for feeding site initiation, while the nematode‐secreted CM was hypothesized to be a vital element in this process (Goverse et al., 2000).

The CM protein produced by H. oryzae contains a unique N‐terminal serine/histidine‐rich domain and a conserved catalytic domain with CM activity. Root‐knot nematodes lack an N‐terminal domain preceding the CM catalytic domain. Cyst nematode CMs do have an N‐terminal domain, but lack homology to the serine/histidine rich domain found in H. oryzae. Vanholme et al. (2009) proved that the N‐terminal sequence was not needed for CM activity in Heterodera schachtii, similar to the here‐obtained results for HoCM. More research is needed to elucidate the function of the N‐terminal region present in cyst nematodes and H. oryzae. Because it is not needed for enzymatic activity, it is possible that this region has another function, for example interaction with host proteins.

A geneencoding ICM was reported in the transcriptome of H. oryzae (Bauters et al., 2014), in R. reniformis (Wubben et al., 2010), and in the genome of the cyst nematode G. rostochiensis (Eves‐van den Akker et al., 2016). As Opperman et al. (2008) suggested while studying the genome of M. hapla, ICM has probably been acquired by horizontal gene transfer from bacteria. The sequence of ICM does not contain a predicted N‐terminal secretion signal, but the fact that the gene is conserved in other PPNs, while being absent in nematodes with a different lifestyle (Bauters et al., 2014), led us to believe it might be involved in plant parasitism. This hypothesis is strengthened by the fact that ICM is also present in the secretome and in the culture supernatants of phytopathogenic fungi, helping to increase susceptibility of the host (Soanes et al., 2008; El‐Bebany et al., 2010; Ismail et al., 2014). ICM is probably secreted by a nonclassical secretory pathway in fungi (Liu et al., 2014), but some fungal ICMs are predicted to have a secretion signal (Ismail et al., 2014). Further research is necessary to provide experimental evidence of ICM secretion by nematodes.

Although it is generally assumed that SA mediates resistance against biotrophic pathogens (Glazebrook, 2005) and H. oryzae kills the cells it is feeding from, in contrast to sedentary nematodes (De Vleesschauwer et al., 2013), SA was shown to be important for an adequate defence response against H. oryzae infection in rice (Nahar et al., 2012). On the other hand, the SA biosynthesis pathway appears changed after H. oryzae infection with up‐regulation of OsPAL1 and down‐regulation of OsICS in shoots (Kyndt et al., 2012). SA biosynthesis in plants can occur through two different pathways, both starting from chorismate as a precursor (Vlot et al., 2009; reviewed in Lefevere et al., 2020). A secreted CM from U. maydis was able to reduce SA accumulation upon pathogen infection (Djamei et al., 2011). ICM, secreted by two fungal pathogens (Verticillium dahliae and Phytophthora sojae), reduces the SA content at least 3‐fold upon transient expression in Nicotiana  benthamiana leaves, rendering the plant more vulnerable to various pathogens (Liu et al., 2014). According to the results presented in this research, HoICM and HoCM do not have an effect on SA. However, upon interpreting these results, one has to keep in mind that the RNA‐Seq and metabolome data presented here came from plants that were not challenged with a pathogen, in contrast to the data of Liu et al. (2014) and Djamei et al. (2011). A GO enrichment analysis performed on RNA‐Seq data from HoICM‐, HoCM_FULL‐, and HoCM_CAT‐expressing rice lines, together with a metabolome analysis, points towards changes in phenylpropanoid content. The GO term “secondary metabolic process” and its child term “phenylpropanoid metabolic process” were reduced in the three different transgenic lines. Both ICM and CM act upon chorismate (or chorismate‐derived metabolites), which is the starting point of the phenylpropanoid pathway (Dempsey et al., 2011). These data were backed up by results extracted from the metabolome analysis on the transgenic lines. Looking at identified compounds, phenylpropanoid content was lower in plants expressing HoCM, whereas results for plants expressing HoICM were less clear. The fact that the PAL genes, considered as the starting point of the phenylpropanoid pathway, are often up‐ and down‐regulated in resistant and susceptible interactions, respectively, shows that phenylpropanoids are important defence molecules. This was shown in several plant–pathogen interaction experiments (Gao et al., 2008; Uehara et al., 2010; Feng et al., 2011; Khanam et al., 2018). In contrast to our observations, Djamei and colleagues (2011) reported an increase in phenylpropanoid biosynthesis products in plants infected with U. maydis, which indicates that chorismate is directed into the phenylpropanoid pathway. CM originating from U. maydis was also able to reduce the SA content of maize upon infection (Djamei et al., 2011). A recent study showed that transient expression of CM originating from the PPN M. incognita was able to reduce the SA content of tobacco leaves after Phytophthora infection (Wang et al., 2018). Our results showed that unchallenged rice plants expressing HoCM or HoICM did not change their SA content in roots or shoots. This observation was strengthened by the fact that no SA‐responsive genes (e.g., OsWRKY45 or OsNPR1) were differentially regulated according to the transcriptome analysis.

Genes involved in “Diterpenoid phytoalexin metabolic process” were inhibited in HoCM_FULL‐ and HoCM_CAT‐expressing plants. Diterpenoid phytoalexins (DPs) have several functions in plants. Their antimicrobial properties have been well documented. DPs accumulate in infection spots where necrotic tissue starts to develop (Huffaker et al., 2011). Transgenic lines deficient in OsCPS4, a phytoalexin biosynthetic gene, showed increased susceptibility to rice blast fungus (Toyomasu et al., 2014). Also in nematodes it was shown that induction of DP biosynthesis increases resistance against nematodes (Verbeek et al., 2019). Some terpenoids might be involved in nematode resistance by accumulation in lesions near the infection site (Veech, 1982), so a reduced terpenoid content could influence susceptibility to nematodes. The mechanism by which expression of HoCM or HoICM influences the diterpenoid phytoalexin biosynthesis pathway remains unclear. So far, it was only suggested that the phenylpropanoid and terpenoid pathway compete for carbon and that carbon flow into these two pathways is tightly regulated (Xie et al., 2008).

The data generated in this paper suggest that the immune system of the host is altered by lowering secondary metabolite content upon secretion of CM and ICM by H. oryzae. Further research is needed to pinpoint the exact metabolites responsible for the drop in immunity and the mechanism by which HoCM and HoICM are able to do this.

EXPERIMENTAL PROCEDURES

Nematode DNA extraction and gene amplification

Genomic DNA was extracted from a batch of mixed stages of H. oryzae or a single hand‐picked nematode as previously described by Bolla et al. (1988), with some minor modifications. Briefly, nematodes were sonicated three times for 10 s to break cell membranes. Samples were incubated for 1 hr at 65°C after which DNA was isolated using phenol and chloroform. Isolated DNA and a cDNA library (Bauters et al., 2014) were used to amplify HoCM and HoICM by PCR with gene‐specific primers (see Table S6). PCR fragments were cloned into pGEM‐T (Promega) according to the standard thymine/adenine protocol and sequenced by the Sanger method (LGC) to verify the insert.

Bioinformatic analysis

Several blast programs (blastn and blastx), available on the NCBI server (http://www.ncbi.nlm.nih.gov/), were used for sequence similarity searches in EST and nonredundant protein databases. Alignments were generated with BioEdit (Hall, 1999). Putative N‐terminal secretion signals were detected by SignalP 4.1 (Petersen et al., 2011). Protein domains were identified by searching the Pfam 27.0 database (Punta et al., 2011). Protein 3D models were constructed with the Protein Homology/anologY Recognition Engine v. 2.0 (Phyre2) (Punta et al., 2011) and visualized by PyMOL 1.6 (http://www.pymol.org/).

Activity assays

Two HoCM fragments (HoCMFull and HoCMCat) were amplified by PCR using the full cDNA sequence as template (primer sequences can be found in Table S6). The first fragment comprised the full sequence of CM without sequence encoding the signal peptide, the second fragment only contained the sequence encoding the catalytic domain (base pair 258 to 747). The complementation assay was performed as described by Vanholme et al. (2009).

HoICM activity was tested using entB‐mutated AN192 E. coli cells. AN192 cells are unable to grow in iron‐limiting conditions. HoICM was amplified by a two‐step PCR starting from cDNA as template to add attB sites to the fragment (primer sequences can be found in Table S6). The fragment was cloned into the destination vector pDEST17 using the standard Gateway technology protocol (Life Technologies). The pDEST17 vector carrying HoICM was transformed via heat shock into E. coli AN192 cells. E. coli AN192 cells transformed with the empty vector were used as negative control. These cells were grown overnight at 37°C in Luria‐Bertani (LB) medium with 100 mg/L carbenicillin. Optical density (OD) was measured at 600 nm and equal amounts of bacteria were added to fresh medium containing different concentrations of 2,2ʹ‐bipyridyl to create iron‐limiting conditions (0, 100, 150, and 200 µM). These different types of media were inoculated with E. coli AN192 cells carrying HoICM or the empty vector. OD was measured at 600 nm at 15, 24, and 48 hr postinoculation. Results were statistically analysed with a Mann–Whitney test. The experiment was performed twice.

Plant transformation

cDNA sequences of HoICM and HoCM were amplified by PCR on a cDNA library (Bauters et al., 2014). Primer sequences are provided in Table S6. Expression vectors (pMBb7Fm21GW‐UBIL, obtained from the Center for Plant Systems Biology [VIB, Belgium]) were created for three different constructs: HoICM (containing the complete open reading frame), HoCM without secretion signal (HoCM_FULL), and only the catalytic domain of HoCM (HoCM_CAT). Vectors were created by standard Gateway cloning. Cloned sequences as well as position in the vector were verified by Sanger sequencing. Rice transformation was done as described by Ji et al. (2015). Leaf samples were taken and DNA was extracted to confirm the presence of the gene of interest by PCR using gene‐specific primers (Table S6). DNA was extracted by heating a leaf sample at 95°C for 15 min in 400 µl of extraction buffer (200 mM Tris.HCl, 250 mM NaCl, 25 mM ethylenediaminetetraacetic acid [EDTA] and 0.5% sodium dodecyl sulphate [SDS]). Afterwards 50 µl chloroform was added. The sample was centrifuged and 200 µl of isopropanol was added to 200 µl of the upper aqueous phase to precipitate the DNA.

RNAextraction, sequencing, and validation

Rice plants were grown under greenhouse conditions at 30°C. In the flowering stage, leaf samples were taken from three plants of a line transformed with the empty vector, and three plants of a line expressing HoCM_FULL or HoCM_CAT. Due to limited seed availability, the line for which most seeds per transformant were available was chosen. Because of technical reasons only two plants were used for lines expressing HoICM. RNA was extracted using the Nucleospin RNA extraction kit (Machery‐Nagel) according to the manufacturer's guidelines with an additional DNase I (Thermo Scientific) treatment. At least 4 µg of RNA of each sample was sent for RNA sequencing (Microsynth). An mRNA library was constructed using the Illumina TruSeq RNA sample preparation kit. The library was sequenced using a NextSeq platform (Illumina). To validate RNA‐Seq results, a set of genes was selected randomly to perform RT‐qPCR. RNA was extracted with the RNeasy Plant Mini kit (Qiagen). cDNA was synthesized and a RT‐qPCR was performed with EIF5c and EXP NARCAI as reference genes to normalize expression results. Primers can be found in Table S6.

Mapping and gene expression profiling

Reads from all libraries were mapped to the O. sativa ssp. japonica reference genome (IRGSP‐1.0) using the Spliced Transcripts Alignment to a Reference (STAR, v2.3.1) software with standard settings (Dobin et al., 2013). The GenomicFeatures (relying on Ensemble's BioMart osativa_eg_gene dataset) and GenomicAlignment R/BioConductor packages were used to obtain an overview of all exons per gene (makeTranscriptDbFromBiomart and exonsBy functions), respectively, to quantify expression per sample for these genes starting from the STAR output (readGAlignmentsFromBam and countOverlaps functions). Additional gene annotation information was obtained from Ensembl (BioMart, canonical identifiers only). Differential gene expression was assessed using the R‐package “EdgeR” version 3.8.5 (Robinson et al., 2010). Data were analysed using the generalized linear model (GLM) approach. Data were normalized using normalization factors calculated by the Trimmed Median of M‐values method with standard settings as implemented in the EdgeR package (Robinson and Oshlack, 2010). The proposed GLM was fitted to the design matrix using the Cox‐Reid profile‐adjusted likelihood method (McCarthy et al., 2012). Differential gene expression levels (Log2 fold change) were calculated and Benjamini and Hochberg adjusted (Benjamini and Hochberg, 1995) false discovery rates (FDR) were estimated to detect DEGs (FDR < 0.05). GO enrichment analysis was performed using agriGO (Du et al., 2010). Gene identifiers with their corresponding Log2FC were used as input in a Parametric Analysis of Gene Set Enrichment (PAGE) (Kim and Volsky, 2005). The Hochberg multitest adjustment method was performed (p < .05).

Phenolic profiling

Roots (approximately 5 mg of dry weight) were homogenized in liquid nitrogen and extracted with 1 ml of methanol. The methanol extract was then evaporated and the pellet dissolved in 200 μl of water/cyclohexane (1:1, vol/vol). Then, 10 μl of the aqueous phase was analysed via reverse‐phase ultrahigh performance liquid chromatography using an Acquity UPLC BEH C18 column. In addition to full MS analysis, a pooled sample was subjected to data‐dependent MS/MS analysis (exclusion duration = 10 s). Integration and alignment of the m/z features were performed via Progenesis QI software v. 2.1 (Waters Corporation). The normalization was set on “external standards” and was based on the dry weight of the samples.

For statistical analysis, the data were split into two sets. A first set comprised the two independent HoICM expression lines (HoICM1 and HoICM4) and the corresponding control line (HoICM set). In the second set, one HoCM_FULL expression line, two independent HoCM_CAT expression lines (HoCM_CAT2 and HoCM_CAT3), and the corresponding control line were included (HoCM set). Both sets were subjected to a one‐way ANOVA using the lm() function in R v. 3.4.2 (R_Core_Team, 2017) with a weighting factor included for the first set. For post hoc testing, the pairwise.t.test(p.adjust = ”bonferroni”) function was applied. The false discovery rate (α = .01) was computed on the ANOVA model significance using the p.adjust(method = ”fdr”) function. Additional filtering on the prerequisites of (a) a significant post hoc test for each of the transgenic lines and (b) abundance changes in the same direction for all transgenic lines as compared to the control lines yielded 54 and 17 differential m/z features (corresponding to 53 and 17 m/z feature groups) for the HoICM and HoCM sets, respectively.

Structural annotation was performed using a retention time window of 1 min, and both precursor ion and MS/MS identity searches. The precursor ion search (10 ppm tolerance) was based on a compound database constructed via instant JChem (ChemAxon), whereas MS/MS identities were obtained by matching against an in‐house mass spectral database (200 ppm fragment tolerance). As none of the MS/MS spectra of the differential m/z features could be identified via database matching, MS/MS spectral elucidation was attempted (see File S2) using gas phase fragmentation rules and in silico MS/MS elucidation software, that is, CSI:FingerID (Böcker and Dührkop, 2016; Dührkop et al., 2015) and CFM‐ID (Allen et al., 2015). A more detailed description of the phenolic profiling protocol can be found in File S2.

Nematode infection assay

The M. graminicola culture, originally isolated in the Philippines, was maintained on O. sativa "Nipponbare" in potting soil at 27°C (16/8 hr light regime). H. oryzae was obtained from infected rice roots sampled from fields in Myanmar. Nematodes were extracted from rice roots using a modified Baerman funnel technique. Rice seeds were first germinated on wet tissue paper at 27°C in the dark for 3 days. Germinated seedlings were transferred to PVC tubes containing sand and absorbent polymer substrate (Reversat et al., 1999). Two‐week old plants (T1 or T2 generation) were infected with 200 (M. graminicola) or 300 (H. oryzae) nematodes/plant. Two weeks after infection, nematode susceptibility was assessed by counting the number of galls (M. graminicola) or nematodes (H. oryzae) per root. Galls/nematodes were visualized by boiling roots in 0.8% acetic acid and 0.013% acid fuchsin for 3 min, after which they were washed under running tap water and destained in acid glycerol. Infection assays with H. oryzae were only conducted on a single line per construct due to difficulties in culturing H. oryzae under laboratory conditions.

FIGURE S1 Validation of RNA‐Seq results. Expression values were calculated compared to expression in the empty vector control (expression level set at 1). Normalization was done with two reference genes (EIF5C and EXP NARCAI). Expression values are plotted on a Log2‐scale. Expression values are a mean of two biological replicates (calculated by taking the average of three technical replicates). Statistical analysis was performed with REST 2009 software (asterisk indicates a significant difference, p < .05). Light grey bars indicate genes for which expression could not be detected in the respective overexpression line, expression values of these genes are equal to zero. Uridine diphosphate (UDP)‐glyctrsfr, UDP‐glycosyltransferase; GlutS, glutathione S‐transferase; ser/thr_kin, serine/threonine‐protein kinase; hist, histone H2B.3; LRR, LRR receptor‐like serine/threonine‐protein kinase; R‐like PK, receptor‐like protein kinase; α‐galac, α‐1,2‐galactosyltransferase; RIK, protein RIK; PM_ATPase, plasma membrane ATPase; gluc6, glucose‐6‐phosphate/phosphate translocator 2. In between brackets the log2FC according to RNA‐Seq is given. Bars represent standard errors, calculated with REST2009 software using Taylor’s series

Click here for additional data file.

FIGURE S2 Root (A) and shoot (B) tissue of 3‐week‐old plants was harvested to quantify SA content on a ultra high performance liquid chromatography (U‐HPLC) system. Plant material was homogenized by grinding in liquid nitrogen (three individual plants per line were used). Extraction of 100 mg of plant material was performed using the modified Bieleski solvent after which it was filtrated and evaporated. Chromatographic separation was performed on a U‐HPLC system equipped with a Nucleodur C18 column (50 × 2 mm, 1.8 μm dp), using a mobile phase gradient consisting of acidified methanol and water. Mass spectrometric analysis was carried out in selected‐ion monitoring (SIM) mode with a Q Exactive™ Orbitrap mass spectrometer (Thermo Scientific), operating in both positive and negative electrospray ionization mode at a resolution of 70,000 full‐width at half maximum. Data were analysed with a nonparametric pairwise Kruskal–Wallis test. No significant differences were detected among the different lines. Ctrl, empty vector line

Click here for additional data file.

FIGURE S3 Representation of the fragmentation pathways of the compound eluting at 3.26 min (id: 3.26_265.0921m/z)

Click here for additional data file.

FIGURE S4 Representation of the fragmentation pathways of the compound eluting at 4.33 min (id: 4.33_309.1185m/z)

Click here for additional data file.

TABLE S1 Table showing the Ct values for two control genes (Exp Narcai and EIF5C) and the transgene (A, HoCM; B, HoICM) that is expressed in the different transgenic lines and controls (wild‐type Nipponbare and empty vector line). Values that are underlined indicate that there was no expression for that gene. NA, not available

Click here for additional data file.

TABLE S2 Overview of the output of RNA‐Seq on the different transgenic lines

Click here for additional data file.

TABLE S3 Differential compounds in the transgenic lines overexpressing HoICM. FC, fold change; FDR, false discovery rate

Click here for additional data file.

TABLE S4 Differential compounds in the transgenic lines overexpressing HoCM. FC, fold change; FDR, false discovery rate

Click here for additional data file.

TABLE S5 MS/MS spectral information for a selection of compounds which are differentially abundant in HoCM or HoICM expressing lines. ID, compound identifier; tR, retention time; m/z, mass to charge ratio; Δppm, deviation of measured mass; MS/MS, relative abundances to the base peak of the MS/MS product ions are indicated between parentheses

Click here for additional data file.

TABLE S6 Table with primer sequences used in this research

Click here for additional data file.

FILE S1 List of differentially expressed genes in transgenic lines compared to the control according to the RNA‐Seq data

Click here for additional data file.

FILE S2 Detailed protocol describing the metabolomics experiment

Click here for additional data file.