Transcriptional profiling of root-knot nematode induced feeding sites in cowpea (Vigna unguiculata L. Walp.) using a soybean genome array

Heterologous microarray platform

In order to elucidate the plant response to root-knot nematodes, infected resistant CB46 (incompatible interaction) was compared with non-infected CB46, and infected susceptible null-Rk (compatible interaction) was compared with non-infected null-Rk. Two time points were chosen for this analysis i.e., 3 and 9 dpi. Nine dpi was selected as a critical time point because sequentially assayed histological sections during 21 days of infection revealed that at 9 dpi the first subtle differences appeared between incompatible and compatible interactions [8]. The 3-dpi samples provided a time point prior to visible differences histologically, between incompatible and compatible interactions.

The average number of soybean probe sets which had "present" call in 9-dpi samples was 10,521, which comprised 28% of the total number of soybean probe sets on the soybean GeneChip. Similarly, the average number of "present" calls in 3-dpi samples was 10,685 (~28.5% of all soybean probe sets). When soybean RNA was used to hybridize the soybean GeneChip, the "present" call percentage ranged from 70-75% [25]. Therefore about 40% (28/70) of the content of the soybean GeneChip is informative for cowpea.

Data quality

Principal component analysis (PCA) was done on triplicate data for each treatment at 9 dpi in order to visualize the overall genome response to nematode infection in resistant and susceptible cowpea genotypes. PCA on conditions (treatments) for the 9-dpi sample are shown in Fig. 1. PCA component 1 (40.32% of total variance) and PCA component 2 (33.62% of total variance) comprised the majority of the described variance (73.94%). The PCA showed a clear separation between the two genotypes when infected with nematodes whereas the two genotypes clustered together when there was no external stimulus (non-infected control). Because we used cowpea near-isogenic lines for this analysis, clustering of the non-infected control samples was expected. This confirmed the robustness of the experimental design. The PCA plot for 3-dpi samples is presented in additional file 1.

Because there were only two biological replicates for the 3-dpi samples, a correlation analysis was performed between the two replicates of each treatment. A representative MA-scatter plot is shown in Fig. 2. Correlation coefficients between replicates ranged from 0.932 to 0.973 at a p value cut-off of 0.05, showing the robustness of the data (correlation coefficients greater than 0.9 are considered to be high [25]).

Gene expression in incompatible and compatible interactions in 9-dpi root samples

At 9 dpi three different comparisons were made. In the first comparison gene expression in the resistant CB46 roots infected with root-knot nematodes (incompatible interaction) was compared with non-infected CB46. Secondly, a comparison in gene expression was made between the infected null-Rk (compatible interaction) and non-infected null-Rk. Finally, a comparison of gene expression was made between the infected Rk and infected null-Rk near-isogenic lines. The final comparison was important because the near-isogenic lines were predicted to show differential gene expression for the genes which are critical for nematode resistance or nematode parasitism.

In the incompatible interaction 552 (~5.3% of total expressed probe sets) genes were significantly differentially expressed between the Rk-infected and non-infected treatments based on the statistical test (see materials and methods). These genes were then passed through a 1.5 fold-change filter. The geometric mean of the normalized expression intensities of all samples under one condition were used to calculate the fold-change ratios. If the ratio of geometric means of the infected sample and the non-infected control for a probe set was ≥ 1.5, then that particular probe set was categorized as 1.5-fold or more up-regulated and if the ratio was ≤ 0.67 then the probe set was categorized as down-regulated by 1.5-fold or more. 141 genes showed 1.5-fold or more up-regulation and 59 genes were down-regulated by 1.5-fold or more in the Rk-infected compared with the Rk-non-infected treatment (Fig. 3a and 3b). In the compatible interaction 1,060 genes passed the statistical filter (~10% of total expressed probe sets). Among these 1,060 genes 218 were 1.5-fold or more up-regulated and 41 genes were 1.5-fold or more down-regulated in the infected null-Rk compared to the non-infected null-Rk treatment (Fig. 3a and 3b). In this context it can be noted that in the current study the number of differentially expressed genes is lower than some of the other microarray studies previously published in the field of plant-microbe interactions. One of the reasons may be due to the use of a heterologous GeneChip. Nevertheless, the information generated will be very valuable as this is the first report on the cowpea root-knot nematode interaction. In the final comparison between the two near-isogenic lines infected with Rk-avirulent root-knot nematodes, 638 genes (~6% of total expressed probe sets) passed the statistical filter. Among the differentially expressed genes only 20 genes were 1.5-fold or more up-regulated in the infected Rk than in the infected null-Rk treatment and 100 genes were 1.5-fold or more down-regulated in the infected Rk than in the infected null-Rk treatment.

