Expression profiling on soybean leaves reveals integration of ER- and osmotic-stress pathways

ER stress-specific response

A total of 306 genes were differentially expressed in response to tunicamycin treatment, with 179 being up-regulated and 127 down-regulated (Figure 1). In response to AZC treatment, 352 genes were differentially expressed, 172 up-regulated and 180 down-regulated. To specifically target the UPR-regulated genes, we considered the overlapping genes that significantly responded to both tunicamycin and AZC treatments. This overlapping transcriptional response consisted of 183 genes, 86 up-regulated and 97 down-regulated (Figure 1, see Additional files 1 and 2).

Among the up-regulated genes, 17% represented genes of unknown functions. In the remaining set of genes for which some functional information was available, those that exhibited an ER-stress response signature predominated (see Additional file 2). More specifically, this up-regulated class of genes included those categorized as having a function in (i) protein folding (ii) ERAD and (iii) translational regulation. These results confirmed the activation of the ER-stress response pathway. In fact, the clones on the array with high homology to the known ER-resident molecular chaperones, BiP and calnexin, or the folding catalyst protein disulfide isomerase (PDI), were strongly up-regulated by both tunicamycin and AZC treatments (Figure 2).

The expression of calnexin, an ER multi-functional protein involved in calcium homeostasis and protein folding, was used as a marker for ER-stress activation in a time-course experiment using real-time RT-PCR. Similar levels of induction of calnexin were observed between treatment with tunicamycin and AZC, possibly representing saturation in the expression (Figure 3A).

We found it particularly interesting that three cDNAs related to the PP2C and PP2A protein phosphatases were up-regulated by ER stress (see Additional file 1, AW471739, AW569267, AW509424), as would be expected if they represented genes involved in UPR signaling in stressed cells. Accordingly, the PP2C-related cDNA is a homolog of the yeast PP2C (AAB64644) that regulates the UPR by dephosphorylation of Ire1 [41], a transmembrane protein kinase/endoribonuclease that triggers the UPR [42, 43].

The identification of wheat MLO (transmembrane domain mildew resistance) allelic variants as endogenous substrates of an ERAD-related quality-control machinery provided direct evidence that an ERAD-like mechanism operates constitutively in plants [44]. As part of the ER-stress response, the activation of this turnover mechanism has been observed in genome-wide analyses of Arabidopsis treated with inducers of ER stress and in a transcript-profiling assay of maize endosperm mutants that display a long-term ER stress response [34, 35, 45]. Here we also observed a tunicamycin and AZC up-regulated repertoire of putative ERAD-related genes in soybean, such as those encoding polyubiquitin, ubiquitin conjugating enzyme, alpha subunit of the proteasome, CDC48 and Derlin (see Additional file 1). These results further confirmed that, like in mammalian cells and in yeast, an ER stress-induced quality control mechanism in plants integrates the cellular response to conditions that alter protein folding in the ER.

During conditions of ER stress in mammals, a dynamic balance between the ER processing capacity and the protein synthesis rate is adaptively achieved through a transient and general down-regulation of protein translation, as a component of the ER-stress response [46]. The results of our microarrays were also effective in identifying a series of up-regulated genes related to ribosomal proteins (60S and 40S subunits; see Additional file 1) that might represent regulatory elements in protein translation. Likewise, we found that a translational inhibitor protein (AW508686), a eukaryotic translation initiation factor 3 subunit 10 (AW317679), and a translation elongation factor 1-gamma (AI960794), which are potentially regulators of protein translation, were also responsive to ER stress. Collectively, the global transcriptional analysis of soybean cDNAs in response to ER stressors clearly unmasked the major branches of the conserved ER-stress response, arguing favorably for a good sampling of the genomic representation on our array and for the biological validation of the global analyses.

Osmotic stress-specific response

In response to osmotic stress caused by PEG treatment, a set of 116 up-regulated and 173 down-regulated genes was observed (Figure 1, see Additional files 3 and 4). Of particular interest were genes in two functional classes: genes directly involved in the stress response as cellular protectants, and regulatory genes involved in signaling events downstream of the osmotic stress response.

