Microarray analysis and scale-free gene networks identify candidate regulators in drought-stressed roots of loblolly pine (P. taeda L.)

Experimental Design

Root samples were collected from 10-month old ramets, (i.e., rooted cuttings) or propagules having root structures essentially indistinguishable from those of pine seedlings. Propagules representing four genotypes (CCLONES 41201, 41369, 44686, and 45226) obtained from full-sib crosses of elite commercial germplasm were grown and harvested by the Forest Biology Research Cooperative (http://fbrc.ifas.ufl.edu) as described [32]. Established ramets were planted in plastic pots (10 cm diameter × 36 cm deep) containing 100% acid-washed fine sand, and maintained in natural daylight in a greenhouse in Gainesville, FL.

The study employed a randomized complete block, split-plot design, with genotype as the whole-plot factor and watering treatment as the split-plot factor. The drought-stress treatment (DS) withheld water until water potential in pre-dawn needles reached -1.75 MPa and needles showed visible wilting, which occurred within seven treatment days. The drought-recovery treatment (DR) withheld water until water potential in pre-dawn needles reached -1.75 MPa, followed by watering to pot capacity on day seven and a subsequent drought recovery period of 48 hours. The control treatment (well-watered, WW) propagules were watered to pot capacity every other day during the treatment period and harvested at the same time as the DS samples. Mean water potential in pre-dawn needles on the WW plants remained at -0.3 MPa ± 0.1. Tissue samples from the DR and WW plants were harvested concomitantly.

Gene expression was evaluated in loblolly roots following one of three treatments (DS, DR, or WW). Four biological replicates were investigated, with each replicate comprised of RNA pooled from three individuals (Figure 1). Labeled targets from each cDNA synthesis were hybridized twice to monitor technical variation. Microarray experiments employed a reference design [39] to identify genes that exhibited treatment-specific differences in loblolly pine roots; the reference sample was comprised of total RNA from P. taeda roots, shoots and needles (Additional File 1).

Quantitative PCR (qPCR)

Primer pairs were validated using a dilution series of P. taeda genomic DNA (80 pg - 50 ng) (Additional File 2). Validation reactions were conducted in duplicate and contained 5 uL genomic DNA, 5 uL primer pairs (0.6 uM final conc.), and 10 uL iQ™ SYBR^® Green Supermix (Bio-Rad Cat#170-8882). Cycling parameters and selection of primers for validation have been published previously [40]. cDNA was synthesized using Invitrogen Super Script III First Strand Synthesis System for RT-PCR (Cat No. 18080-051). Optimal cDNA concentration was determined empirically using primer pairs made for the actin ACT2 gene. The same gene was used to normalize all qPCR data. Reactions were performed in triplicate and averaged. Samples with threshold cycle (Ct) values ≥ 35 or those where the standard deviation across a triplicate was ≥ 0.5 were considered failed reactions and re-run.

The PtGen2 Microarray

cDNA sequences for the PtGen2 microarray were derived from thirty-six cDNA libraries constructed primarily from loblolly root, stem and needle tissues [32]. A total of 25,848 (i.e., 26,496 total spots minus buffer blanks and duplicate spots) cDNA clones were amplified using 0.4 uM pSLFOR2 primer (5'-AAAACGACGGCCAGTGCCAAG-3') and 0.4 uM pSLREV2 primer (5'-GCTTCCGGCTGCTATGTTGTGTGG-3'). Reactions were incubated at 95°C for 3 min, followed by 30 amplification cycles (95°C for 30 sec, 58°C for 1 min, and 72°C for 2.5 min) and then incubated at 72°C for 7 min. All reaction products were verified by agarose gel electrophoresis. cDNAs were re-suspended in spotting buffer (150 mM PO₄ 1.5M betaine, 0.001% SDS) and printed on Corning UltraGAPS™ slides (Corning Inc., Corning, NY) using a BioRobotics Microgrid II spotter (Genomic Solutions, Ann Arbor, MI). The 384 probes in print plate #10 were re-spotted as plate #69, for use in monitoring technical variation across the array. Microarrays were printed by the Vanderbilt Microarray Shared Resource facility.

Loblolly pine 3' EST sequences have been assembled previously [32, 35]. Contig consensus sequence or EST singleton sequence was used as the BlastX query to databases including NCBI nr, UniProtSProt, UniProtTrembl, UniRef, TAIR, Medicago truncutula Database (MT3.0), and Oryza sativa (OSA). Returned annotations with E-values > 1 × 10^-5 were ignored. A master annotation file for identification of PtGen2 probes has been included (Additional File 3). Blast2GO was used to identify gene ontology (GO) terms associated with differentially expressed genes [41].

cDNA Synthesis, Labeling, and Hybridization

Root samples were pulverized under liquid nitrogen using a SPEX model 6850 freezer mill (SPEX, Metuchen, NJ). Total RNA was extracted from sample tissues using a modified protocol of Chang et al. [42, 43] and RNA samples were DNase-treated using the Ambion TURBO DNA-free™ Kit (Applied Biosystems Inc., Foster City, CA). RNA concentration was determined spectrophotometrically and agarose gel electrophoresis was used to ascertain RNA quality and integrity.

Synthesis of amino-allyl and amino-hexyl modified cDNA used the Invitrogen SuperScript™ Indirect cDNA Labeling System kit (Invitrogen Corp., Carlsbad, CA). Reactions contained 20 ug of total RNA, 2 uL anchored oligo(dT)₂₀ primer, 1 uL random hexamer primer, and were carried out at 45°C for 12-14 hr. Modified cDNA was purified on S.N.A.P™ columns (Invitrogen Corp., Carlsbad, CA), precipitated, and re-suspended in warm 2X coupling buffer. Cy-5 and Cy-3 (GE Healthcare Bio-Sciences Corp., Piscataway, NJ) were suspended in 80 uL DMSO and coupled to the modified cDNAs, as per the manufacturer's instructions. Cy-labeled cDNAs were purified using published methods [44], after which synthesis and labeling efficiency were monitored spectrophotometrically.

Microarray processing and hybridization protocols have been described previously [36, 45]. Slides were incubated in a HybChamber™ (Genomic Solutions, Ann Arbor, MI) containing 20 uL 100 mM DTT in the humidity wells. The entire apparatus was wrapped in foil and incubated at 48°C for 14-16 hr.

Data Analysis

Microarrays were scanned using a ProScanArray™ confocal scanner (Perkin Elmer, Waltham, MA) equipped with 532 nm and 635 nm lasers. Raw fluorescence data were processed using ImaGene, ver. 7.5 (Bio-Discovery Inc., El Segundo, CA, USA). Signal means were determined without background correction [46]. Data filtration, log₂ transformation, normalization (print-tip Lowess), and statistical analyses were completed using BRB-Array Tools [47]. Duplicated probes on the array were treated independently during normalization and statistical analyses. Spot signals were filtered as described previously [45, 48] and unigenes with ≥ 50% missing data were not included in the analysis. When necessary, missing values were imputed using the KNNimpute algorithm (K value = 15) [49]. Differentially expressed genes were identified from paired comparisons of the three treatment groups (DS vs. WW, DR vs. WW, DR vs. DS) using a random variance test [47]. The false discovery rate (FDR) was set to 0.01 and a cut-off minimum of 1.5-fold was used. Differentially expressed genes identified in the DS vs. WW comparison were used to extract the normalized log₂ mean ratio data identified from drought-recovery samples. Weighted gene co-expression network analysis was performed as described previously [37, 38] using an R-script modified for this analysis and a network threshold cut off of 0.01 (Additional File 4). Cytoscape Ver. 2.6.3 [50] and NetworkAnalyzer ver. 1.0 [51] were used to visualize the global gene network and calculate topology parameters. Candidate network drivers were identified by ranking nodes based on their degree index and also by comparing to index values determined for closeness, radiality and eccentricity.

Sequence and Microarray Data Accession

Contig and EST sequence data utilized in probe synthesis can be accessed at the Fungen website [35]. Microarray data can be accessed at the NCBI Gene Expression Omnibus (GEO) [52] under accession GPL11184.

Array Performance Validation

Slide-to-slide reproducibility was verified via a hierarchical cluster analysis (centered correlation, average linkage) of all datasets from each treatment group (Figure 2) which showed a clear separation of clustering of the drought-stressed (DS) samples from well-watered (WW) and drought-recovered (DR) samples. Feature-to-feature reproducibility was verified via a correlation analysis of filtered log₂ ratios for 384 probes that are replicated on the array. The average r-values determined for each sample (within treatments) ranged from r = 0.68 to r = 0.83, with an overall mean correlation of r = 0.78 (Additional File 5).

Quantitative PCR (qPCR) was used to validate treatment-specific expression patterns of 16 differentially expressed (DE) genes, eight of which were the most highly up- or down-regulated transcripts identified in this study, and eight of which were moderately expressed transcription factors (Additional File 6). All tested genes exhibited expression patterns that were consistent across the two platforms; however, PtGen2 underestimated the expression ratios in some cases by more than 10-fold when compared to qPCR values. Signal compression of fold-change ratios is a common phenomenon that has been demonstrated with both oligo- and cDNA microarrays [17, 53] and has been attributed to the presence of cross-hybridization signals [54].

Effect of Drought Stress on the P. taeda Root Transcriptome

Comparison of drought-stressed and well-watered roots identified 2445 DE genes (Additional File 7), 1670 of which returned annotation results from BlastX and Blast2GO (E-value cut off ≤ 1 × 10^-5). In addition to Gene Ontology (GO) terms associated with cellular and primary metabolic processes, the third most abundant GO category for the 1848 DE genes that mapped to at least one term was stress response (Figure 3). Three other relevant GO categories were associated with responses to chemical, abiotic, and biotic stresses. Table 1 lists the twenty-five drought stress-induced genes having the highest log₂ expression ratios. Genes thought to encode a Class IV chitinase, late embryogenesis abundant proteins (LEAs), and dehydrins were among the most highly up-regulated genes; while seven genes had no apparent homology to any known proteins. Many pine genes with apparent homology to other genes associated previously with plant stress responses in angiosperms and gymnosperms were up-regulated in drought-stressed roots, including those thought to encode chitinase, dehydrin, defensin, chaperones/HSPs, water-deficit-inducible LP3, galactinol synthase (GolS), glutathione peroxidase, thaumatin, sucrose synthase (SuS), and β-glucosidase (Additional File 7). In particular, 16 DE genes thought to encode LEAs were detected and 11 were up-regulated following drought stress, as were all 15 genes identified as dehydrins. LEAs and dehydrins were also identified as hubs for the predicted global gene network (Table 2). Genes that may be involved in the production of osmolytes were also up-regulated. For example, all DE genes identified as SuS (6 probes), GolS (11 probes), or raffinose synthase (2 probes) were up-regulated after drought-stress, as were genes thought to be involved in proline biosynthesis, e.g. pyrroline-5-carboxylate synthetase (2.2.16.16).

Up-regulated DE genes thought to encode regulatory proteins included ABA-responsive protein orthologs, aminocyclopropane (ACC) oxidases and synthase, and several transcription factors. Transcription factors (TFs) showing increased transcript abundance after drought-stress included genes with notable homology to DREBs (3.2.8.7, 4.4.14.24) and ERFs (7.3.7.6, 8.1.15.4), as well as TFs in the WRKY (7.3.15.5), Dof (2.4.15.5), bZIP (9.1.19.7), MYB (3.1.19.21), and NAC (6.2.15.1) families. Several of these are thought to be central players in signal transduction pathways involved in both drought and abiotic stress responses [1, 3, 55, 56].

Of the 2445 DE genes, 1251 were down-regulated in DS roots. These included genes thought to encode auxin-responsive proteins, and components of cell wall and carbohydrate metabolism (Additional File 7). Numerous genes encoding putative aquaporins were down-regulated after drought stress, as were genes with homology to histones, dirigent-like proteins, cellulose synthases (CesA), expansins, pectate lyases, senescence-associated proteins, subtilisin-like and aspartic proteinases, and xyloglucan endotransglycosidases. For example, of the 18 DE genes thought to encode aquaporins (AQPs), most fell into the plasma membrane intrinsic protein (PIP) subgroup and all were down-regulated, except an Arabidopsis PIP2-6 ortholog (9.3.16.18). Only seventeen TFs were down-regulated and only two of those, AP2/ERF (2.2.22.4) and WRKY (1.2.18.1), have been previously associated with regulatory roles in drought response.

Weighted Gene Correlation Network Analysis (WGCNA) [37, 38] was employed to mine the 2445 DE genes for candidate regulatory genes. As would be expected for a biological system, our data predict a gene co-expression network that displays scale-free topology with inherent modular features (Figure 4A). Eleven co-expression modules (with a minimum of 30 genes each) were identified (Additional File 8). Evaluation of expression patterns from individual network modules revealed groups ranging from 24 to 694 genes that displayed coordinated changes in their transcript populations during drought stress and drought recovery (Figure 4B). As shown in Figure 5, node degree distribution for the pine gene network approximates a power-law distribution, which is a characteristic of scale-free networks [57]. Network statistics demonstrate that the pine root gene network has properties found in other biological networks (Table 2). For example, 1477 nodes were connected by 49,997 edges, representing 2.3% of the total edges possible in the network. Sparsely connected nodes in highly clustered areas are a common characteristic of the hierarchal architecture seen in many biological networks [58]. The average clustering coefficient as well as values determined for the shortest path length and smallest network diameter, are characteristic of network compactness and indicative of "small-world" networks that display high clustering and short path lengths [59].

Transcriptome Responses in Drought-Recovered Roots

Comparison of DR roots with WW roots identified only 76 DE genes (Additional File 7), whereas comparison of DS roots with DR roots identified 1918 DE genes (Additional File 7). More than 70% of the genes in this latter set were also identified in the comparison of DS and WW roots (Figure 6). The data shown in Figure 2 is consistent with the observation that the transcript populations in drought recovered and well-watered roots are far more similar than either is to the transcript population in drought-stressed roots. Taken together, these results suggest that the root transcriptome in P. taeda seedlings can return to its pre-drought state within 48 hours of recovery from the applied stress. Moreover, the drought-recovered plants (at harvest) no longer displayed the wilting phenotype seen in drought stressed plants and DR root tips had changed from the totally brown appearance seen in DS roots to having copious numbers of white root tips characteristic of control roots (data not shown).

Sixty-four of the 76 DE genes identified in the drought recovered samples were up-regulated and 30 of these had little sequence homology to known genes (E-value > 1 × 10^-5). Significant overlap between drought-recovered and drought-stressed expression patterns was seen as 71 of 76 genes identified in DR roots were also identified in DS roots (Figure 6). These included multiple genes encoding apparent LEAs, chitinases, and β-glucosidases (Additional File 7). Possible regulatory genes identified in the drought-recovered sample include a two-component histidine-kinase response regulator protein (11.1.17.15), an ethylene-responsive RNA helicase (8.2.16.16), a gibberellin-regulated protein (10.3.10.14), and TFs belonging to the WRKY (7.3.15.5) and zinc finger (9.3.15.22) classes. With the exception of the gibberellin-regulated protein, putative regulatory genes identified in drought-recovered roots were also found in the drought-stressed roots where they all demonstrated increased expression. Few genes showed reduced expression after drought recovery and only one, an α-pinene synthase, showed expression reduced > 2-fold. Among the genes down-regulated in drought-recovered roots were ones thought to encode pinene synthases, a pectate lyase, a xyloglucan endo-transglycosylase, and a β-glucosidase. Only five DE genes were identified exclusively in drought recovered roots. These appeared to encode a ribosomal protein (11.3.20.24), a xyloglucan endo-transglycosylase (11.2.17.17), a gibberellin-regulated protein (10.3.10.14), and two unknowns.

PtGen2 Microarray

The loblolly pine microarray (PtGen2) that is described here contains 26,496 features, including 25,848 distinct cDNAs, making it the largest publicly available microarray for any conifer species. PtGen2 can also be used for gene expression studies in species other than P. taeda. For example, it has been used to evaluate transcriptome dynamics in staged embryos from P. pinaster[62], and to compare expression profiles from several closely related pine and conifer species (data not shown). The data presented here demonstrate the utility of PtGen2 in the evaluation of transcriptome responses in loblolly pine and identifies classes of genes whose expression profiles are consistent with results from previous drought studies in conifers and angiosperms. The cDNA sequences printed on PtGen2, as well as the contigs to which they belong, can be accessed at the Fungen website [35]. PtGen2 chips, detailed protocols, complete annotation files, and the .gal file can be obtained by contacting the authors.

Molecular Basis of Drought Tolerance in Pine Roots

Root transcriptome responses to drought stress were characterized by 9.6% of the genes interrogated with PtGen2 being differentially expressed by at least 1.5-fold. Forty-nine percent of DE genes were up-regulated and 51% were down-regulated, indicating that down-regulation of transcripts is likely a non-trivial aspect of the drought response in P. taeda roots. A similar percentage of DE genes (12.5%) were also observed in P. taeda needles after cycles of mild and severe drought stress, and an approximately equal number of transcripts were both up- and down-regulated during the early stages of mild and severe drought stress [29]. Analyses of transcriptome responses to water stress in angiosperms, e.g. Arabidopsis[63, 64], rice [65], barley [66, 67], and grape [68] identified 4-16% of the genes interrogated as being differentially expressed. This is in stark contrast to Populus subjected to gradual drought stress where only 1.5% of all genes and < 1% of root-responsive genes were differentially expressed at a ≥ 2-fold level [26]. Using a 2-fold cut-off for our data, the percentage of DE genes in drought-stressed pine roots is reduced to 3.5% which is still significantly higher than the Populus value. A complete explanation for this difference remains to be elucidated.

The complexity of transcriptome responses to drought stress makes comparisons among the number and type of DE genes identified in previous studies difficult as some differences are no doubt due to the duration and severity of stress, the tissue under study, the experimental design, and the statistical treatment of data. Nevertheless, the differential expression patterns seen for specific classes of genes in this study were generally consistent with what has been documented in other plant species subjected to severe water or salt stress. Many of the drought stress responsive genes found in this study, whether identified by the number of probes representing them, their log₂ values, or by gene network analysis, have previously been well documented in angiosperm drought responses. This suggests evolutionary conservation of the molecular mechanisms utilized by both gymnosperms and angiosperms for response to water deficits. In drought-stressed roots, the majority of DE genes displaying increased transcript abundance fall into the tolerance or protective category of drought-responsive genes, probably because the root samples were so severely stressed at the time of harvest that genes involved in early sensory and regulatory responses were underrepresented.

Numerous studies have addressed the physiological recovery from drought stress in pines and conifers [69–72], but very little information is available on transcriptome responses during drought recovery. What is clear from this study is that the transcript population in drought-stressed loblolly roots can return to nearly pre-drought transcript population levels within 48 hours of removing the stress (Figures 2 and 6). Of 76 DE genes identified in DR roots, 93% were found to be expressed at elevated levels in DS roots. One gene that stood out in drought recovery was a gibberellic acid-binding protein. However, we did not see differential expression in any gibberellic acid 2-oxidase genes, even though these genes have frequently been reported as induced by stress conditions, including drought. Gibberellic acid signaling could have relevance for drought tolerance through its effect on plant morphology. Future studies that include multiple time points to more precisely pinpoint when drought recovery begins and how long it takes to complete would help identification of transcriptome changes that occur early in the recovery process.

Candidate Protection Genes

A large gene family commonly seen in drought response studies is the late embryogenesis abundant (LEA) proteins, which includes dehydrins [73]. Their precise function(s) is unknown, but they are believed to be involved in a number of protective processes, including acting as hydration buffers, antioxidants, metal ion binding proteins, and both enzyme and membrane stabilizers (reviewed) [74]. Studies have shown the induction of LEA and dehydrin transcripts, as well as increased production of the encoded proteins, in response to osmotic stress [74, 75]. LEA gene expression was previously shown to be up-regulated in pine shoot and stem tissues responding to drought stress [28, 29]. Twenty-six of the 31 probes on the PtGen2 microarray that appear to encode LEAs or dehydrins showed their genes to be up-regulated in DS roots. Three of these (3.1.11.9, 2.4.16.12, and 1.3.16.12) were highly induced and found to be highly ranked in the gene network analysis. This supports the belief that LEAs play an important role in mitigating the effects of drought stress and also demonstrates the identification of an expected class of drought-responsive genes by PtGen2. These results also suggest that involvement LEA gene function in drought tolerance is conserved among seed plants.

Compatible solutes (osmolytes) are small solutes (primarily sugars and amino acids, such as proline) that are thought to protect cells by maintaining turgor pressure and stabilizing proteins and membranes under osmotic-stress conditions [76, 77]. Sugar-based osmolytes, such as sucrose and raffinose family oligosaccharides (RFOs), are thought to be particularly important for reducing the effects of osmotic stress in plants [77, 78]. DE genes thought to encode sucrose synthase, galactinol synthase, and raffinose synthase were up-regulated in drought stressed roots. Genes thought to be involved in proline synthesis also were up-regulated, while almost all genes identified as "proline-rich protein" were down-regulated.

Desiccation-tolerant plants accumulate high concentrations of sucrose during dehydration which, in conjunction with LEAs and dehydrins, may act to stabilize drying cells through direct interactions with macromolecules and membranes [78]. In maize ear and tassel tissues, Sus transcripts increased during water deficit stress [79]. In A. scabra, a thermotolerant grass species, increased accumulation of sucrose synthase in the roots was proposed to contribute to superior root thermotolerance by regulating sucrose metabolism [80]. Analysis of differentially expressed ESTs has shown increased Sus transcript levels as well as elevated of sucrose levels in water-stressed citrus seedlings [81]. Interestingly, Sus transcript levels were induced in water-stressed Populus leaves, but were consistently repressed in roots during the early stages of water deficit [26]. Up-regulation of Sus gene transcripts in water-stressed pine roots may be related to either the generation of sucrose or to the regulation of hexose flux. For example, both white spruce (P. glauca) and jack pine (P. banksiana) seedlings accumulated higher concentrations of sucrose in roots relative to stems after exposure to elevated salt concentrations [82], while in Arabidopsis, sucrose synthase was shown to control intracellular glucose levels and transitory starch biosynthesis [83]. Similar alterations in sugar metabolism likely occur in pine roots as all genes on the PtGen2 microarray encoding putative invertases and hexokinases were up-regulated under drought stress. One suggested regulatory mechanism for SuS is that its role in sucrose cycling impacts flux through the hexokinase reaction [84] Meanwhile, both hexokinases and cell wall invertases are thought to play key regulatory roles in assimilate partitioning during stress [85, 86].

Up-regulation of SuS genes could also be related to raffinose family oligosaccharide (RFO) biosynthesis. Galactinol synthase catalyzes the first committed step in RFO synthesis, utilizing sucrose and galactinol as substrates and, thus, plays a regulatory role in carbon partitioning [87]. In Arabidopsis, it was shown that GolS1 is responsible for heat stress-dependent synthesis of raffinose [88], while both GolS1 and GolS2 are up-regulated in response to drought treatment. Over-expression of GolS2 in transgenic Arabidopsis leads to improved drought tolerance [89]. It also has been suggested that the osmoprotectant ability of galactinol and raffinose is due to their participation in ROS scavenging, because both can competitively inhibit hydroxyl radical-induced formation of 2,3-dihydroxy-benzoic acid (DHBA) in vitro[90].

An inositol transporter gene (3.2.1.7) was highly ranked in the network analysis, which could suggest a role in osmoprotection. Myo-inositol is a substrate for conjugation with galactose to form galactinol. Accumulation of myo-inositol and its methylated derivatives is correlated with drought tolerance [91], and it has been hypothesized that increased transport of inositols and its derivatives may act to counter osmotic stress [92]. The large number of up-regulated DE genes related to the biosynthesis and regulation of osmolytes in conjunction with an inositol transporter hub gene strongly suggests an important role for osmoprotectants in the response of pine roots to drought stress. Physiological and biochemical studies will be needed to address this hypothesis more thoroughly.

Aquaporins (AQPs) are members of the major intrinsic protein (MIP) family involved in increasing water diffusion across membranes by raising hydraulic or osmotic permeability [93, 94]. Disparate expression patterns for these proteins in response to water deficit is common and this, in conjunction with the large number of AQP genes typically found in plants as well as their highly variable tissue-specific expression patterns, have confounded attempts to determine their specific roles in drought response [28, 95–98]. Transgenic plants that over-express a class of AQPs called plasma membrane intrinsic proteins (PIPs) have shown both increased and decreased drought tolerance [99, 100]. In this study, 17 of 18 probes for pine AQPs responded to drought stress. Most fell into the PIP category, and 16 were down-regulated (one Arabidopsis PIP2-6 ortholog was up-regulated). In Arabidopsis subjected to gradual drought stress, AQPs were generally down-regulated [101], as was the case for PIPs identified in tobacco [102] and Populus roots [26]. Thus, down-regulation of AQPs in drought-stressed pine roots appears to follow a pattern similar to that seen in some angiosperms and may reflect a mechanism for water conservation via reduced membrane permeability that minimizes water efflux into the surrounding soil [95].

Candidate Regulatory Genes

Apparent transcription factors (TFs) from classes known to play a role in drought and abiotic stress responses in plants via ABA-dependent and ABA-independent mechanisms were differentially expressed in water-stressed pine roots. Genes encoding putative DREB1, bZIP, AP2/ERF, MYB, NAC, and WRKY TFs were all up-regulated. An apparent WRKY TF displayed the highest log₂ ratio of any TF identified in DS roots. Interestingly, one of the most common drought response regulatory genes (DREB2) was not seen among the DE genes in either DS or DR roots. On the other hand, DREB1 TFs are more often associated with cold stress in plants [1]. Whether the apparent DREB1 identified in water-stressed pine roots plays a similar role to DREB2 in other plants remains to be determined.

Correlation network analysis did not identify any of the up-regulated TFs in water-stressed pine roots as hub nodes. However, network analysis did identify three up-regulated DE genes that may have a regulatory role in pine root responses to drought stress, cardiolipin synthase/phosphatidyltransferase, 9-cis-epoxycarotenoid dioxygenase, and zeatin o-glucosyltransferase, which could influence phospholipid metabolism, ABA synthesis, and cytokinin availability, respectively. The most likely explanation for why "classic" regulatory TFs, such as DREB2 and bZIP, were not seen as highly ranked network hubs in this study is a matter of experimental design. By the time water-stressed roots were harvested they were so severely stressed that transcription network dynamics may have progressed to a point where genes were more likely to be involved in stress mediation than to response of the onset of water deficiency. This would explain why the majority of up-regulated hubs listed in Table 3 identified genes encoding proteins involved in protection from ROS-mediated cell death (thioredoxin, peroxidase, metallothionein), transport (inositol transporter 2), cell/macromolecule protection (LEA, dehydrin, HSP70), and defense (NBS-LRR, phenylcoumaran benzylic ether reductase).

We identified an apparent WRKY gene up-regulated in both drought stress and drought recovery that is a candidate transcription factor for the regulation of drought response in pine roots. WRKYs comprise a large family of conserved plant TFs that are up-regulated in response to wounding, infection, and abiotic stress and that play multiple roles in developmental and metabolic pathways [103, 104]. In Arabidopsis and rice, co-regulatory networks have been identified for WRKY genes involved in both defense and osmotic stress responses [105]. WRKY TFs induced by drought stress have been identified in rice and pennycress [106–108], while an Arabidopsis WRKY TF has been shown to play a critical role in both ABA response and drought tolerance [109]. Drought tolerance has been conferred on transgenic Arabidopsis expressing rice and soybean WRKY genes [107, 110]. Recently, a WRKY TF was shown to participate in dehydration tolerance by binding to W-box elements of the galactinol synthase promoter in Boea hygrometrica[111]. Given that all galactinol synthase genes showed increased transcript abundance in drought-stressed samples, we can speculate that their expression in pine roots may reflect regulation by a WRKY TF.

ABA is considered something of a universal stress hormone in higher plants and it plays a central role in the regulation of responses to abiotic stress [112]. Although the primary site for ABA synthesis is generally considered to be the roots (from which it is released to the xylem for translocation to the rest of the plant), recent studies have suggested that ABA synthesis may also occur in shoot tissues [113]. Both ABA-dependent and -independent signaling pathways have been described in stress responses to drought, salinity and cold [1]. The gene 9-cis-epoxycarotenoid dioxygenase (NCED) catalyzes the first dedicated step and is the key enzyme in ABA biosynthesis [114]. In transgenic Arabidopsis, an NCED gene was shown to increase endogenous ABA levels resulting in enhanced drought tolerance [115]. In legumes and shrub species undergoing water-deficit induced stress, elevated NCED gene expression has been correlated with concomitant increases in ABA [116, 117]. A gene thought to encode NCED was up-regulated in drought stressed pine roots and ranked high as a network hub. A second hub gene related to ABA signaling (ABA-responsive protein) was also identified from the network analysis. These observations are consistent with the regulatory role that NCEDs are thought to play in response to water stress in angiosperms and makes NCED a logical choice for having a similar function in pine roots.

Other translocated factors, such as cytokinins, may also originate in roots and modulate ABA signaling to impact shoot responses to drought [6, 10]. The importance of cytokinins in drought-induced changes of root morphology has been suggested based on their accumulation in drought-stressed tobacco roots [118]. In this study, a gene encoding an apparent zeatin-O-glucosyl transferase (ZOG) was strongly up-regulated and was identified as a potential network regulator, while a second ZOG and a cytokinin oxidase displayed moderate increases in transcript abundance. Over-expression of ZOG 1 in transgenic tobacco resulted in increased total cytokinin content, although bioactive cytokinins did not change significantly [119]. O-glucosylated forms of trans-zeatin are protected from oxidases and dehydrogenases and may act as cytokinin storage forms that can be converted quickly to bioactive cytokinins by the action of β-glucosidases [120]. Reducing bioactive cytokinin availability via glucosylation may be a regulatory mechanism for lowering the cytokinin:ABA ratio in pine roots to lessen the stimulatory effect of cytokinins on root meristematic tissue.

Correlation network analysis is a powerful method for mining complex "omics" datasets and the evaluation of gene networks often leads to identification of novel genes involved in particular processes or pathways, as well as identification of putative regulatory genes. In plants, this approach has seen increasing use, particularly in model species. For example, a genome-wide co-expression network in Arabidopsis identified 127 functional modules, fourteen of which were associated with stress and defense [121], and correlation networks in Arabidopsis[122] and rice [123] have identified co-expression modules that allude to the involvement of DREB genes in the regulation of trehalose-6-phosphate. However, few studies prior to this have specifically targeted drought responses.

In the current study, a variety of apparent transcription factors of classes known to regulate drought and water-deficit responses in other plants were found among DE genes in pine roots; however, none of these ranked highly in the network analysis. The most highly ranked network hub gene, although it is of unknown function, is of great interest. Other highly ranked hub genes appear related to protective functions rather than regulatory roles. Nevertheless, the apparent NCED and ZOG genes identified in this study have not been reported previously to respond to drought stress in pines. The crucial role these genes likely play in the ABA and cytokinin signaling pathways, respectively, should make them foci for future molecular investigations of water-stress responses in roots.