Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Insights into soybean transcriptome reconfiguration under hypoxic stress: Functional, regulatory, structural, and compositional characterization

  • Thiago J. Nakayama,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

    Affiliation Departamento de Fitotecnia, Universidade Federal de Viçosa, Viçosa, Minas Gerais, Brazil

  • Fabiana A. Rodrigues,

    Roles Conceptualization, Investigation, Methodology

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  • Norman Neumaier,

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Writing – review & editing

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  • Juliana Marcolino-Gomes,

    Roles Investigation

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  • Hugo B. C. Molinari,

    Roles Writing – review & editing

    Affiliation Embrapa Agroenergia, Empresa Brasileira de Pesquisa Agropecuária, Brasília, Distrito Federal, Brazil

  • Thaís R. Santiago,

    Roles Formal analysis, Software

    Affiliation Embrapa Agroenergia, Empresa Brasileira de Pesquisa Agropecuária, Brasília, Distrito Federal, Brazil

  • Eduardo F. Formighieri,

    Roles Formal analysis, Software

    Affiliation Embrapa Agroenergia, Empresa Brasileira de Pesquisa Agropecuária, Brasília, Distrito Federal, Brazil

  • Marcos F. Basso,

    Roles Writing – review & editing

    Affiliation Embrapa Agroenergia, Empresa Brasileira de Pesquisa Agropecuária, Brasília, Distrito Federal, Brazil

  • José R. B. Farias,

    Roles Funding acquisition, Project administration, Resources

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  • Beatriz M. Emygdio,

    Roles Funding acquisition, Project administration, Resources

    Affiliation Embrapa Clima Temperado, Empresa Brasileira de Pesquisa Agropecuária, Pelotas, Rio Grande do Sul, Brazil

  • Ana C. B. de Oliveira,

    Roles Funding acquisition, Project administration, Resources

    Affiliation Embrapa Clima Temperado, Empresa Brasileira de Pesquisa Agropecuária, Pelotas, Rio Grande do Sul, Brazil

  • Ângela D. Campos,

    Roles Investigation, Validation

    Affiliation Embrapa Clima Temperado, Empresa Brasileira de Pesquisa Agropecuária, Pelotas, Rio Grande do Sul, Brazil

  • Aluízio Borém,

    Roles Funding acquisition, Supervision, Writing – review & editing

    Affiliation Departamento de Fitotecnia, Universidade Federal de Viçosa, Viçosa, Minas Gerais, Brazil

  • Frank G. Harmon,

    Roles Writing – review & editing

    Affiliation Department of Plant and Microbial Biology, University of California-Berkeley, Berkeley, California, United States of America

  • Liliane M. Mertz-Henning,

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  •  [ ... ],
  • Alexandre L. Nepomuceno

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision

    alexandre.nepomuceno@embrapa.br

    Affiliation Embrapa Soja, Empresa Brasileira de Pesquisa Agropecuária, Londrina, Paraná, Brazil

  • [ view all ]
  • [ view less ]

Abstract

Soybean (Glycine max) is one of the major crops worldwide and flooding stress affects the production and expansion of cultivated areas. Oxygen is essential for mitochondrial aerobic respiration to supply the energy demand of plant cells. Because oxygen diffusion in water is 10,000 times lower than in air, partial (hypoxic) or total (anoxic) oxygen deficiency is important component of flooding. Even when oxygen is externally available, oxygen deficiency frequently occurs in bulky, dense or metabolically active tissues such as phloem, meristems, seeds, and fruits. In this study, we analyzed conserved and divergent root transcriptional responses between flood-tolerant Embrapa 45 and flood-sensitive BR 4 soybean cultivars under hypoxic stress conditions with RNA-seq. To understand how soybean genes evolve and respond to hypoxia, stable and differentially expressed genes were characterized structurally and compositionally comparing its mechanistic relationship. Between cultivars, Embrapa 45 showed less up- and more down-regulated genes, and stronger induction of phosphoglucomutase (Glyma05g34790), unknown protein related to N-terminal protein myristoylation (Glyma06g03430), protein suppressor of phyA-105 (Glyma06g37080), and fibrillin (Glyma10g32620). RNA-seq and qRT-PCR analysis of non-symbiotic hemoglobin (Glyma11g12980) indicated divergence in gene structure between cultivars. Transcriptional changes for genes in amino acids and derivative metabolic process suggest involvement of amino acids metabolism in tRNA modifications, translation accuracy/efficiency, and endoplasmic reticulum stress in both cultivars under hypoxia. Gene groups differed in promoter TATA box, ABREs (ABA-responsive elements), and CRT/DREs (C-repeat/dehydration-responsive elements) frequency. Gene groups also differed in structure, composition, and codon usage, indicating biological significances. Additional data suggests that cis-acting ABRE elements can mediate gene expression independent of ABA in soybean roots under hypoxia.

Introduction

Regimes of excess water (flooding) influence the distribution and diversity of species in natural ecosystems [1] and lead to yield losses of many farmland crops [2]. The increase in flooding events over the past six decades is associated with climate change, which threatens food security of the growing human population [3]. Among the four major crops (soybean, wheat, maize, and rice), only rice is adapted to soil waterlogging and all are sensitive to submergence [2].

The oxygen diffusion in water is 10,000 times lower than in air [4]. Thus, water surrounding roots (waterlogging) or entire plants (submergence) can cause severe energy crisis once oxygen is required for energy production through mitochondrial respiration [5]. Additionally, even when oxygen is externally available, oxygen deficiency is frequent metabolic status of bulky, dense or metabolically active tissues such as phloem, meristems, seeds, and fruits [6]. It is required in several metabolic pathways, including heme, sterol, and fatty-acid biosynthesis [6].

Orchestrated by complex gene regulatory network, plants and other organisms need to perceive, signal, and promote biochemical, physiological, and morpho-anatomical changes appropriate to survive and thrive under oxygen level fluctuations. The adaptation capacities of crops validated under stress field conditions have shown association with gene duplication events [7]. As an example, the multigenic SNORKEL (SK) and SUBMERGENCE-1 (SUB1) loci are found in deep-water and submergence-tolerant rice varieties, respectively [8,9]. Both loci are members of AP2 (APETALA2)/ERF (Ethylene Responsive Factor) plant-specific transcription factors superfamily. They encode tandem-repeated genes, of which SK2 [8] and SUB1A-1 [9] are majorly responsible for the rapid-growth avoidance escape (SK2) and energy saving quiescence (SUB1A-1) strategies against oxygen deprivation.

Unlike rice, soybean genome does not contain SUB1 and SK2 orthologs [10]. So far, only QTLs (Quantitative Trait Loci) with small effect for waterlogging tolerance have already been predicted in soybean [11,12]. In addition, few soybean studies of transcriptome-wide responses to flooding stress have already been reported, and those that have all examine a single genotype. Chen et al. [13] evaluated the leaf transcriptome from adult soybean seven days after waterlogging. Yin et al. [14] performed transcriptomic analysis of flooding-tolerant line and ABA-treated newly germinated seedlings under hypoxic stress. Others studies analyzed transcriptomes from shoots (cotyledons including hypocotyls) [10] and roots including hypocotyls [15,16] of newly germinated soybean seedlings under hypoxic stress.

In the present study, we analyzed conserved and divergent root transcriptional responses between flood-tolerant Embrapa 45 and flood-sensitive BR 4 soybean cultivars under hypoxic stress by using RNA-seq platform. For progress to understand how soybean genes evolve and respond to hypoxic stress, stable (SGs) and differentially expressed genes (DEGs) were structurally and compositionally characterized, comparing its mechanistic relationship with expression regulation.

Materials and methods

Grain yield

Under field conditions, Embrapa 45 and BR 4 seeds were sown on December 21sh, 2011, at Embrapa Clima Temperado, Pelotas, RS, Brazil (latitude 31°42'S and longitude 52°24'W). The experiment was carried out in randomized blocks, with four replicates (75 plants per plot, at a density of 200,000 plants/ha). The seed emergence occurred six days after sowing. The first waterlogging occurred 18 days after sowing due to heavy rain, lasting five days. We observed mild symptoms of yellowing leaves. On February 15th, 2012, we waterlogged the soil beds, maintaining water level 2 cm above the soil surface for 10 days. The plants developed typical reaction to waterlogging. Harvest was done on May 23th, 2012. All seeds were collected and corrected to 13% moisture content. Data met assumptions of the analysis of variance (ANOVA). Thus, means were compared by the Tukey test 5%.

Plant material for RNA-seq and qRT-PCR analysis

