EST generation and sequence analysis

In order to determine sequence similarities of our contigs and in particular to identify genes that are conserved between different nematode species, we BLAST searched the 6,860 H. glycines contigs versus three databases (Figure 1). About half of the contigs (44%) matched sequences in at least one of these three databases at a threshold value of E = 1e^-20. Examination of the BLAST match distribution revealed that 19% of the contigs that matched all three databases are most likely representing highly conserved genes involved in fundamental housekeeping processes in metazoans, while the 31% of contigs exclusively matching sequences in the cyst nematode database contained genes that likely are important for specific host adaptations of Heterodera spp.

When assessing BLAST hit identities, the cluster that contained the most ESTs (HgAffx.13905.2; 599 ESTs) belonged to a gene coding for a putative cuticular collagen. Identities of other highly represented contigs were actin, tropomyosin and myosin, as well as additional house-keeping genes like ribosomal components, ubiquitin, arginine kinase, synaptobrevin and heat shock proteins. Interestingly, four putative parasitism gene sequences from the esophageal gland cells were among the 40 contigs with the highest EST constituents: three of unknown function (AAO33474.1, AAL78212.1, AAP30835.1) [7, 8] and a β-1,4-endoglucanase-2 precursor.

We further determined which H. glycines genes showed the highest degree of conservation when compared to C. elegans. BLASTX searches of all 6,860 soybean cyst nematode contigs against the Wormpep database revealed that 34.9% matched C. elegans entries at a threshold value of E = 1e^-20. The products of the 25 most conserved genes were heat shock proteins, proteins related to transcription and translation (for example, elongation and splicing factors and RNA polymerase II) and structural proteins, including tubulin and actin, as well as enzymes, including guanylate cyclase (Additional data file 3).

A survey and functional classification of developmentally regulated genes

In order to identify H. glycines genes that are developmentally regulated and to document their expression profiles, we designed a microarray experiment using three complete and independent biological replications (that is, three independent sample series representing three complete life cycles). We identified 6,695 probesets (Additional data file 4) as described in Materials and methods that were differentially expressed with a false discovery rate (FDR) of 5% when observed over the entire life cycle of H. glycines. This group of probesets equals 89% of all H. glycines probesets on the microarray. In other words, the vast majority of H. glycines genes represented on the GeneChip significantly changed expression during the nematode life cycle. We then grouped these 6,695 probesets into 10 clusters based on their expression profiles (Figure 2). As an exemplary gene family, we analyzed the expression pattern of FMRF (Phe-Met-Arg-Phe-NH₂)-related neuropeptide (FaRP)-encoding genes. This group encodes a specific class of neuronally expressed tetrapeptides that are potent myoactive transmitters in nematode neuromusculature [52–56], which are expressed in motor neurons that act on body wall muscle cells [57–59]. Based on our BLAST searches against various databases as detailed above, we identified five probesets for genes encoding FaRPs (HgAffx.23446.1.S1_at, HgAffx.23636.1.S1_at, HgAffx63.1.S1_at, HgAffx20469.1.S1_at, HgAffx.24161.1.S1_at). All five probesets were co-expressed with each other and showed an expression peak in the infective J2 stage (Additional data file 1). These FaRP probesets were differentially expressed when observed over the entire life cycle of the nematode, and, with the exception of HgAffx.20469.1.S1_at, which was found in cluster 4, all probesets were grouped in cluster 7 (Figure 2). The general profile of cluster 4 showed an expression peak in infective J2 and fell steadily in later life stages, while cluster 7 demonstrated the same overall pattern but showed a more pronounced increase from egg to infective J2.

