Soybean (William 82 genotype) plants were grown in 26.5 L pots with a dimension of 30 × 27 × 37 cm (top diameter × bottom diameter × height), filled with a mixture of turface and sand (2:1). Osmocote (Scotts Co., Marysville, OH) with a nutrient content of N:P₂O₅:K₂O = 14:14:14 was added as a nutrient source at a rate of 20 g per pot. The plants were grown in a well-watered conditions [up to vegetative (V4) stage] under a 16 h photoperiod (~10,200 lux) at 28°C. Drought stress was imposed by withdrawing water for 7 days, while flooding stress was imposed by placing the pots into a bigger pot with a trashcan liner filled up to a water level of 4 cm above the soil surface for 7 days. After 7 days of treatment, all the leaves (irrespective of their response to stress) were sampled from the drought and flood treatments and control plants for RNA sequencing and real-time PCR analysis. The soil moisture was measured after the stress imposition in the drought experiment using a PR2 moisture probe (Delta T, UK).

The total RNA was extracted using RNeasy Plant Mini Kit (Qiagen, Valencia, CA) and on-column DNase digestion was performed according to manufacturer's protocol. The concentration and quality of RNA were checked by a Nanodrop 1000 (Grand Island, NY) and by running the product on an agarose gel. The RNA-Seq libraries were prepared using a TruSeq Stranded mRNA LT Sample Prep Kit—Set A (Catalog #: RS-122-2101, Illumina, USA) according to the manufacturer's instructions. Three biological replicates were used for RNA-Seq and real-time PCR analysis.