As representatives of the first class, we detected induction of late embryogenesis abundant (LEA) proteins, heat shock proteins (HSPs), senescence-related proteins, protease inhibitors, enzymes associated with osmolyte distribution, such as sugar transporters, and in osmoprotectant biosynthesis, such as Δ1-pyrroline-5-carboxylate synthase (P5CS), which is involved in proline biosynthesis [47, 48], and sucrose synthase [49, 50] (see Additional file 3). Also in this category, we observed up-regulation of anti-oxidative defense components, such as glutathione peroxidase and glutathione-S-transferase homologs, which are involved in detoxification and protection from reactive oxygen species [51, 52] and lipid transfer proteins that could possibly function in the stress-damaged membrane repair system, or in the regulation of cellular membrane permeability by changing the lipid composition in response to stress [53, 54]. As representatives of the second class, our microarray results detected the presence of two protein phosphatases, PP2C-like and PP1/PP2A PRL1-like coding cDNAs, and a serine/threonine-like protein kinase (see Additional file 3). We also detected PEG-induced members of the NAC family of transcription factors. NAC proteins induced by dehydration have been previously described in Arabidopsis and rice [5, 55].

LEA genes were the most abundant among genes up-regulated by PEG treatment. We examined the osmotic stress response in more detail with real-time RT-PCR analysis of a gene encoding a putative LEA protein (Figure 3A). These data demonstrated that PEG treatment caused a gradual increase in induction of this LEA gene as a 15-fold increase in LEA mRNA level was detected at 4 h post-treatment, reaching 130- and 146-fold-changes at 10 h and 16 h, respectively. In contrast, the LEA gene did not exhibit a significant variation of expression in response to ER-stress inducers and thus could be used as an appropriate monitor of the level of osmotic stress in treated plants.

The antagonistic response to ER stress and osmotic stress has an ER protein-folding signature

Different clones present in the soybean cDNA microarrays that likely encode ER-associated proteins, showed a differential and antagonistic regulation by ER and osmotic stresses. The most dramatic is within the PDI-like sequences where several family members were represented on the array. The protein disulfide isomerases (PDIs) belong to the superfamily of thioredoxin-domain-containing proteins that catalyze the formation of disulfide bonds and play an important role in protein folding. Within the TRX (thioredoxin) superfamily, they make a large gene family, designated PDIL (PDI-like proteins), encompassing disulfide isomerases and oxidoreductases that are associated predominantly with the protein secretory pathway in plants [56]. The biochemically and genetically characterized members of this family were first identified in plants as ER-resident proteins that were induced under ER-stress conditions [57]. Based on sequence comparison, we identified in our soybean arrays four thioredoxin domain-encoding cDNA fragments that were classified as members of the PDIL family. Two clones were highly and specifically induced by both tunicamycin and AZC treatments, whereas two others were only induced by PEG (see Additional files 1 and 3). Analysis of gene expression using real-time RT-PCR in a time-course experiment confirmed the differential regulation of the soybean PDIL gene family by osmotic and ER stresses and revealed that specific members of this family respond inversely to treatment with inducers of ER stress and PEG (Figure 3B). Sequence alignment and the expression pattern suggest to us that both PEG-induced clones refer to the same gene.

Modest overlap of the ER-stress and osmotic-stress transcriptional responses

An overlap of the osmotic-stress and the ER-stress responses is represented by 10 up-regulated genes (see Additional file 5), 75 down-regulated genes (see Additional file 6), 8 ER stress-induced but osmotically repressed genes, and 1 osmotically stress-induced but ER stress-repressed gene (Figure 1, Table 1). Thus, only about 10% of the genes up-regulated by either ER or osmotic stress in our survey population were induced by both treatments. In contrast, a substantial overlap in the genes down-regulated by ER and osmotic stresses was observed, with about 50% and 75% of the genes being in the sets affected by osmotic stress and ER stress, respectively. These results represent a much larger down-regulation of transcripts by ER stress than those reported previously ([34, 35]; see Additional files 2 and 6). Likely, these results reflect substantial differences in experimental design and conditions, including cDNA library origin, plant background, stringency of stress conditions, plant species, setting up and processing of microarrays.