We used InterProScan [60] to conduct functional classification for all 6,860 contigs in all expression clusters of each comparison. The relative abundance of the 25 InterPro domains with the highest representation in each of the 10 clusters for contigs that showed differential expression (FDR 5%) throughout the entire life cycle is summarized in Additional data file 5. While most clusters contained a wide range of genes represented by diverse InterPro domains, collagen domains stood out, in that they accumulated in cluster 2 at a high frequency relative to other InterPro domains. Since it has been suggested that down-regulation of collagens might be a common feature in dauer and infective stages of nematodes [27], we analyzed the expression profiles of H. glycines collagens in more detail. Using a reciprocal BLAST strategy as described in Materials and methods, we identified eight H. glycines probesets representing seven unique contigs orthologous to C. elegans collagens (Table 4). The temporal expression pattern of these seven orthologs was very similar (Figure 3) and congruent with observations in other nematode species [14, 51], which supports the hypothesis that down-regulation of collagen transcription is a conserved characteristic of non-molting infective and dauer-stage nematodes [27].

H. glycines probeset	C. elegans collagen	E-value, score, % identity (BLASTX)**	E-value, score, % identity (TBLASTN)**
HgAffx.10090.1.S1_at*	CE06699	2e-25, 105, 59%	3e-24, 105, 59%
	CE05938	2e-25, 105, 59%	2e-25, 109, 57%
	CE05937	2e-25, 105, 59%	2e-25, 109, 57%
HgAffx.10017.1.S1_at	CE05147	1e-25, 106, 51%	2e-40, 159, 39%
HgAffx.18987.1.S1_at	CE32085	3e-32, 127, 66%	1e-30, 127, 66%
HgAffx.19573.1.S1_at	CE02380	4e-26, 108, 65%	4e-27, 115, 62%
HgAffx.19987.1.S1_at	CE02380	9e-30, 119, 62%	2e-28, 119, 62%
HgAffx.241.1.S1_at	CE29723	3e-32, 127, 56%	3e-32, 132, 52%
HgAffx.241.1.A1_at	CE29723	3e-32, 127, 58%	3e-32, 132, 52%
HgAffx.7962.1.S1_at*	CE04335	5e-82, 293, 69%	2e-87, 318, 85%
	CE04334	5e-82, 293, 69%	2e-87, 318, 85%

*Probeset matched several C. elegans collagens equally well.

**BLASTX of H. glycines nucleotide probesets against C. elegans collagen proteins and TBLASTN of C. elegans collagen proteins against H. glycines nucleotide probesets.

Heterodera glycines orthologs of dauer-enriched C. elegansgenes are more likely to be down-regulated upon transition from infective J2 to parasitic J2 and J3 than other genes

In addition to providing a comprehensive gene characterization and expression resource, we wished to demonstrate the applicability and power of our data by addressing the question of whether the infective J2 stage of H. glycines is biochemically analogous to the C. elegans dauer larva stage. We compiled a list of 1,839 C. elegans genes that were identified by Wang and Kim [61] as so-called dauer-regulated genes by conducting microarray experiments comparing gene expression in C. elegans larvae that were in transition from dauer to non-dauer with that of freshly fed L1 larvae that had been starved. Dauer-regulated genes showed significant expression changes during a dauer exit time course that were not related to the introduction of food [61]. Reciprocal BLAST searches resulted in the identification of 438 H. glycines probesets that could be categorized unambiguously as orthologs of these C. elegans dauer-regulated genes (Additional data file 6). Because of the deliberate redundancy of the Affymetrix GeneChip, these 438 probesets corresponded to 396 unique H. glycines gene predictions. In other words, we identified H. glycines orthologs for 22% of the 1,839 C. elegans dauer-regulated genes (396/1,839).

To test whether the frequency of H. glycines orthologs to C. elegans dauer-regulated and dauer-enriched genes is similar to that of other, randomly chosen genes, we randomly selected 1,000 C. elegans proteins from the Wormpep database (v. 157) and repeated the reciprocal BLAST searches. In these searches, we identified 159 unique H. glycines contigs that fulfilled our criteria (data not shown). In other words, 16% of these randomly selected C. elegans genes have H. glycines orthologs. These analyses showed that C. elegans dauer-regulated genes have a slightly higher frequency (22%) of having orthologs in H. glycines than either dauer-enriched (14%) or random (16%) genes, both having about the same rate.