The single end RNA-Seq reads were generated on the Illumina Genome Analyzer (San Diego, CA) platform. The processing of the initial reads was performed using the Illumina analysis pipeline (in the Fastq format) using custom and default parameters. Additional filtering was performed by removing adaptor sequences and low quality bases. To facilitate the read mapping, the Glycine max reference genome (Gmax1.1version) was indexed by Bowtie (http://www.phytozome.net; Langmead and Salzberg, 2012). The read mapping was performed using the Tophat software package (Trapnell et al., 2009; Kim et al., 2013). The reads were first mapped directly to the genome using indexing and then some of the unmapped reads were resolved by identifying novel splicing events. Two mismatched base pairs were allowed and the multiple position matching was reported up to 40 alignments using the Tophat mapping procedure. The transcriptome raw sequencing data from this study have been submitted on the NCBI (http://www.ncbi.nlm.nih.gov/) database as individual BioProjects: PRJNA324522.

The binary read alignment files were used as input to Cufflinks (Trapnell et al., 2009), which assembled the reads into transfrags (transcripts). The estimated gene abundance was then measured in terms of the fragments per kilobase of transcript per million mapped reads (FPKM). The differentially expressed genes (DEGs) between the two sets of samples were identified using cuffdiff. The significant up-regulated and down-regulated gene lists were obtained for the drought and flood samples, respectively. Only the genes with a log₂ fold change ≥+2 and ≤ −2, but without infinite values and a FDR adjusted p ≤ 0.05 after Benjamini-Hochberg correction for multiple-testing with significance level “yes,” were considered as significantly DEGs.

The DEGs were annotated for gene ontology (GO) terms (Ashburner et al., 2000) and categorized into Molecular Function (MF), Cellular Component (CC), and Biological Process (BP) categories. A gene enrichment test was then performed on each of the gene lists to obtain the significant terms. Fisher's exact test, which is based on the hypergeometric distribution, was used to calculate the p-value. The TreeView (http://jtreeview.sourceforge.net/docs/overview.html) program was used to draw the heatmap of the significant DEGs in response to drought and flood stress.

The drought and flood responsive pathways were identified and plotted using MapMan (Usadel et al., 2009). Multiple biological or metabolic pathways were plotted together with the mapped gene intensity of the fold change (≥ +2 and ≤ −2) by a blue and red schema. The TFs and transcription related genes were also mapped and plotted by MapMan. A more detailed TF family annotation was obtained from plantTFDB (Jin et al., 2014).

Using the collection of the motif sequences from the “Database of Plant Cis-acting Regulatory DNA Elements” (Higo et al., 1999), the 3 kb upstream of the 5′ translation start base was searched for all of the annotated soybean genes (phytozome v9.1). The number of genes was counted for each matched motif (gene count) and the total number of matches to all of the genes by each motif (hit count). Furthermore, counts from the subset of genes in each of the gene sets with all the genes were compared using Fisher's test. Any of the comparisons with a p ≤ 0.05 is presented in Supplementary Tables 3A,B.

Total RNA was extracted from each sample using the Qiagen RNeasy mini kit (Qiagen, CA, USA).

The first strand cDNA from 1 μg of total RNA was synthesized using the EcoDry premix (Clontech, CA, USA), following the manufacturer's instructions. Quantitative PCR (qPCR) was performed using 10-fold diluted cDNA product in a 10 μL reaction volume using the Maxima SYBR Green/ROX qPCR Master Mix (Thermo, Waltham, MA, USA) on ABI7900HT detection system (Foster City, CA, USA). Three biological replicates and two technical replicates were used for analysis. The PCR was performed using two-step cycling protocol as follows: 50°C for 2 min; 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, and 60°C for 1 min (https://www.thermofisher.com/order/catalog/product/K0221). To normalize the gene expression, Actin (Glyma18g52780) was selected as a housekeeping gene. All the primers were designed using QuantPrime software (http://www.quantprime.de) (Supplementary Table 5).

Three biological replications were subjected to well-watered, drought, and flooding conditions at the V4 growth stage. Soil moisture, the major limiting factor used to impose drought stress on growing plants (Else et al., 1995), and was measured before and after the drought treatment. The soil moisture of the well-watered plant was 19% and was reduced to 10% after 7 days of drought treatment; the plants showed symptoms of leaf wilting. The degree of drought and flooding was determined by monitoring the canopy temperature. The canopy temperature (°F) increased from 76.8 ± 1.03 in normal conditions to 85.7 ± 4.42 and 79.4 ± 1.06 under the drought and flooding conditions, respectively.

The RNA samples from the soybean leaves under the control/drought/flooding conditions were used for sequencing by the Illumina Genome Analyzer. Approximately 28 million reads were generated from each sample (Table 1). The RNA-Seq analysis workflow is shown in Supplementary Figure 1 and was utilized for the data analysis. Approximately 72% of the reads were mapped to the soybean reference genome and ~2.5 to 4 million reads were mapped to multiple regions (Table 1).

To quantify and identify the genes, a core set of DEGs under drought and flood stress in soybean were examined. High-throughput RNA-sequencing analysis using Cufflinks pipeline (Trapnell et al., 2012) was performed in the following three combinations: (i) control vs. drought, (ii) control vs. flood, and (iii) drought vs. flood. The DEGs specific to and common between the stresses were also identified (Figures 1,2). A total of 36,968 and 37,434 genes with confident expression in drought and flood conditions, respectively, when compared to the control (Supplementary Tables 1A,B). Out of these genes, 36,239 were found to be common between the two conditions. The top 50 genes with >500 FPKM featured in the drought and flooding conditions were identified (Tables 2, 3). Interestingly, there are 15 genes among the top 50 highly expressed ones enriched for photosystem II (PS II) and dehydrin family proteins featured in both stress conditions. Flooding stress specific genes were found to be enriched for lipoxygenase 2, ethylene forming enzyme, matrixin family proteins, and 12-oxophytodienoate reductase 2, whereas drought specific genes were enriched for S-adenosylmethionine synthetase, calmodulin 5, 1-deoxy-D-xylulose-5-phosphate-reductoisomerase, and carbonic anhydrase (Tables 2, 3).

The distribution trends in terms of fold change ranged from ~8- to 9-fold change for DEGs under drought and flood stress (Supplementary Figure 2). A total of 2724 DEGs were identified under the drought conditions when compared to control, and 1802 genes were up-regulated and 922 genes were down-regulated (Supplementary Table 3A). During flood stress, a total of 3498 DEGs were identified when compared to the control, and 1303 genes were up-regulated and 2195 genes were down-regulated (Supplementary Table 2B). The total number of DEGs was greater in the flood condition compared to the drought. However, a greater number of up-regulated genes were observed under the drought condition compared to the flood condition with an overlap of 166 genes across the three comparisons (Figure 1). A total of 3768 genes were found to be differentially expressed when flood and drought conditions were compared to each other (Supplementary Table 2C). However, among the 2724 and 3498 DEGs identified in the drought and flood conditions, respectively, 823 DEGs were found to be in common (Figure 2A). Among these genes, 461 genes were both up-regulated under the drought and flood conditions, 290 genes were both down-regulated under the two conditions, 62 genes were up-regulated in the drought but down-regulated in the flood conditions, and 10 genes were down-regulated in the drought but up-regulated during flood conditions (Figure 2B).

The profiles of biological processes represented by GO enriched DEGs under drought and flood stresses were studied and represented (Figure 3, Supplementary Table 3). The heat-map revealed a stark difference in biological processes under different GO categories, such as phagocytosis, cell morphogenesis, cell cycle, isoprenoid biosynthesis, and transcriptional regulation represented by the enriched genes in these categories. However, little differences were observed in the profiles of the genes enriched for the GO categories “responsive to stress” and “defense response” under drought and flooding conditions. Biological processes such as “tRNA processing,” “cell wall macromolecule catabolic process,” “ubiquitin-dependent protein catabolic process,” and “carbohydrate transport” were reduced under both conditions (Figure 3). Most of the GO categories showed an opposite and comparable profile under drought and flood conditions, which suggested differential and specific gene regulation in both conditions. The GO categories of “cell redox homeostasis,” “intracellular protein transport,” “phosphatidylinositol phosphorylation,” “phagocytosis,” “embryo development,” and “Golgi vesicle transport” were enriched under flood condition but reduced under the drought conditions. On the contrary, “protein-heme linkage,” “regulation of transcription DNA-dependent,” “nucleoside metabolic process,” and “metabolic process” were enriched under drought but reduced under flooding conditions.

The RNA-Seq expression profiling revealed 289 and 271 differentially expressed TFs under drought and flood stress, respectively. A total of 213 and 199 TFs were up-regulated, whereas 76 and 72 TFs were down-regulated under drought and flood conditions, respectively (Table 4). Interestingly, the majority of differentially expressed TFs, irrespective of a family, were up-regulated under both conditions (Table 4). Genes belonging to the bHLH, ERF, MYB, NAC, and WRKY family represent most of the differentially expressed TFs. The Ethylene Response Factor (ERF) represented the highest number of significantly expressed genes under both drought and flooding conditions followed by bHLB, MYB, NAC, and WRKY (Table 4, Figure 4, Supplementary Figure 5).

The photosynthesis related genes were down-regulated under both drought and flooding conditions. The expression level of three PS-II light harvesting complex (LHC) related genes, Glyma16g28030.1, Glyma16g28070.1, and Glyma10g04630.1, were repressed significantly under drought conditions. Similarly, three photosynthesis light reaction related genes, Glyma07g21150.1, Glyma20g24500.1, and Glyma20g24500.1, and one Calvin cycle related gene, Glyma01g34300.1, were also down-regulated under drought conditions. Only one gene from the photosynthetic system encoding an electron carrier was induced under the drought conditions (Supplementary Figure 6A). However, compared to the drought, more photosynthesis related genes were found to be down-regulated under flooding conditions (Supplementary Figure 6B). The relative expression level of 34 genes encoding PS-II LHC and 45 genes encoding PS-II subunits were repressed (Supplementary Figure 6B). In addition, 41 genes related to the Calvin cycle and eight electron carrier related genes were also repressed. In total, seven genes were down- and one was up-regulated under drought stress. However, 128 genes were down-regulated under flood treatment (Supplementary Figure 6B). The dynamic expression of photosynthesis related genes can be correlated with the chlorophyll content under drought and flood conditions, which concurs with the inhibition of photosynthesis (Supplementary Figures 3A,B).

Changes in the expression levels of a large number of sucrose and starch metabolism pathway genes were observed under the drought and flood conditions. Genes encoding sucrose and starch biosynthesis genes, namely, amylose synthase, ADP-Glucose synthase, and starch and maltose synthase, showed a significant decrease in the transcript level in drought and flooding treatments (Supplementary Figure 4). No change in expression was observed for sucrose phosphate synthase genes (Glyma08g42140.1 and Glyma06g48200.1) in the drought condition, in contrast to the significant repression in expression of these genes under the flood conditions. However, the genes encoding beta-fructofuranosidase (Glyma13g42530.1), fructokinase (Glyma16g32530.1), hexokinase (Glyma11g10130.1), and three maltose synthase genes (Glyma13g28630.1, Glyma11g04210.1, and Glyma15g10480.1) were up-regulated under the drought conditions. Overall, transcript abundance for starch biosynthesis encoding genes was reduced more in the flooding stress compared to the drought stress. This could result in the decrease of starch content under flooding conditions (Supplementary Figure 4B). Only one maltose synthesis gene, Glyma11g04210.1, was induced during flooding (Supplementary Figures 4, 6). We also identified six genes of the soybean sugar-effluxer (SWEET) family that was differentially expressed under the two extreme water regimes. Three genes, Glyma12g36300, Glyma14g27610, and Glyma06g13110, were up-regulated under the drought conditions, and the other three genes, Glyma06g17520, Glyma09g04840, and Glyma13g33950, were down-regulated under flooding conditions. Interestingly, one SWEET gene, Glyma14g27610, was up-regulated in both drought and flooding conditions (Supplementary Tables 1A,B).

The annotated RNA-Seq data was analyzed for the cell wall precursor's related genes. Expression of 11 cell wall precursor's related genes were up-regulated in the drought conditions, whereas not a single gene was observed to be down-regulated. On the contrary, the flooding condition showed a reduced transcript level for 10 cell wall precursor synthase genes and three genes were up-regulated (Supplementary Figure 7). The UDP-Glucosyl Transferase enzyme coding gene under drought stress was found to be highly expressed (FPKM = 1240.53; Table 2). This enzyme is known to be involved in the formation of β-glucans, which are thought to be involved in cell wall formation (Kosegarten et al., 1988).

Genes with similar expression patterns are likely co-regulated and have the same cis-regulatory elements (CREs or known as motifs; Priest et al., 2009). The set of DEGs (up- and down-regulated genes) identified under flood and drought treatments in this study was utilized (Supplementary Figure 2). Using the collection of motif sequences from the “Database of Plant Cis-acting Regulatory DNA Elements,” three kb upstream of the 5′ most translation start base for all of the annotated soy genes was searched (Phytozome v9.1). Sixteen and 76 over-represented motifs in the down- and up-regulated drought response genes, respectively, and 48 and 5 over-represented motifs in the down- and up-regulated flood response genes, respectively were identified (Supplementary Table 4A). No overlapping motifs were found among the drought and flood response up-regulated genes; however, eight overlapped motifs were found between the drought and flooding down-regulated genes. Additionally, no overlapped motifs were found between the drought down-regulated genes and flooding up-regulated genes, and 30 overlapped motifs were found between the up- and down-regulated drought and flood genes (Supplementary Table 4B). Among the 30 overlapped motifs, 24 showed the ACGT sequence, which serves as a binding site for the bZIP class of TFs. It has been reported that the flanking sequence of the ACGT core motif is responsible for differential distribution of the stress responsive regulatory elements in the stress responsive promoters (Suzuki et al., 2005). The nucleotides flanking the ACGT motif determine the binding affinity and specificity of bZIP TFs, which in turn leads to distinct physiological functions in stress responses, light regulation, etc. (Izawa et al., 1993; Foster et al., 1994).

Although drought stress up-regulated genes shared several cis-regulatory elements with down-regulated flood stress genes, few genes were found to be common between drought up-regulated and flooding down-regulated genes. For example, there are 324 up-regulated drought responsive genes and 457 down-regulated flood response genes with only one shared motif, ACGTABREMOTIFA2OSEM. Importantly, this motif acts at the core of ACGT of motif A in the rice ABRE gene (Hattori et al., 2002). However, there were only 16 common genes representing 5% of the drought up-regulated genes and 4% of the flood down-regulated genes.

Drought and flood down-regulated genes had ~15% of their genes in common (Supplementary Table 3D). Interestingly, flood down-regulated genes shared most their over-represented motifs (30/48) with drought up-regulated genes. However, they did not share many genes between the down-regulated flooding response and the up-regulated drought responsive genes (Supplementary Tables 3C,D). The molecular function categories were also examined using the GO enrichment analysis for both of the datasets represented by the ACGTABREMOTIFA2OSEM motif. Interestingly, the top three molecular functions for both of the datasets were “catalytic activity,” “sequence-specific DNA binding TF activity,” and “DNA binding.”

To validate RNA-Seq analysis, qRT-PCR was performed. A few key genes were selected based on increased or reduced transcript abundance under drought and flooding conditions. Although the RNA-Seq values showed slight variations compared with the corresponding values from the qRT-PCR analyses, the expression data and pattern from RNA-Seq were largely consistent with those obtained by qRT-PCR with a high correlation coefficient value (R² = 0.70, 0.76 for drought and flood responsive genes, respectively; Figure 5).

The involvement of TFs was examined in the drought and flood stress responses through both the ABA dependent and independent pathways. The ABA and other upstream signals can trigger the downstream regulation pathway through the TF regulatory network in both drought and flood conditions. In this research, ABA related TFs: AP2, bZIP, MYB, and NAC were up-regulated under both drought and flood conditions (Figure 4). These TF families have been extensively studied in many plants under several different stress conditions (Dubos et al., 2010; Hao et al., 2011; Mizoi et al., 2012; Nakashima et al., 2012,). In several abiotic stress regulation studies, the TFs identified are initially involved in tolerance against a specific stress and later observed to be important for other stresses. For instance, the SUB1A TF was initially identified as a key regulator for the flood responsive pathways in rice, and later transcriptional and functional studies confirmed its important role in drought tolerance (Xu et al., 2006; Fukao et al., 2011). However, the flood response in soybean seems to be independent of the SUB1 pathway since no SUB1A homolog was found to be differentially expressed in our study. On the other hand, recent studies have shown the involvement of SUB1 and core-clock related genes mediated through alternatively spliced forms (Syed et al., 2015).

Genes including TFs involved in the ABA dependent stress response were found to be differentially expressed under both drought and flood stress conditions (Table 4, Supplementary Table 2). Although ABA is the most specifically studied hormone for its regulatory role in different stress conditions, other hormones such as cytokinins, brassinosteroids, auxins, and ethylene also played an important functional and regulatory role (Peleg and Blumwald, 2011; Patil and Nicander, 2013). The highest number of TFs belonging to the ERF family coincided with the higher number of ethylene related genes observed under drought and flooding conditions (Table 4). This suggests the role of ethylene in both drought and flood responsive pathways. Ethylene biosynthesis has been found to be affected by auxin, particularly with the auxin dependent roles of several members of the 1-aminocyclopropane-1-carboxylate synthase (ACS) family genes (Nakatsuka et al., 1998). S-Adenosylmethionine synthetase (SAM synthetase) coding genes were found to be highly enriched and featured among the top 50 highly expressed genes under drought stress. SAM acts as a precursor for polyamines and ethylene biosynthesis, which indicates a high efflux of ethylene hormone under stress to counteract the adverse conditions. Higher ethylene levels are known to induce the ethylene signaling cascade, which in turn involves the ERF TF activation. That is also found to be highly expressed under the stresses imposed in the current study. Higher ethylene levels are manifested as the ethylene triple response (Knight et al., 1910) to prevent root and shoot elongation, radial swelling of the root and shoot, and horizontal growth of the embryonic leaves and meristem instead of vertical growth (a consequence of auxin imbalance across the stem axis). This triple response prevents the plants from damage. It would also be interesting to see the receptor response against elevated ethylene levels as reviewed by Agarwal et al. (2012) because different plants perceive ethylene in their own unique way and the receptors have a very important role in signal perception and the downstream cascade.

The 1-Deoxy-D-xylulose-5-phosphate reductoisomerase (DXR) is an important enzyme involved in the 2-C-methyl-D-erythritol-4-phosphate (MEP) pathway for isoprenoid biosynthesis and was identified among the highly expressed drought responsive genes (Table 2). Xing et al. (2010) reported that dxr mutants with a dwarf phenotype significantly reduced the number of trichomes and abnormal stomatal closure, which was rescued with exogenous application of ABA and gibberellic acids (GAs). Thus, the high expression of DXR under drought stress implies a plant's response to withstand stress and compensate for compromised photosynthetic capacity and abnormal stomatal conductance. Auxin, ethylene, and ABA were previously identified to be involved in the drought responsive pathway in soybean (Le et al., 2012). Additionally, ethylene, ABA, and GA are known to be involved in a regulatory network under flooding stress (Bailey-Serres et al., 2012). Hormonal cross-talk at different levels of the stress response and the resultant synergetic or antagonist interactions have been found to be important for stress tolerance (Chinnusamy et al., 2004; Peleg and Blumwald, 2011). The biosynthesis of jasmonic acid (JA) through the octadecanoic pathway (Vick and Zimmerman, 1984) with α-linolenic acid (LA) as a precursor, is mediated by 12-oxophytodienoic acid reductase 2 (OPDAR2) (Glyma01g44600), which is found to be up-regulated under both drought and flooding stress (Supplementary Tables 1B, 2A) but highly expressed under flooding stress (Table 3). The high expression of OPDRA2 is correlated with higher JA synthesis and other systemic responses. The enzyme has also been shown to detoxify 2,4,6-trinitro toluene (TNT) in polluted water contaminated by explosives as a result of manufacturing and munitions (Beynon et al., 2009). The upsurge of OPDAR2 under flooding stress suggests plants' adaptive response against stress and this could also help in detoxification of contaminants such as TNT if any are present in the water. DEGs involved in different hormonal pathways identified here under drought and flooding stress could be good candidates to manipulate hormonal cross-talk for the effective regulation of stress tolerance in plants.

The genes involved in photosynthesis were found to be preferentially down-regulated under both drought and flood conditions in this study (Supplementary Figure 6). This agrees with the hypothesis that a plant's photosynthetic capacity is affected under stresses such as drought and flooding stresses. Thus, this implies that the reduction in photosynthesis activity is the major effect of drought and flooding in soybean. Similarly, the down-regulation of genes involved in photosynthesis has been previously reported under a water deficit situation in other crops (Wang et al., 2011). Interestingly, Nanjo et al. (2011) reported the up-regulation of photosynthesis related genes under flooding; however, they imposed stress at the 2 day-old seedling stage of soybean and performed transcriptome profiling with hypocotyl tissue. Similarly, during the early stress response in Arabidopsis (low oxygen, 3 h) and gray poplar (flooded 5 h), the up-regulation of photosynthesis related genes was observed (Kreuzwieser et al., 2009). In our study, the down-regulation of photosynthesis related genes was observed after 7 days of flood treatment (Supplementary Figure 6). The down-regulation of photosynthesis related genes has also been reported in the leaves of rape seedlings waterlogged for 3 days (Lee et al., 2014). This suggests that the flooding or low-oxygen stress stimulates the expression of photosynthesis related genes during the initial stress response but it is repressed at later stages. Flooding results in reduced photosynthesis manifested by depleted carbon dioxide (CO₂) concentration in leaves mediated by the stomatal constraint to the exchange gases (Bailey-Serres et al., 2012). Reduced CO₂ levels lead to a depleted electron sink in chloroplasts, which ultimately results in accumulation of photons that are a potential source of ROS. The shrunken electron sink under flooding conditions is also substantiated in our findings and is demonstrated as the repressed expression of eight electron carriers and PS II related genes (Supplementary Figure 6B). The chlorophyll synthesis related genes were found to be repressed under both drought and flooding conditions in this study (Supplementary Figure 3). The lower availability of oxygen (O₂) and CO₂ under flooding caused by the low gas exchange ratio (Bailey-Serres et al., 2012) may result in decreased photosynthesis efficiency and carbohydrate content. This was consistent with the reduced expression level of starch synthesis genes in the present study (Supplementary Figure 4).

Reduced photosynthesis during prolonged abiotic stresses such as drought and flooding can be correlated with the production of ROS (Klok et al., 2002). High levels of ROS can damage the cell through peroxidation of the lipids, oxidation of proteins and other pathways, and finally the death of the cell (Clement et al., 2008; Wrzaczek et al., 2011). Glutathione S-transferase, one of the important enzyme in the regulation of ROS flux, is featured as a highly enriched gene among the top 50 highly expressed genes under flooding stress in this study, which indicates another adaptive response under the stress conditions (Table 3). However, the opposite (high GST expression in response to drought compared to flooding stress) was reported by Oh and Komatsu (2015).

Drought also has severe effects on the stomatal conductance and low CO₂ assimilation. It is often considered to be one of the most prominent reasons for poor photosynthesis. However, it is not the only reason. Another factor is reduced ATP synthesis leading to depleted ribulose 1,5-bisphosphate content (Tezara et al., 1999), which produces a 6C compound when combined with CO₂ in the dark reaction (Calvin cycle) of photosynthesis. The repressed expression levels of genes associated with the Calvin cycle in this study (Supplementary Figure 6A) substantiates the decreased photosynthetic capacity mediated by down-regulated dark reaction genes. However, the observed up-regulated electron carrier related gene in drought stressed plants needs to be justified when the Calvin cycle and LHC related genes are down-regulated. This is because under lesser CO₂ assimilation in drought stress, the photosynthetically accumulated energy becomes excess for the dark reaction, and to counteract this condition, over reduction of the electron transport chain occurs (Biehler and Fock, 1996).

Water stress can change the turgor pressure of the plant cell and it responds accordingly, either by tightening or loosening the wall under the extreme water conditions of drought or flooding conditions (Moore et al., 2008). UDP-glucuronic acid is the major precursor of cell wall polysaccharides and comes from two sources: UDP-glucose and myo-inositol (Siddique et al., 2013). The genes involved in both UDP-glucuronic acid synthesis pathways were up-regulated under drought stress (Supplementary Figure 7A). This can cause the accumulation of cell wall precursors and the increase of cell turgor under drought conditions. However, an adverse result was observed under flooding. The GDP-glucuronic acid synthesis genes were reduced under flooding (Supplementary Figure 4). This implies the loosening of the cell wall and stimulation of shoot elongation under flooding (Voesenek et al., 2006).

Most research on the effect of drought and flooding on sugar metabolism and starch biosynthesis has demonstrated that the carbohydrate levels are altered in leaves (Lemoine et al., 2013). Under stress conditions, the increase in hexose amounts is associated with the induction of starch hydrolysis (Pelleschi et al., 1997). In agreement with these findings, we observed an increased expression of the hexokinase (Glyma11g10130.1) gene in the drought treatment. The gradual slowing down of starch metabolism is a general mechanism under flooding (Fukao et al., 2011) and drought (Burke, 2007) conditions for rice and cotton, respectively, which reserves the energy source for a prolonged energy supply to maintain cell survival. It is conceivable that starch metabolism during flooding and drought response may also be controlled/gated by the up- and down-regulation of carbohydrate metabolism genes in soybean. One of the prominent processes under sugar metabolism that gets affected in water stressed plants is the sugar efflux rate, which is carried out by proteins coded by sugar transporter genes. They ensure sugar transport from the source (leaves) to the sink (roots and seeds) to counteract the depleted sugar levels under severe conditions such as drought and flooding. Utilizing the available resource generated by Patil et al. (2015), six SWEET family genes were featured among the DEGs identified in this study under drought and flood conditions. The findings suggest that SWEET genes in general are induced in drought conditions and repressed in flooding conditions, which indicate that sucrose efflux genes are activated under drought conditions and repressed under flooding conditions. However, one gene (Glyma14g276210) was found to be up-regulated under both of the conditions. Syed et al. (2015) also identified the up-regulation of soybean SWEET genes under drought stress. The activated gene expression under drought stress could be accounted for by sugar transportation due to an uninterrupted energy supply under such adversity. The down-regulation of three SWEET family genes under flooding stress in the current study is contrary to the findings of Lemoine et al. (2013), where it was suggested that the up-regulated sugar effluxer genes may contribute to sugar homeostasis in a flood tolerant genotype. In this study, the repressed genes might be playing the same role by controlling the efflux of sugar under flooding stress and thus putting a check on normal efflux of sugar. However, this hypothesis needs additional substantiation.

Glycolytic enzymes along with inducers of heat shock proteins are key factors in the early response to flooding and drought in soybean (Nanjo et al., 2011; Oh and Komatsu, 2015), which suggests the significance of the glycolytic pathway in the adaptation to flooding and drought conditions. Oh and Komatsu (2015) reported increased expression of all glycolysis related proteins under both drought and flooding responses. The up-regulated expression of genes coding for the glycolytic pathway enzymes, fructokinase (Glyma16g32530.1) and hexokinase (Glyma11g10130.1), under drought stress can be considered an adaptive response of soybean under drought stress. On the contrary, the transcript levels of glycolytic pathway enzymes were not significantly affected under water deficit conditions in chickpea (Khanna et al., 2014). This suggests a legume specific adaptive variability under similar stress.

The analysis of cis regulatory elements showed that although the genes between drought up-regulated and flooding down-regulated DEGs shared a common motif, they might have a similar function with opposite regulation under contrasting stress treatments. It has been reported that proteins with up to 67% amino acid sequence similarity may share similar high-affinity binding sequences, and prefer different low-affinity sites. However, even closely related TFs may have distinct DNA-binding profiles (Badis et al., 2009). The positions of the flanking sequence of motif sites also play an important role in achieving regulatory specificity within the TF families (Ciolkowski et al., 2008; Gordan et al., 2013). Therefore, it is possible that this opposite regulation between drought up-regulated genes and flooding down-regulated genes with the same motif could be due to the different TF binding sites or because of the different site specific cofactors with the same affinity for the same TF binding sites, which leads to activation or repression of the target genes.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.