While 25% of the repressed genes that were ER-stress specific were predicted to be related to the secretory pathway (see Additional file 2), the remaining 75% were co-repressed genes that seem to represent a general response of plants to abiotic stresses, as they reflect an inhibition of photosynthesis and development. In fact, a common effect of osmotic and ER stresses revealed by the microarray analysis was a general decrease in the expression of photosynthesis-related genes, including genes that encode the oxygen-evolving enhancer protein, chlorophyll a/b binding protein, small subunit of RuBisCO, NADPH-protochlorophyllide oxidoreductase (involved in chlorophyll biosysnthesis), a thylakoid membrane phosphoprotein and others. There are at least 32 redundant clones involved in photosynthesis that are down-regulated by all three treatments (Table 1). Consistent with photosynthesic inhibition was our observation of three starch synthase homolog cDNAs, and two clones related to cinnamoyl-CoA reductase, a lignin synthesis-related gene being repressed by both stresses (Table 1). Another class represented in the co-repressed category included members involved in hormone biosynthetic pathways, such as 1-aminocyclopropane-1-carboxylate oxidase (ACC, AW508290), involved in ethylene biosynthesis [58], an IAA-amino acid conjugate hydrolase (AW278733), which regulates the level of the auxin indole-3-acetic acid [IAA; [59]], geranylgeranyl hydrogenase and cytochrome P450 monooxygenase that participate in the gibberellin (GAs) biosynthetic pathway [60, 61].

The overlap in genes that shared the up-regulated response revealed two cDNA fragments that encode putative transcription factors belonging to the NAM family (NAC and ATAF2 homologs), two clones that encode DCD-domain-containing proteins (N-rich proteins), two encoding glutathione-S-transferase, a UBA protein gene, an eukaryotic translation initiation factor 5 (eIF5) gene and two cDNAs whose function is unknown (Table 1). A more precise analysis using real-time RT-PCR confirmed a significant induction of co-up-regulated transcripts by tunicamycin, AZC or PEG treatment, except for glutathione-S-transferase (GST) mRNAs which did not show induction by tunicamycin (Figure 4).

If these genes had a role in regulating this branch of the pathway, they would be predicted to be induced early. We tested this possibility with a time-course experiment. Real-time RT-PCR assays during induction demonstrated that the NAC-containing proteins ATAF2 and NAM exhibited an early kinetics of induction, consistent with their putative role as transcriptional factors (Figure 4). Four-hour treatments were sufficient to saturate their expression, which remained high for the duration of the experiment. A similar kinetic pattern was observed for the N-rich DCD-domain transcripts AI973541 and AW184865, which were strongly induced at 4 h and reached maximum accumulation at 10 h post treatment. Accordingly, the N-rich (AW184865) gene has been shown to be rapidly induced during the hypersensitive response in soybean [62]. In contrast to the early induction of the NAC and N-rich related genes by ER and osmotic stresses, the induction of the remaining co-up-regulated genes occurred with delayed kinetics (Figure 4). The induction of the UBA, eIF5, GST (AW472161) and GST (AW397276) transcripts was initially detected by 10 h post-treatment and continued to increase through the 16 h time point. The kinetic pattern of the co-up-regulated genes clearly defined a class of early response genes that may have regulatory functions and delayed genes that may exhibit protective functions.

To examine directly the interactions of ER and osmotic stress on the co-up-regulated response, we analyzed whether the combination of AZC and PEG treatments promoted an additive increase in expression (Figures 5A and 5B). Five of the nine co-induced genes that we examined (asterisks) were induced by both stimuli in a more than additive fashion. Thus, the ER-stress and osmotic-stress signaling responses are integrated in a synergistic, convergent manner at the gene activation level. In this non-additive response, we also observed integration of the pathways by epistasis as PEG treatment reduced drastically the effect of AZC on the induction of the unknown gene (Figure 5A).

Plant growth and stress treatments