A gene ontology based analysis was carried out to categorize the differentially expressed genes into different functional classes. In the incompatible interaction the most abundant functional class observed for the up-regulated probe sets (Fig. 4a) was genes involved in metabolism (32.8%), followed by proteins with binding function (28.8%) and genes involved in cell rescue and defense (13.7%). In the compatible interaction the most abundant functional classes in up-regulated probe sets (Fig. 4a) were metabolism (30.6%), proteins with binding function (26.9%), and genes involved in protein fate (13.2%). For down-regulated probe sets the most abundant functional classes in the incompatible interaction (Fig. 4b) were proteins with binding function (39.2%), metabolism (21.5%), and interaction with the environment (15.6%), whereas in the compatible interaction the most abundant classes (Fig. 4b) were proteins with binding function (23.6%), cellular transport (18.4%), and systemic interaction with the environment (13.1%). A significant number of probe sets were also categorized under unclassified or unknown proteins in all the above comparisons.

Genes differentially expressed between the two infected NILs were also functionally classified (Fig. 5). In the probe sets up-regulated in infected resistant CB46 over infected susceptible null-Rk, the most abundant classes were cellular transport (36.8%), proteins with binding function (21%), and proteins involved in transcription (15.7%). Among the probe sets which were down-regulated in the resistant genotype over the susceptible genotype, the most abundant functional classes were metabolism (35.1%), proteins with binding function (24.1%), and systemic interaction with the environment (10.9%).

Gene expression in incompatible and compatible interactions in 3-dpi root samples

Similar to the 9-dpi samples, three different comparisons were made for 3-dpi root samples. In the first comparison gene expression in the resistant Rk plants (incompatible interaction) infected with root-knot nematodes was compared with non-infected Rk plants. A comparison in gene expression also was made between the infected null-Rk (compatible interaction) and non-infected null-Rk. Finally, a comparison of gene expression was made between the infected Rk and infected null-Rk near-isogenic lines.

In the incompatible interaction 746 (~6.9% of total expressed probe sets) genes were significantly differentially expressed based on the statistical test. These genes were then passed through a fold-change filter based on log₂ ratios. 65 genes showed 1.5-fold or more up-regulation and 129 genes were down-regulated by 1.5-fold or more in the incompatible interaction (Fig. 6a and 6b). In the compatible interaction 623 genes passed the statistical filter (~5.8% of total expressed probe sets). Among these 623 genes 81 were 1.5-fold or more up-regulated and 148 genes were 1.5-fold or more down-regulated in the compatible interaction (Fig. 6a and 6b). In the final comparison between the two near-isogenic lines infected with Rk-avirulent root-knot nematodes, 197 genes (~1.8% of total expressed probe sets) passed the statistical filter. Among the differentially expressed genes only 4 genes were 1.5-fold or more up-regulated in the resistant CB46 compared to susceptible null-Rk, and 10 genes were 1.5-fold or more down-regulated in the CB46 compared to null-Rk plants.