Previously described [17], using a hydroponic system under greenhouse conditions, the experiment was set in a randomized block design composed of twelve treatments: two cultivars (Embrapa 45 and BR 4), two oxygen conditions [fully aerobic state (normoxy) and hypoxic], and three treatment sampling times (0.5h, 4h, and 28h). Each treatment has three biological replicates (four plantlets per replicate in order to reduce biological variation). At each time point, root tissues were collected and immediately frozen in liquid nitrogen before being stored at -80°C. We compared stressed and unstressed samples at the same time point to remove putative additive effects, such as gene-intrinsic effects (e.g., circadian rhythm [18]), differences in developmental stages among individuals, or any unknown variation between the time points [17].

Total RNA was extracted using Trizol reagent (Invitrogen) and treated with DNase I (Invitrogen) according to manufacturer instructions. RNA concentration and purity were measured using a spectrophotometer (NanoDrop, ND-1000), and the integrity of the molecules was analyzed on 1% agarose gels stained with ethidium bromide.

Transcriptome library construction, deep sequencing, and mapping of reads

For each of twelve treatments, equimolar quantities of purified total RNA from roots of twelve plants were pooled to result one library. Then, the twelve libraries were sent to Fasteris SA (Plan-les-Ouates, Switzerland) for processing and sequencing.

The RNA quality and integrity was checked using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA), of which only samples with a RIN ≥ 8.0 were used. The twelve libraries were processed (poly-A purification, fractionation, cDNA synthesis using random primers, and ligation to bar-coded adapters), fragments of 150–250 bp were isolated and multiplexed, resulting one sequencing library. The sequencing library is a pool of equimolar quantities from twelve initial libraries, each library with a specific barcode for further bioinformatic discrimination. Single end reads were generated by the Illumina HiSeq 2000 (read length 1 × 100 base; one lane of the flow-cell; Illumina, Inc. San Diego, CA). The raw data, deposited in the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5709, was uploaded to the GeneSifter platform (Geospiza, Seattle, WA) for alignment with the soybean genome (assembly Glyma 1.1) [19]. The mapping of reads and transcripts analysis was done as described previously [18].

For structural analysis of non-symbiotic hemoglobin Glyma11g12980.1 transcript, reads from each cultivar were de novo assembled with Trinity [20] (standard parameters with minimum contig length of 400bp) and mapped to the Glyma11g12980.1 reference (Phytozome transcript model) using BWA-MEM (v0.7.5) [21] default settings.

Transcriptomic analysis

For each time point (0.5, 4, and 28h), the expression ratio (fold-change, fc) of genes was performed by dividing transcript abundance values (in RPM = Reads per Mapped Million) from plants under hypoxic and normoxic conditions. The statistical significance of DEGs were obtained by using Bioconductor package edgeR [22], corrected by Benjamini and Hochberg method [23]. We only considered as DEGs those showing fold-change ≥ 2 (up), ≤ -2 (down), adj. p-value ≤ 0.01, and with more than 20 mapped reads (RPM ≥ 9) in at least one of the two compared libraries.