For the microarray experiments, soybean (Glycine max) seeds (cultivar Dare) were germinated in soil (MetroMix-360, Scotts, Marysville, OH) in a growth chamber with a day/night cycle of 9/15 h at 26°C/22°C. The aerial portions of three-week-old plants were excised below the cotyledons and directly placed into 15 ml of 10% (w/v) polyethylene glycol (PEG; MW 8000, Sigma, St. Louis, MO), 10 μg/ml tunicamycin (Sigma) or 50 mM L-azetidine-2-carboxylic acid (AZC, Sigma) solutions. The first trifoliate leaves were harvested after 16 h of PEG or AZC treatment (water control) and after 24 h of tunicamycin treatment (DMSO control, Sigma), then immediately frozen in liquid N₂ and stored at -80°C until use. In all experiments two independent biological replicates were used.

For the real-time RT-PCR experiments, soybean seeds (cultivar Conquista) were germinated in soil and grown in greenhouse conditions (avg. 21°C, max. 31°C, min. 15°C) under natural conditions of light, relative humidity 70%, and approximately equal day and night length. The first trifoliate leaves of three-week-old plants were excised and fed, via the petiole, solutions that induce the osmotic (10% PEG w/v) or ER stress responses (10 μg/ml tunicamycin or 50 mM AZC). After treatments for the times indicated in figure legends, the stressed trifoliate leaves and their untreated counterparts were immediately frozen in liquid N₂ and stored at -80°C until use. To avoid PEG interference in AZC uptake during the combined treatments we pre-treated the plants with AZC for six hours, when PEG was added for an additional 10 hours. Treatment with PEG and AZC simultaneously for 16 hours was found to give similar results. Each stress treatment and RNA extraction were replicated in three independent experiments.

Generation of soybean microarrays