Following the identification of H. glycines orthologs for C. elegans dauer-regulated and dauer-enriched genes, we clustered these genes according to their expression profiles throughout the life cycle. Clustering the 438 probesets for the dauer-regulated orthologs led to their placement into nine groups (Additional data file 2), while clustering of the 74 H. glycines probesets for dauer-enriched gene orthologs resulted in seven distinct groups (Figure 4). It is obvious that not all H. glycines dauer-enriched orthologs were down-regulated from infective J2 to parasitic J2. Indeed, only 41% out of 74 probesets were significantly down-regulated, whereas 22% were up-regulated and 38% did not exhibit a statistically significant change in expression. Similarly, when comparing the infective J2 stage with the J3 stage of H. glycines, only 47% out of 74 probesets were down-regulated. Twenty-two percent were up-regulated, and 31% did not exhibit a statistically significant change in expression. In other words, in both comparisons, the majority of H. glycines genes that are orthologous to C. elegans genes down-regulated upon dauer exit were up-regulated or did not exhibit a statistically significant change in expression.

To determine whether H. glycines genes orthologous to C. elegans dauer-enriched genes behave differently from other H. glycines genes, we compared the dauer-enriched H. glycines orthologs with the entire set of 7,530 H. glycines probesets on the Affymetrix Soybean Genome Array, as well as to the 159 H. glycines genes that we determined to be orthologous to 1,000 randomly chosen C. elegans genes. We found that out of 7,530 H. glycines probesets, 19% were down-regulated when infective J2 are compared to parasitic J2, 21% were up-regulated and 60% did not exhibit a statistically significant change. Similarly, when infective J2 are compared to J3, 26% of all probesets were down-regulated, 20% were up-regulated and 54% did not exhibit a statistically significant change. The 159 H. glycines genes that are orthologous to 1,000 randomly chosen C. elegans genes are represented by 181 probesets on the Affymetrix GeneChip. Of those 181, 27% were down-regulated from infective J2 to parasitic J2, 24% were up-regulated and 50% did not exhibit a statistically significant change. When infective J2 were compared to J3, 26% out of 181 were down-regulated, 20% up-regulated and 54% did not exhibit a statistically significant change. We used a Fisher's exact test [62] to examine whether the proportion of down-regulated genes among the set of dauer-enriched H. glycines orthologs was significantly different from all other genes on the array or from the H. glycines probesets that were orthologous to the 1,000 randomly chosen C. elegans genes, respectively. We found that the observed differences between the proportions of down-regulated probesets between dauer-enriched H. glycines orthologs and the entire set of probesets on the microarray were significant at the 0.05 level in comparisons of both infective J2 versus parasitic J2 (P = 0.000015) and infective J2 versus J3 (P = 0.007160). Similarly, the differences between dauer-enriched H. glycines orthologs and the H. glycines probesets orthologous to random C. elegans genes were significant for comparisons of infective J2 versus parasitic J2 (P = 0.0219190) and for infective J2 versus J3 (P = 0.0094261). If a Bonferroni correction is used to control the overall type I error rate for this family of four tests, all comparisons would remain significant at the 0.05 level except the comparison between dauer-enriched H. glycines orthologs and the H. glycines probesets orthologous to random C. elegans genes for infective J2 versus parasitic J2. In other words, while the majority of H. glycines genes that are orthologous to C. elegans dauer-enriched genes was not down-regulated upon transition to parasitic J2 or J3, the proportion of H. glycines orthologs that were in fact down-regulated was statistically significantly enriched about two times compared to all H. glycines genes on the microarray or to orthologs to random C. elegans genes.