A gene ontology based analysis was carried out to categorize the differentially expressed genes into different functional classes. In the incompatible interaction the most abundant functional class observed for the up-regulated probe sets (Fig. 4a) was genes involved in metabolism (42.8%), followed by proteins with binding function (25.3%), and genes involved in interaction with the environment (15.8%). In the compatible interaction the most abundant functional classes in up-regulated probe sets (Fig. 4a) were proteins with binding function (41.3%), metabolism (29.3%), and genes involved in protein fate (10.6%). For down-regulated probe sets the most abundant functional classes in the incompatible interaction (Fig. 4b) were proteins with binding function (34.7%), metabolism (29.6%), and genes involved in protein fate (20.3%), whereas in the compatible interaction the most abundant classes (Fig. 4b) were proteins with binding function (36%), metabolism (30.4%), and genes involved in transcription (14.4%). A significant number of probe sets were also categorized under unclassified or unknown proteins in all the above comparisons.

Comparison of gene expression between 9-dpi and 3-dpi root samples

A comparison of differentially expressed genes was made between the 9-dpi and 3-dpi samples. In the incompatible interaction 188 genes (137 up-regulated by 1.5-fold or more and 51 down-regulated 1.5-fold or more) were uniquely expressed at 9 dpi and 182 (62 up-regulated by 1.5-fold or more and 120 down-regulated 1.5-fold or more) genes were unique for 3 dpi. 12 genes were differentially expressed at both time points (additional file 8). In the compatible interaction 238 genes (198 up-regulated by 1.5-fold or more and 40 down-regulated 1.5-fold or more) were uniquely expressed at 9 dpi and 208 genes (61 up-regulated by 1.5-fold or more and 147 down-regulated 1.5-fold or more) were unique for 3 dpi. 21 genes were differentially expressed at both time points (additional file 9).

Plant material

Two near-isogenic lines (NIL) differing in presence or absence of gene Rk were used. The two parents used to develop the NIL were M. incognita race 3 resistant cowpea genotype 'CB46' (homozygous resistant, RkRk) and a highly susceptible genotype 'Chinese Red' (homozygous susceptible, rkrk). The F₁ was backcrossed to recurrent parent CB46 (BC₁), homozygous Rk plants were discarded in BC₁F_2, and non-segregating rkrk plants were advanced to the next back-cross (BC₂). Repeated backcrossing and selection was used to recover the rkrk line in the CB46 background. BC₆F₄ progenies were used for all the experiments described here. The rkrk line is referred to as the null-Rk line from here on.

Nematode inoculum

Eggs of M. incognita race 3 (isolate Beltran) cultured on susceptible tomato host plants were extracted from roots using 10% bleach solution [43]. This isolate is avirulent to gene Rk in CB46. Eggs were hatched in an incubator at 28°C and J2 were collected in fresh deionized water. The J2 inoculum was prepared according to the experimental requirements.

Root infections for microarray analysis

Seeds of CB46 and null-Rk cowpea lines were surface-sterilized using 10% (v/v) bleach solution and planted singly in seedling growth pouches. Plants were grown under controlled environmental conditions of 26.7°C ± 0.5°C constant temperature and daily light/dark cycles of 16/8 hours. This temperature was used because it lies within the optimum temperature range of 26 - 28°C for development and reproduction of M. incognita on cowpea in growth pouches [44]. Each of 100 pouches (50 pouches for each genotype) was inoculated with 3000 J2 in 5 ml of deionized water 12 days after planting (dap). This inoculum level was found to be optimum in a previous study (Das and Roberts, unpublished data) and it generated uniform infection throughout the root system. As a result the amount of infected tissue was maximized. An equal number of pouches were mock-inoculated with 5 ml of deionized water to use as non-infected controls. Infected and non-infected plants were arranged in a completely randomized design. Nematode infected root tissue was excised using a sterile scalpel at 3 days post-inoculation (dpi) and 9 dpi, respectively, under a magnifying glass and flash frozen immediately in liquid nitrogen. In previous studies (Das and Roberts, unpublished) the infected root regions showed swelling at 3 dpi which is indicative of initiation of giant cell formation, and at 9 dpi prominent galling was visible on infected roots. The infected tissue was collected based on these visual indicators. For each biological replicate infected tissue was collected from 7 plants picked randomly and pooled together. This was done in order to obtain enough biological material for RNA isolation. Similarly root tissue was collected from equivalent root regions of the control plants (tissue near root tips of secondary and tertiary roots) and flash frozen. Galled tissue was excised by cutting immediately adjacent to the root-gall in order to minimize the amount of non-infected tissue included in the assays. The harvested tissue was stored at -80°C until RNA isolation. A few infected root pieces were stained in acid fuchsin [45] to confirm the nematode infection.

