Figures
Abstract
Switchgrass (Panicum virgatum) has received a lot of attention as a forage and bioenergy crop during the past few years. Gene expression studies are in progress to improve new traits and develop new cultivars. Quantitative real time PCR (qRT-PCR) has emerged as an important technique to study gene expression analysis. For accurate and reliable results, normalization of data with reference genes is essential. In this work, we evaluate the stability of expression of genes to use as reference for qRT-PCR in the grass P. virgatum. Eleven candidate reference genes, including eEF-1α, UBQ6, ACT12, TUB6, eIF-4a, GAPDH, SAMDC, TUA6, CYP5, U2AF, and FTSH4, were validated for qRT-PCR normalization in different plant tissues and under different stress conditions. The expression stability of these genes was verified by the use of two distinct algorithms, geNorm and NormFinder. Differences were observed after comparison of the ranking of the candidate reference genes identified by both programs but eEF-1α, eIF-4a, CYP5 and U2AF are ranked as the most stable genes in the samples sets under study. Both programs discard the use of SAMDC and TUA6 for normalization. Validation of the reference genes proposed by geNorm and NormFinder were performed by normalization of transcript abundance of a group of target genes in different samples. Results show similar expression patterns when the best reference genes selected by both programs were used but differences were detected in the transcript abundance of the target genes. Based on the above research, we recommend the use of different statistical algorithms to identify the best reference genes for expression data normalization. The best genes selected in this study will help to improve the quality of gene expression data in a wide variety of samples in switchgrass.
Citation: Gimeno J, Eattock N, Van Deynze A, Blumwald E (2014) Selection and Validation of Reference Genes for Gene Expression Analysis in Switchgrass (Panicum virgatum) Using Quantitative Real-Time RT-PCR. PLoS ONE 9(3): e91474. https://doi.org/10.1371/journal.pone.0091474
Editor: Christian Schönbach, Nazarbayev University, Kazakhstan
Received: November 20, 2013; Accepted: February 12, 2014; Published: March 12, 2014
Copyright: © 2014 Gimeno et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: Funding came from the University of California Discovery Grants. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The instability of petroleum sources and threat of global climate change have incited a push towards locally produced renewable energy sources. Recent scientific developments have allowed the conversion of renewable biologic materials to liquid or gaseous fuels. Of particular interest is the saccharification of cellulose cell wall, the principal component of plants, to its sugar components that can be fermented into high value liquid ethanol with the potential to displace a substantial portion of petroleum based fuels.
Switchgrass (Panicum virgatum L.) is a perennial warm-season grass native to North America and a productive C4 species. Used initially as a forage crop, switchgrass is currently being pursued as a dedicated cellulosic biofuel feedstock in the Midwestern and Southeastern United Sates. Switchgrass has been shown to produce 540% more renewable energy than energy consumed in its production and has significant environmental benefits related to soil conservation and low greenhouse gas emission [1]. Its perennial nature allows efficient nutrient cycling to below ground organs for the subsequent seasons' growth, a characteristic that allows substantial annual harvests with low nutrient and water requirements [2].
As a result of the increasing interest in this species, numerous studies are under way in switchgrass to understand the basis of water and nutrient- use efficiency [3]–[6] and the mechanisms of cell wall synthesis [7]–[8]. Moreover, an increasing number of genomic resources have been developed including SSR markers [9]–[11], miRNAs [12], ESTs, microarrays [13]–[16] and a transformation system [17].
To verify the function of a candidate gene, the analysis of mRNA expression patterns is one of the most used molecular techniques. Real-time RT-PCR (also referred to as quantitative RT-PCR or qRT-PCR) has become a popular approach for the analysis of gene expression, its main advantages being a high sensitivity, specificity and broad quantification range [18]. The accuracy of qRT-PCR is strongly influenced by the stability of the internal reference genes used for data normalization [19], [20]. It requires the use of constitutively expressed genes as internal controls to correct for variability associated with the different steps of the experimental procedure. The expression levels of these genes should remain constant between different tissues and/or in response to different environmental stimuli. Though there are several genes frequently used as HKGs (Housekeeping genes) (18S sRNA, actin, tubulin), many studies demonstrate that the expression of these genes is not always stable and the selection of alternative genes with stable expression profiles must be identified in each organism across experimental conditions prior to the use of real-time PCR [21]. Several statistical approaches have been proposed for evaluating reference genes such geNorm [22], NormFinder [23] and Best-Keeper [24]. These approaches are based on different statistical algorithms and employ multiple reference genes as the best strategy for normalization of qPCR results.
Until now, the number of studies on gene expression in switchgrass is scant and few genes are used as reference genes for normalization [7], [25]–[28]. A recent switchgrass transcriptome analysis [16] identified a set of stably expressed genes by measuring the transcript levels in all major organ systems at different stages of development. This group of genes could be a source of candidate reference genes but will need to be evaluated as reference genes for each experimental condition.
Given the lack of information on appropriate reference genes in switchgrass, eleven HKGs, commonly used in plant gene expression studies, were selected in order to assess them as internal controls. These genes are elongation factor 1α (eEF-1α), ubiquitin (UBQ6), actin (ACT12), β-tubulin (TUB6), eukaryotic initiation factor (eIF-4a), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), s-adenosyl methionine decarboxylase (SAMDC), α-tubulin (TUA6), cyclophilin (CYP5), splicing factor (U2AF) and FtsH protease 4 (FTSH4). The stability of expression of these candidate reference genes was tested in different switchgrass tissues and under different stress conditions using both geNorm [22] and NormFinder [23] algorithms.
In addition, we provide a detailed expression analysis of seven putative homologues of the cellulose synthase gene family in switchgrass (SGCesAs), as an example for the genes involved in cellulose biosynthesis. Defining the pattern of their expression in different tissues and under different stress conditions may lead to their utilization towards improving cellulose production in energy crops, main sources of cellulosic ethanol. To illustrate the usefulness of the most stable reference genes identified in this work, the best combinations of reference genes were used to perform the expression normalization. Therefore, this study provides useful guidelines and starting point for reference gene selection for expression analysis using qRT-PCR techniques in P.virgatum.
Materials and Methods
Plant material
Switchgrass (Panicum virgatum L. cv. Alamo) seeds on water were cold-treated at 4°C overnight to synchronize germination and planted in the greenhouse in flats filled with UC mix (50% washed sand, 50% sphagnum peat moss). After germination, seedlings of 35–40 cm where transplanted to 7.6-liter pots filled with the UC mix.
Plants were grown in the greenhouse under 14 h photoperiod at 30°C day and 18°C night and the relative humidity was maintained between 40 and 70%. Irrigation was applied in early-morning under a regime of two days fertilization (N:P:K = 261∶30∶326 ppm) followed by one day of deionized water. Plants received 500 ml per day, the amount of water necessary to maintain the field capacity.
Stress treatments and tissue collection
Samples were collected from different tissues and under different stress conditions. Tissue from leaves, stems, rachis and roots was collected from three-month old plants growing under well-watered conditions. Plants to be treated with different abiotic stresses were grown for four weeks under field capacity after transplanting. Plants were subjected to drought treatments by reducing the irrigation to 250 ml per day (50% water) for one week. Salt stress treatment was applied by a progressive increase of the salt concentration in the solution (40, 120, 200 and 300 mM NaCl). Four days of salt treatment were followed by one day of fertilization. Last salt solution concentration was applied for seven days. For cold and heat shock treatments plants were transferred to 4°C and 42°C respectively for 48 h. To apply wounding stress, the upper leaf of a tiller was mechanically wounded 3–5 times by a hemostat perpendicular to the mid-vein located two-thirds from the tip of the leaf. The leaf area wounded was approximate the 10% of the leaf surface. Tissue was collected after 24 hours. Flooding treatment was applied for 30 days by maintaining a water level above the soil surface. Tissue samples were collected from the last fully expanded leaf after each stress treatment. All 10 samples were collected from three replicate plants, giving a total of 30 samples comprised of 12 tissue-specific samples and 18 abiotic stress treatment samples. Samples were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction.
Total RNA isolation and cDNA synthesis
A total of 30 mg of frozen plant material was ground in eppendorf tubes using 3.2 mm stainless steel beads and a shaker SO-10M (Fluid Management, Sassenheim, Netherlands). Total RNA was extracted using the RNeasy Mini kit (Qiagen, Valencia, USA) with the addition of an on-column DNase I digestion according to the manufacturer's instructions. RNA sample concentration and quality was determined using the NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The quality of RNA samples was also assessed by 1% agarose gel electrophoresis.
First-strand cDNA was synthesized from 1 μg of total RNA in a total volume of 20 μl per reaction using the QuantiTect Reverse Transcription kit (Qiagen, Valencia, USA) following the manufacturer's protocol. RNA quantity and quality were verified using 1% agarose gel electrophoresis after the genomic DNA elimination step of this protocol. cDNAs were diluted to a final volume of 400 μl to be used on the quantitative real-time RT-PCR reactions (qRT-PCR).
Selection of candidate reference genes, primer design and validation
A group of genes commonly used as reference genes in different plant species were selected to be evaluated as control genes in P. virgatum. Rice protein sequences from the potential reference genes were used as a query sequence for a TBLASTN search of the switchgrass GenBank EST database. Selected switchgrass ESTs were then used to query the non-redundant protein database using BLASTX for identity verification. The gene name, EST accession number and gene ontology are shown in Table 1. Specific primers for qRT-PCR were designed using the Primer Express 3.0 software (PE Applied Biosystems, Foster City, USA) using the default parameters for a real-time PCR assay. Amplicon lengths varied from 90 to 100 bp, melting temperatures were between 58–60°C and primer lengths between 18–24 bp (Table 2). The primers were screened for hairpins, dimer formation, and target specificity by BLASTN against the nr databank. A pool of cDNA from 30 samples was used to perform the qRT-PCR reactions to determine the primer efficiency. The entire raw fluorescence data was used to test each primer pair amplification with the Excel workbook entitled Data Analysis for Real-Time PCR (DART-PCRv.1.0) [29]. All primer amplification efficiencies were between 97 and 110% (Table 2). Agarose gel electrophoresis and melting curve of the amplification products of each candidate gene were analyzed in order to verify that a single PCR product was amplified for each set of primers.
Development of specific primers for CesA genes in switchgrass
The amino acid sequences of Oryza sativa cellulose synthase gene family were used as query sequences for TBLASTN search of the switchgrass GenBank EST database. TIGR identifiers of the rice genes used are: Os01g54620, Os03g59340, Os07g10770, Os05g08370, Os07g24190, Os03g62090, Os07g14850, Os09g25490, Os10g32980 and Os06g39970. EST sequences obtained (>80% identity) were assembled by CAP3 assembly tool [30] with default parameters. The resulting partial DNA consensus sequences after assembly (File S1) were used to design specific primers on the HVR-II region to discriminate the different members of the family and analyze the expression patterns. Primer design and evaluation (Table S1) and qRT-PCR were performed as described previously.
Real-time RT-PCR analysis
The real-time RT-PCRs were performed using the StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, USA). Each PCR reaction contained 7.5 μl of 2× Fast SYBR Green master mix (Applied Biosystems, Foster City, USA), 2 μl of the diluted cDNA reaction mixture (corresponding to 5 ng of starting amount of RNA) and 200 nM of each primer in a total reaction volume of 15 μl. PCR reactions were performed under the following conditions: 10 min at 95°C and 40 cycles of the one step thermal cycling of 3 s at 95 °C and 30 s at 60°C in a 96-well reaction plate. Specificity of PCR reactions was verified by melting curve analysis of each sample amplified product. Each real-time PCR reaction was performed in triplicate (technical replicates) on three individual plants (biological replicates).
Analysis of reference gene expression stability
The fluorescence raw data, generated by the Applied Biosystems StepOne software 2.0 (PE Applied Biosystems, Foster City, USA), were exported to the DART-PCR v.1.0 [29]. This programmed Excel workbook enables the rapid calculation of threshold cycles, amplification efficiencies and R0 (raw expression values) for every sample. Differences in amplification efficiency were assessed using one-way analysis of variance (ANOVA), based upon the null hypotheses: (i) that amplification efficiency is comparable within sample groups (outlier detection) and (ii) that amplification efficiency is comparable between sample groups (amplification equivalence). Outliers identified were omitted prior to further analysis.
Relative quantities were calculated in Microsoft Excel using the highest expression value as calibrator, and then they were imported to the gene expression stability program geNorm 3.5 [22], [31]. This program estimates an expression stability value (M) for each gene as the average pairwise variation for a particular gene with all the other tested reference genes included in the analysis. The M-value measure relies on the principle that the expression ratio of two ideal internal control genes is identical in all samples, regardless of the experimental conditions or treatments. Genes with the lowest expression stability are removed and a new M value for each of the remaining genes is calculated until only two genes remain. The geNorm program also establishes the minimal number of reference genes required for calculating an accurate normalization factor (NF). The NFn is calculated as the geometric mean of the expression level of the n most stable genes and recalculated adding the next most stable gene. Finally, geNorm calculates the pairwise variation, Vn/Vn+1, between two sequential normalization factors, NFn/NFn+1, reflecting the effect of the addition of a new gene on the normalization factor. A pair-wise variation (V) value higher than 0.15 implies that the added gene has a significant effect on normalization and should be included for calculation of a reliable normalization factor. Additional reference genes should be included to the normalization factor until the added gene has no significant effect. The optimum number of genes is the lowest number of genes with an acceptably low Vn/Vn+1.
NormFinder software [23], [32] was also used to evaluate the candidate normalization genes. The expression values obtained from DART-PCR were used as inputs for calculations. NormFinder has the ability to measure the gene expression stability taking into account intra- and inter-variations in defined samples groups or treatments. The combination of these estimates provides a direct measure of the variation in expression for each gene. This function determinates the best combination of two reference genes for normalization. The use of the two reference genes in combination will be, in most cases, more accurate than just using the most stable gene. The most stable candidate genes within and between the groups are those with the lowest variation values.
Determination of SGCesA expresión profile
The expression profile of the putative SGCesA was calculated in different tissue (leaf, stem, rachis and root) and leaves of plants under different abiotic stress conditions. Normalization factors were calculated as the geometric mean of the expression values (R0) of the reference genes tested. Relative expression levels of the target genes are calculated by dividing the expression value of the target gene by the normalization factor. The measures of expression were made in three biological replicates of three technical repeats per sample.
Results
Selection of candidate reference genes in P. virgatum, primer design and amplification specificity
Based on searches of the literature a group of previously used control genes in plants was selected [33]–[36]. The sequences for 15 candidate reference genes were identified after performing a TBLATN among the switchgrass GenBank ESTdb. From three to four primer pairs per gene were designed for qRT-PCR and finally, primer pairs for 11 genes were selected on the basis of their amplification efficiency and single PCR product. The candidate reference genes encode eEF-1α (elongation factor 1α), UBQ6 (ubiquitin 6), ACT12 (Actin 12), TUB6 (β-tubulin 6), eIF-4a (eukaryotic initiation factor 4a), GAPDH (glyceraldehyde-3-phosphate dehydrogenase C1), SAMDC (s-adenosyl methionine decarboxylase), TUA6 (α-tubulin 6), CYP5 (cyclophilin 5), U2AF (Splicing factor U2af) and FTSH4 (FtsH protease 4) are described in Table 1. The products of these genes are associated with a wide variety of biological functions. Specificity of real-time PCR products was confirmed by the presence of a single peak in the melting curve and the presence of a single band with the expected size (90 to 100 bp) after agarose gel electrophoresis and GelRed nucleic acid staining (Figure 1). The primer sequences and amplification efficiencies are indicated in Table 2.
A) Melting curve of the 11 reference genes showing a single pick (each including three technical replicates of the cDNA pool of the total samples used in this study). B) Agarose gel (1.5%) showing amplification of a specific PCR product of the expected size for each gene tested in this study.
The transcript abundance of these 11 potential reference genes was analyzed by qRT-PCR across a series of three biological repeats of 10 different samples obtained from: leaves from plants under different abiotic stress conditions and leaves, stems, roots and rachis from well watered plants.
Expression profile of the reference genes
The cycle threshold (Ct) values in the qRT-PCR reactions are defined as the cycle at which the fluorescent signal is significantly different from the background and is used to identify the differences in transcript expression levels. The Ct mean values for the 11 reference genes on the 10 samples under study were used to compare the expression rates between genes and within the set of samples. The data analysis showed that candidate reference genes were moderately abundant and exhibited a wide range of expression levels in all tested samples. As can be seen in Figure 2, most of the selected genes presented Ct values that range from 19.6 to 27.1, while majority of the values were distributed between 21 and 24. The expression variation of an individual gene within all samples was smaller in U2AF, FTSH4 and GAPDH and higher in TUB6, TUA6 and SAMDC. The lowest and highest variation values were 0.88 and 5.10 cycles corresponding to U2AF and TUB6, respectively. The most highly expressed genes were UBQ6 and SAMDC having a mean Ct value of 21 cycles. eIF-4a showed the lowest transcript abundance with a Ct of 27 cycles. Although some of these genes showed a low expression variation between the tested samples, it would be necessary to perform further analyses in order to identify the most suitable combination of these genes candidates for normalizing gene expression under certain experimental conditions.
Box-plot graph of Ct value shows the median values as lines across the box. Lower and upper boxes indicating the 25th percentile to the 75th percentile. Whiskers represent the maximum and minimum values.
Expression stability analysis
Two publicly available software tools, geNorm [22] and NormFinder [23] were used to analyze the gene expression stability and the number of reference genes necessary for accurate gene-expression profiling. Expression values were obtained using the DART-PCR analysis excel software which takes in account the amplification efficiency of each gene to obtain relative quantities [29].
We used geNorm to analyze gene expression stability across different sets of samples: 1) all the samples, 2) abiotic stress samples 3) total tissues samples and 4) leaves and stems tissue samples. The genes were ranked accordingly to the average expression M values (Figure 3A). The 11 genes used for the analysis showed high expression stability and presented an M value lower than the cutoff established by Vandesompele et al. [22] (M<1.5), otherwise the most stable reference genes were not identical in all the sets of samples under study. When all the samples were taken together, the average expression stability value (M) for UBQ6 and CYP5 was the lowest, and that of the SAMDC was the highest, indicating that UBQ6 and CYP5 had the most stable expression and that SAMDC was the most variably expressed of the 11 candidate genes. The genes U2AF and FTSH4 had the highest expression stability in samples under abiotic stress conditions and once again, SAMDC was the least stable gene. When only plant tissue samples were considered, UBQ6 and CYP5 were the best candidates for normalization and in the subset of samples composed by samples from leaves and stems the most stable genes were ACT12 and CYP5 while in both cases, TUB6 was the worst gene (Figure 3A).
A. Ranking of the gene expression stability performed in all the samples, abiotic stress samples, tissue samples and leaf and stem samples. The least stable genes are on the left and the most stable genes on the right. B. Pairwise variation (Vn/Vn+1) was analyzed between the normalization factors NFn and NFn+1. Asterisk indicates the optimal number of reference genes required for normalization.
The geNorm program was also used to calculate the optimal number of reference genes required for accurate normalization. Pair-wise variation analysis (Figure 3B) showed that the ideal number of reference genes may be different for distinct set of samples. When evaluating pairwise variation in all plant samples data set, V2/3 value was 0.213 and V3/4 was 0.15. The V3/4 value is equal to the cut-off value of 0.15, indicating that the use of the three most stable reference genes, UBQ6, CYP5 and eIF-4a, is sufficient for accurate normalization. On the other hand, only two control genes are sufficient for accurate normalization in tissues (UBQ6 and CYP5) or tissues under abiotic stress (U2AF and FTHS4), as indicated by the V2/3 values lower or equal than 0.15 in Figure 3B.
Stability of expression was then re-analyzed using the NormFinder algorithm that ranks the best candidate reference genes according to their minimal combined inter- and intra-group variation of expression for normalization factor (NF) calculation. Results are summarized in Table 3. When all the 30 samples data set or tissue samples subset were considered, eIF-4a and eEF-1α were identified as the most stable genes with a stability value between 0.379 to 0.299. The best combination of two genes was eEF-1α and ACT12 which improved the stability value to 0.224 and 0.206, respectively. Although eEF-1α was again identified as the most stable gene for the subset of samples from leaves and stems, a different best combination of two genes was recommended for this set, for which NormFinder proposed eEF-1α and CYP5. Finally, in the samples including vegetative tissue under abiotic stress, NormFinder identified FTHS4 as the most stably expressed gene and the best combination of two genes recommended was CYP5 and FTHS4. Out of the four datasets, SAMDC and TUB6 were the most often identified as the least stable reference genes.
Reference genes validation by quantification of SGCesAs expression profile with different normalization factors
In order to validate the use of the genes selected by geNorm and NormFinder algorithms as reference genes, we examine the effect of using different HKGs on the qRT-PCR-determined expression profile of putative switchgrass homologues to Cellulose Synthase. By assembling the switchgrass ESTs sequences with significant similarity to OsCesA genes, a set of seven contigs were obtained representing potential unique members of the SGCesA family. Further deduced protein sequence analysis of the SGCesA revealed that they contain the characteristic domain structure of CesA proteins [37], [38] that includes four conserved regions (U1, U2, U3 and U4), each containing a D residue or QXXRW sequence (D-D-D-QXXRW motif). The hypervariable region HVR-II located between the U2 and U3 which is found only in plant CesAs is also present (File S2). Gene-specific primers located in the hypervariable region HVR-II were designed to analyze expression patterns of the identified SGCesA genes.
To confirm and compare the results from the SGCesAs expression, different combinations of reference genes were used for normalization. When the set of samples from different tissues was analyzed, normalization was accomplished by using the geometric mean of eEF-1α and ACT12 or UBQ6 and CPY5, which are the best combination of genes determined by NormFinder and geNorm respectively. For the set of samples subjected to abiotic stress, the expression normalization was performed by using CYP5 and FTHS4 (NormFinder) or U2AF and FTHS4 (geNorm). When all of samples under study were analyzed, normalization was calculated by using eEF-1α and ACT12, determined by NormFinder or UBQ6, CYP5 and eIF-4a, the three most stable genes identify by geNorm and recommended by this program after the pairwise variation analysis (Figure 3B). As previously observed in other plant species [39]–[42], the expression levels for the SGCesA family members show variation among different tissues and abiotic stress conditions. The results reveal that the expression level of SGCesA1 and SGCesA7 was the highest in the stems. SGCesA7 was also highly expressed in roots along with SGCesA2. On leaves and rachis a low transcript expression was detected for all the SGCesAs genes (Figure 4a, b, d, e). Among the leaf samples under abiotic stress conditions, expression analysis show that, in most of the tissues, SGCesA5, 6 and 7 transcript levels are higher than the other genes. Especially high is the level of expression of SGCesA5 and SGCesA7 in leaves under flood and wounding stress conditions (Figure 4a, c, d, f). A similar expression pattern was obtained when the best combination of genes identified by geNorm and NormFinder was used. However, a marked increase in the transcript abundance was observed in the samples normalized using the genes proposed by geNorm than the genes proposed by NormFinder (Figure 4a–f). Normalization based on the less stable genes show that the expression pattern varied considerably or the expression level was underestimated (Figure 4g–i).
A. Gene expression normalized with the best combination of reference genes selected by NormFinder for all set of samples under study. B. Gene expression normalized with the best combination of reference genes selected by NormFinder for tissue samples. C. Gene expression normalized with the best combination of reference genes selected by NormFinder for samples under abiotic stress. D. Gene expression normalized with the most stable reference genes selected by geNorm for all samples under study. E. Gene expression normalized with the most stable reference genes selected by geNorm for tissue samples. F. Gene expression normalized with the most stable reference genes selected by geNorm for samples under abiotic stress. G. Gene expression normalized with the less stable genes identify by both algorithms for all samples under study. H. Gene expression normalized with the less stable genes identify by both algorithms for tissue samples. I. Gene expression normalized with the less stable genes identify by both algorithms for samples under abiotic stress.
Genes in the cellulose synthase (CesA) super family have been associated with the primary or secondary cell wall cellulose synthesis. Based on sequence similarities to the known CesA genes from Arabidopsis and other species [13], [43], [44] the identified SGCesA2-7, expressed in most of the tissues, are potentially involved in the primary wall cellulose synthesis. On the other hand, SGCesA1, show similarities to CesA involved in the secondary cell wall synthesis. This gene is highly expressed only in stems representing initiation of vascular development.
Discussion
With the interest in improving switchgrass as a forage or bioenergy crop the need for the development of new tools for genomics has been emphasized [2], [45], [46]. Gene expression analyses could lead to a better understanding of processes involved in commercially-important traits in switchgrass. Quantitative Real Time PCR (qRT-PCR) provides an important tool for gene expression analysis. Although some studies on switchgrass gene expression have used polyubiquitin, actin, eIF-4a, GAPDH or eEF-1α as reference genes [7], [25]–[28] the validation of a suitable set of reference genes is of critical importance. The switchgrass transcriptome analysis and the development of a gene expression atlas (PviGEA) [16] allowed the identification of a group of stably expressed genes in samples of tissues under different stages of development. These resources are useful to select candidate reference genes for gene expression normalization in this specie but a validation of their expression stability is needed across experimental conditions and tissues.
The present work is the first detailed study on the stability of a set of genes that aims to identify a set of control genes for normalization of transcript levels in switchgrass. In this study, we selected a series of candidate reference genes for which sequence information was available in the switchgrass dbEST. After determination of primer amplification efficiencies, 11 candidates were selected for evaluation of their normalization potential along a group of 30 samples from different tissues and experimental conditions.
Since there is no universal accepted method, the comparison of two different computer programs, geNorm [22] and NormFinder [23], was used to select the best internal controls for normalization of gene expression studies in switchgrass. Differences were observed after comparison of the ranking of the candidate reference genes from these two programs (Figure 3, Table 3). This is an expected outcome due to the different mathematical approaches to calculate stability associated with each method [47]–[49]. Nevertheless, the expression stability analysis showed that six of the top seven most stable genes selected in each sample set were common for the two programs. Genes as eEF-1α, e-IF-4a, CYP5 and U2AF are ranked by both algorithms among the most stable genes in the sample sets under study. The eEF-1α gene has been identified as a stable reference gene in other species including other members of the Poaceae family such as Brachiaria [50], lolium [35], perennial ryegrass [51] and rice [33]. Although less frequently, eIF-4a has been also identified as a good reference gene in some tissues as stems in flax [52] and in rice seeds [53]. The cyclophilins have been used in some reference genes studies where they have been identified usually as variable genes [52], [54], [55], however in our results CYP5 is included among the most stable genes. A new reference gene identified in rice, U2AF [34], has been also reported as one of the most stably expressed gene in switchgrass, being an example of information transfer from a model plant to a crop. The ranking of the least stable genes determined by NormFinder was the same as in geNorm, both programs discard the use of SAMDC and TUB6 because they display unacceptable expression stability limiting its use as internal control in P.virgatum. On the other hand, genes used previously as reference genes in switchgrass such as eEF-1α and eIF-4a [7], [27], are classified among the most stable genes by geNorm and NormFinder but other genes, as actin and GAPDH [25]–[28], show variations in the transcript levels in the set of leaf samples under abiotic stress conditions and leaf/stem, respectively. These results suggest that genes used commonly as reference genes should be used with caution requiring a careful evaluation for every set of samples under different experimental conditions before use.
Although the ranking of the different reference genes provided by both programs is similar, our results show that the genes that are most appropriate for use as reference genes can change, in the same set of samples and that each condition require a specific set of reference genes (Figure 3, Table 3). For tissue samples, geNorm recommended the use of UBQ6 and CYP5 and NormFinder recommended eEF-1α and ACT12. For leaves samples of plants under abiotic stress conditions CYP5 and FTHS4 were selected by geNorm and U2AF and FTHS4 for NormFinder. For all set of samples under study, eEF-1α and ACT12, were determined by NormFinder and UBQ6, CYP5 and eIF-4a by geNorm. To determine the effectiveness of selected reference genes, we evaluated their efficacy to normalize seven members of the SGCesA gene family in samples of different tissues and under different abiotic stress conditions (Figure 4). The quantification of SGCesA using different control genes showed that the choice of the internal standard is very important in the normalization of the target gene expression levels. Similar expression patterns were generated when the best combination of genes selected by both programs was used but differences were detected in the transcript abundance of the target genes (Figure 4 a–f). This is due to the highest transcript level of the reference genes selected by NormFinder compare with the genes selected by geNorm.
To quantifying expression, data analysis was performed by a relative approach developed in DART-PCR [29]. The starting fluorescence (R0), which is proportional to the starting template quantity, is calculated taking in account the threshold cycle (Ct), the fluorescence at this cycle (RCt) and the amplification efficiency (E), based on the equation: R0 = RCt×(1+E)−Ct. The R0 of each sample was normalized by the normalization factor obtained as the geometric mean of the reference genes R0 values. Any difference in the Ct, E and RCt values is able to produce differences on the resulting R0 value and consequently, on the normalized expression values. As an example, we compare the normalized samples subjected to flooding conditions. Despite the similar Ct values of the reference genes selected for normalization by NormFinder (CYP5, 23.88 and FTHS4, 22.81) and geNorm (U2AF, 24.16 and FTHS4 22.81) the normalization factors R0 values were 1.24×10−8 and 8.50×10−9, respectively. The use of these NFs results in a difference of approx. 30% between the relative gene expressions of the genes (Figure 4 c, f). We consider that the genes selected by both programs are valid for normalization but when the transcript abundance of the gene under study is low, reference genes with also low abundance will improve the measurement accuracy. Next we used the genes identified by both programs with the lowest expression stability for normalization (Figure 4g–i). Expression patterns varied considerably in comparison with the result from the use of the most stable genes. For example, the low expression level of TUB6 in leaf samples leads to an overestimation of SGCesA5 (Figure 4h) and the high transcript level of SAMDC gives an underestimation of the SGCesAs in samples under abiotic stress were the transcript abundance is undetectable in some samples (Figure 4i). It suggested that not only the stability but also the abundance of a reference gene affect the normalized results and demonstrate the need to evaluate the stability of the reference genes before to be used in a set of samples.
Conclusions
To the best of our knowledge, this work is the first study aimed to validate a set of candidate reference genes for gene expression normalization using qRT-PCR in switchgrass. In the present research 11 candidate reference genes were identified in the EST database of switchgrass and their expression stability was evaluated across a set of 30 samples using the computer programs geNorm and NormFinder. The use of the reference genes selected by both programs to determine the expression profile of the target genes allowed us to determine the most adequate combination of control genes. These results may constitute a starting point to select reference genes in the future to more accurate normalization in other tissues and under other experimental conditions in switchgrass.
Supporting Information
File S1.
SGCesA partial sequences. Consensus sequences after Switchgrass ESTs assembly. In gray gene specific primer in the HVR-II of the SGCesAs.
https://doi.org/10.1371/journal.pone.0091474.s001
(PDF)
File S2.
Alignment of the partial-length deduced amino acid sequence of putative P. virgatum cellulose synthase proteins. Amino acid sequence alignment of the highly conserved U1, U2 and U3/4 regions of the predicted SGCesAs. Domains are indicated below the sequence as follows: shaded amino acid with light gray are U domains, the D, D, D, QXXRW motif is indicated in red, the conserved region is underlined with light gray shading within a box, and the variable region (HVR-II) with dark gray shading without a box.
https://doi.org/10.1371/journal.pone.0091474.s002
(PDF)
Table S1.
Primer sequence, amplicon and efficiency of the SGCesA genes evaluated in this study.
https://doi.org/10.1371/journal.pone.0091474.s003
(PDF)
Acknowledgments
We would like to thank Dr. Peetambar Dahal for early discursion and comments and for his helpful advice on qRT-PCR experiments and data analysis.
Author Contributions
Conceived and designed the experiments: JG. Performed the experiments: JG. Analyzed the data: JG. Contributed reagents/materials/analysis tools: AV EB. Wrote the paper: JG. Responsible for the biological assays and carrying out the tissue collections and sample preparation: JG NE. Supervised the study and critically revised the manuscript: AV EB.
References
- 1. Schmer M, Vogel K, Mitchell R, Perrin R (2008) Net energy of cellulosic ethanol from switchgrass. Proceedings of the National Academy of Sciences of the United States of America 105: 464–469.
- 2.
Parrish D, Fike J (2005) The biology and agronomy of switchgrass for biofuels. Critical Reviews in Plant Sciences: 423–459.
- 3. Muir JP, Sanderson MA, Ocumpaugh WR, Jones RM, Reed RL (2001) Biomass production of ‘Alamo’ switchgrass in response to nitrogen, phosphorus, and row spacing. Agronomy Journal 93: 896–901.
- 4. Heaton E, Voigt T, Long SP (2004) A quantitative review comparing the yields of two candidate C-4 perennial biomass crops in relation to nitrogen, temperature and water. Biomass & Bioenergy 27: 21–30.
- 5. Vogel KP, Brejda JJ, Walters DT, Buxton DR (2002) Switchgrass biomass production in the Midwest USA: Harvest and nitrogen management. Agronomy Journal 94: 413–420.
- 6. Guretzky JA, Biermacher JT, Cook BJ, Kering MK, Mosali J (2011) Switchgrass for forage and bioenergy: harvest and nitrogen rate effects on biomass yields and nutrient composition. Plant and Soil 339: 69–81.
- 7. Fu CX, Mielenz JR, Xiao XR, Ge YX, Hamilton CY, et al. (2011) Genetic manipulation of lignin reduces recalcitrance and improves ethanol production from switchgrass. Proceedings of the National Academy of Sciences of the United States of America 108: 3803–3808.
- 8. DeMartini JD, Pattathil S, Miller JS, Li HJ, Hahn MG, et al. (2013) Investigating plant cell wall components that affect biomass recalcitrance in poplar and switchgrass. Energy & Environmental Science 6: 898–909.
- 9. Okada M, Lanzatella C, Saha MC, Bouton J, Wu RL, et al. (2010) Complete Switchgrass Genetic Maps Reveal Subgenome Collinearity, Preferential Pairing and Multilocus Interactions. Genetics 185: 745–760.
- 10. Tobias CM, Hayden DM, Twigg P, Sarath G (2006) Genic microsatellite markers derived from EST sequences of switchgrass (Panicum virgatum L.). Molecular Ecology Notes 6: 185–187.
- 11. Wang YW, Samuels TD, Wu YQ (2011) Development of 1,030 genomic SSR markers in switchgrass. Theoretical and Applied Genetics 122: 677–686.
- 12. Matts J, Jagadeeswaran G, Roe BA, Sunkar R (2010) Identification of microRNAs and their targets in switchgrass, a model biofuel plant species. Journal of Plant Physiology 167: 896–904.
- 13. Tobias CM, Sarath G, Twigg P, Lindquist E, Pangilinan J, et al. (2008) Comparative Genomics in Switchgrass Using 61,585 High-Quality Expressed Sequence Tags. Plant Genome 1: 111–124.
- 14. Palmer NA, Saathoff AJ, Kim J, Benson A, Tobias CM, et al. (2012) Next-Generation Sequencing of Crown and Rhizome Transcriptome from an Upland, Tetraploid Switchgrass. Bioenergy Research 5: 649–661.
- 15.
Wang YX, Zeng X, Iyer NJ, Bryant DW, Mockler TC, et al.. (2012) Exploring the Switchgrass Transcriptome Using Second-Generation Sequencing Technology. Plos One 7..
- 16. Zhang JY, Lee YC, Torres-Jerez I, Wang MY, Yin YB, et al. (2013) Development of an integrated transcript sequence database and a gene expression atlas for gene discovery and analysis in switchgrass (Panicum virgatum L.). Plant Journal 74: 160–173.
- 17. Somleva MN, Tomaszewski Z, Conger BV (2002) Agrobacterium-mediated genetic transformation of switchgrass. Crop Science 42: 2080–2087.
- 18.
Bustin S (2000) Absolute quantification of mRNA using real-time reverse transcription polymerase chain reaction assays. Journal of Molecular Endocrinology: 169–193.
- 19.
Bustin SA (2002) Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems. Journal of Molecular Endocrinology 29..
- 20.
Huggett J, Dheda K, Bustin S, Zumla A (2005) Real-time RT-PCR normalisation; strategies and considerations. Genes and Immunity: 279–284.
- 21. Gutierrez L, Mauriat M, Guenin S, Pelloux J, Lefebvre JF, et al. (2008) The lack of a systematic validation of reference genes: a serious pitfall undervalued in reverse transcription-polymerase chain reaction (RT-PCR) analysis in plants. Plant Biotechnology Journal 6: 609–618.
- 22. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, et al. (2002) Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biology 3: RESEARCH0034.
- 23.
Andersen C, Jensen J, Orntoft T (2004) 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 Research: 5245–5250.
- 24. Pfaffl MW (2001) A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Research 29: 6.
- 25. Somleva MN, Snell KD, Beaulieu JJ, Peoples OP, Garrison BR, et al. (2008) Production of polyhydroxybutyrate in switchgrass, a value-added co-product in an important lignocellulosic biomass crop. Plant Biotechnology Journal 6: 663–678.
- 26. Escamilla-Trevino LL, Shen H, Uppalapati SR, Ray T, Tang YH, et al. (2010) Switchgrass (Panicum virgatum) possesses a divergent family of cinnamoyl CoA reductases with distinct biochemical properties. New Phytologist 185: 143–155.
- 27. Saathoff AJ, Tobias CM, Sattler SE, Haas EJ, Twigg P, et al. (2011) Switchgrass contains two cinnamyl alcohol deshydrogenases involved in lignin formation. BioEnergy Research 4: 120–133.
- 28.
Xu B, Sathitsuksanoh N, Tang YH, Udvardi MK, Zhang JY, et al.. (2012) Overexpression of AtLOV1 in Switchgrass Alters Plant Architecture, Lignin Content, and Flowering Time. Plos One 7..
- 29. Peirson S, Butler J, Foster R (2003) Experimental validation of novel and conventional approaches to quantitative real-time PCR data analysis. Nucleic Acids Research 31: e73.
- 30. Huang XQ, Madan A (1999) CAP3: A DNA sequence assembly program. Genome Research 9: 868–877.
- 31.
geNorm software website. Available: http://medgen.ugent.be/~jvdesomp/genorm/. Accesed 2011 Jan, 11.
- 32.
NormFinder software website. Available: http://www.mdl.dk/publicationsnormfinder.htm. Accesed 2011 Jan 11.
- 33. Jain M, Nijhawan A, Tyagi AK, Khurana JP (2006) Validation of housekeeping genes as internal control for studying gene expression in rice by quantitative real-time PCR. Biochemical and Biophysical Research Communications 345: 646–651.
- 34. Narsai R, Ivanova A, Ng S, Whelan J (2010) Defining reference genes in Oryza sativa using organ, development, biotic and abiotic transcriptome datasets. Bmc Plant Biology 10: 56–69.
- 35.
Dombrowski J, Martin R (2009) Evaluation of reference genes for quantitative RT-PCR in Lolium temulentum under abiotic stress. Plant Science: 390–396.
- 36. Hong S, Seo P, Yang M, Xiang F, Park C (2008) Exploring valid reference genes for gene expression studies in Brachypodium distachyon by real-time PCR. Bmc Plant Biology 8: 112–123.
- 37. Delmer DP (1999) Cellulose biosynthesis: Exciting times for a difficult field of study. Annual Review of Plant Physiology and Plant Molecular Biology 50: 245–276.
- 38. Richmond TA, Somerville CR (2000) The cellulose synthase superfamily. Plant Physiology 124: 495–498.
- 39. Holland N, Holland D, Helentjaris T, Dhugga KS, Xoconostle-Cazares B, et al. (2000) A comparative analysis of the plant cellulose synthase (CesA) gene family. Plant Physiology 123: 1313–1323.
- 40. Tanaka K, Murata K, Yamazaki M, Onosato K, Miyao A, et al. (2003) Three distinct rice cellulose synthase catalytic subunit genes required for cellulose synthesis in the secondary wall. Plant Physiology 133: 73–83.
- 41. Kalluri UC, Joshi CP (2004) Differential expression patterns of two cellulose synthase genes are associated with primary and secondary cell wall development in aspen trees. Planta 220: 47–55.
- 42. Ranik M, Myburg AA (2006) Six new cellulose synthase genes from Eucalyptus are associated with primary and secondary cell wall biosynthesis. Tree Physiology 26: 545–556.
- 43. Appenzeller L, Doblin M, Barreiro R, Wang H, Niu X (2004) Cellulose synthesis in maize: isolation and expression analysis of the cellulose synthase (CesA) gene family. Cellulose 11: 287–299.
- 44. Carroll A, Mansoori N, Li S, Lei L, Vernhettes S, et al. (2012) Complexes with Mixed Primary and Secondary Cellulose Synthases Are Functional in Arabidopsis Plants. Plant Physiology 160: 726–737.
- 45.
McLaughlin S, Kszos L (2005) Development of switchgrass (Panicum virgatum) as a bioenergy feedstock in the United States. Biomass & Bioenergy: 515–535.
- 46.
Bouton J (2007) The economic benefits of forage improvement in the United States. Euphytica: 263–270.
- 47.
Cruz F, Kalaoun S, Nobile P, Colombo C, Almeida J, et al.. (2009) Evaluation of coffee reference genes for relative expression studies by quantitative real-time RT-PCR. Molecular Breeding: 607–616.
- 48. Paolacci AR, Tanzarella OA, Porceddu E, Ciaffi M (2009) Identification and validation of reference genes for quantitative RT-PCR normalization in wheat. Bmc Molecular Biology 10: 27.
- 49. Artico S, Nardeli SM, Brilhante O, Grossi-de-Sa MF, Alves-Ferreira M (2010) Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real-time quantitative RT-PCR data. Bmc Plant Biology 10: 49–61.
- 50. Silveira ED, Alves-Ferreira M, Guimaraes LA, da Silva FR, Carneiro VTD (2009) Selection of reference genes for quantitative real-time PCR expression studies in the apomictic and sexual grass Brachiaria brizantha. Bmc Plant Biology 9: 10.
- 51. Lee JM, Roche JR, Donaghy DJ, Thrush A, Sathish P (2010) Validation of reference genes for quantitative RT-PCR studies of gene expression in perennial ryegrass (Lolium perenne L.). Bmc Molecular Biology 11: 14.
- 52. Huis R, Hawkins S, Neutelings G (2010) Selection of reference genes for quantitative gene expression normalization in flax (Linum usitatissimum L.). Bmc Plant Biology 10: 14.
- 53. Li Q, Sun S, Yuan D, Yu H, Gu M, et al. (2010) Validation of Candidate Reference Genes for the Accurate Normalization of Real-Time Quantitative RT-PCR Data in Rice During Seed Development. Plant Molecular Biology Reporter 28: 49–57.
- 54. Tong ZG, Gao ZH, Wang F, Zhou J, Zhang Z (2009) Selection of reliable reference genes for gene expression studies in peach using real-time PCR. Bmc Molecular Biology 10: 13.
- 55. Reid KE, Olsson N, Schlosser J, Peng F, Lund ST (2006) An optimized grapevine RNA isolation procedure and statistical determination of reference genes for real-time RT-PCR during berry development. Bmc Plant Biology 6: 11.