The identities of H. glycines genes that followed the expression pattern of their dauer-enriched C. elegans orthologs (that is, they were down-regulated upon transition to infective J2 or J3) reflect a wide range of effector functions and biochemical pathways, including peptidases, epoxide and glycoside hydrolases, phosphate transporters and neuropeptide-like proteins. H. glycines genes that did not follow the C. elegans pattern of expression (that is, they were not down-regulated) span an equally diverse group of genes and include carbohydrate kinase, catalase and glutathione peroxidase (Table 5).

Metabolism in C. elegans dauer larvae and H. glycinesinfective J2 is dissimilar

Validation of microarray results by quantitative RT-PCR

Quantitative real-time reverse transcription PCR (qRT-PCR) was used to validate selected microarray results. We analyzed the expression patterns of six genes representing different expression patterns for each of the five consecutive pairs of life stages (egg/infective J2, infective J2/parasitic J2, parasitic J2/J3, J3/J4, J4/female), giving a total of 30 different genes. The qRT-PCR template was the same biological material used for hybridization of the microarrays, and the reactions were performed in technical triplicates. Additional data file 7 shows that there is qualitative agreement between our microarray approach and qRT-PCR for 28 out of 30 genes.

Nematode cultivation, cDNA library generation, sequencing and clustering

H. glycines strain OP-50 [67] was cultivated under greenhouse conditions and isolated as described previously [68]. Unidirectional Uni-Zap lambda libraries (Stratagene, La Jolla, CA, USA) were generated for H. glycines strain OP-50 eggs, infective J2, J3, J4 and virgin females. The mRNA concentration used ranged between 1.0 μg (J3) and 4.6 μg (eggs). The libraries were sequenced at the Washington University Medical School Genome Sequencing Center (St Louis, MO, USA), and the sequencing results were deposited in GenBank. Contigs were formed by Affymetrix for the design of the H. glycines group of probesets of the Affymetrix Soybean Genome Array GeneChip and all consensus sequences and contig size details can be accessed at Affymetrix [69].

Nematode cultivation for microarray experiments

For each replication, 40 pots with 10 seeds each of Kenwood 94 soybeans were planted in a 2:1 sand:soil mixture in the greenhouse. Two weeks after planting, each pot, containing an average of 7 to 8 germinated seedlings, was inoculated with 15,000 to 20,000 H. glycines strain OP-50 [67] infective J2. The inoculum was collected by setting up two hatch chambers, each containing about two million H. glycines OP-50 eggs, and allowing the eggs to hatch for four days. From the same batch of eggs used in the hatch chamber, 50,000 eggs were collected and flash frozen in liquid nitrogen for use as the egg stage of the replication. After 4 days, the hatched infective J2 were collected and counted, and an aliquot of 50,000 larvae was flash frozen in liquid nitrogen for use as the infective J2 stage of the replication. The remainder of hatched infective J2 larvae was divided up among the 40 pots for seedling inoculation. Four days after infection, 12 pots were collected, and the soil was washed away from the root systems of these pots to isolate the parasitic J2 stage. Eight days after infection, another 12 pots were harvested for collection of J3 juveniles, and, 14 days after infection, a further 10 pots were used to isolate J4 juveniles. Finally, 21 days after infection, the final 6 pots were harvested for collection of adult females. These stages were isolated as published previously [68]. All stages were divided in about 30 mg aliquots in 1.5 ml screw cap tubes, flash frozen in liquid nitrogen and stored at -80°C.

RNA extraction for microarray experiments and GeneChip procedures

Frozen (-20°C) 1.0 mm zirconia beads (BioSpec Products, Bartlesville, OK, USA) were added to frozen nematode tissue in a 1.5 ml screw cap tube in a 1:1 tissue:bead ratio. RNA was isolated using the Versagene kit from Gentra Systems (Minneapolis, MN, USA). On ice, 800 μl lysis buffer was mixed with 8 μl tri 2-carboxyethyl phosphine (TCEP), and 50 μl was added to one screw cap tube with about 30 mg nematode tissue. Tissue was homogenized in a beadbeater (BioSpec Products) at setting 48 for 10 s and chilled on ice for 60 s. This step was repeated three times. Two more tubes for the same life stage were treated in the same way. The suspension was collected by pipetting and transferred into a new tube. Beads were rinsed with remaining buffer, and the standard Versagene protocol was followed from then on. For all stages and all repetitions, we obtained 113 to 320 mg frozen tissue and 7.85 to 173.38 μg total RNA. The RNA concentration and quality of each sample were determined by a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and by RNA Nanochip on a 2100 Bioanalyzer (Agilent Technologies Inc., Palo Alto, CA, USA). RNA was submitted to the Iowa State University GeneChip Facility, where standard procedures recommended by Affymetrix were followed for reverse transcription and labeling of the probes and for hybridization and scanning of the GeneChips.

