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

Transcriptome-Wide Identification of Reference Genes for Expression Analysis of Soybean Responses to Drought Stress along the Day

  • Juliana Marcolino-Gomes,

    Affiliations Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil, Department of Biology, State University of Londrina, Londrina, Paraná, Brazil

  • Fabiana Aparecida Rodrigues,

    Affiliation Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil

  • Renata Fuganti-Pagliarini,

    Affiliation Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil

  • Thiago Jonas Nakayama,

    Affiliations Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil, Department of Crop Science, Federal University of Viçosa, Viçosa, Minas Gerais, Brazil

  • Rafaela Ribeiro Reis,

    Affiliations Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil, Department of Biology, State University of Londrina, Londrina, Paraná, Brazil

  • Jose Renato Bouças Farias,

    Affiliation Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil

  • Frank G. Harmon,

    Affiliations Plant Gene Expression Center, ARS/USDA, Albany, California, United States of America, Department of Plant and Microbial Biology, University of California-Berkeley, Berkeley, California, United States of America

  • Hugo Bruno Correa Molinari,

    Affiliation Embrapa Agroenergy, Brazilian Agricultural Research Corporation, Brasília, Brazil

  • Mayla Daiane Correa Molinari,

    Affiliations Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil, Department of Biology, State University of Londrina, Londrina, Paraná, Brazil

  • Alexandre Nepomuceno

    alexandre.nepomuceno@embrapa.br

    Affiliation Embrapa Soybean, Brazilian Agricultural Research Corporation, Londrina, Paraná, Brazil

Abstract

The soybean transcriptome displays strong variation along the day in optimal growth conditions and also in response to adverse circumstances, like drought stress. However, no study conducted to date has presented suitable reference genes, with stable expression along the day, for relative gene expression quantification in combined studies on drought stress and diurnal oscillations. Recently, water deficit responses have been associated with circadian clock oscillations at the transcription level, revealing the existence of hitherto unknown processes and increasing the demand for studies on plant responses to drought stress and its oscillation during the day. We performed data mining from a transcriptome-wide background using microarrays and RNA-seq databases to select an unpublished set of candidate reference genes, specifically chosen for the normalization of gene expression in studies on soybean under both drought stress and diurnal oscillations. Experimental validation and stability analysis in soybean plants submitted to drought stress and sampled during a 24 h timecourse showed that four of these newer reference genes (FYVE, NUDIX, Golgin-84 and CYST) indeed exhibited greater expression stability than the conventionally used housekeeping genes (ELF1-β and β-actin) under these conditions. We also demonstrated the effect of using reference candidate genes with different stability values to normalize the relative expression data from a drought-inducible soybean gene (DREB5) evaluated in different periods of the day.

Introduction

As sessile organisms, plants must endure environmental changes during the day and across seasons. These environmental oscillations strongly affect light, temperature, nutrient and water availability, acting as a powerful selective pressure that have shaped adaptive mechanisms in plants during their evolutionary history. As a result, these organisms have developed a complex molecular network that confers adaptive advantages by coordinating their metabolism with predictable daily and seasonal changes, known as the circadian clock [1].

The circadian clock is composed of a core of interconnected transcriptional–translational feedback loops, which are entrained by signals such as light and temperature to adjust metabolism to the environment. In plants, the clock controls a number of physiological and developmental processes. For example, the expression of chlorophyll biosynthesis genes is regulated by the circadian clock to peak at the end of the night, which is an important mechanism to ensure photosynthesis in subsequent light periods of the day, whereas the products of photosynthesis modulate the rhythm [2]. The circadian clock also allows plants to coordinate flowering with favorable seasons to increase their fitness [1], as well as it controls the rate of starch degradation [3] and nitrogen assimilation and utilization pathways [4].

In addition to normal day/night variations, plants are subject to other environmental variations via biotic and abiotic stresses. Among the abiotic stresses, drought stands out as the factor with the greatest impact on yield of important crops worldwide, including soybean. Different mechanisms are employed by plants to protect themselves against water deficits, including changes in stomatal conductance [5], osmotic adjustment [6], the accumulation of osmoprotectant molecules [7], and the activity of antioxidant proteins [8]. Because the circadian clock is known to improve organism fitness according to environmental conditions, a significant number of studies addressing the relationships between water deficit stress and the circadian clock have been conducted, providing consistent evidence of this interaction [913].

The metabolic and physiological adjustments performed in response to drought stresses usually involve the reconfiguration of the transcriptome [12,13], and therefore the analysis of gene expression in response to water deficits during the day is an interesting strategy [10]. One of the most sensitive methods for the quantification of gene expression is the fluorescence-based quantitative real-time PCR (RT-qPCR), which is increasingly being used. The advantages of this technique include its practical simplicity combined with the possibility of measuring small amounts of RNA in a wide range of samples, rapidly and with high specificity.

Thus, RT-qPCR is an important tool that allows the relative quantification of transcript abundance and can therefore be used to evaluate gene expression responses to environmental changes, such as diurnal oscillations and abiotic stresses, including drought. However, because most of the quantitative RNA data obtained are not absolute, but relative, accurate quantification of gene expression relies on the use of appropriate reference genes. These genes should be stably expressed, showing a transcript abundance that is strongly correlated with the total mRNA present in the samples to allow the normalization of gene expression data [14]. Normalization is a key step in RT-qPCR analysis, as it reduces/eliminates variations due to variations in RNA extraction, reverse transcription yields or amplification efficiency, allowing comparisons of mRNA concentrations across different samples, playing a critical role in the accurate quantification of relative gene expression [14]. Although several genes have been indicated as good references, it is known that even housekeeping genes may exhibit altered expression in response to experimental treatments, sampling times and the life cycle [1518].