Gene set enrichment analysis was performed using Singular Enrichment Analysis (SEA) provided by agriGO (http://bioinfo.cau.edu.cn/agriGO/) [24]. We chose hypergeometric test, corrected by Hochberg FDR method, plant GO slim database. Soybean pathways of DEGs related to amino acids and its derivatives were analyzed using KEGG Mapper—Search&Color Pathway [25]. The Relative Synonymous Codon Usage (RSCU) was calculated with CodonW 1.4.4 (http://mobyle.pasteur.fr/cgi-bin/portal.py?#forms::CodonW). The copy numbers of soybean nuclear tRNA genes was extracted from the genomic tRNA database (http://gtrnadb.ucsc.edu/) [26]. Samples sizes used for structural and compositional analysis among groups of genes are shown in the Table A in S1 Supporting Information.

The Microsoft Excel (Microsoft, Redmond, WA) was used for TATA-box searching in promoter regions and structural/compositional characterization of genes. For hypothesis testing on binary data, we used Microsoft Excel add-in Real Statistics Resource Pack (ver. 3.3.1, http://www.real-statistics.com/).

qRT-PCR analysis

For quantitative real-time PCR (qRT-PCR) analysis, the synthesis of cDNA, design of primers, and expression analysis of genes used to verify reliability of RNA-seq expression data were done as described previously [17] (Table D in S1 Supporting Information).

Results and discussion

In the present study, we analyzed root RNA-seq data from flood-tolerant Embrapa 45 and flood-sensitive BR 4 soybean cultivars that showed contrasting grain yields when cultivated in waterlogged soil (Fig 1). In order to evaluate pairwise RNA-seq data, the relative expression of six common hypoxia-responsive genes [Trehalose-6-Phosphate Synthase (Glyma17g07530), Ascorbate Peroxidase (Glyma12g07780), Sucrose Synthase (Glyma13g17420), Alternative Oxidase (Glyma04g14800), non-symbiotic Hemoglobin (nsHB; Glyma11g12980), and Nitrate Reductase (Glyma06g11430)] were determined by qRT-PCR. The non-log-transformed qRT-PCR and RNA-seq expression data were consistent for all these genes (Fig 2) showing a strong positive Pearson correlation (r = 0.95; P < 0.001), indicating reliability in our transcriptome analysis.

thumbnail
Fig 1. Grain yields of two soybean cultivars (flood-tolerant Embrapa 45 and flood-sensitive BR 4) under two moisture regimes in the soil.

Control indicates well-watered conditions (70% available water in the soil). n = 4, except for BR 4 under waterlogging (n = 3). Each biological repetition consisted of 75 plants, from which grain yield were converted to 200,000 plants per hectare. Means values (± S.E.M.) followed by different capital letters between cultivars under same soil condition, and lowercase letters between soil conditions for same cultivar, significantly differ according to Tukey test 5%.

https://doi.org/10.1371/journal.pone.0187920.g001

thumbnail
Fig 2. Validation of pairwise RNA-seq data through qRT-PCR.

Six hypoxia-responsive genes (TPS: Trehalose-6-Phosphate Synthase, Glyma17g07530; APX: Ascorbate Peroxidase, Glyma12g07780; SUSY: Sucrose Synthase, Glyma13g17420; AOX: Alternative Oxidase, Glyma04g14800; nsHB: non-symbiotic Hemoglobin, Glyma11g12980; and NR: Nitrate Reductase, Glyma06g11430) were analyzed in flood-tolerant Embrapa 45 (Blue chart) and flood-sensitive BR 4 (Green chart) soybean cultivars. The transcripts abundance of the target genes from plants subjected to hypoxic conditions for different periods of time were compared with the respective controls (normoxic condition). Differential gene expression statistically significant: *p < 0.05, **p < 0.01, and ***p < 0.001. From qRT-PCR, raw data was normalized using the ELF1B and ACTB reference genes [17]. The no qRT-PCR amplification of Glyma11g12980 is clarified in Fig B in S1 Supporting Information.

https://doi.org/10.1371/journal.pone.0187920.g002

Transcriptome reconfiguration

Induction of signaling genes and down-regulation of genes related with energy-consuming processes under hypoxia.

From 54,174 predicted protein-coding genes in the soybean genome (assembly Glyma 1.1) [19], 2,656 up-regulated (URGs) and 4,970 down-regulated genes (DRGs) were found in both cultivars under hypoxic stress. Of these total, after 0.5, 4, and 28h of root hypoxia, 1,144, 5,687, and 3,761 genes were differentially expressed, of which 89, 28, and 46% were URGs, respectively (Fig 3). In Arabidopsis thaliana, another flood-sensitive species, more URGs were also observed in roots after 0.5h of hypoxic stress, and URGs were more prevalent deeper into stress conditions [27]. In contrast to our findings in soybean data, where DRGs number decreased after 28h, DRGs remained after 168h of waterlogging in roots of flood-tolerant gray poplar (Populus × canescens) [28]. Even so, Embrapa 45 showed fewer URGs and more DRGs than BR 4, of which the greatest difference in DEGs between cultivars was after 28h where Embrapa 45 (4,216 DRGs) had more DRGs than BR 4 (2,582 DRGs) (Fig 3).

thumbnail
Fig 3. Number and GO enrichment analysis of up-, down-, and stable-regulated genes in flood-tolerant Embrapa 45 (E), flood-sensitive BR 4 (B), and in both soybean cultivars (C).

Differentially expressed genes: fold-change ≥ 2 (up), ≤ -2 (down); adj. p ≤ 0.01; RPM ≥ 9 (control or stress datasets). Stable-regulated genes: fold-change ≤ 1.1 or ≥ -1.1; RPM ≥ 9 (control and stress datasets). RPM ≥ 9: at least 20 reads in all datasets. Red GO names are cited in the text.

https://doi.org/10.1371/journal.pone.0187920.g003

Gene function was determined by identifying Gene Ontology (GO) categories for DEGs and SGs. The most enrichment for GO categories was found after 0.5h in URGs from both cultivars and after 4-28h in DRGs mainly from Embrapa 45 (Fig 3). Noteworthy GO categories in URGs were gene expression regulation (GO:0010468), more specifically for transcription (GO:0006350) and protein modification (GO:0006464) involving transcription factors (GO:0003700) and kinases (GO:0016301) activities. Overrepresented in DRGs were energy-demanding processes, including transport (GO:0005215) and biosynthesis (GO:0009058), as well as translation (GO:0006412), most of which encode ribosomal proteins (GO:0005198; structural molecule activity). Further, transcription factors (GO:0003700), kinases (GO:0016301) and transporters (GO:0005215) were enriched in DEGs, while more general functions such as binding (GO:0005488) in SGs. Our results show that hypoxia induces controlling/signaling genes and suppresses genes related with energy-consuming processes in soybean. Therefore, both induction and repression of genes may be important for flooding tolerance.

Hypoxic soybean roots experience changes in amino acids and amino acid-related metabolism.

After 4h of hypoxia, URGs and DRGs were enriched for reorganization of cellular amino acid and derivative metabolic processes (GO:0006519) (Fig 3). Alterations in amino acid metabolism have been previously observed in hypoxic roots of soybean [29], Lotus japonicus [30], and gray poplar [28], including high accumulation of alanine and GABA (Gamma-Amino Butyric Acid) as well as reduction of aspartate level. In agreement, we observed up-regulation of alanine aminotransferase (Glyma07g05130) and glutamate decarboxylase (Glyma08g09670) at all-time points in both cultivars (S2 File).

Glutamate is directly related to alanine and aspartate metabolism via transamination, and to GABA via decarboxylation [31]. Interestingly, we observed induction of the genes NADH-dependent glutamate synthase (Glyma04g41540 and Glyma19g16486) and aspartate aminotransferase (Glyma14g13480, Glyma17g33050, Glyma06g08670, and Glyma01g32360), as well as repression of ferredoxin-dependent glutamate synthase (Glyma03g28410 and Glyma19g31120), ATP-dependent asparagine synthase (Glyma11g27480, Glyma11g27720, Glyma14g37440, and Glyma18g06840) (ArrayExpress database, accession number E-MTAB-5709). These responses suggests roots change NADH oxidation to save ATP for glutamate synthesis under hypoxia.

While expression of aspartate kinase, aspartate semialdehyde dehydrogenase and homocysteine S-methyltransferase increases under hypoxia, we observed repression of polyamines and phenylpropanoids biosynthesis-related genes, including upstream genes from shikimate pathway (Fig A in S1 Supporting Information, S2 File). Among DRGs was found S-adenosyl-L-methionine (SAM) synthase (EC 2.5.1.6). SAM connects to ethylene, polyamines, and phenylpropanoid-derived lignin pathways (Fig A in S1 Supporting Information) as well as histone and nucleic acid methylation for gene expression regulation [32,33]. Studies involving exogenous application or endogenous production of polyamines via genetic manipulation have shown increased tolerance to a broad spectrum of abiotic stresses [34], opening opportunities for improvement of soybean flooding tolerance by genetic engineering approaches.

Same gene, different regulation between cultivars: Identification of candidate genes for flooding-tolerance.

The phosphoglucomutase (Glyma05g34790) gene was up-regulated in Embrapa 45 and down-regulated in BR 4 soybean cultivar after 4h of hypoxic stress (Fig 4). Its up-regulation in the flood-tolerant cultivar is in accordance with a shift in sucrose catabolism from ATP-dependent invertase-hexokinase to energy-saving SuSy-UGPase pathway [5]. In both cultivars, SuSy (Glyma19g40041, Glyma09g08550, and Glyma15g20180) and UGPase (Glyma11g33160) genes were observed up-regulated, while genes down-regulation were invertase (Glyma10g35890) and hexokinase (Glyma05g35890, Glyma07g12190, and Glyma17g37720) genes.

thumbnail
Fig 4. Relative expression (in fold-change) and expression level (in RPM) in flood-tolerant Embrapa 45 and flood-sensitive BR 4 soybean cultivars along three time points.

The analyzed genes were phosphoglucomutase (Glyma05g34790), N-terminal protein myristoylation (Glyma06g03430), fibrillin (Glyma10g32620), suppressor of phyA-105 (Glyma06g37080), and non-symbiotic hemoglobin (nsHB; Glyma11g12980). Differentially expressed genes: fold-change ≥ 2 (up), ≤ -2 (down); adj. p ≤ 0.01; RPM ≥ 9 (control or stress datasets). RPM ≥ 9: at least 20 reads in all datasets.

https://doi.org/10.1371/journal.pone.0187920.g004

Higher induction of an unknown protein related to N-terminal protein myristoylation (Glyma06g03430), suppressor of phyA-105 (Glyma06g37080), and fibrillin (Glyma10g32620) were observed in flood-tolerant Embrapa 45. The A. thaliana ortholog of Glyma10g32620 (AT3G23400) is required for resistance to multiple stresses [35]. Moreover, differential expression of fibrillin genes correspond to plastoglobule number in leaves of contrasting soybean genotypes under drought and waterlogging stresses [36].

Although the nsHB gene (Glyma11g12980) exhibited similar expression ratio (in fold-change) between the two cultivars, the expression level (in RPM) in Embrapa 45 was lower under normoxic and hypoxic conditions (Fig 4). The last 284 nucleotides of the nsHB 3’UTR (3’Untranslated Region) are absent only in Embrapa 45 (Fig B in S1 Supporting Information). This explains the absence of qRT-PCR amplification in Embrapa 45 (Fig 2) and the similar transcriptional profile (qRT-PCR and RNA-seq) in both cultivars (Fig 2). Considering the important role of nsHB to protect plants under hypoxic stress [37], further study is required to understand the alternate 3’UTR structures influence transcription, transcript stability, and protein abundance.

Function, structure, and composition of gene groups: Comparing its mechanistic relationship with expression regulation

To further understanding how soybean genes evolve and respond to hypoxia, the top 40 transcripts (all time points in both cultivars) and top 500 transcripts (at least one time point in both cultivars, except for S500 keeping all time points as criterion) stable and DEGs were structurally/compositionally characterized, comparing its mechanistic relationship with expression regulation.

The top gene groups differ functionally.

As noted above, the genes aspartate aminotransferase (Glyma01g32360) and SAM decarboxylase (Glyma17g07830) were ranked in top 40 up-regulated (U40) and down-regulated (D40) groups, respectively. Likewise, phenylpropanoid/flavonoid related (D500) and SuSy-UGPase (U500) genes as well as three paralogs of Glyma06g03430 (N-terminal protein myristoylation) (U500) belonged to top 500 groups.

Genes associated with ethylene biosynthesis (ACC oxidase: Glyma05g36310), glycolysis (Pyruvate kinase: Glyma10g34490), ethanol fermentation (pyruvate decarboxylase: Glyma01g29190 and Glyma03g07380; alcohol dehydrogenase: Glyma04g39190 and Glyma14g27940), biotic stress defense (kunitz trypsin protease inhibitor: Glyma16g33770; polygalacturonase inhibiting protein: Glyma05g25370 and Glyma08g08360), and flooding governing acclimation (ERFVII related to A. thaliana HRE2: Glyma19g40070) were found in U40.

In the D40 were genes involved in antimicrobial defense (cysteine-rich secretory proteins: Glyma13g32500, Glyma13g32510, Glyma13g32530, and Glyma13g32540), gene regulation (bZIP: Glyma05g22860 and Glyma17g17100; MYB: Glyma13g20510; NAC: Glyma08g18470), oxygen consuming (2-oxoglutarate oxygenase: Glyma03g23770, Glyma07g12210, and Glyma08g18000; cytochrome P450: Glyma03g02410 and Glyma04g03790), and transport of dicarboxylate (Glyma07g39580), sulfite (Glyma07g30570), manganese (2+) and iron (2+) (Glyma08g08090, Glyma16g28340, and Glyma09g21920). These ions can accumulate to toxic levels in waterlogged plants [3840].

The top 40 stable group (S40) has genes associated with protein (Glyma09g08830, Glyma06g01120, Glyma18g10060, Glyma18g32830, and Glyma04g23560) and nucleic acid binding (Glyma03g01920, Glyma18g11950, Glyma18g32190, Glyma06g48070, Glyma09g06750, Glyma12g29270, Glyma04g04880, and Glyma01g44460). For flooding stress, these genes are promising candidate reference genes for qRT-PCR normalization, given their higher stability compared to stable genes commonly used in the literature [17].

Top gene groups have different TATA box, ABRE, and DRE motif usage in promoter sequences. Does ABRE mediate gene expression independent of ABA in hypoxic soybean roots?

Transcription of protein-coding genes in eukaryotes requires numerous protein factors to recognize specific DNA loci. The core promoter region, proximal to the transcription start-site (TSS), recruits general transcription factors involved in basal transcription [41] and cis-regulatory elements from extended promoter recruits DNA-bound transcription factors (activators or repressors) to fine-tune the transcriptional control [42].

The general regulator TATA-box binding protein (TBP) is required for transcription initiation by all three eukaryotic RNA polymerases [43]. TBP can bind to various DNA sequences but has higher affinity for the consensus TATA-box [44]. Based on previous work [4547], we scanned for the core sequence TATA extending 4 bp in the 3' direction within the 50 bp region upstream of the predicted TSS (between -50 and -1). We found more DEGs (from 21% in U500 to 40% in D40) with the consensus TATA box sequence TATA(T/A)ATA than SGs (at most 6% in S500) (Fig 5, Table B in S1 Supporting Information). Our results are in agreement with hexamer sequences TATA(T/A)A over-represented in A. thaliana, Oryza sativa, and Glycine max genomes [48]. Similarly, Tirosh et al. [45] observed a correlation between consensus TATA-containing genes being differentially expressed and TATA less-containing genes stable expressed in yeast, metazoans, and plants.

thumbnail
Fig 5. Percentage of genes with which putative promoter region contain consensus TATA box, ABRE, and/or DRE/CRT motifs.

Statistical significances from pairwise comparisons are provided in Table B in S1 Supporting Information.

https://doi.org/10.1371/journal.pone.0187920.g005

Higher TBP turnover at consensus TATA- compared to TATA-less promoters is associated with specific coactivators [46,49]. Coactivators are multisubunit complexes represented by SAGA (Spt-Ada-Gcn5-Acetyltransferase), TFIID (transcription Factor II D), related with TBP binding on TATA and TATA-less promoters, respectively [49], and Mediator [50]. The latter is organized into head, middle, tail, and Cdk8 kinase modules to converge and transmit signals from sequence-specific transcription factors to RNA polymerase [51,52]. Here, mediator components Glyma13g16910 (head MED20a) and Glyma13g31480 (tail MED16) were at least three times down- and nine times up-regulated after 4 and 28h of hypoxia, respectively, in both cultivars. The head module MED20a subunit participates in transcription regulation of miRNA and protein-coding genes involved with plant development, time flowering, and fruit size in A. thaliana [53]. In contrast, the tail module MED16 component regulates several biotic [54,55] and abiotic [5658] stress responses in plants. Med16 is required for transcriptional activation of cold- and dehydration-inducible genes that have C-repeat/dehydration-responsive elements (CRT/DRE) promoter motifs controlled by CRT/DRE-binding transcription factors (CBF/DREB) [57,58].

The CRT/DRE and CBF/DREB, the cis-acting ABA-responsive elements (ABRE) and corresponding trans-acting factors ABRE-binding proteins/ABRE-binding factors (AREB/ABF) play important roles in abiotic stress tolerance in plants [59]. Based on the genomatix database (http://www.genomatix.de/) [60], we found that cis-acting ABRE and DRE/CRT motifs were most frequent in the putative promoter region (1000 bp up- and 100 bp down-stream the TSS) of top up- compared to top down-regulated genes (Fig 5, Table B in S1 Supporting Information). The AREB/ABF genes Glyma06g04353 and Glyma19g30230 increased their mRNA abundance (~2 fc) after 28h of hypoxia. The higher differential expression of GmDREB1B;1 (Glyma20g29410 [61]) and of GmDREB2A;2 (Glyma14g06080 [62]) was observed after 0.5h (>20 fc) and 28h (>24 fc) of hypoxia, respectively. Both these genes are also induced by heat, cold, drought, and salinity stress, and up-regulate ABA receptor GmPYL21 (Glyma13g30210 [63]) [61,62]. GmPYL14 (Glyma14g06100, up-regulated by GmDREB2A;2 [62]), and GmPYL21 were up-regulated under hypoxia (>9 and >3 fc after 4h, respectively). These receptors interact with the phosphatase GmPP2C1 (Glyma13g16640, at most 3 fc under hypoxia), inhibiting it in an ABA independent manner [63]. Interestingly, Kidokoro et al. [61] observed that transcriptional activation of GmPYL21 by GmDREB1B;1 can enhance ABRE-mediated gene expression in an ABA-independent manner under cold stress, although ABA levels are not increased under such condition. We propose that ABRE-mediated gene expression is ABA-independent in hypoxic soybean roots. Although exogenous application of ABA improves flooding tolerance of plants [6466], endogenous ABA decreases in roots under flooding stress [6769]. We observed down-regulation of genes related to ABA biosynthesis (Glyma19g06540, Glyma06g08944, Glyma13g27220, Glyma11g21160, and Glyma14g04950) and up-regulation of ABA inactivation genes GmCYP707As (Glyma16g20490, Glyma01g35660, and Glyma09g35250 [70]). In addition to GmCYP707As, hypoxia changes expression of soybean orthologs of A. thaliana genes involved in ABA signaling (ABA receptors Pyl4-6; phosphatases ABI1-2, HAB1-2, HAI1-3, and AHG3; AFP1-4; MAPKKK18) dependent on SRK2D/E/I (SNF1-related kinases SRK2D/SnRK2.2, SRK2E/SnRK2.6, and SRK2I/SnRK2.3) (Genevestigator microarray database [71]) (S3 File). SRK2D/E/I are key activators of AREB/ABF proteins [72]. Glyma01g39020, which is most similar in amino acid sequence to SRK2D/E/I, is up-regulated 2- and 3-fold after 4 and 28h of hypoxia. Moreover, orthologous genes of CIPK-6 and CIPK-25 (Glyma17g08270, Glyma04g06520, and Glyma06g06550) were up-regulated at all-time points (from 4 to 40 fc) under hypoxia. These are kinases involved in Ca+2-mediated expression of DREB1-2 and ABA signaling [73,74]. Strikingly, although ABI1-2, PP2CA, and HAI1-2 genes from A. thaliana are down-regulated under hypoxic stress, we observed strong up-regulation of soybean orthologs of Glyma01g43460 (HAI3) transcripts at all-time points (S3 File). These A. thaliana genes are up-regulated under drought, osmotic, and salinity stresses (S3 File) and soybean genes under drought stress [18], stresses associated with ABA production. In agreement with a role for AREB in tolerance of diverse abiotic stresses, our results indicate although ABA level decreases in soybean roots under hypoxia, ABRE-mediated gene expression may occur. In this context, AREB/ABF are powerful candidates to improve flooding tolerance.

Top gene groups differed in structure, composition, and codon usage. How may hypoxia alter translation?

Compared to SGs, DEGs are smaller, with shorter CDS (coding sequence) length and fewer introns (Fig 6, Table C in S1 Supporting Information). Similar results were observed in A. thaliana [75], yeasts, and mice [76] for genes responsive to other stresses. Shorter genes with fewer introns demand less energy [77] and can have faster expression dynamics [78]. Interestingly, SGs also seem to have improvement of energetic and time costs. We compared soybean SGs with cognate proteins [79], and analyzed A. thaliana data sets of immunopurified polysomal mRNAs [80] and translationally inactive mRNAs [81]. The results suggest that energy is saved from translation by down-regulation of cognate soybean proteins under hypoxia (Fig C in S1 Supporting Information). In this context, stable mRNAs from A. thaliana are sequestered into stress granules and poorly associated with translating ribosomes under hypoxia. Upon reoxygenation, they are rapidly released from stress granules forming new polyribosomes, minimizing the need for de novo transcription.

thumbnail
Fig 6. Structural features among gene groups.

The groups are formed by top 40 and 500 ranked up- (U), down-regulated (D), and stable-genes (S) to hypoxia, and genome. The box is determined by the 25th and 75th percentiles with a line as the median and a black square as the mean of the data. Error bars extend 1.5 times the interquartile range from the 25th and 75th percentiles. Statistical significances from pairwise comparisons are provided in Tables B and C in S1 Supporting Information.

https://doi.org/10.1371/journal.pone.0187920.g006

The higher intron number and number of splicing variants in soybean SGs (Fig 6, Tables B and C in S1 Supporting Information) are in agreement with the higher gene body (i.e., transcript region) methylation in SGs [75,82,83], involved in splicing regulation [8486]. Here, whereas SGs exhibited steeper 5’ to 3’ negative G+C and CpG gradient, no decrement from start to middle region of DEGs were observed (Fig D and Table C in S1 Supporting Information). The low strength but diverged G+C and CpG patterns among gene groups can be influenced by gene structure, recombination, and DNA methylation [87]. Moreover, C3pG1 (C at third position of a codon binding G at first position of a neighbor codon) and G1 content were higher among di- and mono-nucleotides, respectively, and C3 was the mono-nucleotide with more diverged content among gene groups, mainly at CDS middle region.

The relative synonymous codon usage (RSCU) analysis for all 59 synonymous codons showed that C3 content divergences highlighted in 2-fold degenerate pyrimidine ending codons between CDS edges as well as among gene groups at CDS middle region (Fig 7). It is noteworthy that 2-fold degenerate codons ending in pyrimidines seem to be in majority [maybe exclusively, including in soybean (Fig E in S1 Supporting Information)] decoded by G34 tRNAs (G at position 34 of the tRNA, the first anticodon position) in all three domains of life (archaea, bacteria, and eukarya) [88,89]. This suggests strong positive selection to discriminate correct cognate C3 and wobble U3 codons from the incorrect near-cognate codons G3 and A3 (e.g., CAC/U histidine versus CAG/A glutamine). C3 over U3 2-fold degenerate ending codon bias also occurs at evolutionarily conserved amino acids sites across 12 fly drosophilid species and correspond to higher levels of G34-to-Q34 substitution [90]. This opens a question if Q34 tRNA influences the compositional divergence among soybean gene groups. Remarkably, queuine (q), the free base of Q (queuosine), is only synthesized by bacteria and salvaged by most eukaryotes [91]. Example is the legume model Medicago truncatula, in which rhizobium Q synthesis is required for effective nitrogen-fixing symbiosis [92]. In contrast, this is not observed in Brassicaceae, including A. thaliana, given the absence of genes encoding transglycosylases that catalyze q insertion in target G34U35N36 tRNAs (N = any base), found in Medicago [91] and soybean (Glyma04g10706, Glyma01g40041, Glyma13g10281, Glyma11g05250, Glyma06g10555, and Glyma08g48310). Besides Q, wybutosine (yW) is another hypermodified nucleoside derived from G, but yW is found exclusively at position 37 (neighboring anticodon sequence) of tRNAs that decode phenylalanine (codons UUU and UUC; biased here among soybean gene groups), important to reduce polyuridine translational frameshift errors [93].

thumbnail
Fig 7. Heat map of relative synonymous codon usage (RSCU) for start, middle, and end 30 CDS (coding sequence) codons ending in pyrimidine (C and T) from different soybean gene groups.

The groups are formed by top 40 and 500 ranked up- (U), down-regulated (D), and stable-genes (S) to hypoxia, and genome (G). Asp, Cys, Tyr, His, Phe, and Asn are coding by 2-fold degenerate pyrimidine ending codons. Heat map of RSCU for all 59 synonymous codons including for whole CDS extension is provided in Fig E in S1 Supporting Information.

https://doi.org/10.1371/journal.pone.0187920.g007

Here, DEGs related to modification of inosine, queuosine, wybutosine, pseudouridine, and methylation (tRNA and DNA methyltransferases), as well as aminoacylation of tRNAs (amino acid charging) were observed under hypoxia (S4 File). In addition, our results indicate that hypoxia impairs biosynthesis of methionine, mainly SAM, and consequently DNA and tRNA methylation (Fig A in S1 Supporting Information, S2 File). Although S-methylmethionine cycle seems to alleviate methionine decrement by transcriptional up-regulation of homocysteine S-methyltransferase (ec: 2.1.1.10), genes encoding folate-dependent methionine synthase (ec: 2.1.1.14) as well as SAM synthase (EC 2.5.1.6) and SAH hydrolase (EC 3.3.1.1) were down-regulated under hypoxic stress (Fig A in S1 Supporting Information, S2 File). This can also occur at protein level [79].

Differential expression of aminoacyl-tRNA synthetases (aaRS) also occur in soybean leaf under drought stress, with more aaRS DEGs in wild-type than in transgenic lines overexpressing BiP chaperone [94]. The BiP (Binding Protein) major regulates the endoplasmic reticulum (ER) stress [95]. Here, many ER stress related genes [96] were differentially expressed (S5 File), including down-regulation of BiP (Glyma05g36600, Glyma05g36620, Glyma08g02940, and Glyma08g02960; from -2 to -10 fc after 4 and 28h in both cultivars under hypoxia). Among these aaRS differentially expressed, Glyma14g11711 (aspartyl-tRNA synthetase) was up-regulated at all-time points in the two cultivars (from 5 to 16 fc). Curiously, transgenic plants overexpressing an aspartyl tRNA synthase (AspRS) orthologue (At4g31180) improve tolerance to biotic stress [97]. The At4g31180 is target of a synthetic isomer of GABA, called BABA (β-Amino Butyric Acid) [97]. BABA primes plants to enhance tolerance against broad-spectrum of biotic and abiotic stresses [98]. Inhibition of AspRS activity by BABA accumulates uncharged tRNA Asp, which as others uncharged tRNAs is signaling molecule to attenuate translation by Gcn2-eIF2α system, and subsequently alleviate ER stress [99]. Based on this, further studies may help to elucidate involvement of amino acids metabolism in tRNA modifications, translation accuracy/efficiency, and ER stress under hypoxia. In addition, BiP and aaRS are candidates for biotechnology applications for improvement in flooding tolerance.

Supporting information

S1 Supporting Information.

File containing all supporting Tables (A-D) and Figures (A-E). Table A in S1 Supporting Information. Samples sizes used for structural and compositional analysis among groups of genes. Table B in S1 Supporting Information. Fisher’s pairwise comparisons between gene groups. Table C in S1 Supporting Information. Mann-Whitney pairwise comparisons. Table D in S1 Supporting Information. Primer pairs used in the study. Fig A in S1 Supporting Information. KEGG pathways showing transcriptional changes for amino acids and derivative metabolic process in both Embrapa 45 and BR 4 soybean cultivars. Fig B in S1 Supporting Information. Structural divergence of non-symbiotic hemoglobin Glyma11g12980.1 transcript between BR 4 and Embrapa 45 cultivars by mapping (A) and de novo assembly (B) of reads. Fig C in S1 Supporting Information. Steady-state stable mRNAs tend to be translationally down-regulated (A), decreasing their association with ribosomes (of which include RPL18 component) by their sequestration into stress granules (ribonucleoprotein complexes including UBP1C) under hypoxia, whereas upon reoxygenation they are rapidly released from stress granules forming new polyribosomes (B). Fig D in S1 Supporting Information. Compositional features for full (normalized to 90 nucleotides), start, middle, and end 90 nucleotides length from coding sequences from different soybean gene groups. Fig E in S1 Supporting Information. Heat map of RSCU (Relative Synonymous Codon Usage) for full, first, middle, and last 30 CDS (coding sequence) codons length from different soybean gene groups.

https://doi.org/10.1371/journal.pone.0187920.s001

(PDF)

S1 File. Tab-delimited file showing raw data for gene expression from flood-tolerant Embrapa 45 (E) and flood-sensitive BR 4 (B) soybean cultivars under normoxia (N) and hypoxia (H) along 0.5, 4, and 28h.

Top 40 and 500 ranked up- (U), down-regulated (D), and stable-genes (S) are shown in separate excel sheets named accordingly.

https://doi.org/10.1371/journal.pone.0187920.s002

(XLSX)

S2 File. Differentially expressed genes (DEGs) related with amino acids and derivative metabolic process in hypoxic roots of flood-tolerant Embrapa 45 (E) and flood-sensitive BR 4 (B) soybean cultivars along 0.5, 4, and 28h.

DEGs: fold-change ≥ 2 (up), ≤ -2 (down); adj. p ≤ 0.01; RPM ≥ 9 (control or stress datasets). RPM ≥ 9: at least 20 reads in all datasets.

https://doi.org/10.1371/journal.pone.0187920.s003

(XLSX)

S3 File. Differentially expressed genes (DEGs) from hypoxic roots of flood-tolerant Embrapa 45 (E) and flood-sensitive BR 4 (B) soybean cultivars along 0.5, 4, and 28h, and Arabidopsis thaliana orthologs related with ABA metabolism and signaling dependent of SRK2D/E/I responding to diverse treatments (from Genevestigator microarray database) [71].

https://doi.org/10.1371/journal.pone.0187920.s004

(XLSX)

S4 File. Differentially expressed genes (DEGs) related with modification of I34, Q34, yW37, pseudouridine, and methylation (tRNA and DNA methyltransferases), as well as aminoacylation of tRNAs (amino acid charging) in hypoxic roots of flood-tolerant Embrapa 45 (E) and flood-sensitive BR 4 (B) soybean cultivars along 0.5, 4, and 28h.

https://doi.org/10.1371/journal.pone.0187920.s005

(XLSX)

S5 File. Differentially expressed genes (DEGs) related with endoplasmic reticulum stress [unfolded protein response (UPR) and ER stress-induced plant-specific cell death signaling pathways] [96] analyzed in hypoxic soybean roots of flood-tolerant Embrapa 45 (E) and flood-sensitive BR 4 (B) soybean cultivars along 0.5, 4, and 28h.

DEGs: fold-change ≥ 2 (up), ≤ -2 (down); adj. p ≤ 0.01; RPM ≥ 9 (control or stress datasets). RPM ≥ 9: at least 20 reads in all datasets.

https://doi.org/10.1371/journal.pone.0187920.s006

(XLSX)

Acknowledgments

We thank Andrei Stecca Steindorff for plotting BWA-MEM results and Maria Cristina Neves de Oliveira for statistical discussions.

References

  1. 1. Silvertown J, Dodd M, Gowing DJ, Mountford JO. Hydrologically defined niches reveal a basis for species richness in plant communities. Nature. 1999;400: 61–63. Available: http://oro.open.ac.uk/27574/
  2. 2. Bailey-Serres J, Lee SC, Brinton E. Waterproofing Crops: Effective Flooding Survival Strategies. Plant Physiol. 2012;160: 1698–1709. pmid:23093359
  3. 3. Bailey-Serres J, Fukao T, Gibbs DJ, Holdsworth MJ, Lee SC, Licausi F, et al. Making sense of low oxygen sensing. Trends Plant Sci. Elsevier Ltd; 2012;17: 129–138. pmid:22280796
  4. 4. Armstrong W. Aeration in Higher Plants. Adv Bot Res. 1980;7: 225–332.
  5. 5. Bailey-Serres J, Voesenek LACJ. Flooding Stress: Acclimations and Genetic Diversity. Annu Rev Plant Biol. 2008;59: 313–339. pmid:18444902
  6. 6. Geigenberger P. Response of plant metabolism to too little oxygen. Curr Opin Plant Biol. 2003;6: 247–256. pmid:12753974
  7. 7. Mickelbart M V, Hasegawa PM, Bailey-Serres J. Genetic mechanisms of abiotic stress tolerance that translate to crop yield stability. Nat Rev Genet. Nature Publishing Group; 2015;16: 237–251. pmid:25752530
  8. 8. Hattori Y, Nagai K, Furukawa S, Song X, Kawano R, Sakakibara H, et al. The ethylene response factors SNORKEL1 and SNORKEL2 allow rice to adapt to deep water. Nature. 2009;460: 1026–1030. pmid:19693083
  9. 9. Xu K, Xu X, Fukao T, Canlas P, Maghirang-Rodriguez R, Heuer S, et al. Sub1A is an ethylene-response-factor-like gene that confers submergence tolerance to rice. Nature. 2006;442: 705–708. pmid:16900200
  10. 10. Tamang BG, Magliozzi JO, Maroof MAS, Fukao T. Physiological and transcriptomic characterization of submergence and reoxygenation responses in soybean seedlings. Plant, Cell Environ. 2014;37: 2350–2365. pmid:24433575
  11. 11. Vantoai TT, Martin SKS, Chase K, Boru G, Schnipke V, Schmitthenner AF, et al. Identification of a QTL associated with tolerance of soybean to soil waterlogging. Crop Sci. 2001;41: 1247–1252.
  12. 12. Nguyen VT, Vuong TD, VanToai T, Lee JD, Wu X, Rouf Mian MA, et al. Mapping of quantitative trait loci associated with resistance to Phytophthora sojae and flooding tolerance in soybean. Crop Sci. 2012;52: 2481–2493.
  13. 13. Chen W, Yao Q, Patil GB, Agarwal G, Deshmukh RK, Lin L, et al. Identification and Comparative Analysis of Differential Gene Expression in Soybean Leaf Tissue under Drought and Flooding Stress Revealed by RNA-Seq. Front Plant Sci. 2016;7: 1–19.
  14. 14. Yin X, Hiraga S, Hajika M, Nishimura M, Komatsu S. Transcriptomic analysis reveals the flooding tolerant mechanism in flooding tolerant line and abscisic acid treated soybean. Plant Mol Biol. Springer Netherlands; 2016;0: 0. pmid:28012053
  15. 15. Nanjo Y, Maruyama K, Yasue H, Yamaguchi-Shinozaki K, Shinozaki K, Komatsu S. Transcriptional responses to flooding stress in roots including hypocotyl of soybean seedlings. Plant Mol Biol. 2011;77: 129–144. pmid:21656040
  16. 16. Komatsu S, Yamamoto R, Nanjo Y, Mikami Y, Yunokawa H, Sakata K. A Comprehensive Analysis of the Soybean Genes and Proteins Expressed under Flooding Stress using Transcriptome and Proteome Techniques. J Proteome Res. 2009;8: 4766–4778. pmid:19658438
  17. 17. Nakayama TJ, Rodrigues FA, Neumaier N, Marcelino-Guimarães FC, Farias JRB, de Oliveira MCN, et al. Reference genes for quantitative real-time polymerase chain reaction studies in soybean plants under hypoxic conditions. Genet Mol Res. 2014;13: 860–871. pmid:24615050
  18. 18. Rodrigues FA, Fuganti-Pagliarini R, Marcolino-Gomes J, Nakayama TJ, Molinari HBC, Lobo FP, et al. Daytime soybean transcriptome fluctuations during water deficit stress. BMC Genomics. BMC Genomics; 2015;16: 505. pmid:26149272
  19. 19. Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. Nature Publishing Group; 2010;463: 178–183. pmid:20075913
  20. 20. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29: 644–52. pmid:21572440
  21. 21. Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010;26: 589–595. pmid:20080505
  22. 22. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26: 139–140. pmid:19910308
  23. 23. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. 1995. pp. 289–300.
  24. 24. Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38: W64–W70. pmid:20435677
  25. 25. Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012;40: D109–D114. pmid:22080510
  26. 26. Chan PP, Lowe TM. GtRNAdb: a database of transfer RNA genes detected in genomic sequence. Nucleic Acids Res. 2009;37: D93–D97. pmid:18984615
  27. 27. Van Dongen JT, Fröhlich A, Ramírez-Aguilar SJ, Schauer N, Fernie AR, Erban A, et al. Transcript and metabolite profiling of the adaptive response to mild decreases in oxygen concentration in the roots of arabidopsis plants. Ann Bot. 2009;103: 269–280. pmid:18660497
  28. 28. Kreuzwieser J, Hauberg J, Howell KA, Carroll A, Rennenberg H, Millar AH, et al. Differential response of gray poplar leaves and roots underpins stress adaptation during hypoxia. Plant Physiol. 2009;149: 461–73. pmid:19005089
  29. 29. Oliveira HC, Sodek L. Effect of oxygen deficiency on nitrogen assimilation and amino acid metabolism of soybean root segments. Amino Acids. 2013;44: 743–755. pmid:22990842
  30. 30. Rocha M, Licausi F, Araújo WL, Nunes-Nesi A, Sodek L, Fernie AR, et al. Glycolysis and the tricarboxylic acid cycle are linked by alanine aminotransferase during hypoxia induced by waterlogging of Lotus japonicus. Plant Physiol. 2010;152: 1501–1513. pmid:20089769
  31. 31. Morot-Gaudry J-F, Job D, Lea PJ. Amino Acid Metabolism. In: Morot-Gaudry J-F, Lea PJ, editors. Plant Nitrogen. Berlin, Heidelberg: Springer Berlin Heidelberg; 2001. pp. 167–211. https://doi.org/10.1007/978-3-662-04064-5_7
  32. 32. Li W, Han Y, Tao F, Chong K. Knockdown of SAMS genes encoding S-adenosyl-l-methionine synthetases causes methylation alterations of DNAs and histones and leads to late flowering in rice. J Plant Physiol. 2011;168: 1837–1843. pmid:21757254
  33. 33. Hori H. Methylated nucleosides in tRNA and tRNA methyltransferases. Front Genet. 2014;5: 1–26.
  34. 34. Minocha R, Majumdar R, Minocha SC. Polyamines and abiotic stress in plants: a complex relationship. Front Plant Sci. 2014;5: 175 (11–17). pmid:24847338
  35. 35. Singh DK, Maximova SN, Jensen PJ, Lehman BL, Ngugi HK, McNellis TW. FIBRILLIN4 Is Required for Plastoglobule Development and Stress Resistance in Apple and Arabidopsis. Plant Physiol. 2010;154: 1281–1293. pmid:20813909
  36. 36. Mutava RN, Prince SJK, Syed NH, Song L, Valliyodan B, Chen W, et al. Understanding abiotic stress tolerance mechanisms in soybean: A comparative evaluation of soybean response to drought and flooding stress. Plant Physiol Biochem. Elsevier Masson SAS; 2015;86: 109–120. pmid:25438143
  37. 37. Gupta KJ, Igamberdiev AU. Reactive Nitrogen Species in Mitochondria and Their Implications in Plant Energy Status and Hypoxic Stress Tolerance. Front Plant Sci. 2016;7: 1–6.
  38. 38. Lamers LPM, Govers LL, Janssen ICJM, Geurts JJM, Van der Welle MEW, Van Katwijk MM, et al. Sulfide as a soil phytotoxin—a review. Front Plant Sci. 2013;4: 1–14.
  39. 39. Laanbroek HJ. Bacterial cycling of minerals that affect plant growth in waterlogged soils: a review. Aquat Bot. 1990;38: 109–125.
  40. 40. Shabala S, Shabala L, Barcelo J, Poschenrieder C. Membrane transporters mediating root signalling and adaptive responses to oxygen deprivation and soil flooding. Plant Cell Environ. 2014;37: 2216–33. pmid:24689809
  41. 41. Nikolov DB, Burley SK. RNA polymerase II transcription initiation: A structural view. Proc Natl Acad Sci. 1997;94: 15–22. pmid:8990153
  42. 42. Fessele S, Maier H, Zischek C, Nelson PJ, Werner T. Regulatory context is a crucial part of gene function. Trends Genet. 2002;18: 60–63. pmid:11818130
  43. 43. Vannini A, Cramer P. Conservation between the RNA Polymerase I, II, and III Transcription Initiation Machineries. Mol Cell. Elsevier Inc.; 2012;45: 439–446. pmid:22365827
  44. 44. Bonham AJ, Neumann T, Tirrell M, Reich NO. Tracking transcription factor complexes on DNA using total internal reflectance fluorescence protein binding microarrays. Nucleic Acids Res. 2009;37: e94–e94. pmid:19487241
  45. 45. Tirosh I, Weinberger A, Carmi M, Barkai N. A genetic signature of interspecies variations in gene expression. Nat Genet. 2006;38: 830–834. pmid:16783381
  46. 46. Tora L, Timmers HTM. The TATA box regulates TATA-binding protein (TBP) dynamics in vivo. Trends Biochem Sci. Elsevier Ltd; 2010;35: 309–314. pmid:20176488
  47. 47. Kim JL, Burley SK. 1.9 Å resolution refined structure of TBP recognizing the minor groove of TATAAAAG. Nat Struct Biol. 1994;1: 638–653. pmid:7634103
  48. 48. Maruyama K, Todaka D, Mizoi J, Yoshida T, Kidokoro S, Matsukura S, et al. Identification of Cis-Acting Promoter Elements in Cold- and Dehydration-Induced Transcriptional Pathways in Arabidopsis, Rice, and Soybean. DNA Res. 2012;19: 37–49. pmid:22184637
  49. 49. van Werven FJ, van Teeffelen HAAM, Holstege FCP, Timmers HTM. Distinct promoter dynamics of the basal transcription factor TBP across the yeast genome. Nat Struct Mol Biol. Nature Publishing Group; 2009;16: 1043–8. pmid:19767748
  50. 50. Paul E, Zhu ZI, Landsman D, Morse RH. Genome-wide association of mediator and RNA polymerase II in wild-type and mediator mutant yeast. Mol Cell Biol. 2015;35: 331–42. pmid:25368384
  51. 51. Eyboulet F, Wydau-Dematteis S, Eychenne T, Alibert O, Neil H, Boschiero C, et al. Mediator independently orchestrates multiple steps of preinitiation complex assembly in vivo. Nucleic Acids Res. 2015;43: 9214–9231. pmid:26240385
  52. 52. Yang Y, Li L, Qu L-J. Plant Mediator complex and its critical functions in transcription regulation. J Integr Plant Biol. 2016;58: 106–118. pmid:26172375
  53. 53. Kim YJ, Zheng B, Yu Y, Won SY, Mo B, Chen X. The role of Mediator in small and long noncoding RNA production in Arabidopsis thaliana. EMBO J. 2011;30: 814–822. pmid:21252857
  54. 54. Zhang X, Wang C, Zhang Y, Sun Y, Mou Z. The Arabidopsis Mediator Complex Subunit16 Positively Regulates Salicylate-Mediated Systemic Acquired Resistance and Jasmonate/Ethylene-Induced Defense Pathways. Plant Cell. 2012;24: 4294–4309. pmid:23064320
  55. 55. Wang C, Yao J, Du X, Zhang Y, Sun Y, Rollins JA, et al. The Arabidopsis Mediator Complex Subunit16 Is a Key Component of Basal Resistance against the Necrotrophic Fungal Pathogen Sclerotinia sclerotiorum. Plant Physiol. 2015;169: 856–872. pmid:26143252
  56. 56. Wathugala DL, Hemsley PA, Moffat CS, Cremelie P, Knight MR, Knight H. The Mediator subunit SFR6/MED16 controls defence gene expression mediated by salicylic acid and jasmonate responsive pathways. New Phytol. 2012;195: 217–230. pmid:22494141
  57. 57. Boyce JM, Knight H, Deyholos M, Openshaw MR, Galbraith DW, Warren G, et al. The sfr6 mutant of Arabidopsis is defective in transcriptional activation via CBF/DREB1 and DREB2 and shows sensitivity to osmotic stress. Plant J. 2003;34: 395–406. pmid:12753580
  58. 58. Knight H, Mugford SG, Ülker B, Gao D, Thorlby G, Knight MR. Identification of SFR6, a key component in cold acclimation acting post-translationally on CBF function. Plant J. 2009;58: 97–108. pmid:19067974
  59. 59. Nakashima K, Yamaguchi-Shinozaki K, Shinozaki K. The transcriptional regulatory network in the drought response and its crosstalk in abiotic stress responses including drought, cold, and heat. Front Plant Sci. 2014;5: 170. pmid:24904597
  60. 60. Cartharius K, Frech K, Grote K, Klocke B, Haltmeier M, Klingenhoff A, et al. MatInspector and beyond: Promoter analysis based on transcription factor binding sites. Bioinformatics. 2005;21: 2933–2942. pmid:15860560
  61. 61. Kidokoro S, Watanabe K, Ohori T, Moriwaki T, Maruyama K, Mizoi J, et al. Soybean DREB1/CBF-type transcription factors function in heat and drought as well as cold stress-responsive gene expression. Plant J. 2015;81: 505–518. pmid:25495120
  62. 62. Mizoi J, Ohori T, Moriwaki T, Kidokoro S, Todaka D, Maruyama K, et al. GmDREB2A;2, a Canonical DEHYDRATION-RESPONSIVE ELEMENT-BINDING PROTEIN2-Type Transcription Factor in Soybean, Is Posttranslationally Regulated and Mediates Dehydration-Responsive Element-Dependent Gene Expression. Plant Physiol. 2013;161: 346–361. pmid:23151346
  63. 63. Bai G, Yang D-H, Zhao Y, Ha S, Yang F, Ma J, et al. Interactions between soybean ABA receptors and type 2C protein phosphatases. Plant Mol Biol. 2013;83: 651–664. pmid:23934343
  64. 64. Ellis MH, Dennis ES, Peacock WJ. Arabidopsis roots and shoots have different mechanisms for hypoxic stress tolerance. Plant Physiol. 1999;119: 57–64. pmid:9880346
  65. 65. Hwang S, Vantoai TT. Abscisic Acid Induces Anaerobiosis Tolerance in Corn. Plant Physiol. 1991;97: 593–597. Available: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1081048/# pmid:16668440
  66. 66. Komatsu S, Han C, Nanjo Y, Altaf-Un-Nahar M, Wang K, He D, et al. Label-Free Quantitative Proteomic Analysis of Abscisic Acid Effect in Early-Stage Soybean under Flooding. J Proteome Res. 2013;12: 4769–4784. pmid:23808807
  67. 67. Argamasilla R, Gómez-Cadenas A, Arbona V. Metabolic and Regulatory Responses in Citrus Rootstocks in Response to Adverse Environmental Conditions. J Plant Growth Regul. 2014;33: 169–180.
  68. 68. Arbona V, Gómez-Cadenas A. Hormonal Modulation of Citrus Responses to Flooding. J Plant Growth Regul. 2008;27: 241–250.
  69. 69. Hsu F-C, Chou M-Y, Peng H-P, Chou S-J, Shih M-C. Insights into hypoxic systemic responses based on analyses of transcriptional regulation in Arabidopsis. Shiu S-H, editor. PLoS One. 2011;6: e28888. pmid:22194941
  70. 70. Zheng Y, Huang Y, Xian W, Wang J, Liao H. Identification and expression analysis of the Glycine max CYP707A gene family in response to drought and salt stresses. Ann Bot. 2012;110: 743–756. pmid:22751653
  71. 71. Hruz T, Laule O, Szabo G, Wessendorp F, Bleuler S, Oertle L, et al. Genevestigator V3: A Reference Expression Database for the Meta-Analysis of Transcriptomes. Adv Bioinformatics. 2008;2008: 1–5. pmid:19956698
  72. 72. Fujita Y, Yoshida T, Yamaguchi-Shinozaki K. Pivotal role of the AREB/ABF-SnRK2 pathway in ABRE-mediated transcription in response to osmotic stress in plants. Physiol Plant. 2013;147: 15–27. pmid:22519646
  73. 73. He L, Yang X, Wang L, Zhu L, Zhou T, Deng J, et al. Molecular cloning and functional characterization of a novel cotton CBL-interacting protein kinase gene (GhCIPK6) reveals its involvement in multiple abiotic stress tolerance in transgenic plants. Biochem Biophys Res Commun. Elsevier Inc.; 2013;435: 209–215. pmid:23660187
  74. 74. Meena MK, Ghawana S, Dwivedi V, Roy A, Chattopadhyay D. Expression of chickpea CIPK25 enhances root growth and tolerance to dehydration and salt stress in transgenic tobacco. Front Plant Sci. 2015;6: 1–11.
  75. 75. Aceituno FF, Moseyko N, Rhee SY, Gutiérrez RA. The rules of gene expression in plants: Organ identity and gene body methylation are key factors for regulation of gene expression in Arabidopsis thaliana. BMC Genomics. 2008;9: 438. pmid:18811951
  76. 76. Jeffares DC, Penkett CJ, Bähler J. Rapidly regulated genes are intron poor [Internet]. Trends in Genetics. 2008. pp. 375–378. pmid:18586348
  77. 77. Lynch M, Marinov GK. The bioenergetic costs of a gene. Proc Natl Acad Sci. 2015;112: 201514974. pmid:26575626
  78. 78. Bentley DL. Coupling mRNA processing with transcription in time and space. Nat Rev Genet. Nature Publishing Group; 2014;15: 163–175. pmid:24514444
  79. 79. Wang X, Oh M, Sakata K, Komatsu S. Gel-free/label-free proteomic analysis of root tip of soybean over time under flooding and drought stresses. J Proteomics. Elsevier B.V.; 2016;130: 42–55. pmid:26376099
  80. 80. Branco-Price C, Kaiser KA, Jang CJH, Larive CK, Bailey-Serres J. Selective mRNA translation coordinates energetic and metabolic adjustments to cellular oxygen deprivation and reoxygenation in Arabidopsis thaliana. Plant J. 2008;56: 743–755. pmid:18665916
  81. 81. Sorenson R, Bailey-Serres J. Selective mRNA sequestration by OLIGOURIDYLATE-BINDING PROTEIN 1 contributes to translational control during hypoxia in Arabidopsis. Proc Natl Acad Sci. 2014;111: 2373–2378. pmid:24469793
  82. 82. Elhaik E, Pellegrini M, Tatarinova T V. Gene expression and nucleotide composition are associated with genic methylation level in Oryza sativa. BMC Bioinformatics. 2014;15: 23. pmid:24447369
  83. 83. Tatarinova T, Elhaik E, Pellegrini M. Cross-Species Analysis of Genic GC3 Content and DNA Methylation Patterns. Genome Biol Evol. 2013;5: 1443–1456. pmid:23833164
  84. 84. Shukla S, Oberdoerffer S. Co-transcriptional regulation of alternative pre-mRNA splicing. Biochim Biophys Acta—Gene Regul Mech. Elsevier B.V.; 2012;1819: 673–683. pmid:22326677
  85. 85. Luco RF, Allo M, Schor IE, Kornblihtt AR, Misteli T. Epigenetics in Alternative Pre-mRNA Splicing. Cell. 2011;144: 16–26. pmid:21215366
  86. 86. Takuno S, Gaut BS. Gene body methylation is conserved between plant orthologs and is of evolutionary consequence. Proc Natl Acad Sci. 2013;110: 1797–1802. pmid:23319627
  87. 87. Glémin S, Clément Y, David J, Ressayre A. GC content evolution in coding regions of angiosperm genomes: a unifying hypothesis. Trends Genet. 2014;30: 263–270. pmid:24916172
  88. 88. Novoa EM, Pavon-Eternod M, Pan T, Ribas de Pouplana L. A Role for tRNA Modifications in Genome Structure and Codon Usage. Cell. 2012;149: 202–213. pmid:22464330
  89. 89. Grosjean H, de Crécy-Lagard V, Marck C. Deciphering synonymous codons in the three domains of life: Co-evolution with specific tRNA modification enzymes. FEBS Lett. Federation of European Biochemical Societies; 2010;584: 252–264. pmid:19931533
  90. 90. Zaborske JM, Bauer DuMont VL, Wallace EWJ, Pan T, Aquadro CF, Drummond DA. A Nutrient-Driven tRNA Modification Alters Translational Fidelity and Genome-wide Protein Coding across an Animal Genus. Malik HS, editor. PLoS Biol. 2014;12: e1002015. pmid:25489848
  91. 91. Zallot R, Brochier-Armanet C, Gaston KW, Forouhar F, Limbach P a, Hunt JF, et al. Plant, Animal, and Fungal Micronutrient Queuosine Is Salvaged by Members of the DUF2419 Protein Family. ACS Chem Biol. 2014;9: 1812–1825. pmid:24911101
  92. 92. Marchetti M, Capela D, Poincloux R, Benmeradi N, Auriac M-C, Le Ru A, et al. Queuosine Biosynthesis Is Required for Sinorhizobium meliloti-Induced Cytoskeletal Modifications on HeLa Cells and Symbiosis with Medicago truncatula. Cascales E, editor. PLoS One. 2013;8: e56043. pmid:23409119
  93. 93. Jackman JE, Alfonzo JD. Transfer RNA modifications: nature’s combinatorial chemistry playground. Wiley Interdiscip Rev RNA. 2013;4: 35–48. pmid:23139145
  94. 94. Carvalho HH, Brustolini OJB, Pimenta MR, Mendes GC, Gouveia BC, Silva PA, et al. The Molecular Chaperone Binding Protein BiP Prevents Leaf Dehydration-Induced Cellular Homeostasis Disruption. Tran L-SP, editor. PLoS One. 2014;9: e86661. pmid:24489761
  95. 95. Reis PAA, Rosado GL, Silva LAC, Oliveira LC, Oliveira LB, Costa MDL, et al. The binding protein BiP attenuates stress-induced cell death in soybean via modulation of the N-rich protein-mediated signaling pathway. Plant Physiol. 2011;157: 1853–65. pmid:22007022
  96. 96. Silva PA, Silva JCF, Caetano HDN, Machado JPB, Mendes GC, Reis PAB, et al. Comprehensive analysis of the endoplasmic reticulum stress response in the soybean genome: conserved and plant-specific features. BMC Genomics. BMC Genomics; 2015;16: 783. pmid:26466891
  97. 97. Luna E, van Hulten M, Zhang Y, Berkowitz O, López A, Pétriacq P, et al. Plant perception of β-aminobutyric acid is mediated by an aspartyl-tRNA synthetase. Nat Chem Biol. 2014;10: 450–456. pmid:24776930
  98. 98. Cohen YR. β-Aminobutyric Acid-Induced Resistance Against Plant Pathogens. Plant Dis. 2002;86: 448–457.
  99. 99. Kørner C, Du X, Vollmer M, Pajerowska-Mukhtar K. Endoplasmic Reticulum Stress Signaling in Plant Immunity—At the Crossroad of Life and Death. Int J Mol Sci. 2015;16: 26582–26598. pmid:26556351