Experimental design of microarray experiments and GeneChip data analysis

Expression was measured using a total of 18 Affymetrix Soybean Genome Array GeneChips (three replications × six life stages) using a randomized complete block design with replications as blocks. Prior to performing the analysis, the Affymetrix signal data were transformed using the natural log (ln) and normalized by median centering (that is, the median ln signal from each particular chip was subtracted from all ln signals on the chip). The normalized data for each gene were analyzed separately using a standard linear model with fixed effects for replications and stages. F tests, resulting in p-values, were performed to test for a difference in expression between the life stages for each probeset. A q-value was computed for each p-value using the method described by Storey and Tibshirani [70]. These q-values can be used to identify differential expression while maintaining approximate control of the FDR. For example, FDR is controlled at approximately 5% if the sets of tests with q-values at or below 0.05 are declared significant. See [71] for a discussion of linear modeling and FDR control in the context of plant microarray experimentation.

Clustering was used to organize and summarize the observed expression patterns of differentially expressed genes. For each probeset, the mean normalized expression level was estimated for all six life stages. The six estimated values were standardized to have mean 0 and standard deviation 1 within each probeset. The Euclidian distance between any pair of standardized expression profiles was used as a measure of dissimilarity in all clustering algorithms. This approach considers genes with similar expression patterns to be close in six-dimensional space and is equivalent to using (1-r)^0.5 as the measure of dissimilarity, where r is the Pearson correlation coefficient between non-standardized expression profiles.

K-medoids clustering [72] was performed on those probesets with the 1,000 lowest p-values (< 0.0000143) for the overall test for expression change across the life cycle. Using the Gap statistic [73], we determined the number of clusters to be 10. Probesets with q-values less than 0.05 were added to these 10 base clusters to produce the clustering depicted in Figure 2. All other clusters were produced using hierarchical agglomerative clustering, using average linkage to measure the dissimilarity between clusters.

All clusters and related figures were generated using the free open-source statistical software package R [74]. Hierarchical clustering was carried out using the R function hclust; K-medoids clustering was carried out using the function pam from the R cluster library.

BLAST searches of H. glycinescontigs

We built three nematode-specific nucleotide databases as follows: cyst nematodes (13,643 sequences) - all sequences from the GenBank query "Globodera [ORGN] or Heterodera [ORGN] not Heterodera glycines [ORGN]"; non-cyst nematodes (933,882 sequences) - results of the GenBank query "Nematoda [ORGN] not Globodera [ORGN] not Heterodera [ORGN]"; and non-nematodes (3,607,410 sequences) - GenBank 'nt' database (nucleotide version of 'nr') minus all Nematoda [ORGN] sequences. The BLASTN parameters were set to expect a 75% target frequency as follows: M = 1, N = -1, Q = 3, R = 3, B = 10, V = 10, lcmask, golmax = 10, topcomboN = 1, filter = seg. The following parameters were used for BLASTX searches of all 6,860 contigs against non-redundant GenBank (downloaded 29 November 2005) and Wormpep v. 152: filter = seg, lcfilter, W = 4, T = 20, E = 100, B = 25, V = 25, topcomboN = 1, golmax = 10.

InterProScan

