Analysis of Gene expression in soybean (Glycine max) roots in response to the root knot nematode Meloidogyne incognita using microarrays and KEGG pathways

Plant and nematode procurement

Glycine max cv Williams 82 and M. incognita population LESREC (Lower Eastern Shore Research and Education Center) were grown in a greenhouse at the United State Department of Agriculture Soybean Genomics and Improvement Laboratory, Beltsville, MD, USA. M. incognita eggs were harvested from roots of G. max cv Williams 82 2-4 months after inoculation using a method modified from those previously described in Meyer et al.[22] and Nitao [23]. Soybean seedlings were grown in Promix (Premier Horticulture INC., Quakertown, PA, USA) for one week in 20 × 20 × 10 cm flats, then moved to sand (The Stone Store, Hanover, MD). Three thousands eggs were used to inoculate roots of 7 day old soybean seedlings (cv. Williams 82). Soybean roots at 12 dai, 10 wai, and control uninfected plants were washed with sterile water, flash frozen in liquid nitrogen, ground to a fine powder and frozen at -80°C until use. The infected roots were collected at 12 days after infection.

Nematodes were stained in infected roots using a modified protocol of Byrd et al.[24] and Mahalingam et al.[25]. Briefly, roots were washed in gently flowing tap water to remove soil and debris, cut to 2 cm segments, and placed in a small beaker, then soaked in 20-30 ml of 10% commercial Clorox (chlorine bleach, 5.25% NaOC1) for 3 min. The roots were rinsed in tap water and then transferred into a 50 ml glass bottle containing 20 ml of distilled water and left to boil in a microwave with loosened caps. A 500 μl of acid fuschin (0.15 g/10 ml H2O) and 500 μl of glacial acetic acid were added to the root samples and heated to boiling in a microwave twice. The roots were left to cool to room temperature before removing the excess stain with running tap water using Miracloth on the top of the bottle. A 20 ml of clearing reagent (1/3 lactic acid+ 1/3 glycerol + 1/3 distilled water) were added to roots and roots were left to destain for two hours to overnight. The nematodes were stained red as observed in the roots under a dissecting microscope. General chemical reagents were obtained from Sigma Chemical Co (St. Louis, MO).

RNA extraction and microarray analyses

RNA was extracted from 100 mg each of the three different root samples using the Ultra Clean Plant RNA Isolation Kit (MOBIO, Carlsbad, CA). Gene expression analysis was performed using the GeneChip^® Soybean Genome Array (Affymetrix^®, Santa Clara, CA, USA) containing more than 37,500 probe sets as described in Klink et a l. [9, 10]. In this GeneChip technology, each high density spot (gene) is represented by 11-probe pairs (11 μm feature size), which allows multiple independent measurement for each transcript. GeneChip^® Soybean Genome Array details are available at the Affymetrix^® website [26]. The microarrays were hybridized and scanned at the Laboratory of Molecular Technology, SAIC-Frederick, National Cancer Institute at Frederick, Fredrick, MD, USA.

