- Research article
- Open Access
Infection and genotype remodel the entire soybean transcriptome
BMC Genomics volume 10, Article number: 49 (2009)
High throughput methods, such as high density oligonucleotide microarray measurements of mRNA levels, are popular and critical to genome scale analysis and systems biology. However understanding the results of these analyses and in particular understanding the very wide range of levels of transcriptional changes observed is still a significant challenge. Many researchers still use an arbitrary cut off such as two-fold in order to identify changes that may be biologically significant. We have used a very large-scale microarray experiment involving 72 biological replicates to analyze the response of soybean plants to infection by the pathogen Phytophthora sojae and to analyze transcriptional modulation as a result of genotypic variation.
With the unprecedented level of statistical sensitivity provided by the high degree of replication, we show unambiguously that almost the entire plant genome (97 to 99% of all detectable genes) undergoes transcriptional modulation in response to infection and genetic variation. The majority of the transcriptional differences are less than two-fold in magnitude. We show that low amplitude modulation of gene expression (less than two-fold changes) is highly statistically significant and consistent across biological replicates, even for modulations of less than 20%. Our results are consistent through two different normalization methods and two different statistical analysis procedures.
Our findings demonstrate that the entire plant genome undergoes transcriptional modulation in response to infection and genetic variation. The pervasive low-magnitude remodeling of the transcriptome may be an integral component of physiological adaptation in soybean, and in all eukaryotes.
How many genes are truly involved in the response of organism to a challenge such as pathogen infection, and what are the roles of those genes? Global assays of gene expression, for example by microarray analysis, are typically conducted to test the hypothesis that a small, defined set of genes are responsible for an organism's response to some challenge. Gene expression changes below a certain threshold (commonly 2 fold) are often disregarded as being irrelevant and/or unreliable. A major challenge in evaluating the importance of low magnitude transcriptional changes is that the level of replication used in a typical microarray experiment is insufficient to detect small changes given the technical and biological variability in the system. Although several methods appear to be promising for precise quantification of gene expression, it remains uncertain what constitutes a significant change in response to treatments [1, 2].
High-density oligonucleotide arrays such as Affymetrix GeneChips can detect up to 90% of all the mRNAs in a transcriptome [3–5]. For example, nearly 90% of all yeast mRNAs could be detected in cells grown under both rich and minimal media growth conditions, with approximately 50% being present at average levels between 0.1 and 1 copy per cell . Of the 31,000 genes on Affymetrix Rat Genomic 230 2.0 GeneChip microarrays, 18,200 (58.7%) could be detected in growing rat bone . In a study with human abdominal aortic aneurysms, of the 18,057 genes common to Affymetrix and Illumina arrays, 11,542 (64%) were expressed in either aneurysmal or normal abdominal aorta . Approximately 26,500 of the soybean genes (70%) on the Affymetrix GeneChip could be detected in soybean cyst nematode (SCN)-colonized root pieces.
Markedly varied numbers of genes, from only a few up to several thousands, have been reported to be differentially expressed in response to diverse challenges, depending on the system and the statistical methodology. For instance, of the approximately 6,200 protein-encoding genes in the Saccharomyces cerevisiae (yeast) genome, over 1,000 showed significant changes in mRNA levels during sporulation . In rat, 8,002 out of 18,200 expressed genes (44.0%) had a significant change in gene expression during growth, about half up-regulated and half down-regulated . In Arabidopsis thaliana, 939 out of approximately 24,000 genes showed a statistically significant response to cold stress, with 655 up-regulated and 284 down-regulated .
One of the most profound challenges an organism can suffer is pathogen infection. In a meta-analysis of 32 studies involving 785 transcriptomic experiments with 77 different host-pathogen interactions , 5042 human genes showed transcriptional changes in response to at least one challenge, and a cluster of 511 co-regulated genes was identified as representing a common infection response. During infection of the plant Arabidopsis by the bacterial pathogen Pseudomonas syringae, approximately 2,000 of the approximately 8,000 genes monitored showed significant expression level changes . In soybean, the Affymetrix GeneChip has been used to profile gene expression during infection with soybean rust fungus and soybean cyst nematode (SCN) [4, 11–14]. During nematode infection, 429 of 35611 soybean transcripts (which account for 1.2%), while 1850 out of 7430 SCN genes (24.9%) showed expression changes .
To identify genes involved in the responses of several soybean genotypes to infection by the oomycete pathogen Phytophthora sojae, we conducted a very large-scale microarray experiment using Affymetrix GeneChips. Three soybean genotypes (V71-370, Sloan and VPRIL9) were included within each of the 29 experimental blocks. Replicates of each set of the three genotypes, incubated in the same growth chamber, were harvested at three different times (9 am, 10:30 am, and 12 pm). For each soybean line, approximately 30 seedlings were inoculated on the roots with P. sojae and after 5 days, 7.5 mm root segments were collected from immediately above and below the upper margin of the visible lesion (referred to from here as the "Upper" and "Lower" infection courts). From approximately 20 mock-inoculated plants, 15 mm root segments ("Mock") were collected from positions on the roots matching the lesions on the pathogen-inoculated plants. RNA was extracted from each of the samples, labeled and hybridized to Affymetrix soybean GeneChips. This paper reports our analyses of a total of 648 GeneChips with the samples obtained from 24 blocks from which a full set of 27 successful microarray assays were obtained.
The entire soybean transcriptome is significantly remodeled by genotypic differences and Phytophthora infection
The primary GeneChip data were pre-processed using GC-RMA background correction, quantile normalization and median polish summarization. To identify factors that significantly affected expression of individual genes, Linear Mixed-Model Analysis (LMMA) was used, and three methods to control the False Discovery Rate (FDR) were employed.
LMMA indicated that 89% of the detectable genes (i.e. those with a significant hybridization signal determined by the MAS5 presence-absence test) were significantly affected by infection in both the upper and lower infection courts, even when using the most conservative FDR method (BH-FDR) at a level of 0.0001 (Table 1). With the least conservative method (q-value) at a level of 0.05, 98.3% of the detectable genes were identified as significantly affected. The analysis also identified very large numbers of genes whose transcription level was significantly affected by genotype. For instance, at moderate stringency (TST-adjusted p ≤ 0.01), 25,923 genes (91.4% of all detectable genes) showed significant expression differences among the three soybean lines. Furthermore, 24,669 (87.0%) showed a genotype-sensitive response to infection (significant genotype × infection interaction; TST-adjusted p ≤ 0.01). Genes showing significant responses to infection and/or infection × genotype interaction constituted 98.2% to 99.4% of the detectable transcriptome (27,831 or 28,187 genes with TST-adjusted p ≤ 0.01 or 0.05 respectively). Genes showing significant responses to genotype and/or infection × genotype interaction constituted 97.7% to 99.4% of the detectable transcriptome (27,710 or 28,172 genes with TST-adjusted p ≤ 0.01 or 0.05 respectively). The experiment was also able to identify a large number of genes that responded significantly to the subtle differences in the times of harvest (9 am, 10:30 am, 12 noon) (4,453 genes = 15.7% at TST-adjusted p ≤ 0.01), although as expected this number was smaller than the overall response to infection and genotype.
To verify that our statistical analysis was effectively controlling false discovery, we separated the 24 blocks into two "sub-groups" consisting of odd-numbered and even-numbered blocks respectively. We then included the effect of sub-group in our linear mixed model and used LMMA to identify genes showing significant differences between the sub-groups. Very few genes showed differences under even the most relaxed criteria (45 at TST-FDR adjusted p ≤ 0.05; Table 2), while at moderate criteria (TST-FDR adjusted p ≤ 0.01) none showed differences. Similarly, when the 24 blocks were randomized into 6 sub-groups, each containing 4 blocks, only 20 genes showed significant differences among sub-groups (TST-FDR adjusted p ≤ 0.01). These small numbers of genes that show differences among sub-groups may reflect genes that are unusually sensitive to small variations in conditions used to assay the different biological replicates.
To verify that our results did not depend on the preprocessing methods (background subtraction and data normalization), we preprocessed the GeneChip data with the Bioconductor MAS5 algorithm and re-analyzed the results by LMMA. For the fixed factors, genotype, treatment, and genotype × treatment interaction, the agreement between the results of the two methods was between 92% and 98% at all levels of significance (Table 3). For the time factor the agreement was less, but still at least 77.0%.
To independently verify the results of the LMMA analysis, we employed the non-parametric Wilcoxon signed rank test to analyze the significance of infection responses; data from each of the infection chips were paired with the corresponding data from mock inoculation. The genes identified as significant by this method matched those identified by LMMA with an agreement of > 99%, with either method of pre-processing (GC-RMA or MAS5) (Table 4).
Transcriptional changes display both consistency and fluctuation
To assess to what extent the transcriptional changes we observed were biologically reproducible, we compared our results to those obtained from a pilot experiment carried out two years previously with cultivars V71-370 and Sloan. That pilot experiment had a near identical design except that all plants were harvested at one time during the morning. Also, the upper and lower infection courts were harvested together. To enable comparison with the 72 sample experiment, in which the upper and lower infection courts were harvested separately, the results from the upper and lower courts were averaged. The overall correlation in the expression changes measured for all genes by microarray analysis in the two experiments was very high (R2 = 0.85) indicating a high level of reproducibility. Even genes showing less than two-fold perturbations by infection showed excellent consistency between the two experiments. For example, in V71-370, Gma.3473.1.S1_at (small heat shock protein) was down-regulated 1.2 fold and 1.7 fold in the two experiments, and in Sloan down-regulated 1.7 fold and 2.2 fold. As another example, Gma.6640.1.S1_at (secretory peroxidase) was up-regulated 1.46 fold and 1.07 fold in V71-370, and down-regulated 1.08 fold and 1.16 fold in Sloan. In a third example, GmaAffx.60616.1.S1_at (Hcr9-OR2A) was down-regulated 1.5 fold and 1.35 fold in V71-370, and down-regulated 1.47 fold and 1.8 fold in Sloan. Only one measurement showed a statistically significant difference between the two experiments (GmaAffx.18905.1.S1 cell death associated protein: in Sloan; FDR adjusted p < 0.05). This agreement is remarkable given that this is a two-organism interaction and that two experiments were conducted two years apart with different batches of soybean seed.
During the pilot experiment, transcript levels of 22 genes, including 4 housekeeping genes, were also measured by quantitative real-time PCR (qRT-PCR), using the same RNA preparations as for the microarrays. For the 22 selected genes, the correlation between the microarray measurements from the two experiments was 0.90 (R2) and the correlation between the microarray and qRT-PCR measurements was 0.88 (R2) (Table 5) [see also Additional file 1]. In most cases, the agreement between the qRT-PCR results and the microarray results was very good. In only two cases (GmaAffx.64261.1.A1_at quinone oxidoreductase and GmaAffx.75675.1.A1_at phenylalanine ammonia-lyase, both in V71-370) were statistically significant differences (FDR-adjusted p < 0.05) observed between the qRT-PCR and the array results; possibly in these cases the two assay methods differed in their sensitivity to the presence of mRNAs from close paralogs and homeologs of the assayed genes.
To further explore the biological reproducibility of the measured gene expression differences, we divided the 24 blocks into 6 sub-groups, each consisting of 4 consecutive blocks (12 measurements) and compared the results for each of the sub-groups for the same set of 22 genes (Table 5). Since the 24 blocks each constitute independent biological replicates spread over 8 months, each of the 6 sub-groups also constitute independent biological replicates. In most cases, only minor fluctuations were observed. Two genes showed significant fluctuations (Gma.17993.1.S1_at P21 protein and GmaAffx.32612.1.S1_at HI4'OMT; least significant difference, FDR adjusted p < 0.05).
Comparing the power of large and small experiments
To compare the power of our experiment with experiments of a size comparable to those more commonly found in the literature, we again divided the 24 blocks into 6 sub-groups, each consisting of 4 consecutive blocks, and then reanalyzed the data by LMMA for each sub-group separately. As expected, each individual data set had less power than the overall data set. An average of 40.4%, 65.7% and 33.8% genes were detected as significant for genotype, treatment and genotype × treatment interaction, respectively, with TST-FDR control at a level of 0.01 [see Additional file 2]. However, different sub-experiments detected different subsets of the genes detected in the overall data set (Figure 1). Thus genes significant in all the 6 sub-groups comprised only 48.9% of those detected in the overall experiment [see Additional file 3]. This lack of agreement among the gene lists became worse at more stringent (lower) FDR levels. As is well-known, increased stringency produces more precise results (fewer false positives) but at the expense of a sometimes major decrease in power, i.e. an increase in the number of false negatives. However, the union of all genes that were significant in any of the six sub-groups comprised 81.0% of the genes with significant changes detected in the overall experiment [see Additional file 3]. These results demonstrate that the most common error made with small experiments is false negatives, rather than false positives, and the common practice of combining results of multiple experiments by only considering the intersection of the gene lists is unnecessarily extremely conservative. A simple union, the techniques of meta-analysis, or where possible a joint analysis, may be most appropriate for combining data from multiple independent experiments.
Low-magnitude transcriptomic remodeling differs among functional categories
Many researchers use two-fold change as an arbitrary cutoff. Our results show that 13.1 to 23.5% of the significantly changed genes in the upper infection court and 29.8% to 36.8% in the lower infection court changed by two-fold or greater in the three tested genotypes (Figure 2) [see Additional file 4]. However, 8.2 to 15.0% of the genes with significantly changes showed only subtle changes (fold change of less than 1.2-fold). To examine the functional significance of gene expression changes of various magnitudes, we plotted the distribution of changes for genes in six functional categories relevant to infection, namely defense and disease, signal transduction, transcription, intracellular trafficking, cell structure and metabolism (Figure 3) [see Additional files 5 to 7]. In the lower infection court of the susceptible cultivar, Sloan, which had the greatest level of infection, 23,556 genes showed significant changes in response to infection. Compared to the entire set of modulated genes, the distribution of changes was significantly different for all six categories (p ≤ 0.01, Kolmogorov-Smirnov test for two samples) (Figure 3). The disease and defense category and to a lesser extent metabolism showed a bias towards strongly up-regulated genes, whereas transcription and to a lesser extent signal transduction and intracellular trafficking, showed a bias towards down-regulation of genes. Similar patterns were observed in the upper infection court of Sloan and in both courts of the resistant cultivar V71-370 [see Additional files 5 to 8]. To test if significant changes in distribution were present among genes showing low magnitude modulation, the distribution comparison was restricted to genes showing less than two-fold or less than 1.5-fold modulation. Five of the categories (all except disease and defense), and three categories (transcription, cell structure and metabolism), respectively, showed significantly different distributions from the overall gene set (p ≤ 0.01). Even among genes with less than a 1.2-fold change, those in the transcription and metabolism categories still showed significantly different distributions (p ≤ 0.01) (Table 6) [see Additional file 8].
Infection profoundly affects the physiology of host cells, including the levels of mRNAs within each cell. The ability to measure changes in host mRNA levels during infection is greatly limited by the intrinsic variability of the tissue and the difficulty of accurately reproducing infection across multiple replicates. As result, typical studies report significant changes in only a few percentage of the genes being assayed [3–5, 7, 8, 10]. Similar difficulties limit the power of many other experiments that measure transcriptional responses to genetic and physiological perturbations. Here we have shown, by means of an experiment of unprecedented statistical power, that significant changes in the levels of nearly all host transcripts can be measured during infection of soybean with the eukaryotic pathogen Phytophthora sojae. Up to 36.8% of the changes were greater than 2-fold in magnitude but the majority were of lower magnitude. Up to 15% of the changes were less than 1.2-fold in magnitude. Similar widespread reprogramming of the transcriptome was also detected in response to genotypic differences among soybean cultivars, and even in response to different times of the day at which plants were harvested.
An experiment such as this, combining 72 biological replicates, each comprised of 20–30 plants, is very successful at identifying overall trends, in this case identifying overall responses to pathogen infection. In many of the individual genes we examined, the overall measurement also was a very good representation of a consistent response to infection observed across many independent experiments. On the other hand, some genes showed significant fluctuations from experiment to experiment. A major challenge in the future will be to design experiments that can dissect these changes and ask whether the fluctuations simply result from technical variations in the assay, whether the fluctuations result from small uncontrolled variations in the conditions in which the biological material is grown, or the more interesting possibility that fluctuations in mRNA are a normal physiological event. It is possible for example that organisms have evolved the ability to tolerate significant changes in the levels of mRNAs or other cellular components, and there is no evolutionary advantage to imposing extremely ranges on mRNA levels.
By re-analyzing subsets of our data that represented the scale of replication more commonly found in microarray experiments (4 replicates), we showed that approximately half of the transcriptional changes that we detected could not be observed with the small-scale experiments. Taking the intersection of six independent subsets greatly reduced the number of detected changes further, showing that this common practice for combining the results of multiple experiments is extremely and unnecessarily conservative. On the other hand, taking the union of multiple gene lists more closely approximated the results from the full-scale experiment.
A recent expression QTL study in Saccharomyces cerevisiae showed variations in the levels of 79%, 69% and 47% of the detectable transcripts in response to treatment, genotype and treatment × genotype interactions respectively . Those findings, together with ours, suggest that widespread, low-magnitude transcriptional remodeling may be a normal process during physiological adaptation in eukaryotes, but one that is missed by conventional experimental designs. The extensive treatment × genotype interactions observed in both these studies suggest that the transcriptional reprogramming is genetically controlled.
The overall physiological response of an organism or cell to a stimulus may require coordinated changes in a wide array of cellular components. Those changes in turn may require compensating or reinforcing changes in an even wider array of functionally-connected cellular components. As a result, understanding how a specific set of transcriptional changes relates causally to a physiological change at a whole system level, is a major challenge. Our analysis suggests that low magnitude expression changes may be of functional significance. The observation that the patterns of low-magnitude transcriptional remodeling were significantly different among most functional categories, in some cases even among genes with less than 20% perturbation, is consistent with the hypothesis that low amplitude remodeling has functional significance. However, we cannot rule out that some percentage of genes may show low magnitude transcriptional modulation that has no functional significance, that is, they represent a low uncontrollable level of transcriptional variation that soybean (and other organisms) may have evolved to tolerate. Interestingly, most of the transcriptional remodeling of defense and disease response genes was of high magnitude, while low-magnitude remodeling was widespread in the other functional categories. We speculate that strong transcriptional modulation of disease and defense-related genes is required for the host to directly engage the pathogen, while the numerous other genes that function in other categories are coordinately modulated to support or adjust to the direct response.
Biological gene regulatory networks are highly interconnected systems. Non-linear, synergistic interactions are common. Large numbers of genes with low magnitude transcriptional modulation could potentially be just as important in conferring phenotypes and mediating physiological adaptation as the small numbers of genes that show large magnitude modulations. However, understanding the role of pervasive low magnitude remodeling may require using computational modeling approaches at a systems level, as well as improved technologies for accurately and cheaply measuring those changes. Systems approaches will also be needed to develop a deeper understanding of how consistent small magnitude changes and stochastic fluctuations are integrated to produce phenotypes.
Our findings, consistent through two different normalization methods and two different statistical analysis procedures, suggest that almost the entire soybean genome undergoes significant transcriptional remodeling in response to pathogen infection and genetic variation. The majority of the transcriptional differences are less than two-fold in magnitude. The low amplitude modulation of gene expression (less than two-fold changes) is highly statistically significant, even for modulations of less than 20%. These findings demonstrate that low amplitude transcriptional modulation forms an integral component of physiological adaptation in soybean, and, we speculate, in all eukaryotes.
Plant materials and pathogen isolate
Three soybean genotypes with varying degrees of resistance to Phytophthora sojae were used for this study. V71-370 is a resistant genotype, VPRIL9 is moderately susceptible and Sloan and is highly susceptible[16, 17]. P. sojae isolate, PT2004C2.S1, isolated from Ohio in 2004, and virulent to soybean differentials with Rps1a, Rps1b, Rps1k, Rps2, Rps3a, Rps3c, Rps4, Rps5, Rps6, or Rps7, was used in this study.
Inoculation assays were performed using the slant board technique . Briefly, 7-day-old soybean seedlings, grown in vermiculite, were thoroughly rinsed under tap water and placed in inoculation trays (20–30 plants per sample). The plants were wounded at 2 cm below the beginning of the root zone by scraping the epidermis with a scalpel and then were inoculated with a mycelial agar slurry from a 7-day-old culture or agar alone. Root tissue samples 7.5 mm long were collected at 5 days post inoculation (dpi) from immediately below (the treatment "Lower") and above (the treatment "Upper") the upper lesion margin from each seedling with characteristic lesions. In each case, the root tissue samples from like treatments and harvests were pooled. Thus each set of 30 pathogen-inoculated plants yielded two pooled samples, one pooled from 30 sections from the lower infection courts and one pooled from 30 sections from the upper infection courts. For the mock-inoculated plants (the treatment "Mock"), 15 mm tissue sections were taken spanning from 7.5 mm below to 7.5 mm above the position corresponding to the average lesion length measured from the inoculated samples. Thus each set of 20 mock-inoculated plants yielded one pooled mock sample. All the plants for each block were grown in the same growth chamber (Environmental Growth Chambers, Chagrin Falls, Ohio; Model M-48 with TC2 microcontroller unit) with day and night temperatures settings of 27°C and 21°C, respectively and relative humidity averaging 75 to 90%.
Overall experimental design
The data for this study comes from control inoculations that formed part of a large study of a recombinant inbred line population that was organized into 29 blocks. All plants for a block were raised, inoculated and incubated together in the same growth chamber. Each block included the same three control genotypes (V71-370, Sloan and VPRIL9). Three replicates of the three genotypes were inoculated with P. sojae (generally 30 plants) or mock-inoculated (generally 20 plants), then harvested 5 days later at three different times during the morning of harvest (9 am, 10:30 am, and 12 pm). Each block yielded 27 sets of samples, which were the complete combinations of the 3 genotypes, 3 treatments and 3 time windows. Thus the experiment was balanced with respect to the factors, soybean genotypes, treatment ("Upper", "Lower", and "Mock") and inoculation time window. Since the differences between the samples harvested at 9 am, 10:30 am and 12 pm were small, in this study we refer to the 3 × 29 sets of samples obtained from the 29 blocks as biological replicates. The 24 blocks in which all the control samples passed quality control were used for statistical analyses.
RNA extraction and microarray assay
Total RNA was isolated from each of the pooled plant samples using the QIAGEN RNeasy® Plant Mini Kit, following the manufacturer's protocol for total RNA extraction from plant and filamentous fungal tissue with the following minor modifications. 250 mg of fine tissue powder was used with appropriate scaled-up buffer volumes, and the Buffer RLT-suspended samples were thoroughly vortexed for 2 minutes and then incubated in 56°C water-bath for 3 min and kept in ice for 1 min before the first spin. The quantity and quality of total RNA were checked with both a NanoDrop ND-1000 Spectrophotometer and an Agilent 2100 Bioanalyzer.
Microarray procedures were performed at the Core Laboratory Facility (CLF) of the Virginia Bioinformatics Institute. The standard eukaryotic gene expression assay protocols were followed as described in the Affymetrix GeneChip® Expression Analysis Technical Manual. Briefly, as the first step, biotin-labeled cRNA was generated using the One-Cycle Target Labeling and Control Reagents (Affymetrix) and 1 μg of total RNA. Second, 20 μg of biotin-labeled cRNA was fragmented in Fragmentation Buffer. Third, fragmented cRNA was hybridized to a soybean GeneChip at 45°C for 16 h in an Affymetrix hybridization oven (model 640). Fourth, the GeneChips were washed and stained with streptavidin-phycoerythrin using the fluidics protocol EukGE-WS2v5-450 in the Affymetrix® GeneChip® Fluidics Station 450. Finally, the stained chips were scanned with an Affymetrix GCS3000 7G Scanner. The Affymetrix GeneChip® Operating Software (GCOS, v1.4.0.036) was used to provide instrument control, first-level data analysis, and data management for the entire GeneChip System. Samples within experimental blocks were randomized before being provided to the CLF, to reduce possible block effects arising from processing samples together during the microarray analysis work flow.
Low level data pre-processing
Quality control was performed using a variety of tools from the MAS5 and AffyPLM quality control toolsets. These included analysis of percentage of probes called present, the scaling factor and average background of each chip. In addition, ratios of 3' to 5' probes for β-Actin were used to check for RNA degradation. Calculation of the interquantile ranges of the normalized unscaled standard error and of the relative log expression was used to remove outliers. All chips in the 24 blocks selected for this study passed these QC tests. Low-level analysis of the raw GeneChip data began with filtering the non-expressed genes based on the Affymetrix Microarray Suite version 5 (MAS5) algorithm for calls (as implemented in the Bioconductor package affy, http://bioconductor.org/packages/2.0/bioc/html/affy.html) taking into account the experimental design. The default parameter τ (0.015) in the MAS5 present call was used, and a probe set (which approximately represents a gene, as 95% of the soybean probe sets are unique) was declared to be detectable (or present) if it had a present call in 40% of the chips.
The principal data pre-processing method used in our data analyses was GC-RMA, which included a GC-RMA background correction, quantile normalization, and computation of gene summary values from the corrected probe-level data. Background correction was performed with the model-based procedure  using sequence information as implemented in the Bioconductor package gcrma http://bioconductor.org/packages/2.0/bioc/html/gcrma.html. Quantile normalization  and Tukey's median-polish algorithm based gene summary values were computed using a modified C program capable of processing several thousand GeneChips simultaneously (L. Bao, unpublished results).
As an alternative method for data pre-processing, we used the Affymetrix MAS5 default algorithms available in the affy package. Global scaling was implemented to assure that all the chips had the same trimmed mean signal intensity using all soybean probe sets.
Our finding that most detectable transcripts show significant changes raises important questions about the assumptions underlying normalization procedures. MAS5 assumes that the sum of the expression levels of all transcripts remains constant. Quantile normalization assumes this also, and assumes in addition that the distribution of probe signals remains constant across all arrays. Invariant set approaches assume that a pre-determined set of genes show no changes. The assumption that the sum of the expression levels of all transcripts remains constant is a highly robust one because the amount of labeled cRNA added to a microarray hybridization is always normalized. Thus even if the amount of mRNA per cell in the source tissue varied (which should be reflected in variation in the sum of the expression levels of all transcripts) this variation is cancelled by the experimental procedure. Accurately capturing changes in the amount of mRNA per cell in source tissue is currently technically infeasible.
To estimate whether normalization had systematically biased the measured transcriptional changes, we determined the mode of the changes. Under the assumption that the most common transcriptional changes will be the smallest in magnitude, the mode should lie close to 1.0 (i.e. no change). In fact the mode in each contrast studied ranged between 0.99 and 1.01.
Linear mixed model analysis of gene expression levels
A total of 648 soybean GeneChips from 24 complete experimental blocks were selected for the high-level data analyses. The high-level analyses were performed separately for each soybean probe set, using Linear Mixed Model Analysis (LMMA) in the PROC MIXED procedure of SAS (SAS 9.1 for Windows, SAS Institute Inc., Cary, NC, USA) . The fixed factors of the linear mixed model included genotype, treatment (Upper, Lower and Mock), and time (9 am, 10:30 am, and 12 pm), and the genotype by treatment interaction. Random factors included in the model were: block, the interactions block by genotype, block by treatment, block by time and block by genotype by treatment, and residual. Variance components were estimated by the method of Residual Maximum Likelihood, and F-tests were performed for all fixed factors.
All F-tests across all genes and all the fixed factors were considered as one family of tests, and the resultant collection of p-values was used to compute adjusted p-values based on two False Discovery Rate (FDR) controlling methods, the linear step-up method of Benjamini and Hochberg  and the two-stage linear step-up method (TST) of Benjamini et al . In addition, the p-values were used to compute the q-values for the positive FDR criterion as described by Storey and Tibshirani . Because the differences between the FDR methods were small, we mainly report results obtained with the TST method.
We also estimated and tested individual contrasts of interest, most of which were linear functions of the cultivar by treatment interactions. Because the F-test for the genotype by treatment interaction was found to be significant for the majority of probe sets (24,669 out of 28,351), we employed a one-step approach for testing the significance of individual contrasts, i.e. we tested all contrasts of interest for all genes rather than for only those genes with significant F-test for the genotype by treatment interaction. We were interested in identifying the differentially expressed genes for the contrast "infection response", which compares gene expression between pathogen and mock inoculation in a given cultivar and infection court, and these are used to determine which genes are up- or down-regulated within a given infection court as a result of the pathogen attack in a given cultivar. The infection responses were split into two subtypes, i.e., Upper vs. Mock, and Lower vs. Mock. The collection of the t-tests for all contrasts across all genes was considered as a family of tests, and significant contrasts were determined with the same FDR methods as those used for the F-tests.
Wilcoxon signed rank tests of differential gene expressions
As an alternative high-level analysis to LMMA, we used the non-parametric Wilcoxon Signed Rank Test using the wilcox.test function (R package stats version 2.6.0). For each gene and within each cultivar, we formed two pairs of samples, Upper and Mock, and Lower and Mock, respectively. Each consisted of the GC-RMA log2 scale signal intensity data for the paired treatments (e.g. Upper and Mock) (sample size = 3 × 24 = 72). The collection of p-values across all genes was considered as a family of tests, and significant contrasts were determined with the same FDR methods as those used earlier. The same procedure was used to perform the Wilcoxon signed rank test and FDR control with the preprocessed data generated by MAS5.
Kolmogorov-Smirnov test for two distribution patterns
We used the Kolmogorov-Smirnov (KS) test to test for significant differences in the distributions of gene expression changes for specific categories of genes. The expression changes in response to pathogen inoculation were considered for genes in which the change was significant by LMMA contrast analysis of infection responses using SAS Proc Mixed. The distributions of expression changes were examined for genes in six functional categories related to infection: "defense and disease", "signal transduction", "transcription", "intracellular trafficking", "cell structure" and "metabolism". For this purpose we used the functional category of each gene drawn from annotation of the Affymetrix soybean GeneChip by the Goldberg group at the University of California, Los Angeles http://estdb.biology.ucla.edu/seed; we have reviewed this annotation extensively and found it very reliable. In each case, the distribution within the functional category was compared to the distribution for the total gene set, including the genes within the functional category, to test the null hypothesis that the observed distribution of expression changes affected all genes irrespective of functional category. This is a more conservative test than testing the distribution of a subset of genes against the remainder of the genes. To test for distribution differences within specific ranges of gene expression changes, three ranges were chosen: all genes with less than two-fold regulation (up or down); all genes with less than 1.5-fold regulation (up or down); and all genes with less than 1.2-fold regulation (up or down). Then the distribution of changes was compared between a functional category and all genes within the chosen expression range. For all these tests we used the ks.test function in the R package stats (version 2.7.0). Because we compared the distribution of the genes in a category to that of all genes, and because the gene expression profiles have complex correlation structures and thus the values for individual genes may not be independent , we computed an empirical p value for each test by the following procedure. For any gene category with n genes, we randomly sampled n genes from the list of all significant genes and recomputed the KS test statistic. This process was repeated 10,000 times, and the collection of ordered KS statistics was then used to compute an empirical p-value for the category. Finally, FDR control was performed on the set of p-values from all KS tests using the TST method .
Quantitative real time RT-PCR (qRT-PCR) assay
To provide a comparison between the microarray assays and another commonly used measure of gene expression differences, qRT-PCR, gene expression differences (infected versus mock) for a set of selected genes were measured by qRT-PCR during a pilot experiment. The pilot experiment included four biological replicates and was conducted using the same methods as described above, with the following modifications: for each of the two cultivars (V71-370 and Sloan) four replicates of 30 plants (grown together in the same growth chamber) were harvested 5 days after inoculation, and a single segment of infected tissue (treatment "Inoculated") comprising the upper and lower infection courts was excised from each plant. RNA was extracted from each pool of 30, then equal amounts of the RNAs were pooled from four biological replicates for microarray analysis. In parallel, an equal number of mock-inoculated plants were harvested and RNA extracted pooled in the same way.
For qRT-PCR assays, 22 soybean genes of interest (including 4 housekeeping (HK) genes) with varied levels of gene expression were selected (Table 5). Equal amounts of RNA from each of the four biological replicates were pooled for the qRT-PCR assays. The primers [see Additional file 9] were designed using the Beacon Designer 4.0 (Premier Biosoft International, Palo Alto, CA) and synthesized by Integrated DNA Technologies, Inc.(Coralville, IA). The amplicons for development of standard curves were prepared first using the purified total RNA as a template and the SuperScript™ III First-Strand Synthesis System for RT-PCR (Invitrogen™ Life Technologies) and then using conventional PCR with the synthesized cDNA as template. The amplicons were purified using a QIAquick PCR Purification Kit (Qiagen) and checked for quality in the Bioanalyser.
qRT-PCR assays were carried out by the Virginia Bioinformatics Institute Core Laboratory Facility. Briefly, total RNA (1 μg) was transcribed to cDNA using first strand cDNA synthesis reagents (Invitrogen Corp., Carlsbad, CA) in a total volume of 20 μl. Standard curves were produced with serial 10-fold dilutions of cDNA products starting from 10 pg/μl. Each 25 μl qRT-PCR reaction consisted of 300 nM sense and anti-sense primers, 1 μl 10× diluted cDNA, and 12.5 μl SYBR Green I PCR Mastermix (Applied Biosystems, Foster City, CA). Each reaction was run in triplicate for both the standard and samples. PCR reactions were performed on a Bio-Rad i-cycler (BioRad, Hercules, CA) under the following conditions: 95°C for 3 min, 40 cycles of 95°C for 10 s, 56°C for 45 s to calculate cycle threshold (CT) values, followed by 95°C for 1 min, 55°C for 1 min, and 80 times of 55°C for 10 s, increasing temperature by 0.5°C each cycle to obtain melt curves. The BioRad iCycler IQ 3.1 Optical System Software was used and PCR efficiency (E) was estimated using the equation (1+E) = 10(-1/slope) . Pathogen-inoculated vs Mock expression ratio of a gene of interest was calculated from the equation:
Where E denotes PCR efficiency, CT denotes cycle threshold value, ΔCT is the cycle threshold difference between a mock sample and its corresponding pathogen-inoculated sample, gene denotes the gene of interest, NF denotes the normalization factor calculated from the geometric mean of the raw ratios (i.e., (1 + E h ) ΔCTh, here h denotes a house keeping gene) of the three selected housekeeping genes . Four housekeeping genes (the first four rows of Table 5, including actin) were initially evaluated, based on the microarray data from the pilot experiment. The best three (marked [HK] in Table 5) were selected using the gene-stability measure and ranking method .
The microarray data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus  and are accessible through GEO Series accession number GSE11611 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11611.
Aimone JB, Gage FH: Unbiased characterization of high-density oligonucleotide microarrays using probe-level statistics. J Neurosci Methods. 2004, 135: 27-33. 10.1016/j.jneumeth.2003.11.016.
Seo J, Hoffman EP: Probe set algorithms: is there a rational best bet?. BMC Bioinformatics. 2006, 7: 395-10.1186/1471-2105-7-395.
Wodicka L, Dong H, Mittmann M, Ho MH, Lockhart DJ: Genome-wide expression monitoring in Saccharomyces cerevisiae. Nat Biotechnol. 1997, 15: 1359-1367. 10.1038/nbt1297-1359.
Puthoff DP, Ehrenfried ML, Vinyard BT, Tucker ML: GeneChip profiling of transcriptional responses to soybean cyst nematode, Heterodera glycines, colonization of soybean roots. J Exp Bot. 2007, 58: 3407-3418. 10.1093/jxb/erm211.
Meyer RA, Meyer MH, Ashraf N, Frick S: Changes in mRNA gene expression during growth in the femoral head of the young rat. Bone. 2007, 40: 1554-1564. 10.1016/j.bone.2007.01.013.
Lenk GM, Tromp G, Weinsheimer S, Gatalica Z, Berguer R, Kuivaniemi H: Whole genome expression profiling reveals a significant role for immune function in human abdominal aortic aneurysms. BMC Genomics. 2007, 8: 237-10.1186/1471-2164-8-237.
Chu S, DeRisi J, Eisen M, Mulholland J, Botstein D, Brown PO, Herskowitz I: The transcriptional program of sporulation in budding yeast. Science. 1998, 282: 699-705. 10.1126/science.282.5389.699.
Lee BH, Henderson DA, Zhu JK: The Arabidopsis cold-responsive transcriptome and its regulation by ICE1. Plant Cell. 2005, 17: 3155-3175. 10.1105/tpc.105.035568.
Jenner RG, Young RA: Insights into host responses against pathogens from transcriptional profiling. Nat Rev Microbiol. 2005, 3: 281-294. 10.1038/nrmicro1126.
Tao Y, Xie Z, Chen W, Glazebrook J, Chang HS, Han B, Zhu T, Zou G, Katagiri F: Quantitative nature of Arabidopsis responses during compatible and incompatible interactions with the bacterial pathogen Pseudomonas syringae. Plant Cell. 2003, 15: 317-330. 10.1105/tpc.007591.
Ithal N, Recknor J, Nettleton D, Hearne L, Maier T, Baum TJ, Mitchum MG: Parallel genome-wide expression profiling of host and pathogen during soybean cyst nematode infection of soybean. Mol Plant Microbe Interact. 2007, 20: 293-305. 10.1094/MPMI-20-3-0293.
Klink VP, Overall CC, Alkharouf NW, MacDonald MH, Matthews BF: Laser capture microdissection (LCM) and comparative microarray expression analysis of syncytial cells isolated from incompatible and compatible soybean (Glycine max) roots infected by the soybean cyst nematode (Heterodera glycines). Planta. 2007, 226: 1389-1409. 10.1007/s00425-007-0578-z.
Klink VP, Overall CC, Alkharouf NW, MacDonald MH, Matthews BF: A time-course comparative microarray analysis of an incompatible and compatible response by Glycine max (soybean) to Heterodera glycines (soybean cyst nematode) infection. Planta. 2007, 226: 1423-1447. 10.1007/s00425-007-0581-4.
Tucker ML, Burke A, Murphy CA, Thai VK, Ehrenfried ML: Gene expression profiles for cell wall-modifying proteins associated with soybean cyst nematode infection, petiole abscission, root tips, flowers, apical buds, and leaves. J Exp Bot. 2007, 58: 3395-3406. 10.1093/jxb/erm188.
Smith EN, Kruglyak L: Gene-environment interaction in yeast gene expression. PLoS Biol. 2008, 6: e83-10.1371/journal.pbio.0060083.
Burnham KD, Dorrance AE, VanToai TT, St. Martin SK: Quantitative trait loci for partial resistance to Phytophthora sojae in soybean. Crop Science. 2003, 43: 1610-1617.
Cicek MS, Chen P, Saghai Maroof MA, Buss GR: Interrelationships among agronomic and seed quality traits in an interspecific soybean recombinant inbred population. Crop Sci. 2006, 46: 1253-1259. 10.2135/cropsci2005.06-0162.
Wu Z, Irizarry RA, Gentleman R, Murillo FM, Spencer F: A model based background adjustment for oligonucleotide expression arrays. J Amer Stat Assoc. 2004, 99: 909-917. 10.1198/016214504000000683.
Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003, 19: 185-193. 10.1093/bioinformatics/19.2.185.
Caldo RA, Nettleton D, Peng J, Wise RP: Stage-specific suppression of basal defense discriminates barley plants containing fast- and delayed-acting Mla powdery mildew resistance alleles. Mol Plant Microbe Interact. 2006, 19: 939-947. 10.1094/MPMI-19-0939.
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc, Ser B. 1995, 57: 289-300.
Benjamini Y, Krieger AM, Yekutieli D: Adaptive linear step-up procedures that control the false discovery rate. Biometrika. 2006, 93: 491-507. 10.1093/biomet/93.3.491.
Storey JD, Tibshirani R: Statistical methods for identifying differentially expressed genes in DNA microarrays. Methods Mol Biol. 2003, 224: 149-157.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005, 102: 15545-15550. 10.1073/pnas.0506580102.
Pfaffl MW: A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29: e45-10.1093/nar/29.9.e45.
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3: RESEARCH0034-10.1186/gb-2002-3-7-research0034.
Edgar R, Domrachev M, Lash AE: Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002, 30: 207-210. 10.1093/nar/30.1.207.
We thank J. Mills, S. A. Berry, H. Wang, F. Cruz, W. Pipatpongpinyo, K. Broders, C. Cruz, M. A. Ortega, S. Gordon, K. Gearhart, E. Denardo, S. Anderson, P. Roongsattham, F. Tierney, E. Helliwell, A. Johnston, L. Leedy, B. James, S. Kale, B. Durham, R. Presler, G. Torto-Alalibo, L. Waller, B. Smith, L. Falin, A. Ferreira, D. Dou, M. C. Chibucos, N. Galloway, L. Fakunmoju, K. Cooper, and B. Howard for technical assistance, and T. Jewell for manuscript preparation. This work was supported by the NSF Plant Genome Research Program (Award DBI-0211863).
Writing team: BMT, LZ, IH, AED; Project management: BMT, MASM., IH., AED, SKSM, LZ, ST; Inoculation assays and sampling: SXM, AED; RNA preparations and microarray assays: LZ, RH, FDA, KK, CE, AJ; Database management: ST, LZ, KK; Statistical analyses: LZ, LB, BMT, IH. All authors have read and approved the final manuscript.
Electronic supplementary material
Additional file 1: Correlation among gene expression changes in response to P. sojae infection measured in independent experiments by qRT-PCR and microarrays. The data provided compare gene expression changes in response to pathogen infection measured in two experiments 2 years apart, and the correlation between qRT-PCR and microarray measurements. (PDF 207 KB)
Additional file 2: Supplemental Table 1. Intersection of six non-random sub-groups with the whole experiment. The data describe the intersections among six independent sub-groups of the whole experiment at different FDR (TST) levels. (PDF 64 KB)
Additional file 3: Supplemental Table 2. Comparisons among the intersection or union of the sub-groups and the whole experiment. The data demonstrate the differences in the power between smaller experiments and the complete experiment, and different ways of combining the results of the smaller experiments. (PDF 58 KB)
Additional file 4: Supplemental Figure 2. Distribution of gene expression changes during infection of genotype VPRIL9. The data describe the fractions of genes showing different levels of gene expression changes during infection of moderately susceptible genotype VPRIL9. (PDF 88 KB)
Additional file 5: Supplemental Figure 3. Distributions of expression changes in the Sloan upper infection court among six functional categories. The data describe the fractions of genes with different expression changes in the Sloan upper infection court within specific functional categories. (PDF 135 KB)
Additional file 6: Supplemental Figure 4. Distributions of expression changes in the V71-370 lower court among six functional categories. The data describe the fractions of genes with different expression changes in the V71-370 lower court within specific functional categories. (PDF 147 KB)
Additional file 7: Supplemental Figure 5. Distributions of expression changes in the V71-370 upper court among six functional categories. The data describe the fractions of genes with different expression changes in the V71-370 upper court within specific functional categories. (PDF 147 KB)
Additional file 8: Supplemental Table 3. Distributions of expression differences of genes in functional categories compared with the overall differences distribution. Detailed data summarizing numbers of genes with different expression differences in six functional categories compared with the overall distribution, in two cultivars and both infection courts. (PDF 81 KB)
Additional file 9: Supplemental Table 4. Genes and primers used for quantitative real time RT-PCR assay. The Affymetrix ID, functional annotation and primer sequences of the 22 genes used in the qRT-PCR assay. (PDF 62 KB)
About this article
Cite this article
Zhou, L., Mideros, S.X., Bao, L. et al. Infection and genotype remodel the entire soybean transcriptome. BMC Genomics 10, 49 (2009). https://doi.org/10.1186/1471-2164-10-49
- False Discovery Rate
- Transcriptional Change
- Soybean Cyst Nematode
- Soybean Genotype
- Transcriptional Modulation