InterProScan was run using InterPro data files [75] dated November 2005 (iprscan_PTHR_DATA_12.0.tar). InterProScan translated all 6,860 unique contig sequences in six frames and then ran its suite of domain-finding tools. We required a minimum translation length of 20 amino acids to be considered by InterProScan, and we used the EGC.0 translation table. Due to the six-frame translation, each contig typically had several alignments amongst the significant open reading frame (ORF) found in the translation. We kept, as representative of each contig, the single longest aligning ORF that contained an InterPro domain, even though InterProScan may have found several ORFs for each contig with alignments to some domain or motif. Results were parsed into files representing expression clusters.

Identification of C. elegansdauer and collagen orthologs

A previous microarray study identified 1,984 so-called dauer-regulated and 540 dauer-enriched genes in C. elegans [61]. We manually compiled a list of 1,839 dauer-regulated and 488 dauer-enriched C. elegans genes using these data (the remaining genes were either deleted or merged with other genes in Wormbase). To isolate H. glycines orthologs, we conducted BLAST searches between these C. elegans genes and all 7,530 H. glycines probeset nucleotide sequences (threshold values 1^e-15, bit score at least 50, identity at least 30%). We then conducted an additional BLAST search between the H. glycines probesets that passed those criteria and the entire Wormpep (v. 159) database using the same cutoff criteria as above. Only those ortholog pairs in which the C. elegans dauer-regulated or dauer-enriched hit was either the best one compared to C. elegans hits of the entire Wormpep database or within 5% of the best hit's bit score were kept. These searches generated a list of 438 H. glycines probesets orthologous to C. elegans dauer-regulated genes and 74 H. glycines probesets orthologous to C. elegans dauer-enriched genes.

To identify collagen orthologs, we downloaded the sequences for collagens listed at the Sanger Institute web site [76] and followed basically the same BLAST strategy, with the exception that we required a higher stringency with a bit score of at least 100.

qRT-PCR

The transcript abundance of 30 differentially expressed cDNA clones was analyzed by qRT-PCR to confirm microarray results. Gene-specific primers were designed, and the sequences are shown in Additional data file 7. DNase-treated RNA (10 ng) was used for cDNA synthesis and PCR amplification using an iScript One-Step RT-PCR kit (BIO-RAD, Hercules, CA, USA) according to manufacturer's protocol. The PCR reactions were performed using an iCycler (BIO-RAD) under the following conditions: 50°C for 10 minutes, 95°C for 5 minutes and 40 cycles of 95°C for 30 s and 60°C for 30 s. Following PCR amplification, the reactions were subjected to temperature ramp to create the dissociation curve, measured as changes in fluorescence as a function of temperature, by which the non-specific products can be detected. The dissociation program was 95°C for 1 minute, 55°C for 10 s, followed by a slow ramp from 55°C to 95°C. Three replicates of each reaction were performed, and constitutively expressed Actin 1 (AF318603) was used as internal control to normalize gene expression levels. Quantifying the relative changes in gene expression was performed using the 2^-ΔΔCTmethod as described by Livak and Schmittgen [77].

Data

The Affymetrix Soybean Genome Array GeneChip raw and normalized data files were deposited in the ArrayExpress database [78] under accession number E-MEXP-1110. This database is MIAME-compliant.

GenBank accession numbers for EST generated in this study: BF013452-BF014867, BF249436-BF249525, BG310659-BG310919, BI396450-BI397002, BI451481-BI451765, BI704104-BI704170, BI749028-BI749665, BI773548-BI773559, CA939105-CA940989, CB238504-CB238733, CB238735-CB238737, CB238740, CB238743-CB238746, CB238750-CB238751, CB238755-CB238759, CB238763, CB238766-CB238769, CB238772-CB238774, CB238776-CB238778, CB238780-CB238782, CB238784, CB238788, CB238790-CB238793, CB278389-CB280465, CB281104-CB281884, CB299031-CB299936, CB373779-CB376363, CB377726-CB380382, CB824093-CB826589, CB934836-CB935636, CD747803-CD749135.