Affymetrix^© soybean Genechip data was analyzed as described in Klink et al.[9, 10, 15] with additional analysis using PAICE [18]. Briefly, microarray gene expression data was imported into MATLAB (Natick, MA) Bioinformatics Toolbox. Normalization of the probe sets was performed using RMA (Robust Multiarray Analysis). The resultant calculated output was the log base 2 of the expression values, enabling scaling of the dataset. Volcano plots were produced, which graphically illustrate gene expression fold-change with respect to statistical significance. The plots were produced using fold-changes >= |2.0| and p-values < = 0.05 with respect to the control (untreated roots of cv. Williams 82). The t-test was used in calculating p-values. False Discovery Rate analysis was further utilized against significantly expressed genes [16]. The False Discovery Rate tool 'Significance Analysis of Microarrays' (SAM v3.0; [8] was used with a false discovery rate of 10%.

RT-PCR for confirmation of differential gene expression in infected soybean roots at 12-dai and 10-wai

Pathway Analysis

Biochemical pathway analysis was conducted using PAICE (P athway A nalysis and I ntegrated C oloring of E xperiments) [18]. This software program maps expression levels of genes encoding enzymes found in the KEGG biochemical pathways database. Gene expression levels are denoted using color codes displayed at the pathway nodes depicted by enzyme EC numbers. Besides the pathway mapping feature, PAICE colors EC accessions using gradients of green and red to represent induced and suppressed gene expression respectively. The colors can be tailored to the analysis through the software package.

PAICE has the unique ability to color in yellow the expression of genes having multiple family members that lack a consensus in gene expression, i.e. some members are over-expressed and others are under-expressed.

Histological examination of RNK infection

At 12 dai, galls can be identified as small swellings along the soybean root (Figure 1a). Within the gall the nematode has started feeding and can be visualized by staining with acid fuchsin [24, 25] to monitor nematode invasion and development inside the roots (Figure 1b). Mature galls are present on soybean roots at 10 wai (Figure 1c). Within the gall, mature female M. incognita can be identified easily by staining (Figure 1d).

Transcript profiling of galls formed by M. incognita infection

A comparison of gene expression at12 dai compared to control led to the identification of 1867 genes with greater than 1.5-fold change in expression (Figure 2). Of these, 1278 genes increased and 589 genes decreased in expression. Transcripts encoding leghemaglobin C1 (AI973819) increased the most at 386-fold. The most down-regulated gene was BF070134 with homology to a putative senescence protein 12 and to ERD7; its transcripts were 77-fold lower than in the control. There were 2108 genes with altered expression in galls at 10 wai. Of these, 1460 genes increased in expression and 648 genes decreased in expression. The transcript of the gene encoding pathogenesis related protein PR1a (CK605838) increased the most at 258-fold. As in the 12 dai experiment, the most down-regulated gene was BF070134 with transcripts 172-fold lower than the control. When gene expression at 10 wai was compared directly to 12 dai, 827 genes were up-regulated, while 535 genes were down-regulated. In this case, transcripts of the gene encoding the cysteine-rich plant defense protein, defensin (PDF2.3; BI321308), increased the most at 63-fold, while the transcripts of the gene encoding xylene serine peptidase 1, subtilase (AW309540) decreased the most at 126-fold.

Mitosis and cell division

Cell wall modifications

Numerous genes involved in cell wall remodeling were shown to be differentially expressed after infection with M. incognita (Table 4). Genes encoding four members of the expansin enzyme family (expansins A1, A5, A16, and 45) were up-regulated at both time points (Table 4). Also, genes encoding endo-1,3-beta-glucanase [EC 3.2.1.60] family members were differentially expressed (Table 4). Most of them were up-regulated, while one member (AW831345) was down-regulated 27-fold in the 12 dai and 44-fold in the 10 wai. A gene encoding the cell wall-modifying xyloglucan endotransglycosylase/hydrolase [EC 2.4.1.207] was down-regulated at both time points, while the gene encoding endoxyloglucan transferase A2 was up-regulated at both time points (Table 4). Expression of a gene encoding pectin esterase [EC 3.1.1.11] increased 24-fold and 47.5-fold at 12 dai and 12 wai, respectively. In addition, a gene encoding cellulose synthase [EC 2.41.12] was down-regulated by 9.3-fold and 4.5-fold at 12 dai and 10 wai, respectively (Figure 3).

Probe set	Genbank accession No.	12 dai_Fold Change	10 wai_Fold Change	Gene Family
Gma.13110.1.S1_at	CD394837	24.68157	63.82636	Expansin A1
Gma.7784.1.A1_at	BM269985	9.148853	14.15569	Expansin A1
GmaAffx.90097.1.S1_at	CF805822	6.665491	7.075862	Expansin A16
GmaAffx.67477.2.S1_at	AW509184	6.36759	3.798536	Expansin A5
Gma.14216.1.A1_at	CD417846	No change	1.854507	Expansin 45
GmaAffx.86629.2.S1_at	BQ742947	10.92815	33.72573	Endo-1,3-beta-glucanase
Gma.6574.1.S1_at	BQ628332	5.983501	4.524413	Beta-1,3-glucanase
GmaAffx.71901.1.S1_s_at	BE022356	No change	4.844002	Beta-1,3-glucanase
Gma.2801.1.S1_at	BI785739	No change	4.21076	Beta-glucosidase
GmaAffx.91055.1.S1_at	CF806780	4.618147	3.953466	Beta-1,3-glucanase
Gma.17299.2.S1_at	BG237465	No change	3.947161	Beta-1,3-glucanase-like protein
GmaAffx.87027.1.S1_at	AW831345	-27.19943	-43.94289	Beta-1,3-glucanase 2
GmaAffx.89396.1.S1_x_at	CK605572	-6.914317	-9.101676	Xyloglucan
GmaAffx.34697.1.S1_at	AW433263	-16.99092	-30.42019	Endoxyloglucan transferase A4

Also, in phenylpropanoid biosynthesis, genes encoding a family of 21 extensin peroxidases [EC 1.11.1.7] that participate in lignin biosynthesis were differentially regulated (Figure 4). The extensin gene with the highest increase in expression (AW309606) increased by 95-fold while the extension gene with the largest decrease in expression (BI970840) decreased by -16-fold.

Carbon and energy metabolism

The expression of numerous genes encoding enzymes in glycolysis, the tricarboxylic acid cycle and in amino sugar synthesis was altered. For example, the gene encoding UDP-glucuronate 4-epimerase [EC 5.1.3.6] was greatly down-regulated with a 21-fold decrease in transcript abundance. Also, the gene encoding GDP mannose 4,6-dehydratase [EC 4.2.1.47] was down-regulated 20.5- and 5.3-fold at 12 dai and 10 wai (data not shown). In the glycolysis pathway, many genes were up-regulated during infection, including genes encoding glucose-6-phosphate isomerase [EC 5.3.1.9], transcripts of which were increased by 14-fold at 12 dai (Figure 5)

Defense-related genes

There are multiple changes in expression of genes encoding enzymes of the alpha-linolenic acid pathway leading to several important defense-related compounds, including jasmonic acid (Figure 6). At 12 dai, the change in expression of genes encoding enzymes leading to jasmonic acid seems conflicting. For example, the gene encoding palmitoyl-CoA hydrolase [EC 3.1.2.2], an enzyme that is needed for the biosynthesis of the alpha-linolenic acid, is suppressed. Linolenic acid can be a substrate for lipoxygenase [EC 1.13.11.12]. Genes encoding six lipoxygenase family members were differentially expressed. Mostly of them were up-regulated with highest induction of 22-fold for one member (AW705829). One gene family member (BI426746) was down-regulated by 3.8-fold. In addition, some gene family members encoding allene oxide synthase (hydroperoxide dehydratase [EC 4.2.1.92]), leading to jasmonic acid synthesis, were up-regulated, while others were down-regulated. The abundance of transcripts of the gene encoding the next enzyme in the pathway, allene oxide cyclase [EC 5.3.99.6] was also decreased. Yet, a gene encoding OPDA reductase [EC 1.3.1.42] was up-regulated. In contrast, most of the genes encoding enzymes involve in jasmonic acid biosynthesis were down-regulated at 10 wai.

Genes encoding transcription factors and other proteins

Changes in gene transcripts were accompanied by changes in expression of transcription factors, especially those in the WRKY family of transcription factors. Our microarray results indicated that genes encoding several family members of WRKY genes were down-regulated at 12 dai, including genes encoding WRKY6, 15, and 22 (-10.5, -14.3, and -8.5 FC, respectively). In contrast, at 10 wai, genes encoding WRKY 21 and 70 were up-regulated at 117 and 42 FC, respectively.

Several pathogenesis-related proteins are induced in plants during infection with any pathogen or by wounding, including nematode infection, and induction of many of these is affected by salicylic acid, jasmonic acid or ethylene [29]. In our microarray data, genes encoding pathogen related proteins (PR) such as PR3 were down-regulated at 12 dai and genes encoding PR3 at 10 wai showed a mixed response; some were up-regulated while others were down-regulated. The three copies of the pathogen related protein PR1 gene were over-expressed by 78.23, 97.56, and 138.50 fold, respectively (Data not shown).

Confirmation of differential gene expression by quantitative PCR

Quantitative PCR was conducted to confirm gene expression patterns revealed by microarray analysis. We measured transcript abundance of 14 genes that showed increased or decreased transcript abundance by microarray analysis (Table 5; Figure 7). The trends in up- or down-regulation of gene transcripts were consistent between microarrays and quantitative PCR results except for expression of the gene encoding lipoxygenase family member LOX1 at 10 wai. However, we did observe differences in levels of expression between methods. Differences in fold-change in gene expression as measured by microarray and qRT-PCR have been reported in previous studies [15, 18, 28].

Gene	Fold Change comparison for 12 dai		Fold Change comparison for 10 wai
Jasmonic acid	RT-PCR	Microarray	RT-PCR	Microarray
LOX1	65.3	22.6	42.2	-9.9
AOS	0.0	-59.3	0.0	-77.4
AOC	0.1	-6.3	0.3	0
OPR2	9.1	5
OPR3			0.1	-8.8
Cell Wall
PECT.	49.4	47.5
CELL.	1.0	-4.6
CCOA-OMT12	25.0	2.6
CCOA-OMT10			1.5	1.9
C3H12	2.7	2.9
C3H10	0	0	9.4	8.4
F5H	78.6	15.4	88.9	20.5
Cell Cycle
CDKB2	16.7	5.2	5.8	3.1
CYCD3	1.0	0	0.2	-2.9

Note: Gray boxes refer to untested genes in a specific time point.

Regulators of the cell cycle and cell division

The cell cycle is regulated by two types of cyclin-dependent kinases. CDKA is required for cells to enter the S and M phases. CDKB1 and CDKB2 are expressed during the G2 and M phases and are responsible for the G2-M transition [30]. Our microarray results indicate that genes encoding some members of the cyclin-dependent kinases family were differentially expressed at 12 dai and 10 wai. Over-expression of the gene encoding CKB2 at 12 dai (5.2 and 4.6 FC) correlates with the increase in plant nuclear division that occurs at the infection site due to M. incognita infection and feeding. Cells selected by M. incognita for feeding become multinucleate giant cells. This increase in nuclear division continues until 14 dai [3–5].

Cell cycle activation in giant cells has also been observed by Engler et al.[31]. In that study, the transition from S to G2 and G2 to M phase was reported after the over-expression of a GUS gene driven by the cycB2 or cycA2 promoters at one to nine days after infection with M. incognita. Expression of the CDKB2 gene at 12 dai was higher than at 10 wai, i.e., 5.2- versus 3.1-fold, respectively (Table 1). Ramsay et al. [32] found that cyc D3 is essential to stimulate the G1 phase of the cell cycle in root knot nematode infected giant cells. In this investigation, the two types of CycD3 (CycD3.2 and CycD3.3) were shown to be relatively more strongly expressed as compared to that of LeCycA1.1, LeCycB1.1 and LeCycD3.1 in giant cells induced by Meloidogyne spp. compared with other cyclin-dependent kinases. They observed PCR amplification of CyD3.2 and CycD3.3, while no amplification of cycA.1, CycB1.1 and CycD3.1 was observed. Our data showed a suppression of gene expression of the gene encoding cycD3 which is important for the regulation of the G1-S transition. In addition, at 10 wai we found an increase in gene expression of CKS1 (CDK-Subunit 1), a protein that prevents CDK from driving the cell cycle into S phase. This result suggests that at the earlier time point, the giant cells reach maturity and then the genes required for nuclear division are turned off.

Cell wall modification and remodeling

Due to multiple nuclear divisions of selected cells with no coincident cell division, the giant cells sometimes reach more than 400-times the size of a normal cell and may contain more than one hundred nuclei [4]. This expansion of the giant cell is accompanied by extensive cell wall modification. Our microarray data indicate altered expression of many genes involved in cell wall extension and remodeling (Table 4). We found that the genes encoding a cell wall-modifying xyloglucan endotransglycosylase/hydrolase [EC 2.4.1.207] and endoxyloglucan transferase A2 are differentially expressed in soybean roots after infection with M. incognita (12 dai and 10 wai). These enzymes play a role in softening and breaking down the cell wall [33]. Genes encoding many endo-1,4-glucanase family members were up-regulated at both time points (12 dai and 10 wai). Endo-1,4-glucanase is involved in cell wall remodeling and expansion. LCM was used to isolate giant cells formed in tomato by M. javanica to examine gene expression [34]. Numerous transcripts of genes involved in cell wall remodeling were also identified in the cDNA library of giant cells 4 dai, including transcripts of genes encoding pectin methylesterase and pectinesterase. Goellner et al. [35] identified genes encoding endo-1,4-glucanases that were up-regulated in feeding cells formed by M. incognita and cyst nematode in tobacco plants. Also, Mitchum et al.[36] found that the promoter of an endo-1,4-β-glucanase gene was strongly activated in feeding cells formed by Meloidogyne incognita as indicated by strong promoter-driven GUS expression.

The increase in expression of the gene encoding expansin A in our results is consistent with other investigations, wherein the expansin (LeEXPA5) genes in A. thaliana and tomato were shown to be up-regulated in developing giant cells after infection with Meloidogyne[37, 38]. Moreover, down-regulation of genes encoding cellulose synthase and over-expression of genes encoding pectin esterase that degrades pectin to pectate coincide with a breakdown of the cell wall during the early time points of infection with Meloidogyne. Our results are consistent with those of Jammes et al.[38], who found genes encoding pectin esterases and pectate lyases were activated in Arabidopsis thaliana roots after infection with Meloidogyne incognita and the cell wall loosening process was activated during the development of the giant cell as well. In giant cells formed in tomato by M. javanica, there is an 8-fold and 7.3-fold increase in expression of the gene encoding pectinesterase U1 precursor [34]. Giant cells formed by M. javanica in roots of Arabidopsis were collected by LCM and analyzed by Barcala et al. [39]. Genes encoding cellulose synthase, expansin, pectate lyase, endoxyloglucan transferase also were all up-regulated in these cells coinciding with cytoskeleton rearrangements that occur during giant cell development.

Nutrients supply for M. incognita

The nematode uses a large amount of plant resources to develop and reproduce. This demand for energy and carbon is reflected in the numerous genes involved in glycolysis and gluconeogenesis that are up-regulated in the soybean root (Figure 5). For example, we found many genes encoding enzymes in the glycolysis pathway and amino sugar synthesis to be up-regulated. Mostly, the changes in gene expression occurred early in infection (12 dai). In addition to their roles in pathways that provide energy and carbon for the nematode, some of these genes have an essential role in the soybean-M. incognita interaction. For example, at 12 dai, the gene encoding UDP-glucuronate 4-epimerase (EC 5.1.3.6) is highly down-regulated (-21 fold). In Arabidopsis, a mutation in this gene resulted in hyper-sensitivity to the cyst nematode, Heterodera schachtii[40]. In Azospirillum brasilense, this enzyme is important for lipopolysaccharide synthesis which is important in the bacterium-plant root interaction. A mutation in this gene resulted in the failure of the bacteria to respond to several stresses and antimicrobial compounds. It also affected the ability of the bacteria to form biofilms [41]. This enzyme may be important in allowing the host to respond to M. incognita invasion [42, 43].

In the glycolysis and gluconeogenesis pathways, many genes were shown to be up-regulated, including the gene encoding glucose-6-phosphate isomerase [EC: 5.3.1.9]. The gene encoding this enzyme was also shown to be up-regulated in cucumber plants after treatment with Trichoderma asperrellum T34 [44]. The enzyme is essential in salinity tolerance in the alga Dunaliella salina[45].

Not only do nematodes require large quantities of carbon and energy from its host, they also use starch during juvenile development. Starch is stored in syncytia formed by Heterodera schachtii in roots of Arabidopsis[46]. The authors postulate that the starch is also used to compensate for fluctuating levels of sugar during the course of nematode development and feeding. The high metabolic rate of cells was suggested by the increased expression of ribosomal genes in giant cells induced by M. javanica in tomato roots [47].

Degradation of the cell walls could result in release of sugar which in turn will be channeled to glycolysis as reflected in the activation of the genes encoding enzymes in the glycolytic pathway. Also, since some enzymes in glycolysis participate in the biosynthesis of pentose, purines and pyrimidines, there could be an increase in production of nucleotides required for DNA replication.

Plant defense system

When a nematode invades a plant root, it must repress or control the plant defense response so it can successfully establish its permanent feeding site [4–6]. Our microarray data showed significant changes in expression of genes related to the defense response against pathogens. The pathway leading to jasmonic acid biosynthesis is one of the pathways associated with pathogen resistance and genes in this pathway were significantly affected by Meloidogyne incognita infestation at both time points (12 dai and 10 wai; Figure 6). At 12 dai, six of seven members of the lipoxygenase gene family were up-regulated. Lipoxygenase (LOX) is essential to oxylipin biosynthesis and has an important function in the plant defense response against wounding and pathogen attack [48]. Reduction of LOX activity resulted in an increase in susceptibility of transgenic potato plants to insect attack [49]. Over-expression of the gene encoding lipoxygenase could mean that more 9-HPOTrE would be produced. This is one of the major products of lipoxygenase (Figure 6). Interestingly, 9-HPOTrE is involved in the activation of the plant defense response directly or through its metabolites. In potato plants, 9-HPOTrE is produced in response to injury or infection. The role of 9-HPOTrE in the plant defense response suggests that there may be another pathway of LOX-mediated defense response [-50]. 9-HPOTrE could also be a substrate for allene oxide synthase (Figure 6) to produce OPDA, the precursor for jasmonic acid. At 10 wai, the abundance of the lipoxygenase transcript was much lower than the 12 dai time point. Three of seven gene family members encoding lipoxygenase were down-regulated. Also, all of the allene oxide synthase family members were greatly down-regulated in addition to some other genes encoding enzymes in the jasmonic acid biosynthetic pathway (Figure 6). This indicates that at 12 dai the plant defense system is still struggling to fight the infestation, but after prolonged infection (10 wai) most of the genes that are responsible for the production of one major defense hormone, jasmonic acid, were turned off. Genes in this pathway could be targets for testing whether resistance to nematode infestation can be increased in transformed plants by over-expression of these genes.

We found a number of genes encoding PR proteins that were differentially expressed in soybean roots 12 dai with M. incognita. The genes encoding PR-1, PR-2 and PR-5 protein families increased in expression. Genes encoding PR-2 and PR-5 are expressed in sid mutants of Arabidopsis that do not accumulate SA. However, genes encoding PR1 are known to be induced by salicylic acid (SA) [50, 51]. This suggests that salicylic acid or its derivatives may be synthesized at 12 dai and 10 wai by M. incognita infection. Interestingly, there are two different possible routes to salicylic acid production [52]. The pathway that has the most scientific support involves isochorismate synthase [53] and its genes are not represented on the microarray. The other pathway leading to SA production involves phenylalanine. In this pathway, we found a large increase in the expression of the gene encoding tyrosine aminotransferase [TAT; EC:2.6.1.5]. If the abundance of this enzyme is increased, then this could lead to increased phenylalanine production. Genes encoding phenylalanine ammonia-lyase [EC:4.3.1.24]; and salicylate 1-monooxygenase [1.14.13.-] were over-expressed 6.9 and 2.9 FC, respectively. SA induces the expression of PR-1 [54]. In giant cells of tomato formed by M. incognita 4dai, several genes involved in the phenylpropanoid pathway were detected, notably phenylalanine ammonia-lyase and cinnamyl alcohol dehydrogenase (35). Transgenic tobacco over-expressing PR-1 was more resistant to blue mold and black shank caused by Peronospora tabacina and Phytophthora parasitica f. sp. nicotianae, respectively [29]. Although SA treatment of tomato plants that were inoculated with Meloidogyne incognita did not completely eliminate nematode infection, it enhanced the synthesis of PR-1, which resulted in a significant increase in resistance to the nematode [55]. In contrast, Portillo et al. [56] reported a decrease in transcripts of the PR-1 precursor in giant cells collected from tomato infected with M. javanica by LCM as measured using qRT-PCR. Genes encoding PR-3 and PR-4 family proteins are reported to be up-regulated by jasmonic acid and ethylene [29]. Also, PR-4 showed ribonuclease activity against fungal protein in wheat [29]. Furthermore, a gene encoding a pathogenesis-related protein was reported as up-regulated in the interaction of M. incognita with Arabidopsis at 21 dai [57]. They also reported the increased expression of genes encoding several proteinase inhibitor proteins and leucine-rich repeat family proteins.