The microarray slides consisted of 5,760 amplified cDNA fragments from soybean libraries prepared from RNA of developing seeds [86]. ESTs from these libraries were placed into contigs to identify unigenes [86], therefore a low redundancy in our set of clones is expected. The cloned cDNA fragments were amplified with M13 primers, purified using PCR Cleanup Filter Plates (Millipore, Bedford, MA), and eluted in water (according to the manufacturer's protocol). An aliquot of each amplified fragment reaction was separated through a 1% (w/v) agarose gel and visualized with ethidium bromide to assess size, quality and quantity. The purified PCR products were transferred to 384-well plates, and diluted with an equal volume of DMSO (Sigma). Finally, the PCR products were arrayed onto UltraGAPS slides (Corning, Corning, NY) using a 417 TM Arrayer (Affymetrix, Santa Clara, CA), cross-linked by exposure to UV light at 250 mJ and baked at 75°C for 2 h.

Isoforms of the ER stress-related molecular chaperone BiP were amplified with gene-specific primers (see Additional file 7) and included in the arrays. They consisted of three soybean isoforms A, C and D [39, 40].

RNA extraction and labeled cDNA preparation

Total RNA was extracted from frozen leaves with TRIzol (Invitrogen, Carlsbad, CA) according to the instructions from the manufacturer, and further purified through silica columns. The quality and integrity of the RNA was monitored by spectrophotometry and agarose gel electrophoresis, respectively. For the microarray hybridizations, 10 μg of total RNA were reverse-transcribed with the SuperScriptTM Indirect cDNA Labeling System (Invitrogen) in the presence of cyanine-3-dUTP (Cy3-dUTP) or cyanine-5-dUTP (Cy5-dUTP; Amersham Biosciences, Piscataway, NJ) according to the manufacturer's instructions. RNA samples from each biological replicate were labeled twice, once with each dye, to control for dye-specific effects on the hybridizations.

For the real-time RT-PCR, 2 μg of total RNA were treated with DNase (Promega, Madison, WI) and fractionated through RNA purification columns (Qiagen, Valencia, CA). Reverse transcription was carried out using M-MLV reverse transcriptase (Invitrogen) and oligo-dT (18, IDT, Coralville, IA) primers (according to the protocol of the manufacturer). Prior to the real-time RT-PCR assays, the quality of the cDNA was assessed by PCR with gene-specific primers for ubiquitin associated protein (UBA; AW508375) to test for genomic DNA contamination, as these primers amplify a larger fragment size from genomic DNA than from cDNA.

Microarray hybridization, scanning and data analysis

Soybean cDNA microrrays were subjected to a similar hybridization protocol as described by [87]. Briefly, the microarray slides were incubated for 45 min in 50 ml of a pre-hybridization buffer 5× SSC, 0.1% (w/v) SDS and 1% (w/v) of BSA (all from Sigma), washed sequentially in ultra-pure water and iso(2)-Propanol (Fisher, Waltham, MA), and air-dried. Slides were then incubated with Cy3- and Cy5-labeled cDNA (20 μl) from treated and control samples for 20 h at 42°C in a water bath protected from light. The hybridization buffer consisted of 0.5% SDS (w/v), 5× SSC, 5× Denhardt's, 50% (v/v) formamide, 0.5 μg/μl denatured calf thymus DNA (all from Sigma) and 0.5 μg/μl polyA RNA (Amersham Biosciences). Following incubation, slides were washed sequentially in three steps in the following solutions: (1) 1× SSC, 0.2% (w/v) SDS, (2) 0.1× SSC, 0.2% (w/v) SDS, and (3) 0.1× SSC.

Microarray slides were scanned at the Cy3 (530 nm) and Cy5 (650 nm) wavelengths with a ScanArray 4000 laser fluorescent scanner (Packard Bioscience, Perkin Elmer, Wellesley, MA), at a laser power of 100% and photomultiplier tube (PMT) gain of 75%. The image analysis and calculation of mean background-subtracted intensity of the spots were performed using QuantArray software version 2.2 (PerkinElmer). Normalization based on the LOWESS algorithm [88] and data analysis were performed using Genespring software version 7.2 (Agilent Technologies, Santa Clara, CA).

Genes were considered differentially expressed if they met both of two criteria. The first was an average fold change greater than 1.63 for PEG, 1.50 for TUN and 1.46 for AZC. These values were calculated based on the mean of gene expression ratio for each treatment plus two standard deviations (expression ratio values above 2 were not included for this calculation). This filtering criterion based on cut-off values of fold-change provided a pool of candidate genes differentially expressed. The second criterion, based on the t-test, was used to determine the statistical significance of the differences observed for the selected genes. The null hypothesis of the t-test [(mean of log2ratio)/SE] was rejected at 5% of probability. Additionally, genes with greater differences in expression (above 2 fold) in both biological replicates (but not necessarily in both technical replicates) were considered differentially expressed. A low stringency in our microarray data analysis was applied because relative expression of genes co-regulated by stress treatments was more accurately measured by real-time RT-PCR.

Annotations and Arabidopsis homologs of the soybean clones were assigned based on the top BlastX predictions against the GenBank [89] and The Arabidopsis Genome Initiative databases [90]. For most of the corresponding proteins there is not a demonstrated biochemical function; therefore, we refer to them as "like" proteins.

Note: The data discussed in this publication have been deposited in NCBIs Gene Expression Omnibus [91] and are accessible through GEO Series accession number GSE8992.

Real-time RT-PCR data analysis

For the quantitative RT-PCR assays, sequences of cDNA and primers are listed in Additional file 7. Analysis of expression of calnexin, an ER-stress responsive gene, and seed maturation protein PM30, a drought-induced gene, were used as positive controls for the respective stress treatments.

To select an endogenous control gene for data normalization in real-time RT-PCR analysis, we analyzed three genes encoding histone H2A, 60S ribosomal protein L30, and RNA helicase, which had been chosen because they had low and consistent expression ratios in the microarray results. The RNA helicase was used to normalize all values in the real-time RT-PCR assays, because it exhibited the lowest variation in expression values among treatments.

Real-time RT-PCR reactions were performed on an ABI7500 instrument (Applied Biosystems, Foster City, CA), using SYBR^® Green PCR Master Mix (Applied Biosystems). The amplification reactions were performed as follows: 2 min at 50°C, 10 min at 95°C, and 40 cycles of 94°C for 15 sec and 60°C for 1 min. To confirm quality and primer specificity, we verified the size of amplification products after electrophoresis through a 1.5% agarose gel, and analyzed the Tm (melting temperature) of amplification products in a dissociation curve, performed by the ABI7500 instrument.

Fold variation in gene expression was quantified using the comparative Ct method: 2^-(ΔCtTreatment – ΔCtControl), which is based on the comparison of expression of the target gene (normalized to the endogenous control) between experimental and control samples.