RNA isolation

RNA from nematode infected and non-infected root tissue was isolated using RNeasy plant mini kit (QIAGEN Inc., Valencia, CA, USA) according to the manufacturer's protocol. One volume of Plant RNA Isolation Aid (Ambion, Austin, TX, USA) per unit mass of frozen tissue (ml/g) was added before the tissue homogenization step for removal of common contaminants such as polysaccharides and polyphenolics. RNA was treated with RNase-Free DNase set (QIAGEN Inc., Valencia, CA, USA) to digest any genomic DNA which might be present. RNA was quantified using a UV-spectrophotometer. RNA quality and integrity was examined using RNA Lab-On-A-Chip (Caliper Technologies Corp., Mountain View, CA, USA) evaluated on an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA, USA).

Soybean genome array

Phylogenetic relationships based on the conserved sequences within Papilionoideae legumes imply that Vigna (cowpea) is closely related to soybean [46]. Since a commercial cowpea genome array was not available, a soybean genome array (Affymetrix Inc., Santa Clara, CA, USA) was used for transcriptome profiling in cowpea. The soybean genome array contains 37,500 probe sets derived from soybean (Glycine max L.) unigenes. This represents 61% of the total probe sets on the chip, with the remainder targeting two pathogens important for soybean research, of which 15,800 (26%) probe sets target Phytophthora sojae (a water mold) and 7,500 (12%) probe sets target Heterodera glycines (soybean cyst nematode). This array uses probe sets composed of 11 probe pairs to measure the expression of each gene. Each probe pair consists of a perfect match (PM) probe and a mismatch (MM) probe (see also http://www.affymetrix.com/products_services/arrays/specific/soybean.affx).

Array hybridization

Double-stranded complementary deoxyribonucleic acid (cDNA) was synthesized using SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen) and T7-oligo (dT) promoter primers. The IVT Labeling Kit (Affymetrix) was then used to synthesize biotin-labeled complementary RNA (cRNA) from template cDNA by in vitro transcription. Twelve to 16 μg labeled cRNA was fragmented by metal-induced hydrolysis to 35-200 base fragments following Affymetrix protocols. 10 μg labeled, fragmented cRNA was then hybridized at 45°C with rotation for 16 h in an Affymetrix microarray Hybridization Oven 320 on Affymetrix soybean genome arrays. The arrays were washed and stained using streptavidin phycoerythrin on an Affymetrix Fluidics Station 450. The arrays were scanned on a Hewlett-Packard GeneArray scanner. cRNA synthesis and array hybridizations were performed in the Genomics Core Facility http://www.genomics.ucr.edu at the University of California, Riverside.

Data analysis

For 9-dpi samples three biological replicates were used for each of the four treatments (Rk infected and non-infected, and Rk-null infected and non-infected), requiring 12 soybean GeneChips. For the 3-dpi samples two biological replicates were used for each treatment requiring 8 GeneChips. The data from all 20 chips (CEL and CHP files) are publicly available in Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/, platform GPL 4592, series GS13631). Expression signals were first analyzed in GeneChip operating software 1.3 (GCOS, Affymetrix Inc.) to determine the "present" probe set list. To detect "present" calls using GCOS software we used all the probe pairs in the probe set as "stat pairs". The definition of the term "stat pairs" is the number of probe pairs per probe set used in the analysis. During data analysis a specified subset of probe pairs can be selected by a probe mask file, but if the default settings are used no probe mask file will be applied and all the probe pairs will be used as "stat pairs" (11 probe pairs in the case of the soybean GeneChip). The detection algorithm uses probe pair intensities to generate a detection p-value and assign a "present", "marginal", or "absent" call. Each probe pair in a probe set has a potential vote in determining whether the measured transcript is or is not "present". The vote is described by the discrimination score (R), which is calculated for each probe pair and compared to a predefined threshold, Tau. Probe pairs with R higher than Tau vote "present" and the voting result is summarized as a p-value. The greater the number of discrimination scores (R) that are above Tau, the smaller the p-value and the more likely the given transcript is truly present in the sample. Only probe sets with a "present" call in all three replicates of at least one treatment were considered to be "expressed".