In this context, a reference gene must be experimentally validated for specific tissues, genotypes and experimental designs. The soybean genes TUA (Glyma08g12140), TUB (Glyma03g27970), ELF1-β (Glyma13g04050), β-actin (Glyma15g05570) and GAPDH (Glyma06g01850) have been widely used as references in gene expression studies on drought responses [15,19]. On the other hand, isopentenyl diphosphate (IPP2), actin and ubiquitin are the most commonly used reference genes in studies investigating circadian/diurnal oscillations [2025]. Thus, no study conducted to date has evaluated the expression stability of reference genes for the study of both water deficit stress and circadian oscillations in soybean. Hence, in this study, after evaluating gene expression in response to drought during the day, we present a novel set of reference genes suitable for the normalization of relative expression data from combined studies on water deficit and diurnal oscillations.

Material and Methods

Selection of reference genes using the RefGenes tool

To evaluate the stability of genes expressed in response to drought during the day, we used the RefGenes tool from the Genevestigator platform [26], available at [https://www.genevestigator.com/gv/plant.jsp]. The Genevestigator platform provides a database of normalized and well-annotated microarray experiments, allowing asses the transcriptome of several organisms; the RefGenes tool enables searching for genes with minimal expression variance across a chosen set of arrays at the Genevestigator platform. For the purposes of this study we performed analysis of gene expression variance in 59 microarray libraries from soybean subjected to drought, heat and distinct light periods. To select the candidate reference genes presenting range of expression levels detected by RT-qPCR, we uploaded a list of the Gene Models (“Glyma”) from reference genes commonly used for gene expression normalization in soybean, previously described by Hu and colleagues (2009) [19], presented in Table 1.

thumbnail
Table 1. Commonly used RT-qPCR reference genes from soybean, according to Hu and colleagues (2009).

https://doi.org/10.1371/journal.pone.0139051.t001

Selection of reference genes using RNA-seq libraries

Under a second approach, we evaluated the expression stability of genes from 36 RNA-seq libraries. These libraries were synthesized from leaves of the drought-sensitive soybean genotype BR16, subjected to moderate drought stress on V2 developmental stage, sampled over a 24 h timecourse, with 4h intervals [10]. These data are deposited in the NCBI’s Gene Expression Omnibus [GEO; http://www.ncbi.nlm.nih.gov/geo/] repository and are accessible through GEO Series accession number GSE69469 (geospiza.com/Products/AnalysisEdition.shtml). To compare gene expression between different times and conditions, we log2-transformed the normalized reads per mapped million (RPM) value.

In this analysis, we selected genes that exhibited minimal expression variance across the libraries, presenting Coefficient of variation lower than 5%, with a range of expression similar to that of commonly used soybean reference genes[19], presented in Table 1.

Primer design

Primers for the new candidate reference genes were designed based on soybean Gene Model sequences [http://www.phytozome.net/search.php?method=Org_Gmax] using the program Primer3 Plus [27], available at [http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi]. The primer sequences were determined for the 3’ end of each gene whenever possible, and the amplicons spanned up to 150 base pairs (bp). The primer sequences were subjected to BLAST searches against the soybean genome [http://www.phytozome.net/search.php?method=Org_Gmax] to verify the specificity of each primer, as recommended by the Minimum Information for Publication of Quantitative Real-Time PCR Experiments guideline (MIQE) [14]. The primers for the commonly used reference genes ELF1-β and β-actin were selected from [28] and [29], respectively. The primers for the target gene (GmDREB5) were selected from [30]. Standard curves were produced from serial dilutions of a cDNA pool to estimate the efficiency of the PCR amplification with each pair of primers. Information on the primers may be visualized in Table 2.

Plant material and treatment application

Plant material was obtained from experiments performed as described by Marcolino-Gomes and colleagues [10]. Briefly, seeds from the soybean BR16 genotype, which exhibits drought-sensitive characteristics [31], were cultivated in peat pots (Jiffy) with Supersoil® (Scotts Miracle-Gro Company, Marysville, Ohio, USA) under optimal growth conditions in controlled growth chambers until reaching the V2 developmental stage [32], when water was withheld to induce a moderate water deficit. Control plants were maintained near field capacity for the unstressed treatment. The soil moisture was calculated by the gravimetric humidity (GH), which corresponds to the percentage of water in the soil in relation to the dry weight of the soil. The volume of irrigation was adjusted to 70% (GH) (near field capacity) for the unstressed treatment, 30% GH for the moderate stress treatment. The pots were weighed twice a day, and water was added to maintain the treatments at the desired GH values. The middle leaflet of the first trifoliate leaf was collected from plants in each treatment group at 4 h intervals from the time the lights came on (8:00 a.m. = Zeitgeiber Time (ZT) 0), during a 24 h timecourse (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20), and were immediately frozen in liquid nitrogen and stored at −80°C until further use. All of the experiments were conducted with three biological replicates, with each replicate consisting of two plants, whose tissues were collected together and pooled.

From RNA extraction to cDNA synthesis

Each replicate tissue set was ground to a fine powder in liquid nitrogen, and total RNA was isolated using the TRIzol reagent (Invitrogen) according to the manufacturer’s instructions. The obtained RNA concentration and purity were measured using a spectrophotometer (NanoDrop, ND-1000), and contaminating DNA in the total RNA was removed using the Turbo DNA-free kit, according to the manufacturer’s instructions (Life Technologies, Grand Island, NY, USA). After DNAse treatment, the integrity of the molecules was analyzed on 1% agarose gels stained with ethidium bromide, and high-quality total RNA was used to synthesize cDNA strands (Superscript III First Strand Synthesis, Invitrogen/Life Technologies, Grand Island, NY, USA). The quality of the cDNA and contamination with genomic DNA were examined using a standard PCR assay with primers that spanned an intronic region of the β-actin soybean gene. High-quality cDNA was used to analyze the transcripts in each treatment.

RT-qPCR analyses

Standard curves were produced from serial dilutions of a cDNA pool to estimate the efficiency of the PCR amplification with each pair of primers. RT-qPCR amplifications were performed in a 7300 RT-qPCR Thermocycler (Applied Biosystems/Life Technologies, Grand Island, NY, USA) with the following cycling parameters: 50°C for 2 min, 95°C for 10 min and 45 cycles at 95°C for 2 min, 60°C for 30 seconds and 72°C for 30 seconds. Each amplification reaction contained 2 μL of cDNA from serial dilutions, 60–200 nM each forward and reverse primer (Table 2), 500 nM ROX (passive reference), 6.5 μL of Platinum® SYBR® Green qPCR SuperMix (Invitrogen/ Grand Island, NY, USA), and ultrapure water to a final volume of 12.5 μL. Data were collected during the extension phase, and dissociation curves were generated by heating each reaction from 60 to 95°C and taking readings at one-degree intervals to verify the specificity of the primers. A control sample, obtained via performing RT-qPCR with no template, was also assayed to confirm that the samples were not contaminated. The primer concentrations were adjusted to achieve efficiency rates higher than 85%, as detailed in Table 2.

After carrying out the efficiency analysis, the expression levels of the candidate reference genes were analyzed separately at 6 time points (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20) in plants under control vs. drought stress condition in order to assess their expression stability along the day. The expression of the target gene (GmDREB5-like; Glyma12g33020) [30] was also measured under the experimental conditions described above. The reactions were performed in triplicate with cycling parameters similar to those described above for the amplification efficiency analysis.

Stability analysis

To validate and compare the suitability of the candidate reference genes for use in normalization, we evaluated their expression stability in response to drought along the day under the experimental conditions described above. For this purpose, Cycle threshold values (Ct) were transformed into non-normalized relative quantities (Q; linear scale). Here, Q = EΔCt, where E is the amplification efficiency, and ΔCt is the lowest Ct from the data set minus the sample Ct. The non-normalized relative quantities were analyzed using NormFinder [33] and geNorm [34] software to assess the expression stability of the reference genes.

Normalization of target gene expression

The relative expression level of the drought-responsive gene GmDREB5-like [30] was measured in leaf samples from BR16 plants subjected to moderate drought, sampled over a 24 h timecourse, under the experimental conditions described above. For each time point (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20), three biological replicates, with three technical replicates each, were analyzed. Target expression was normalized using a combination of 2 reference genes with high (FYVE and GOL84), intermediate (ELF1-β and β -actin) and low (DNAJ and NCL1) expression stabilities. Plants grown under normal water conditions (control plants) were used to calibrate relative expression.

The gene expression analysis was performed using the Rest2009 software package [35], which allows the input of different amplification efficiencies for the reference and target genes and provides the statistical significance of expression levels through randomization (Pair Wise Fixed Reallocation Randomisation Test©), with 10,000 interactions and bootstrapping of the data. At the randomization tests, the observed values were repeatedly and randomly reallocated to the two groups and the apparent effect (expression ratio in our case) was noted each time. The proportion of these effects which are as great as that actually observed in the experiment gave us the P-value of the test. In the applied Pair Wise Fixed Reallocation Randomisation Test©, for each sample, the CP values for reference and target genes were jointly reallocated to control and sample groups (= pairwise fixed reallocation), and the expression ratios was calculated on the basis of the mean values. Hypothesis testing was conducted to determine whether the differences between the control and treatment conditions were significant [35].

Results and Discussion

Screening of candidate reference genes

To date, most of the studies on reference genes have focused on validating a subset of commonly used reference genes for specific contexts [15,16,18]. Although these studies have their merits, they attempt to identify the best candidates from a small set of genes. A recent analysis demonstrated that reference genes are preferably selected by adopting a complete genome strategy, rather than from a handful of commonly used reference genes [26]. In this context, we searched reference genes showing high expression stability in 59 microarray libraries from soybean subjected to drought stress, heat and different photoperiods [26].

This tool allowed us to perform in silico identification of genes showing high expression stability in 59 microarray libraries from soybean subjected to drought, heat and distinct light periods. The candidate reference genes obtained in this analysis were pre-validated by checking their expression across all microarrays available on the Genevestigator platform (3458 arrays) (data not shown). The expression profiles of the five most stable new reference genes in response to drought and to diurnal oscillations was compared with the expression of the commonly used soybean reference genes (Table 1), as shown in Fig 1.

thumbnail
Fig 1. Expression of commonly used reference genes and new candidates from microarray databases.

The in silico analyses were performed with Genevestigator software, using data from soybean subjected to drought stress, heat and different photoperiods[27]. The box plots represent the variation in the signal intensity (log2 scale). Interquartile range (IQR) values are shown for each gene across the dataset, represented by the middle boxes, which encompass the middle 50% of scores for the group. The upper and lower whiskers represent scores outside the middle 50%. Outliers are plotted separately as asterisks on the chart.

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

Then, in a second approach, we selected genes that exhibited minimal expression variance across 36 cDNA libraries synthesized from drought-stressed soybean plants sampled over a 24 h timecourse [10,36]. The expression profiles of the two most stable new reference genes and the commonly used references are shown in Fig 2.

thumbnail
Fig 2. Variation in the expression of commonly used reference genes and new candidates from the RNA-Seq database.

The in silico gene expression analyses were performed using data from soybean grown under drought stress across a 24 h timecourse. (A) The box plots represent the variation in the signal intensity (log2 scale); the middle boxes represent the middle 50% of scores for the group; the upper and lower whiskers represent scores outside the middle 50%. (B) Coefficient of variation (CV%).

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

In summary, the data mining approaches using microarray and RNA-seq databases allowed us to select a new set of candidate reference genes for validation via RT-qPCR, composed of seven soybean genes: Glyma13g24060, Glyma01g40510, Glyma13g17500, Glyma08g05790, Glyma11g38000, Glyma08g41240 and Glyma10g44020.

The majority of the selected candidate genes are related to the plant’s primary metabolism. For example, Glyma01g40510 encodes a cysteine desulfurase (CYST) similar to nitrogen fixation S (NIFS)-like 1 from Arabidopsis; Glyma08g05790 encodes a protein that participates in Golgi vesicles transport (Golgin-84) [37,38]; Glyma11g38000 produces an RNA (cytosine-5)-methyltransferase (NCL1) involved in epigenetic modifications of tRNA [39,40]; and Glyma13g17500 produces an FYVE domain protein, present in kinases and lipases in Arabidopsis, that recognizes phosphoinositide signals [41].

The differential expression of some of the selected candidate genes has been reported during biotic and abiotic stress responses. The gene Glyma13g24060, for example, encodes a protein similar to a NUDIX hydrolase protein from Arabidopsis. The NUDIX hydrolase family is widespread, from eukaryotes to Archaea, and consists of pyrophosphohydrolases that act upon substrates with a general nucleoside diphosphate structure, including (deoxy)ribonucleoside diphosphates and triphosphates, nucleotide sugars, coenzymes and RNA caps [42,43]. Members of the NUDIX family have been reported to be induced by salt, drought, heat, and cold in Chrysanthemum lavandulifolium [44]. Similarly, Glyma10g44020 encodes a protein from the DnaJ/Hsp40 cysteine-rich domain superfamily, which is described as being involved in diverse cellular processes (protein folding, translocation, and degradation) [45], including biotic and abiotic stress responses [4649]. However, analysis of the soybean NUDIX (Fig 1) and DNAJ (Fig 2) genes in microarray and RNA-Seq databases showed high expression stability in response to drought and diurnal oscillations, suggesting that these soybean genes could be reliable candidate reference genes for drought studies.

Furthermore, we identified a gene (Glyma08g41240) that encodes an RNA-dependent RNA polymerase from a mitovirus (Fig 2). A recent study on the soybean mitochondrial genome revealed the presence of a 0.5 kb insertion (at rps10 intron) that is 57.4% identical to a mitovirus RNA polymerase gene, which might have been horizontally transferred during recent evolution. Although the effect of this insertion remains unknown, analysis of the insert’s position suggests that it might affect the function of the mitochondrial rps10 gene, which encodes the ribosomal protein S10 [50].

The preliminary analysis of this set of genes (Figs 1 and 2) revealed that in the majority of cases, the candidate reference genes presented less variation than the commonly used reference genes selected from the literature [15,19]. Additionally, in silico pre-validation of the NUDIX, CYST, FYVE, Golgin-84, and NCL1 genes across 3,458 microarrays using GeneVestigator platform revealed that these genes are unresponsive to a wide variety of conditions, including abiotic stresses, such as heat, salinity, and cold, and show little variation between developmental stage and genotypes, being responsive only to infections by Phytophthora sojae, Phakopsora pachyrhizi and Bradyrhizobium japonicum (data not shown).

Stability analysis

To experimentally validate this new set of candidate reference genes, we examined their individual properties and compared the stability of their expression with the commonly used reference genes ELF1-β and β-actin (Table 1) using RT-qPCR assays. In previous studies, ELF1-β and β-actin showed high expression stability across different levels of drought [15,28,30]. However, no study conducted to date has investigated the expression stability of these genes during the diurnal cycle. Validation was carried out on soybean leaves from plants subjected to a moderate water deficit (30% of gravimetric humidity (GH)), sampled across a 24 h timecourse, with 4 h intervals, from 8:00 a.m. to 4:00 a.m., corresponding to Zeitgeber Time (ZT) 0 to ZT20. The expression level of each reference gene was evaluated separately in drought-treated and control plants at all sampling times (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20).

To assess the stability values for the reference genes, we performed analysis using the NormFinder [33] and geNorm softwares [34]. The NormFinder software employs a variance estimation approach to rank genes according to combined inter- and intragroup expression variation across a given set of experimental conditions, which were drought stress and diurnal oscillations in this case. The NormFinder [33] is known to perform in a more robust manner and with less sensitivity toward the co-regulation of candidate genes compared with other software. On the other hand, the NormFinder software may be less robust with small sample sizes compared to the geNorm algorithm [33]. Thus, we performed additional stability analysis using geNorm software [34] to compare the performance of the candidate genes using stability analysis tools with different algorithms. In general, the results of geNorm analysis were similar to those from NormFinder, with slight differences regarding ranking positions (Fig 3A and 3B). This good consistence between both outcomes strengthens the robustness of the results obtained. Small differences in rank position among the two software are expected because the statistical algorithms used are distinct: the geNorm detects the two reference genes whose expression ratios show the least variation from those of the other tested genes [34], whereas the NormFinder takes intra- and intergroup variation into account for calculations [33].

thumbnail
Fig 3. Stability analysis of candidate reference genes.

Gene expression stability was measured in leaf tissues of soybean subjected to drought stress and control conditions at different times of day (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20). The analysis was performed with NormFinder and geNorm softwares. Genes were ranked according to their stability values and M values from (A) NormFinder and (B) geNorm, respectively. The genes were plotted on the x-axis from the most to the least stable.

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

According these stability analysis, many of the newer reference genes indeed exhibited greater expression stability than the conventionally used reference genes (ELF1-β and β -actin) (Fig 3A and 3B). The FYVE, NUDIX and Golgin-84 genes were the most stable, suggesting that they are the most suitable for normalizing expression data from combined studies addressing drought treatment and diurnal oscillation (Fig 3A and 3B).

Although most of the candidate reference genes performed well in response to the applied experimental conditions, the NCL1 and DNAJ genes showed lower stability (Fig 3A and 3B). The NCL1 gene produces an RNA (cytosine-5)-methyltransferase involved in epigenetic modifications of tRNA [39,40]. The analysis of intragroup variation performed on NormFinder for NCL1 revealed that under control conditions this was one of the least stable genes (Fig 4A), whereas in stressed plants it was among the most stable (Fig 4B), indicating that NCL1 is not a suitable reference gene for studies on gene expression oscillation during the day in plants under control conditions. Previous studies show that the oscillation of genes that cycle daily in normal conditions may be altered due the imposition of abiotic stresses, like cold and drought stress. The mRNA levels of some chestnut genes, like TOC1 and LHY cycle daily in seedlings and adult plants, as expected. However, during chilling stress (4°C), mRNA levels of these genes were higher and did not oscillate [51]. Similar events were observed in soybean plants subjected to severe drought stress, where a general reduction in the amplitude of the daily oscillation was observed for most clock genes, including the GmPRR3-like, GmPRR7-like, GmPRR9- like, GmGI-like, GmZTL-like, and GmCHE-like genes [10].

thumbnail
Fig 4. Intra- and intergroup variation of gene expression.

The intragroup variation within the (A) control (non-stressed plants) and (B) stressed (plants under drought stress) and (C) the intergroup variation are presented.

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

Furthermore, these results illustrate that gene expression stability during the day may vary in response to stressful conditions, like changes in the plant’s water status (e.g. normal hydration versus water deficit conditions).

In contrast, the expression of the DNAJ gene was unstable in both control and stressed conditions (Fig 4). Genes from the DNAJ family have been reported to be drought responsive in many species, and its lack of stability may therefore be explained in part by its possible involvement in drought responses in soybean, as studies have reported that drought-responsive genes oscillate in response to diurnal oscillations [10,12,13,30,52]. Additionally, the intragroup analysis shows that Golgin-84 and FYVE (Fig 4A) are the most reliable reference genes for gene expression normalization when studying diurnal oscillations in plants under normal water conditions, whereas NUDIX and NCL1 (Fig 4B) are the most strongly indicated for use in studies under drought conditions. Furthermore, the intergroup variation analysis allowed us to identify the genes NUDIX and DNAJ as the most and least stable genes, respectively, for data normalization in studies comparing gene expression under control and drought stress, without considering the time of day (Fig 4C).

Normalizing the expression of a target gene

A previous study showed that the normalization of soybean genes under circadian regulation using unstable reference genes may lead to erroneous data interpretation [19]. To demonstrate the effect of data normalization using reference genes with different stability values in response to drought stress during the day, we evaluated the relative expression of the drought-responsive gene GmDREB5-like [30] using a combination of 2 reference genes with high (FYVE and GOL84), intermediate (ELF1-β and β -actin) and low (DNAJ and NCL1) expression stabilities (Fig 5). The GmDREB5 gene expression was previously investigated in response to drought in short-term stress conditions [30], and its expression was normalized by ELF1-β and β –actin. As expected, our results show that the expression of GmDREB5-like in response to drought oscillates during the day, as previously described for this gene [30] and for other genes of the DREB subfamily [10,5358] (Fig 5).

thumbnail
Fig 5. Normalization of the target gene GmDREB5-like.

Gene expression was measured in leaf tissues of soybean subjected to drought stress at different times of day (ZT0, ZT4, ZT8, ZT12, ZT16 and ZT20). The raw data were normalized using a combination of 2 candidate reference genes with high (FYVE and GOL84), intermediate (ELF1-β and β -actin) and low (DNAJ and NCL1) expression stabilities. The relative expression of GmDREB5-like under drought stress was determined after comparison with the control sample (non-stressed plants). Asterisks indicate statistically significant changes in gene expression between drought stressed and non-stressed plants. The analyses were performed using REST 2009 software.

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

However, our results demonstrated that choosing reference genes with diverse stability of expression can lead to differences on gene expression data interpretation when evaluating combined studies on water deficit and diurnal oscillations. As shown in Fig 5, at ZT16 the normalization of GmDREB5-like gene using the least stable genes (DNAJ and NCL1) resulted in higher expression levels than observed for normalization using genes with high (FYVE and GOL84) and intermediate expression stability (ELF1-β and β –actin). Additionally, slight changes in gene expression, such as the down-regulation of the target gene at ZT0 and ZT20, were detected only by normalization using genes with high (FYVE and GOL84) and intermediate expression stability (ELF1-β and β –actin) (Fig 5A). These results emphasize the importance of selecting reference genes with stable expression for accurate gene expression analysis on drought responses and diurnal oscillations.

Conclusions

Here, by analyzing experiments involving both drought and diurnal oscillations, we demonstrated the importance of selecting reference genes under the specific studied conditions. From a transcriptome-wide background, we selected a new set of candidate reference genes for the normalization of data obtained in studies on drought and diurnal oscillations.

The experimental validation of this new set of candidate reference genes revealed that FYVE, NUDIX and Golgin-84 were the most stably expressed genes in soybean plants under control and drought conditions along the day, and are therefore considered the best reference genes for the studied conditions. Our results highlight that the selection of reference genes is crucial for the proper quantification of relative expression data obtained under these specific experimental conditions.

Acknowledgments

We thank the National Center for Soybean Research—Embrapa Soybean—for the use of the greenhouse and laboratory facilities.

Author Contributions

Conceived and designed the experiments: JMG FGH ALN. Performed the experiments: JMG TJN RRR. Analyzed the data: JMG TJN FGH ALN. Contributed reagents/materials/analysis tools: JRBF. Wrote the paper: JMG FGH RFP FAR ALN HBCM MDCM.

References

  1. 1. Yerushalmi S, Green RM. Evidence for the adaptive significance of circadian rhythms. Ecol Lett. 2009;12: 970–81. pmid:19566794
  2. 2. Haydon MJ, Mielczarek O, Robertson FC, Hubbard KE, Webb A a R. Photosynthetic entrainment of the Arabidopsis thaliana circadian clock. Nature. Nature Publishing Group; 2013;502: 689–92. pmid:24153186
  3. 3. Graf A, Smith AM. Starch and the clock: the dark side of plant productivity. Trends Plant Sci. Elsevier Ltd; 2011;16: 169–75. pmid:21216654
  4. 4. Gutiérrez R, Stokes T. Systems approach identifies an organic nitrogen-responsive gene network that is regulated by the master clock control gene CCA1. Proc Natl Acad Sci. 2008;105: 4939–4944. pmid:18344319
  5. 5. Mahdieh M, Mostajeran A, Horie T, Katsuhara M. Drought stress alters water relations and expression of PIP-type aquaporin genes in Nicotiana tabacum plants. Plant Cell Physiol. 2008;49: 801–13. pmid:18385163
  6. 6. Babita M, Maheswari M, Rao LM, Shanker AK, Rao DG. Osmotic adjustment, drought tolerance and yield in castor (Ricinus communis L.) hybrids. Environ Exp Bot. Elsevier B.V.; 2010;69: 243–249.
  7. 7. Molinari HBC, Marur CJ, Daros E, de Campos MKF, de Carvalho JFRP, Filho JCB, et al. Evaluation of the stress-inducible production of proline in transgenic sugarcane (Saccharum spp.): osmotic adjustment, chlorophyll fluorescence and oxidative stress. Physiol Plant. 2007;130: 218–229.
  8. 8. Carvalho MHC. Drought stress and reactive oxygen species. Plant Signal Behav. 2008;3: 156–165. pmid:19513210
  9. 9. Legnaioli T, Cuevas J, Mas P. TOC1 functions as a molecular switch connecting the circadian clock with plant responses to drought. EMBO J. Nature Publishing Group; 2009;28: 3745–57.
  10. 10. Marcolino-Gomes J, Rodrigues FA, Fuganti-Pagliarini R, Bendix C, Nakayama TJ, Celaya B, et al. Diurnal oscillations of soybean circadian clock and drought responsive genes. PLoS One. 2014;9: e86402. pmid:24475115
  11. 11. Sanchez A, Shin J, Davis SJ. Abiotic stress and the plant circadian clock. Plant Signal Behav. 2011;6: 223–231. pmid:21325898
  12. 12. Wilkins O, Bräutigam K, Campbell MM. Time of day shapes Arabidopsis drought transcriptomes. Plant J. 2010;63: 715–727. pmid:20553421
  13. 13. Wilkins O, Waldron L, Nahal H, Provart NJ, Campbell MM. Genotype and time of day shape the Populus drought response. Plant J. 2009;60: 703–15. pmid:19682285
  14. 14. Bustin S a, Benes V, Garson J a, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55: 611–22. pmid:19246619
  15. 15. Ma S, Niu H, Liu C, Zhang J, Hou C, Wang D. Expression stabilities of candidate reference genes for RT-qPCR under different stress conditions in soybean. PLoS One. 2013;8: e75271. pmid:24124481
  16. 16. Miranda VDJ, Coelho RR, Viana AAB, de Oliveira Neto OB, Carneiro RMDG, Rocha TL, et al. Validation of reference genes aiming accurate normalization of qPCR data in soybean upon nematode parasitism and insect attack. BMC Res Notes. 2013;6: 196. pmid:23668315
  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–71. pmid:24615050
  18. 18. Le DT, Aldrich DL, Valliyodan B, Watanabe Y, Ha C Van, Nishiyama R, et al. Evaluation of candidate reference genes for normalization of quantitative RT-PCR in soybean tissues under various abiotic stress conditions. PLoS One. 2012;7: e46487. pmid:23029532
  19. 19. Hu R, Fan C, Li H, Zhang Q, Fu Y-F. Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR. BMC Mol Biol. 2009;10: 93. pmid:19785741
  20. 20. Song YH, Lee I, Lee SY, Imaizumi T, Hong JC. CONSTANS and ASYMMETRIC LEAVES 1 complex is involved in the induction of FLOWERING LOCUS T in photoperiodic flowering in Arabidopsis. Plant J. 2012;69: 332–42. pmid:21950734
  21. 21. Huang W, Pérez-García P, Pokhilko A, Millar AJ, Antoshechkin I, Riechmann JL, et al. Mapping the core of the Arabidopsis circadian clock defines the network structure of the oscillator. Science (80-). 2012;336: 75–9. pmid:22403178
  22. 22. Jiang Y, Han YZ, Zhang XM. Expression profiles of a CONSTANS homologue GmCOL11 in Glycine max. Russ J Plant Physiol. 2011;58: 928–935.
  23. 23. Galeou A, Prombona A. Light at night resynchronizes the evening-phased rhythms of TOC1 and ELF4 in Phaseolus vulgaris. Plant Sci. Elsevier Ireland Ltd; 2012;184: 141–147. pmid:22284718
  24. 24. Hudson KA. The Circadian Clock-controlled Transcriptome of Developing Soybean Seeds. Plant Genome J. 2010;3: 3.
  25. 25. Li G, Siddiqui H, Teng Y, Lin R, Wan X, Li J, et al. Coordinated transcriptional regulation underlying the circadian clock in Arabidopsis. Nat Cell Biol. Nature Publishing Group; 2011;13: 616–22. pmid:21499259
  26. 26. Hruz T, Wyss M, Docquier M, Pfaffl MW, Masanetz S, Borghi L, et al. RefGenes: identification of reliable and condition specific reference genes for RT-qPCR data normalization. BMC Genomics. BioMed Central Ltd; 2011;12: 156. pmid:21418615
  27. 27. Rozen S, Skaletsky H. Primer3 on the WWW for General Users and for Biologist Programmers. In: Misener S, Krawetz SA, editors. Bioinformatics Methods and Protocols Methods in Molecular Biology. Humana Press; 1999. pp. 365–386.
  28. 28. Jian B, Liu B, Bi Y, Hou W, Wu C, Han T. Validation of internal control for gene expression study in soybean by quantitative real-time PCR. BMC Mol Biol. 2008;9: 59. pmid:18573215
  29. 29. Stolf-Moreira R, Lemos EG de M, Abdelnoor RV, Beneventi MA, Rolla AAP, Pereira S dos S, et al. Identification of reference genes for expression analysis by real-time quantitative PCR in drought-stressed soybean. Pesqui Agropecu Bras. Pesquisa agorpecuaria brasileira; 2011;46: 58–65.
  30. 30. Marcolino-Gomes J, Rodrigues FA, Oliveira MCN, Farias JRB, Neumaier N, Abdelnoor RV, et al. Expression Patterns of GmAP2/EREB-Like Transcription Factors Involved in Soybean Responses to Water Deficit. PLoS One. 2013;8: e62294. pmid:23667465
  31. 31. Oya T, Nepomuceno AL, Neumaier N, Renato J, Farias B, Tobita S, et al. Drought Tolerance Characteristics of Brazilian Soybean Cultivars: Evaluation and characterization of drought tolerance of various Brazilian soybean cultivars in the field. Plant Prod Sci. 2004;7: 129–137.
  32. 32. Fehr WR, Caviness CE, Burmood DT, Pennington JS. Stage of Development Descriptions for Soybeans, Glycine Max (L.) Merrill1. Crop Sci. 1971;11: 929.
  33. 33. Andersen CL, Jensen JL, Ørntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64: 5245–50. pmid:15289330
  34. 34. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3: RESEARCH0034. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=126239&tool=pmcentrez&rendertype=abstract
  35. 35. Pfaffl MW, Horgan GW, Dempfle L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30: e36. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=113859&tool=pmcentrez&rendertype=abstract pmid:11972351
  36. 36. 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
  37. 37. Munro S. The golgin coiled-coil proteins of the Golgi apparatus. Cold Spring Harb Perspect Biol. 2011;3. pmid:21436057
  38. 38. Osterrieder a. Tales of tethers and tentacles: golgins in plants. J Microsc. 2012;247: 68–77. pmid:22591132
  39. 39. Pavlopoulou A, Kossida S. Phylogenetic analysis of the eukaryotic RNA (cytosine-5)-methyltransferases. Genomics. Elsevier Inc.; 2009;93: 350–7. pmid:19135144
  40. 40. Squires JE, Preiss T. 5-Methylcytosine as a modification in RNA. In: Craig J, Wong NC, editors. Epigenetics: A Reference Manual. 1st ed. Norfolk: Horizon Scientific Press; 2011. p. 449.
  41. 41. Jensen RB, Cour TLA, Albrethsen J, Nielsen M, Skriver K. identification of PtdIns3 P-binding residues by comparison of classic and variant FYVE domains. 2001;173: 165–173.
  42. 42. Ogawa T, Yoshimura K, Miyake H, Ishikawa K, Ito D, Tanabe N, et al. Molecular characterization of organelle-type Nudix hydrolases in Arabidopsis. Plant Physiol. 2008;148: 1412–24. pmid:18815383
  43. 43. Goyer A, Hasnain G, Frelin O, Ralat M a, Gregory JF, Hanson AD. A cross-kingdom Nudix enzyme that pre-empts damage in thiamin metabolism. Biochem J. 2013;454: 533–42. pmid:23834287
  44. 44. Huang H, Cao H, Niu Y, Dai S. Expression Analysis of Nudix Hydrolase Genes in Chrysanthemum lavandulifolium. Plant Mol Biol Report. 2012;30: 973–982.
  45. 45. Rajan VB V, D’Silva P. Arabidopsis thaliana J-class heat shock proteins: cellular stress sensors. Funct Integr Genomics. 2009;9: 433–46. pmid:19633874
  46. 46. Wang G, Cai G, Kong F, Deng Y, Ma N, Meng Q. Overexpression of tomato chloroplast-targeted DnaJ protein enhances tolerance to drought stress and resistance to Pseudomonas solanacearum in transgenic tobacco. Plant Physiol Biochem. 2014;82: 95–104. pmid:24929777
  47. 47. So H-A, Chung E, Lee J-H. Molecular characterization of soybean GmDjp1 encoding a type III J-protein induced by abiotic stress. Genes Genomics. 2013;35: 247–256.
  48. 48. Liu JZ, Whitham S a. Overexpression of a soybean nuclear localized type-III DnaJ domain-containing HSP40 reveals its roles in cell death and disease resistance. Plant J. 2013;74: 110–121. pmid:23289813
  49. 49. Bekh-Ochir D, Shimada S, Yamagami A, Kanda S, Ogawa K, Nakazawa M, et al. A novel mitochondrial DnaJ/Hsp40 family protein BIL2 promotes plant growth and resistance against environmental stress in brassinosteroid signaling. Planta. 2013;237: 1509–1525. pmid:23494613
  50. 50. Chang S, Wang Y, Lu J, Gai J, Li J, Chu P, et al. The mitochondrial genome of soybean reveals complex genome structures and gene evolution at intercellular and phylogenetic levels. PLoS One. 2013;8: e56502. pmid:23431381
  51. 51. Ramos A, Pérez-Solís E, Ibáñez C, Casado R, Collada C, Gómez L, et al. Winter disruption of the circadian clock in chestnut. Proc Natl Acad Sci U S A. 2005;102: 7037–42. pmid:15860586
  52. 52. Mizuno T, Yamashino T. Comparative transcriptome of diurnally oscillating genes and hormone-responsive genes in Arabidopsis thaliana: insight into circadian clock-controlled daily responses to common ambient stresses in plants. Plant Cell Physiol. 2008;49: 481–7. pmid:18202002
  53. 53. Kidokoro S, Maruyama K, Nakashima K, Imura Y, Narusaka Y, Shinwari ZK, et al. The phytochrome-interacting factor PIF7 negatively regulates DREB1 expression under circadian control in Arabidopsis. Plant Physiol. 2009;151: 2046–57. pmid:19837816
  54. 54. Nakamichi N, Kusano M, Fukushima A, Kita M, Ito S, Yamashino T, et al. Transcript profiling of an Arabidopsis PSEUDO RESPONSE REGULATOR arrhythmic triple mutant reveals a role for the circadian clock in cold stress response. Plant Cell Physiol. 2009;50: 447–62. pmid:19131357
  55. 55. Fowler SG, Cook D, Thomashow MF. Low temperature induction of Arabidopsis CBF1, 2, and 3 is gated by the circadian clock. Plant Physiol. 2005;137: 961–8. pmid:15728337
  56. 56. Dong MA, Farré EM, Thomashow MF. Circadian clock-associated 1 and late elongated hypocotyl regulate expression of the C-repeat binding factor (CBF) pathway in Arabidopsis. Proc Natl Acad Sci U S A. 2011;108: 7241–6. pmid:21471455
  57. 57. Mikkelsen MD, Thomashow MF. A role for circadian evening elements in cold-regulated gene expression in Arabidopsis. Plant J. 2009;60: 328–339. pmid:19566593
  58. 58. Bieniawska Z, Espinoza C, Schlereth A, Sulpice R, Hincha DK, Hannah M a. Disruption of the Arabidopsis circadian clock is responsible for extensive variation in the cold-responsive transcriptome. Plant Physiol. 2008;147: 263–279. pmid:18375597