Data normalization and further analysis was carried out in GeneSpring GX 7.3 (Agilent Technologies, Palo Alto, CA, USA). Robust Multiarray Average (RMA, [47, 48]) normalization was performed. Each chip was normalized to the 50th percentile and each gene was normalized to the median. As we were only interested in plant response to nematode infection, all the probe set data from P. sojae and H. glycines were excluded from any further analysis.

Principal component analysis (PCA) is often used to reduce multidimensional data sets to lower dimensions for summarizing the most important part of the data while simultaneously filtering out the background errors. PCA involves the calculation of the eigenvalue decomposition of a data covariance matrix or singular value decomposition of a data matrix, usually after mean centering the data for each attribute. The results of PCA are usually discussed in terms of component scores and loadings [49]. PCA on conditions (treatments) based on all genes which were present in at least one chip in the 9-dpi and 3-dpi samples were carried out to visualize the overall genome response to nematode infection in the resistant and susceptible cowpea genotypes.

For the 3-dpi samples, with only two biological replicates available, a Pearson correlation coefficient was calculated for normalized values of all probe sets between the two replicates of each treatment to determine the robustness of the data. This analysis was carried out in dChip software [50].

Differentially expressed genes were identified using a one-way analysis of variance (ANOVA) with a p-value cut-off of 0.05. A multiple testing correction was performed using the Bonferroni error correction model [51]. False discovery rate (FDR) was set at 5.0%. Subsequently, differentially expressed genes were filtered for 1.5-fold change in expression level between the control and nematode infected treatment for both genotypes and also between the nematode infected treatments of the resistant and susceptible genotypes.

Validation of the use of Affymetrix soybean genome array for cowpea transcriptional profiling

In our recent related work [12] the Affymetrix soybean genome array was used successfully to identify single feature polymorphisms in cowpea and the statistical data were validated using PCR amplicon sequencing. Thus, we were able to correctly identify polymorphisms between two cowpea genotypes at a resolution as high as the single nucleotide level. Also, we conducted a small analysis to look at the sequence homology between sequence information files (SIF) of 30 probe sets selected to carry out PCR validation of predicted SFPs and their corresponding cowpea sequences. The homology ranged from 87% to as high as 94.5%. Though this analysis is not exhaustive, it provided a good indication that the homology between cowpea and soybean genomes is quite high at least in the SIF regions from where the probe sets were designed. For this work RNA was used as surrogate for genomic DNA. These data established that the soybean probe sets faithfully measure cowpea transcripts, validating the general reliability of the soybean-based platform for cowpea.

Annotations and functional classification of genes

The soybean genome array unigene sequences were used to query (using blastx) Arabidopsis translated gene models (version 7.0) from The Arabidopsis Information Resource (TAIR, http://www.arabidopsis.org) and Medicago truncatula 2.0 assembly release http://www.medicago.org. Annotations for the Affymetrix soybean probe sets were compiled into a browser called HarvEST:SoyChip which can be accessed online http://www.harvest-web.org or downloaded for Windows installation http://harvest.ucr.edu/. The E value cut-off for the gene annotations was equal to or less than E-10, E0 being a near perfect match.

Gene ontology based classification was obtained by transferring the corresponding Arabidopsis gene models to Munich Information Center for Protein Sequences Arabidopsis thaliana FunCat database (MIPS, http://mips.gsf.de/proj/funcatDB/search_main_frame.html). Arabidopsis gene models were taken from HarvEST:SoyChip.