Physiologia Plantarum 152: 529–545. 2014

© 2014 Scandinavian Plant Physiology Society, ISSN 0031-9317

Identification and validation of reference genes for Populus euphratica gene expression analysis during abiotic stresses by quantitative real-time PCR Hou-Ling Wanga,b , Jinhuan Chena , Qianqian Tiana , Shu Wanga , Xinli Xiaa,∗ and Weilun Yina,b,∗ a b

National Engineering Laboratory for Tree Breeding, College of Biological Sciences and Technology, Beijing Forestry University, Beijing 100083, China Key Laboratory for Silviculture and Conservation, College of Forestry, Beijing Forestry University, Beijing 100083, China

Correspondence *Corresponding authors, e-mail: [email protected]; [email protected] Received 26 January 2014; revised 5 March 2014 doi:10.1111/ppl.12206

Populus euphratica is the only arboreal species that is established in the world’s largest shifting-sand desert in China and is well-adapted to the extreme desert environment, so it is widely considered a model system for researching into abiotic stress resistance of woody plants. However, few P. euphratica reference genes (RGs) have been identified for quantitative real-time polymerase chain reaction (qRT-PCR) until now. Validation of suitable RGs is essential for gene expression normalization research. In this study, we screened 16 endogenous candidate RGs in P. euphratica leaves in six abiotic stress treatments, including abscisic acid (ABA), cold, dehydration, drought, short-duration salt (SS) and long-duration salt (LS) treatments, each with 6 treatment gradients. After calculation of PCR efficiencies, three different software tools, NormFinder, geNorm and BestKeeper, were employed to analyze the qRT-PCR data systematically, and the outputs were merged by means of a non-weighted unsupervised rank aggregation method. The genes selected as optimal for gene expression analysis of the six treatments were RPL17 (ribosomal protein L17) in ABA, EF1𝛼 (elongation factor-1 alpha) in cold, HIS (histone superfamily protein H3) in dehydration, GII𝛼 in drought and SS, and TUB (tubulin) in LS. The expression of 60S (the 60S ribosomal protein) varied the least during all treatments. To illustrate the suitability of these RGs, the relative quantifications of three stress-inducible genes, PePYL1, PeSCOF-1 and PeSCL7 were investigated with different RGs. The results, calculated using qBasePlus software, showed that compared with the least-appropriate RGs, the expression profiles normalized by the recommended RGs were closer to expectations. Our study provided an important RG application guideline for P. euphratica gene expression characterization.

Introduction Perennial woody plants dominate lots of natural land ecosystems and play an important role in fixing carbon,

producing oxygen and protecting the environment. Different with herbaceous plants they usually have long life cycles, spanning decades of years, several centuries or even longer. Therefore, they suffer more challenges in

Abbreviations – Cpn60𝛽, 60-kDa chaperonin β-subunit; CT, cycle threshold; CV, coefficient of variance; EF1𝛼, elongation factor-1 alpha; HIS, histone superfamily protein H3; LIPL, lipocalin-like gene; LR, low-ranked reference gene; LS, long-duration salt; LTP, lipid transfer protein; miRNA, microRNA; NF, normalization factor; OD, optical density; PV, pairwise variation; qRT-PCR, quantitative real-time polymerase chain reaction; RA, rubisco activase; RG, reference gene; RP, ribosomal L27e protein; RPL17, ribosomal protein L17; RQ, relative quantification; RSMC, relative soil moisture content; SS, short-duration salt; TR, top-ranked reference gene; TUB, tubulin; UBQ, ubiquitin.

Physiol. Plant. 152, 2014

529

complicated natural environments during the entire year, including water deficiency, daily temperature swings, variational intensity of the sunlight and differential soil nutrient conditions, etc. Compared with biotic stresses, these abiotic stresses are the main factors that limit their distribution and affect their growth. Studies on the resistance mechanisms of woody plants to abiotic stresses will provide the possibility of cultivating fine varieties. Populus euphratica, which is characterized by a wide temperature range, as well as salinity and aridity tolerance, is the only arboreal species that is established in the world’s largest shifting-sand desert, the Taklimakan Desert in China (Gries et al. 2003) and is well-adapted to the extreme desert environment, so P. euphratica is increasingly considered as a model system for woody plant research into abiotic stress resistance (Ottow et al. 2005). However, its stress resistance mechanisms are far from understood, and many stress-related genes and microRNAs (miRNAs) need to be identified. Studies of P. euphratica have mainly involved analyses of single gene functions (Ottow et al. 2005, Chen et al. 2009, Ma et al. 2010), transcriptome sequencing (Tang et al. 2013) or miRNA identification (Li et al. 2011). In these studies the gene expression levels upregulated or downregulated by stresses are important and need to be calculated using appropriate normalization methods. Quantitative real-time polymerase chain reaction (qRT-PCR), which shows high sensitivity, specificity and reliability, has become one of the most commonly used gene expression profiles detection methods. It can detect low-abundance mRNAs and is capable of detecting slight variations in gene expression in different tissues or different treatments in the same plant (Gachon et al. 2004). Accurate and reliable gene expression data analyzed by qRT-PCR relies on moderate amounts of starting material, high quality RNA and efficient enzymes, appropriate normalization techniques and suitable reference genes (RGs) used as endogenous controls (Bustin 2002, Wan et al. 2010). Among these approaches, the use of suitable RGs for internal standardization of target gene expression is a basic prerequisite for any qRT-PCR (Pfaffl et al. 2004). RGs are usually used without validating their stability, as it is assumed that they are expressed at constant levels throughout different plant tissues, physiological conditions or developmental stages and are not influenced by exogenous treatments, but extensive transcriptomic data mining and RG validation in model plants have shown that the reliability of these endogenous controls can be influenced by the plant species, growth conditions and organs/tissues examined, and they exhibit either upregulation or downregulation under different 530

experimental conditions (Pfaffl et al. 2004, Huis et al. 2008). In qRT-PCR, inappropriate RGs can result in biased target gene expression profiles, so it is important to identify the best RGs to use in each biological system before performing gene expression normalization. In P. euphratica, no systematic study on validating RGs has been conducted as in Arabidopsis (Czechowski et al. 2005), rice (Jain et al. 2006), soybean (Jian et al. 2008), Oryza (Narsai et al. 2010), Gossypium (Artico et al. 2010), flax (Huis et al. 2010), petunia (Mallona et al. 2010), Eucalyptus (de Oliveira et al. 2012) and sesame (Wei et al. 2013). According to these previous studies, the most suitable RGs are those associated with basic cellular processes and which are constitutively expressed, such as those related to the ubiquitin degradation process, polyubiquitin, ubiquitin-conjugating enzymes and ubiquitin ligases (Czechowski et al. 2005, Jain et al. 2006, Hong et al. 2008, Silveira et al. 2009, Artico et al. 2010), and elongation factors (Jain et al. 2006, Jian et al. 2008, Silveira et al. 2009, Mallona et al. 2010). Nevertheless, different species, developmental stages, tissues and treatments tend to require different RGs, and few RGs have been validated that show stable expression across all experimental conditions. Candidate control genes need to be further evaluated with more comprehensive and accurate analytical techniques to determine the most appropriate one. To this end, programs such as geNorm (Vandesompele et al. 2002), NormFinder (Andersen et al. 2004), BestKeeper (Pfaffl et al. 2004) and qBasePlus (Hellemans et al. 2007) are commonly used for statistical analyses. In this work, we evaluated 16 potential RGs in P. euphratica leaves under six different abiotic stresses, including abscisic acid (ABA), cold, dehydration, drought, short-duration salt (SS) and long-duration salt (LS) treatments, each with six treatment gradients. The qRT-PCR data were then analyzed using above-mentioned four programs. The RGs determined as best or worst were used to investigate the expression of PePYL1 (homolog of AtPYL1/AtPYR1/AtRCAR12, E-value 2.2e-82; ABA receptor in Arabidopsis, responds to ABA stress) (Fujii et al. 2009, Ma et al. 2009, Park et al. 2009), PeSCOF-1 (homolog of GmSCOF-1, E-value 1.1e-78; cold-inducible zinc finger protein in soybean, responds to cold stress) (Kim et al. 2001) and the plant-specific GRAS/SCL transcription factor PeSCL7 (homolog of AtSCL7 identified in P. euphratica; can be induced by drought and salt stress and enhances salt and drought tolerance in transgenic Arabidopsis thaliana) (Ma et al. 2010). As a consequence PeSCL7 can be taken as a positive control for P. euphratica genes expression normalization during abiotic stresses.

Physiol. Plant. 152, 2014

Materials and methods Plant materials One-year-old seedlings of P. euphratica, obtained from the Xinjiang Uygur Autonomous Region of China, were planted in individual pots containing mainly sandy soil and placed in a greenhouse at Beijing Forestry University. The temperature in the greenhouse was 20–25∘ C with a 16-h photoperiod (4:00–20:00 h). Each pot contained three to five individuals, and three pots were used for each stress treatment gradient. For expression analysis, uniformly developed plants (60–80 cm high, with 70–80 leaves) were subjected to various stress treatments. For ABA treatment, the leaves of young plants were sprayed with a 200 μM ABA solution. For cold stress, plants were transferred to a room held at 4–6∘ C. Dehydration was performed by removing the plants from the plots and exposing them on filter paper to air at 70% relative humidity (RH) and 25∘ C. SS stress was performed by carefully taking plants from the soil and planting in water for 3 days, then the water was replaced with 350 mM NaCl solution and the salt water was changed every 6 h to avoid hypoxia stress. In these four treatments, mature leaves from the same position were collected at 1, 3, 6, 9 and 12 h after the treatment. Each treatment gradient was performed within three pots with 9–15 individuals and non-treated plants were used as controls. In the drought-stress treatment, P. euphratica plants were subjected to soil water deficiency at six relative soil moisture content (RSMC) levels for 1 month, based on previous research (Hsiao 1973). Plants grown in soil that was sufficiently irrigated each day were used as controls; the RSMC was maintained at 70–75% by transpiration. For serious drought stress, the RSMC was 10–15%, and the plants were divided into groups from the highest to the lowest RSMC (groups A–F). LS stress was conducted by watering young plants in soil with 200 mM NaCl solution, and leaves were gathered after 5, 10, 15, 20 and 25 days stress. The leaf water potential (WP) was measured using a PsyPro WP data logger (Wescor, Logan, UT). The net photosynthetic rate, light response curve, water conductance, intercellular CO2 concentration and transpiration rate were measured using the Li-6400 Photosynthesis System (Li-Cor, Lincoln, NE). All samples were immediately frozen in liquid nitrogen and stored at −80∘ C. Total RNA isolation, quality control and cDNA synthesis Frozen samples were ground in liquid nitrogen to a fine powder with a pestle and a mortar. Total RNA was extracted using the cetyltrimethyl ammonium bromide

Physiol. Plant. 152, 2014

(CTAB) method for plants (Chang et al. 1993), dithiothreitol (DTT) was used in place of 𝛽-mercaptoethanol. Potentially contaminating DNA was eliminated by treatment with DNase I in the last step. The quality and quantity of RNA were assessed with a NanoDrop2000 spectrophotometer (Thermo, West Palm Beach, FL) by determining the optical density (OD)260 /OD280 and OD260 /OD230 ratios, respectively. The RNA samples were also examined by 1% agarose gel electrophoresis for 15 min and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) according to the manufacturer’s instructions; 1.8 μg of total RNA was reverse-transcribed using the TIANGEN FastQuant RT Kit (with gDNase) (Qiagen, Düsseldorf, Germany) according to the manufacturer’s protocol; 20 μl of cDNA was diluted 1:10 with nuclease-free water. Selection of candidate RGs and primers design RGs assumed to have stable expression are usually used to normalize target gene expression in many species. We searched for RGs already identified in other species by performing a BLAST search of the Phytozome network on the Populus trichocarpa genome (http://www. phytozome.net/search.php?method=Org_Ptrichocarpa) and we also chose genes with low expression variation in previous high-throughput approaches (Gu et al. 2004, Bogeat-Triboulot et al. 2007, Brinker et al. 2010, Qiu et al. 2011). Three candidate genes already used as RGs for Populus were selected, including the glucosidase II 𝛼-subunit gene (Bogeat-Triboulot et al. 2007), the 18S ribosomal RNA gene (Brentner et al. 2008) and Actin (Gu et al. 2004, Chen et al. 2009). In addition, 13 novel candidate RGs were also selected, including 60S (the 60S ribosomal protein), Cpn60𝛽 (60-kDa chaperonin β-subunit), EF1𝛼 (elongation factor-1 alpha), eIF-5A (eukaryotic initiation factor 5A), GAPDH (glyceraldehyde-3-phosphate dehydrogenase), HIS (histone superfamily protein H3), LIPL (lipocalin-like gene), LTP (seed storage/lipid transfer protein), RA (rubisco activase), RP (ribosomal L27e protein), RPL17 (ribosomal protein L17), TUB (tubulin) and UBQ (ubiquitin). The products of these genes are involved in a wide variety of biological functions. For primers design, Primer6 software upgraded from Primer3 (Rozen and Skaletsky 2000) was used (amplicon length 60–200 bp, primer length 18–30 bp, optimal Tm (annealing temperature) at 60 ± 2∘ C, % GC between 40 and 60, avoid adenine at the 3′ end). Then, these primers were screened for hairpins, target specificity and dimer formation by BLASTN searches (http://www.ncbi.nlm.nih.gov/BLAST) against the P. trichocarpa databank. The selected primer pairs were also 531

tested for specificity by qRT-PCR, followed by dissociation curve observation and agarose gel electrophoresis. RT-PCR and qRT-PCR conditions The same cDNA samples and primers were used for both RT-PCR and qRT-PCR analyses. The PCR products were examined on a 2% agarose gel stained with GoldenView. Real-time PCR was performed using SYBR Green detection chemistry and running in triplicate in 96-well plates with the ABI StepOnePlus instrument (ABI, Foster City, CA). Reactions were prepared in a total volume of 20 μl containing 10 μl 2× TIANGEN SuperReal PreMix Plus, 2 μl ROX Reference Dye (Qiagen), 1 μl single stranded cDNA [single-stranded circular DNA (sscDNA), corresponding to 9 ng of total RNA] and 0.3 μM of each primer. Blank controls were run in triplicate for each master mix. The cycling parameters were: 95∘ C for 15 min, 45 cycles of 20 s at 95∘ C and 60 s at 60∘ C. Dissociation curves for each amplicon were then analyzed by heating the amplicon from 40 to 100∘ C and reading at each ∘ C to verify the specificity of each amplification reaction. Bioinformatics and statistical analysis of gene expression stability by the three tools The data analysis strategy is described in detail according to the flow chart. The PCR efficiency (E) was calculated by the equation below (Pfaffl 2001). ( ) CT = k lg X0 + b E=

−1 10 ∕K

(1) (2)

In Eqn 1 X0 is the absolute content of cDNA, which was diluted 4-, 8-, 16-, 32- and 64-fold, and in Eqn 1 and Eqn 2 k stands for the slope. In all 36 samples, the expression levels of the 16 RGs were determined by the cycle threshold (CT) values. The software packages used in this step including SigmaPlot v. 12.5 (http://www.sigmaplot. com/), geNorm v. 3.5 (http://medgen.ugent.be/∼jvde somp/genorm/), the Excel add-in of NormFinder v. 0.953 (http://www.mdl.dk/publicationsnormfinder.htm), BestKeeper v. 1 (http://www.gene-quantification.de/bestkee per.html) and qBasePlus v. 2.5 (http://www.biogazelle. com/). Gene suitability rankings merged by RankAggreg The RankAggreg v. 0.4-3 (http://cran.r-project.org/web/ packages/RankAggreg/) package of R program v. 3.0.1 (http://www.r-project.org/) was used to merge 532

the stability measurements obtained from the three Excel-based tools and establish a consensus rank of RGs (Pihur and Datta 2009). The rankings from four terms ‘Stab’, ‘CV’, ‘r2 ’ and ‘M’ were used as input. The R script file example is available in Appendix S1 (see Supporting Information) or in page 5 of the PDF (http://cran.rproject.org/web/packages/RankAggreg/RankAggreg.pdf). According to the size of the ranking list (16) we used Monte Carlo algorithm to calculate and visually present the rank aggregation by line chart. The rank aggregation aims at finding an aggregated ranking that minimizes the distance to each of the ranked lists in the input set. The distance among ordered lists is calculated using the Spearman foot rule function. Because geNorm yielded the same M stability value in the chart of the two most stable genes, we distinguished these two genes’ order by the initial M value when calculating the normalization factor (NF) value. Calculate relative quantification values using a single RG or a combination of multiple RGs The Pfaffl mathematical model was employed to calculate the relative quantification of a target gene in comparison to a single RG (Pfaffl 2001). The relative expression ratio (target gene)/(RG) is calculated based on E and CT values of an unknown sample vs a endogenous control. ( Etarget

Ratio = ( Ereference

)ΔCT

target

-

gene(control sample)

gene

)ΔCT

reference

-

gene(control sample)

(3)

gene

( )ΔCTt Et Ratio = ( )ΔCT r Er

(4)

Eqn 3 is abbreviated to Eqn 4, in Eqn 3 and 4 Etarget gene or Et is the real-time PCR efficiency of target gene transcript; Ereference gene or Er is the real-time efficiency of RG transcript; 𝛥CT target gene or 𝛥CT t is the CT deviation of control – sample of target gene transcript and 𝛥CT reference gene or 𝛥CT r is the CT deviation of control – sample of RG transcript. The RG is what we validated as the best or worst in this study. Real-time PCR efficiencies were calculated according to Eqn 1 and Eqn 2. When using a combination of multiple RGs the equation was transformed from Eqn 3. ( )ΔCT target gene(min−sample) Etarget gene Ratio = (5) Normalization Factor Ratio =

( )ΔCTt Et NF

(6)

Physiol. Plant. 152, 2014

Eqn 5 is abbreviated to Eqn 6, in Eqn 5 and 6 Etarget gene , Et , 𝛥CT target gene and 𝛥CT t had the same meanings as those in Eqn 3 or 4. NF is calculated using qBasePlus or geNorm software (geNorm is part of Biogazelle’s qBasePlus software now and is available for 299 EUR). In Eqn 5 ‘min’ is the minimum CT value of all samples. The CT values are transformed to quantities, and the highest relative quantities for each gene are set to 1. These raw not-yet normalized RG quantities are the required input. After calculation of a NF, the relative expression ratio [or named relative quantification (RQ)] can be calculated according to the Eqn 5. The RQ values were obtained by comparing with non-stressed plants. RQ values calculated using Eqn 3 for the non-stressed plants were 1 while using Eqn 5 the values were not convenient for comparing, after we got the values through Eqn 5 we divided the RQ values of stressed plants (sample) by the RQ values of non-stressed plants (control) and we got the final RQ values. Expression analysis of abiotic stress inducible genes In order to validate the applicability of the identified RGs, the relative abundance of three genes that are induced by abiotic stress were calculated. The genes used were three novel genes, PePYL1/PePYR1/PeRCAR12, PeSCOF-1 and PeSCL7 which were previously identified in P. euphratica (Ma et al. 2010). CT values and amplification efficiency values were processed with qBasePlus. The single RG refers to the ranks based on RankAggreg results, whereas the combinations were recommended by geNorm (the number) and RankAggreg (the RGs). The RQ values of these three genes in comparison to a single RG or a combination of multiple RGs are calculated based on Eqns 3 and 5.

efficiency, followed by GII𝛼 and 60S. The PCR efficiencies are listed in Table 1. For PCR specificity detection the presence of a single peak in the melting curve was obtained after amplification (Fig S2). The amplified products were also further analyzed by 2% agarose gel electrophoresis and GoldView staining and only one band of the expected size was detected in each experiment (Fig. 1). Expression profiles of the RGs To obtain the expression rates of the 16 selected genes, all PCR assays were performed in triplicate. The CT values for the 16 genes were pooled so that the expression level of a gene could be compared with each other’s. All data were analyzed according to the flow chart plotted in Fig. 2. The average CT values ranged from 12.05 to 34.89 in each treatment (Fig. 3A); 18S had the highest expression level, whereas TUB was the lowest. Except for 18S and TUB, most of the genes’ CT values were concentrated (between 20 and 30). The data distribution is shown as a box-plot in Fig. 3B; in all 36 tested samples, the 18S data was the most concentrated while the CT values of GAPDH and TUB were scattered. The greatest expression variation was displayed by Actin, GAPDH, RA and TUB (more than six cycles), whereas 18S, 60S, Cpn60𝛽, LIPL and UBQ (four or less cycles) were more stable. Analysis of the data by different treatments showed that the cold, SS had a relatively slight influence on the expression of the 16 RGs, and dehydration had a serious impact on the CT values [see standard deviation (SD) values in Fig. 3A or Table S1], with expression showing varying degrees of decline after the stress. These results indicated that the selected genes showed variable expression levels during the six abiotic stresses, so it was necessary to screen out the best RG for target gene expression normalization.

Results RNA quality, PCR specificity and amplification efficiency

Statistical analysis of qRT-PCR data using three bioinformatics programs

The quality and quantity of RNA were assessed with a NanoDrop2000 spectrophotometer (Thermo) by determining the OD260 /OD280 ratios which were between 1.8 and 2.0, and the OD260 /OD230 ratios were larger than 1.5. When detecting by gel electrophoresis and microfluidic devices, the integrities of RNA samples were good and little evidence of degradation was noted (Fig. S1). PCR specificity was confirmed by the use of primers designed to be specific based on gene searches (NCBI). The PCR efficiencies of these 16 candidate RGs varied from 1.896 to 2.097, among which EF1𝛼 had the lowest amplification efficiency and 18S showed the optimal

As the expression levels of the 16 RGs varied in the different abiotic stress treatments, we ranked their stabilities to determine the optimal RG. Data from each treatment were analyzed separately, and then added together. The data were transferred into three Excel-based programs (NormFinder, geNorm and BestKeeper) according to the flow chart shown in Fig. 2. CT values were first log-transformed and used as input for the NormFinder tool. This algorithm requires data from a minimum of three genes and a minimum of two samples per group. It can take into account intra- and inter-group variations for NFs. The software calculates

Physiol. Plant. 152, 2014

533

534

Physiol. Plant. 152, 2014

AY652861.1 (GenBank No.) Potri.007G093700.1 Potri.006G192700.1 AJ773373 (GenBank No.) Potri.006G130900.1 Potri.018G107300.1 Potri.010G055400.1

DQ388455.1 (GenBank No.) Potri.005G072300.1 AJ778489 (GenBank No.) AJ779446 (GenBank No.) AJ780799 (GenBank No.) Potri.001G342500.1 AJ777362 (GenBank No.) Potri.003G126800.1 Potri.014G115100.1 Potri.001G142500.1

Potri.002G119300.1 KJ511863 (GenBank No.)

GII𝛼 HIS LIPL LTP RA RP RPL17 TUB UBQ PePYL1

PeSCOF-1 PeSCL7

P. trichocarpa transcript name

18S 60S Actin Cpn60β EF1𝛼 eIF-5A GAPDH

Gene abbreviation

AT1G27730 AT3G50650

AT5G63840 AT4G40040 AT5G58070 AT5G05960 AT2G39730 AT3G22230 AT1G04270 AT5G23860 AT5G03240 AT5G46790

AT3G41768 AT5G65220 AT3G12110 AT5G56500 AT1G07940 AT1G13950 AT3G04120

Arabidopsis ortholog locus 18S ribosomal RNA 60S ribosomal protein L35 Actin 3 60-kDa chaperonin β-subunit Elongation factor-1 alpha Eukaryotic initiation factor 5A Glyceraldehyde-3-phosphate dehydrogenase Glucosidase II a-subunit Histone superfamily protein H3 Lipocalin-like gene Seed storage/lipid transfer protein Rubisco activase Ribosomal L27e protein family Ribosomal protein L17 Tubulin beta chain Ubiquitin family START domain protein/ABA receptor Zinc finger protein Plant-specific GRAS/SCL transcription factor

Tentative annotation

AGAGGTGTCACTACGAAGGCATC/CTGGCAAGGCTGGAACATTAAGGT CGAGAATAAAGAGCAATGGAGGG/ATGACGAAACTGTG AGTAGTGG

CTCTCATTGAGCCGGCAAAT/CCCCCCTTCAAGCATAAGG ACTGCTCGTAAGTCTACTGGAGG/GCGGTAACGGTGAGGCTTCTTC CGCCATTCTTGCCCATCATTCCT/TCCCTGTTCTTTCGCCTTCTCC GATCGTCCCGTGGGCTACAAGT/GATGTGAACTCCAGCAGGTTGATAGG GTGGGTCTCAGGTGTTGGTGTTG/CTCTTCACATTCTCTTGCTCCTTGAC GTTACACGCTGGATGTGGACTTG/AACCACCTGTTCTTGCCTGTCTT GCATAAGGGAGGGTCAAATGGAAGT/TCATTCTCTGCTCTTGCCTTTACG GGAGGTGGAACTGGATCAGGAATG/GGCATTGTAAGGCTCAACCACTGT AGACCTACACCAAGCCCAAGAAGAT/CCAGCACCGCACTCAGCATTAG GCCACAGACCTACAAGCACTTCAT/GACCTGCCGTTCATCGTCCAAT

TCAACTTTCGATGGTAGGATAGTG/CCGTGTCAGGATTGGGTAATTT AGGTGAACTCTTGATGCTTCGTCTT/CTTCTCTTCCATTGCCTGTCCAACT AAGATTCCGTTGTCCAGAGGTCCT/GAACATAGTAGAGCCACCACTGAGAAC TGCTCGCAAGTGACAATCCTAAGTT/GCACCGATTCAGGCTCCTTAATGT TCCGTCTTCCACTTCAGGATGTCT/GTCACGACCATACCAGGCTTCAG TCGGACGAGGAGCACCACTT/TGCAAGGACGGTTCTTGATGACTAT ATGAAGGACTGGAGAGGTGGAAGG/CACAGTAGGAACACGGAAGGACATT

Primer sequences (forward/reverse)

Table 1. Gene description, primer sequences, amplicon length and PCR efficiency for Populus euphratica RGs selection.

147 197

65 111 179 112 150 136 163 138 139 154

145 177 160 189 100 118 138

1.886 2.043

2.009 1.968 1.930 1.911 1.902 1.935 1.911 1.914 2.097 2.083

2.008 1.972 2.027 1.922 1.896 1.965 1.904

Amplicon PCR length (bp) efficiency

Fig. 1. Quantitative real-time PCR amplification specificity and size. Amplified fragments obtained after qRT-PCR were separated by 2% agarose gel electrophoresis. Specific amplification primers were designed and used for the 16 candidate RGs.

CT

NormFinder

efficiency

BestKeeper

geNorm

scatter-plot Stab

CV

r2

M

box-plot RankAggreg

gene RQ

PV

qBasePlus

Fig. 2. Data analysis flow chart. CT values were calculated using different algorithms. Each candidate RG has one efficiency value. The CT data were checked for distribution by scatter diagram and box plot. Circles indicate statistical results to be merged with RankAggreg. Stab, NormFinder stability value; r2 , determination coefficient – regression from BestKeeper; M, classical stability value.

and ranks the stability values, and the lowest stability value indicates the most stably expressed gene. The results of the NormFinder analysis are shown in Table 2. In ABA treatment, the best RG was RPL17 (stability value 0.166). In cold treatment, the best RG was GII𝛼, GII𝛼 was also detected as the most stably expressed gene during drought and SS stresses, and in dehydration and LS treatments, the most optimal RGs were HIS (0.058) and 18S (0.037), respectively. Among all treatments, EF1𝛼 varied the least and TUB varied the most. The stability value of the total was higher than those of the separate treatments. The average expression stability M values of all genes were calculated using geNorm. M values are defined as the mean variation of a gene in relation to all others. The lower the M value, the more stable the gene. To identify RGs with stable expression, this software recommends an M value below a threshold of 1.5. Some authors (Coll et al. 2010, de Almeida et al. 2010, Mallona et al. 2010, Taylor et al. 2010) are accustomed to using 0.5 as the threshold limit. As shown in Fig. 4, in the SS treatment, 13 of 16 genes showed highly stable expression, with the thresholds below 0.5, which was similar to the results in

Physiol. Plant. 152, 2014

the cold treatment: 9 of 16 gene thresholds were below 0.5. In contrast, in the dehydration treatment, the stability values were very variable; the M values for 14 genes were greater than 0.5. According to the geNorm analysis, the two best RGs for qRT-PCR normalization were HIS and UBQ for ABA, EF1𝛼 and RA for cold, 60S and EF1𝛼 for dehydration, Actin and GII𝛼 for drought, GII𝛼 and RPL17 for SS, and RP and TUB for LS. Overall, 60S and UBQ were the most stable genes. BestKeeper ranks the stability values using the coefficient of variance (CV) and the SD. RGs are identified as the most stable genes when they exhibit the lowest CV and SD (CV ± SD). Another function is intended to establish the best-suited standards out of the candidates and to merge them in a NF called the BestKeeper index, which is designed to determine a reliable NF but not to evaluate the goodness of each candidate gene independently, called the coefficient of determination (r2 ). The larger r2 the better (closer to 1). In order to obtain more credible results, we used both the CV and the r2 to rank the stability values of the 16 candidate RGs. The ranks are shown in Table 3, and the detailed calculated process and values can be seen in the Table S2. None of these genes had the top rank by both the CV and the r2 . For instance, with ABA treatment, considering only the CV, 60S was the most stably expressed gene, but ranked by r2 , RPL17 was the best. According to the CV values, BestKeeper calculated the greatest stability for RP in cold, TUB in dehydration and LS, EF1𝛼 in drought, HIS in SS and 18S in total. In contrast, ranked by the r2 , GII𝛼 was the most stable in drought, short-duration and total, LTP was the most stable during cold, HIS was the most stable in dehydration treatment and RP was the most stable in the LS treatment. Determination of the number of RGs to use for normalization After extensive studies to ensure stability, a single RG is generally used for target gene normalization in most qRT-PCR research. However, quantification of gene expression relative to multiple RGs has increasingly 535

A

α

β

α

α

β

α

B

Fig. 3. qRT-PCR CT values for the 16 candidate RGs. (A) Mean values in different treatments. (B) Box-whisker plot showing the CT variation among 36 test samples. A line across the box depicts the median. In each box, the upper and lower edges indicate the 25th and 75th percentiles. Whiskers go from the minimal to maximal value or, if the distance from the first quartile to the minimum CT value is more than 1.5 times the interquartile range (IQR), from the smallest value included within the IQR to the first quartile. Dots indicate outliers (values smaller than the lower quartile by 1.5 times the IQR or larger than the upper quartile by 1.5 times the IQR). Table 2. Populus euphratica RGs for normalization and their expression stability values calculated using the NormFinder program. Rank position

Abscisic acid Gene

Stab

Cold Gene

Dehydration Stab

Gene

Stab

Drought Gene

Stab

1 RPL17 0.166 GII𝛼 0.048 HIS 0.058 GII𝛼 0.037 2 RP 0.187 60S 0.057 eIF-5A 0.092 Actin 0.075 3 Cpn60𝛽 0.255 EF1𝛼 0.091 Actin 0.111 60S 0.101 4 LIPL 0.267 HIS 0.110 UBQ 0.119 RA 0.120 5 18S 0.268 Actin 0.117 GII𝛼 0.121 GAPDH 0.123 6 UBQ 0.269 RA 0.122 EF1𝛼 0.148 UBQ 0.124 7 GII𝛼 0.277 RP 0.139 60S 0.190 Cpn60𝛽 0.128 8 EF1𝛼 0.288 eIF-5A 0.149 LIPL 0.248 HIS 0.134 9 HIS 0.309 UBQ 0.152 RP 0.255 RP 0.138 10 GAPDH 0.319 GAPDH 0.164 GAPDH 0.287 eIF-5A 0.157 11 eIF-5A 0.334 18S 0.246 LTP 0.309 RPL17 0.159 12 60S 0.341 TUB 0.253 RPL17 0.356 EF1𝛼 0.183 13 TUB 0.346 Cpn60𝛽 0.286 RA 0.363 TUB 0.223 14 RA 0.376 LIPL 0.341 Cpn60𝛽 0.457 LTP 0.293 15 Actin 0.403 RPL17 0.497 TUB 0.570 18S 0.382 16 LTP 0.439 LTP 0.527 18S 0.928 LIPL 0.456 Best combination RP and RPL17 60S and GII𝛼 eIF-5A and HIS Actin and GII𝛼 0.125 0.038 0.055 0.043

been accepted. This computational method requires the calculation of a NF that merges data from several internal genes. The minimum number of RGs that should be used is estimated by computing the pairwise variation (PV) of two sequential NFs (Vn/n + 1) as the standard deviation of the logarithmically transformed NFn/NFn + 1 ratios, reflecting the effect of including an additional gene (Vandesompele et al. 2002). Once the PV value for n genes is below a cutoff of 0.15, additional genes are considered not to improve normalization (Mallona et al. 2010). The results in Fig. 5 also suggested that when normalizing gene expression, at least two RGs should be used except in the case of dehydration treatment, in 536

Short-duration salt

Long-duration salt

Gene

Gene

Stab

GII𝛼 0.035 HIS 0.042 RPL17 0.047 GAPDH 0.054 RP 0.065 18S 0.090 Actin 0.094 LIPL 0.109 Cpn60𝛽 0.125 EF1𝛼 0.133 eIF-5A 0.138 60S 0.150 UBQ 0.165 TUB 0.209 RA 0.343 LTP 0.432 GII𝛼 and HIS 0.027

Stab

18S 0.037 RP 0.040 TUB 0.061 60S 0.075 UBQ 0.088 HIS 0.090 eIF-5A 0.124 EF1𝛼 0.226 GAPDH 0.407 Cpn60𝛽 0.437 RPL17 0.468 RA 0.477 GII𝛼 0.523 Actin 0.542 LTP 0.601 LIPL 0.698 18S and RP 0.028

Total Gene

Stab

EF1𝛼 0.251 60S 0.284 HIS 0.305 Cpn60𝛽 0.307 UBQ 0.312 RP 0.318 eIF-5A 0.342 GII𝛼 0.363 Actin 0.451 LIPL 0.458 18S 0.459 RPL17 0.560 LTP 0.631 RA 0.649 GAPDH 0.720 TUB 0.817 60S and EF1𝛼 0.195

which three was the minimum number. The PV values showed the same trend as seen for stability measurements: the treatments with the lowest average PV were cold, SS. In contrast, gene expression in the dehydration stress, including in total, showed greater variability with higher PV values (table in Fig. 5). Consensus list and validation of the stability of P. euphratica RGs In determining the suitability of RGs for normalization of target gene expression, the three different software programs gave significantly different results in ranking

Physiol. Plant. 152, 2014

Fig. 4. Average expression stability values (M) and ranking of the candidate RGs calculated using geNorm in all 36 cDNA samples. A lower value of the average expression stability indicates more stable expression.

patterns. Because we arranged the internal genes into four lists according to the rank positions generated by each of the four statistical approaches: geNorm M values, NormFinder stability values, BestKeeper CV and r2 values, the discrepancies are probably related to differences in the algorithm each program uses in gene ranking. To provide a consensus result, RankAggreg with non-weighted unsupervised rank aggregation method was used to merge the four outputs with the aim of obtaining an optimal list of genes for each abiotic stress treatment. The results of the merged data revealed that the two best RGs for normalization were RPL17 and RP for ABA treatment, EF1𝛼 and 60S for cold treatment, HIS

Physiol. Plant. 152, 2014

and 60S for dehydration treatment, GII𝛼 and Actin for drought treatment, GII𝛼 and RPL17 for SS, and TUB and RP for LS. When considering all the treatments, 60S and UBQ were the best choices (Fig. 6A–G). On the basis of the number of RGs to use suggested by geNorm and the ranking list suggested by RankAggreg, the best combination of RGs in each treatment was shown in Table 4. Expression analysis of PePYL1, PeSCOF-1 and PeSCL7 In order to confirm the validity of the procedure for selection of the control genes detailed above, the RQs of three 537

Table 3. Rankings of candidate RGs in order of their expression stability as calculated by BestKeeper. Abscisic acid Rank position 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Cold

r2

CV

Dehydration r2

CV

r2

CV

Drought CV

r2

Short-duration salt CV

60S RPL17 RP LTP TUB HIS EF1𝛼 GII𝛼 HIS RP GAPDH EF1𝛼 RPL17 18S GAPDH Actin UBQ RP HIS Cpn60𝛽 eIF-5A 60S LTP RA TUB eIF-5A GII𝛼 UBQ RA RA GII𝛼 RP LIPL 18S GAPDH EF1𝛼 EF1𝛼 LIPL GII𝛼 Cpn60𝛽 UBQ RP GII𝛼 Actin LIPL eIF-5A LTP 60S 18S Actin 60S 60S 60S GAPDH TUB GII𝛼 UBQ Actin HIS RPL17 RP RPL17 RPL17 18S Actin HIS HIS GII𝛼 EF1𝛼 UBQ RA Actin GII𝛼 60S GAPDH EF1𝛼 eIF-5A GII𝛼 Cpn60𝛽 HIS Cpn60𝛽 Actin RP Actin eIF-5A 60S eIF-5A HIS Cpn60𝛽 eIF-5A RPL17 HIS TUB RP EF1𝛼 Cpn60𝛽 RA RP 60S Cpn60𝛽 18S LIPL LIPL LIPL UBQ eIF-5A EF1𝛼 TUB LIPL UBQ 18S TUB GAPDH Actin GAPDH TUB UBQ LTP TUB Cpn60𝛽 GAPDH RPL17 LTP LTP LTP 18S GAPDH eIF-5A LTP UBQ Cpn60𝛽 TUB RPL17 LIPL RA RA EF1𝛼 RPL17 RA RA 18S LIPL 18S LTP

Long-duration salt

r2

CV

r2

GII𝛼 RPL17 GAPDH LTP HIS 18S Actin RA RP Cpn60𝛽 60S UBQ TUB LIPL eIF-5A EF1𝛼

TUB RP HIS 60S eIF-5A 18S UBQ EF1𝛼 GAPDH Cpn60𝛽 LTP GII𝛼 RPL17 Actin RA LIPL

RP LTP UBQ RA 60S RPL17 TUB HIS 18S eIF-5A LIPL GAPDH GII𝛼 Actin EF1𝛼 Cpn60𝛽

Total CV

r2

18S GII𝛼 UBQ Cpn60𝛽 RP UBQ LTP RP LIPL 60S 60S 18S EF1𝛼 LIPL HIS EF1𝛼 GII𝛼 HIS RPL17 eIF-5A eIF-5A Actin Cpn60𝛽 TUB Actin RPL17 TUB GAPDH GAPDH RA RA LTP

Pairwise variation (PV)

0.25

0.20

*

0.15

*

* *

0.10

* * 0.05

*

0.00 Abscisic acid Cold Dehydration Drought Short-duration salt Long-duration salt Total

V2/3 0.120 0.082 0.164 0.110 0.041 0.078 0.209

V3/4 0.098 0.100 0.134 0.107 0.081 0.075 0.194

V4/5 0.099 0.079 0.116 0.078 0.069 0.098 0.151

V5/6 0.113 0.063 0.114 0.073 0.055 0.074 0.145

V6/7 0.102 0.061 0.088 0.063 0.046 0.074 0.168

V7/8 0.093 0.051 0.090 0.058 0.045 0.096 0.159

V8/9 0.162 0.052 0.085 0.056 0.046 0.152 0.170

V9/10 V10/11 V11/12 V12/13 V13/14 0.153 0.127 0.116 0.107 0.109 0.057 0.063 0.073 0.091 0.086 0.091 0.080 0.094 0.113 0.101 0.053 0.050 0.058 0.053 0.066 0.047 0.047 0.045 0.043 0.044 0.131 0.133 0.116 0.164 0.143 0.148 0.137 0.152 0.147 0.151

V14/15 V15/16 0.101 0.098 0.125 0.116 0.093 0.144 0.087 0.089 0.083 0.096 0.150 0.155 0.140 0.157

Fig. 5. PV to determine the optimal number of control genes for accurate normalization. The pairwise variation (Vn/Vn + 1) was analyzed between the NFs NFn and NFn+1 using geNorm software, where n is the number of genes involved in the NF. Asterisk and red figure indicate the optimal number of genes for normalization.

P. euphratica genes were examined. These included the already identified gene PeSCL7, a positive control used in P. euphratica gene abiotic stress-regulation studies; a START domain protein known as the ABA receptor PePYL1/PePYR1/PeRCAR12, and the cold-inducible zinc finger protein PeSCOF-1. Both of them have similar functions in P. euphratica as that in Arabidopsis. This was also confirmed by semi-quantitative RT-PCR (Fig. S3). To confirm whether the use of different RGs [low-ranked reference gene (LR); or top-ranked reference gene (TR)] affected the expression profiles of the three target genes, we compared the expression levels using 538

either just one LR, just one TR or a combination of multiple TRs (Fig. 7). The results indicated that PePYL1 was induced by ABA stress and PeSCOF-1 was induced by cold stress, and both of their RQs increased at first and then decreased. However, the entire or maximum values were quite different when normalization was performed using different RGs. For instance, for the gene PePYL1 in ABA treatment, when LR LTP was used for normalization, the maximum RQ was 6.43 at 12 h (Fig. 7A), but when the TR RPL17 was used, the maximum RQ was 11.61 at 3 h (Fig. 7B, green bar), which was closer to the maximum RQ of 9.88 at 3 h calculated with a combination

Physiol. Plant. 152, 2014

A

C

α 6 eI 0S F5A R P G IIα R A A ct in G HI A S P D H 18 S TU B C LIP pn L 6 R 0β P L1 7 LT P U B Q

E

G I A Iα ct in G R A A P D H U B Q 60 S C R pn P 60 β eI HIS F5 E A F1 T α R UB P L1 7 LT P 18 S LI P L

R A TU B 18 S

H IS 60 A S ct i U n B Q G eI IIα F5 E A F1 α G LIP A L P D H R RP P C L1 pn 7 60 LT β P

D

E

R P 18 S 6 eI 0S F5A H I U S B E Q G F1 A α P C DH pn 6 R 0β P L1 7 R A G II A α ct i LT n P LI P L

B

G R IIα G PL1 A 7 P D H H IS R P 18 A S ct L in C IP pn L 60 E β F1 α 60 S U B T Q eI UB F5A R A LT P

F

TU

R

P

L1

F1

7 R P H IS U B E Q F1 α 18 S G IIα C 60 pn S 60 L β eI IPL F5A G TU A B P D H R A A ct i LT n P

B

60 U S B Q H E IS F1 α R P C 18 pn S 60 G β II L α eI IPL F5 A A R ctin P L1 7 LT P T G U A B P D H R A

G

Fig. 6. Rank aggregation of gene lists using the Monte Carlo algorithm. Visual representation of rank aggregation using Monte Carlo algorithm with the Spearman footrule distances. Panels A–G represent different treatments. (A) ABA, (B) cold, (C) dehydration, (D) drought, (E) SS, (F) long-duration salt and (G) total. The solution of the rank aggregation is shown in a plot in which genes are ordered based on their rank position according to each stability measurement (gray lines). The mean rank position of each gene is shown in black, as well as the model computed by the Monte Carlo algorithm (red line).

of the TRs RPL17 and RP (Fig. 7B, purple bar). Under SS treatment when using a TR, PeSCL7 expression appeared to be induced quickly with the level of the product peaking at 1 h (Fig. 7J), which was similar to the results of semi-quantitative RT-PCR and previous study (Ma et al. 2010). Meanwhile, if a LR was used, the RQs rose insignificantly. In cold and drought treatments, the maximum RQs of the three different normalization patterns also varied (Fig. 7C, D, K, L). Taking this one step further, in dehydration and LS, the differences in RQs were extremely large. Normalized using a LR, the PeSCL7 expression levels were opposite and not upregulated but downregulated by the stresses, especially in dehydration (Fig. 7E). For PeSCL7 under these two stresses, the use of the different RGs (LR or TR) fundamentally changed the confirmed PeSCL7 expression profiles. Therefore, our results underline the fact that the use of particular RGs can influence the expression value of the gene of interest, and the RGs defined as stable (TR) might be suitable for quantitative gene expression normalization.

Physiol. Plant. 152, 2014

Table 4. Best combination of RGs based on the geNorm and RankAggreg programs. Experimental treatments ABA

Cold

Dehydration

Drought

SS

LS

Total

RPL17 RP

EF1𝛼 60S

HIS 60S Actin

GII𝛼 Actin

GII𝛼 RPL17

TUB RP

60S UBQ HIS EF1𝛼 RP

Discussion Identification of normalization genes for P. euphratica qRT-PCR is broadly used for the accurate and sensitive quantification of gene transcript levels, even for those genes whose transcript levels are low. To obtain valid qRT-PCR results, it is necessary to choose an appropriate internal control. The optimal internal control 539

gene should have similar expression levels in different cell types, developmental stages or sample treatments in experimental conditions. However, few of these endogenous controls show stable expression in all conditions (Pfaffl et al. 2004, Huis et al. 2010). Therefore, it is important to validate the expression stability of a control gene under specific experimental conditions prior to use in qRT-PCR normalization. Generally, most gene expression studies use a single internal control for normalization (Suzuki et al. 2000, Deng et al. 2013, Yang et al. 2013, Zachová et al. 2013), and the validity of the conclusions depends highly on the control applied. However, in order to avoid erroneous data and misinterpretation of experimental results that may be triggered by the use of a single RG, normalization with multiple RGs is becoming the standard (Vandesompele et al. 2002, Artico et al. 2010, Liu et al. 2013). The identifications of RGs on woody plants are limited to species as pine (Gonçalves et al. 2005), longan tree (Lin and Lai 2010), conifers (de Vega-Bartol et al. 2011), Platycladus (Chang et al. 2012), Eucalyptus (de Oliveira et al. 2012) and two species of Populus (Brunner et al. 2004, Xu et al. 2011). One of the most important woody species for fixing carbon, Populus, grows in abundance in many countries. The earliest RG selection for Populus was for the poplar P. trichocarpa × Populus deltoides (cottonwood hybrid). UBQ was the most stable gene among 10 potential control genes in a variety of poplar tissues when detecting changes associated with seasonal development and tree aging (Brunner et al. 2004). Another Populus RG study was in hybrid poplar (Populus × euramericana cv. Nanlin895), which used nine genes to examine macroscopic developmental time in the adventitious rooting of hardwood cuttings. EF1𝛼 and 18S were regarded as the most appropriate RGs for qRT-PCR analysis in this complex developmental process (Xu et al. 2011). However, no RG selection for stress resistance in woody plants has been reported, and there is also very little in stress resistant herbaceous plants. As cold, salinity, drought or other abiotic stresses adversely affect plant growth and productivity, it is important to develop stress-tolerant plants (Mahajan and Tuteja 2005). This study was conducted in another important Populus group, P. euphratica, which is considered a model system for woody plants in abiotic stress resistance research. The expression consistency of the 16 genes of this species when used as internal controls was validated among six abiotic stress conditions. Influence of stress on gene expression patterns Stress signals are perceived through diverse sensors, including many second messengers, phytohormones 540

(such as ABA), signal transducers (such as protein kinases and phosphatases) and transcription factors, resulting in the activation of a large number of stress-related genes and the synthesis of diverse functional proteins (Zhu 2002, Tang et al. 2012). These processes have the potential to greatly affect the transcriptional levels of all 16 selected RGs. The six stresses were found to have minimum impact on 18S expression and a large effect on TUB and GAPDH expressions; these effects may be associated with their overall levels, with 18S being the highest and TUB the lowest (Fig. 3). Evaluation of the six stresses individually revealed that dehydration had the largest impacts, followed by ABA and LS (Table 2, Fig. 4, stability and M values are large). The reason might be that dehydration leads to severe water deficit and seriously affect cellular activities, while water provides the medium within which most cellular functions take place. Data mining strategies and consensus list of genes for normalization As a single software package may introduce bias, we employed three tools in our analysis. As discussed by other authors, NormFinder and BestKeeper examine primarily CT values, whereas geNorm evaluates relative quantification, a consequence of which is that dissimilarities in PCR efficiency can affect stability measurements (Jian et al. 2008). In addition, the programs are intrinsically biased in some of their algorithms because they assume that the data are normally distributed. For example, BestKeeper is based on Pearson correlation analysis, which requires normally distributed and variance homogeneous data (Pfaffl et al. 2004). We opted to use the two different evaluation indices CV and r2 calculated by BestKeeper, and our analysis merged different statistics, some of which are CT-based and others RQ-based, with the aim of counteracting this biasing influence. To summarize the results of our dataset analysis, RankAggreg was used to merge the data. The comprehensive results recommended the use of a combination of RPL17 and RP for ABA treatment, and a combination of HIS, 60S and Actin for dehydration (Table 4, Fig. 6A, C, D). Except cold and SS treatments UBQ ranks well for use with the other four treatments. One of the most commonly used RGs, UBQ, was selected as the best RG, as determined previously for A. thaliana (Czechowski et al. 2005), rice (Narsai et al. 2010), Brachiaria brizantha (Silveira et al. 2009) and Gossypium hirsutum (Artico et al. 2010). Another commonly used RG is EF1𝛼, which has been confirmed as one of the best for rice (Jain et al. 2006), soybean (Jian et al. 2008), B. brizantha (Silveira et al.

Physiol. Plant. 152, 2014

Physiol. Plant. 152, 2014

541

18S LTP

E

I

green purple green purple yellow green purple

RPL17 RPL17 and RP HIS HIS and 60S HIS, 60S and Actin GIIα GIIα and RPL17

J

F

J

I

B

F

E

UBQ LIPL LIPL

C G K

L

H

D

GIIα and Actin TUB TUB and RP

green purple

EF1 EF1 and 60S GII purple

green purple green

K

G

C

L

H

D

Fig. 7. Relative quantifications of PePYL1 (A and B), PeSCOF-1 (C and D) and PeSCL7 (E–L) in the different abiotic stress treatments with different RGs after normalization. Samples A, C, E, G, I and K were normalized with a single low-ranked RG (LR) while samples B, D, F, H, J and L were normalized with a single top-ranked RG (TR) and a combination of multiple top-ranked RGs (TRs). The RGs in use were displayed in the table. The y-axis is consistent in the same treatment. Values shown are relative expression levels ±SE (n = 3 experiments).

LTP

A

B

A

2009) and Petunia hybrida (Mallona et al. 2010). In this study, EF1𝛼 was one of the best choices for use with ABA, cold and dehydration treatments (Fig. 7A–C). This may suggest that the RGs validated in P. euphratica basically correspond to the above-mentioned species. Simultaneously, some novel RGs appropriate for qRT-PCR normalization were also identified in P. euphratica, for instance RP or RPL17 in ABA, cold, SS and LS treatments (Fig. 7A, B, E, F). In particular, 60S varied the least (Fig. 7G). Taken together, these results suggest that RGs exhibit similar expression patterns in different plant species, because they are implicated in basic cell metabolism and cellular functions. However, in specific species some RGs are regulated differently, and each species may have its own stably expressed genes that can be used for qRT-PCR normalization. P. euphratica gene expression For P. euphratica stress studies, we recommend that at least two RGs be used for reliable quantification, except for dehydration treatment, which requires three (Fig. 5). Because dehydration had the greatest impact on the stability of all candidate RGs, in order to ensure accurate quantification, it is better to use more RGs. The first value smaller than 0.15 is V3/V4 (V3/V4 = 0.134 < 0.15, n = 3; Fig. 5, Table 4, and yellow bar in Fig. 7F). Our results revealed that the choice of RGs modified the determined expression profiles of the target gene, PeSCL7. Although the expression levels normalized by a single top-ranked RG or a combination of multiple top-ranked RGs were not significantly different, the use of a single low-ranked RG had a great impact on the results. Our results showed important differences in the PeSCL7 expression levels obtained using different normalization methods, and suggest that care should be taken when choosing endogenous controls to interpret gene expression. Different normalization procedures are needed to verify the data. In addition, in this study, the PeSCL7 expression levels normalized by top-ranked RGs were close to those of the previous study, probably meaning that the top-ranked RGs identified in this study are appropriate for use in quantitative gene expression studies in P. euphratica.

Conclusion In summary, each program with different algorithm identified different genes as those best suitable for performing normalization. RankAggreg (Pihur and Datta 2009), an R package for non-weighted rank aggregation, was employed to merge the data in an unbiased way and give identical weight to the output of the different programs. 542

Our composite results showed that RPL17 was the best RG for use in ABA treatment, EF1𝛼 was the best for cold treatments, HIS was the best for dehydration treatment, GII𝛼 was the best RG in drought and SS treatments, and TUB was the best for LS treatment. 60S varied the least during the different treatments. In each treatment, the lower-ranked genes should be avoided as RGs, whereas the recommended ones may be suitable for gene normalization. In addition, we recommend that at least two RGs should be used for reliable quantification, except in dehydration treatment, which requires three (Fig. 5). The applicability of the RGs already validated as suitable ones was confirmed by comparing the RQs of PePYL1, PeSCOF-1 and PeSCL7 using only one RG or a combination of multiple RGs (Fig. 7). The results demonstrated that the use of unsuitable RGs may result in biased expression levels. Therefore, we recommend that a RG stability test should be done before performing gene expression studies in P. euphratica. Acknowledgements – We would like to thank Jing Hu (University of Alabama, Birmingham, USA) for providing us the QBASE software and some practical guidances; and Sha Tang, Xian Wu, Peng Shuai, Tao Pang, Xiao Han, Wei Han and Chao Yuan for their help and insightful comments on the experiments. This work was supported by the Hi-Tech Research and Development Program of China (2013AA102701), the National Natural Science Foundation of China (31270656), Joint Programs of the Scientific Research and Graduate Training from BMEC (Stress Resistance Mechanism of Poplar) and 111 Project of Beijing Forestry University (B13007).

References de Almeida MR, Ruedell CM, Ricachenevsky FK, Sperotto RA, Pasquali G, Fett-Neto AG (2010) Reference gene selection for quantitative reverse transcriptionpolymerase chain reaction normalization during in vitro adventitious rooting in Eucalyptus globulus Labill. BMC Mol Biol 11: 73 Andersen CL, Jensen JL, Orntoft TF (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 Res 64: 5245–5250 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 Biol 10: 49 Bogeat-Triboulot MB, Brosché M, Renaut J, Jouve L, Le Thiec D, Fayyaz P, Vinocur B, Witters E, Laukens K,

Physiol. Plant. 152, 2014

Teichmann T, Altman A, Hausman JF, Polle A, Kangasjärvi J, Dreyer E (2007) Gradual soil water depletion results in reversible changes of gene expression, protein profiles, ecophysiology, and growth performance in Populus euphratica, a poplar growing in arid regions. Plant Physiol 143: 876–892 Brentner LB, Mukherji ST, Merchie KM, Yoon JM, Schnoor JL, Aken BV (2008) Expression of glutathione S-transferases in poplar trees (Populus trichocarpa) exposed to 2,4,6-trinitrotoluene (TNT). Chemosphere 73: 657–662 Brinker M, Brosché M, Vinocur B, Abo-Ogiala A, Fayyaz P, Janz D, Ottow EA, Cullmann AD, Saborowski J, Kangasjärvi J, Altman A, Polle A (2010) Linking the salt transcriptome with physiological responses of a salt-resistant Populus species as a strategy to identify genes important for stress acclimation. Plant Physiol 154: 1697–1709 Brunner AM, Yakovlev IA, Strauss SH (2004) Validating internal controls for quantitative plant gene expression studies. BMC Plant Biol 4: 14 Bustin SA (2002) Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems. J Mol Endocrinol 29: 23–39 Chang S, Puryear J, Cairney J (1993) A simple and efficient method for isolating RNA from pine trees. Plant Mol Biol Rep 11: 113–116 Chang E, Shi S, Liu J, Cheng T, Xue L, Yang X, Yang W, Lan Q, Jiang Z (2012) Selection of reference genes for quantitative gene expression studies in Platycladus orientalis (Cupressaceae) using real-time PCR. PLoS One 7: e33278 Chen J, Xia X, Yin W (2009) Expression profiling and functional characterization of a DREB2-type gene from Populus euphratica. Biochem Biophys Res Commun 378: 483–487 Coll A, Nadal A, Collado R, Capellades G, Kubista M, Messeguer J, Pla M (2010) Natural variation explains most transcriptomic changes among maize plants of MON810 and comparable non-GM varieties subjected to two N-fertilization farming practices. Plant Mol Biol 73: 349–362 Czechowski T, Stitt M, Altmann T, Udvardi MK, Scheible WR (2005) Genome-wide identification and testing of superior reference genes for transcript normalization in Arabidopsis. Plant Physiol 139: 5–17 Deng X, Zhou S, Hu W, Feng J, Zhang F, Chen L, Huang C, Luo Q, He Y, Yang G, He G (2013) Ectopic expression of wheat TaCIPK14, encoding a calcineurin B-like protein-interacting protein kinase, confers salinity and cold tolerance in tobacco. Physiol Plant 149: 367–377 Fujii H, Chinnusamy V, Rodrigues A, Rubio S, Antoni R, Park SY, Cutler SR, Sheen J, Rodriguez PL, Zhu JK (2009) In vitro reconstitution of an abscisic acid signalling pathway. Nature 462: 660–664

Physiol. Plant. 152, 2014

Gachon C, Mingam A, Charrier B (2004) Real-time PCR: what relevance to plant studies? J Exp Bot 55: 1445–1454 Gonçalves S, Cairney J, Maroco J, Oliveira MM, Miguel C (2005) Evaluation of control transcripts in real-time RT-PCR expression analysis during maritime pine embryogenesis. Planta 222: 556–563 Gries D, Zeng F, Foetzki A, Arndt SK, Bruelheide H, Thomas FM, Zhang X, Runge M (2003) Growth and water relations of Tamarix ramosissima and Populus euphratica on Taklamakan desert dunes in relation to depth to a permanent water table. Plant Cell Environ 26: 725–736 Gu R, Fonseca S, Puskás LG, Hackler L, Zvara Á, Dudits D, Pais MS (2004) Transcript identification and profiling during salt stress and recovery of Populus euphratica. Tree Physiol 24: 265–276 Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J (2007) qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol 8: R19 Hong SY, Seo PJ, Yang MS, Xiang F, Park CM (2008) Exploring valid reference genes for gene expression studies in Brachypodium distachyon by real-time PCR. BMC Plant Biol 8: 112 Hsiao TC (1973) Plant responses to water stress. Annu Rev Plant Physiol Plant Mol Biol 24: 519–570 Huis R, Hawkins S, Neutelings G (2010) Selection of reference genes for quantitative gene expression normalization in flax (Linum usitatissimum L.). BMC Plant Biol 10: 71 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. Biochem Biophys Res Commun 345: 646–651 Jian B, Liu B, Bi Y, Hou W, Wu C, Han T (2008) Validation of internal control for gene expression study in soybean by quantitative real-time PCR. BMC Mol Biol 9: 59 Kim JC, Lee SH, Cheong YH, Yoo CM, Lee SI, Chun HJ, Yun DJ, Hong JC, Lee SY, Lim CO (2001) A novel cold-inducible zinc finger protein from soybean, SCOF-1, enhances cold tolerance in transgenic plants. Plant J 25: 247–259 Li B, Qin Y, Duan H, Yin W, Xia X (2011) Genome-wide characterization of new and drought stress responsive microRNAs in Populus euphratica. J Exp Bot 62: 3765–3779 Lin YL, Lai ZX (2010) Reference gene selection for qPCR analysis during somatic embryogenesis in longan tree. Plant Sci 178: 359–365 Liu M, Shi J, Lu C (2013) Identification of stress-responsive genes in Ammopiptanthus mongolicus using ESTs generated from cold- and drought-stressed seedlings. BMC Plant Biol 13: 88

543

Ma Y, Szostkiewicz I, Korte A, Moes D, Yang Y, Christmann A, Grill E (2009) Regulators of PP2C phosphatase activity function as abscisic acid sensors. Science 324: 1064–1068 Ma HS, Liang D, Shuai P, Xia XL, Yin WL (2010) The saltand drought-inducible poplar GRAS protein SCL7 confers salt and drought tolerance in Arabidopsis thaliana. J Exp Bot 61: 4011–4019 Mahajan S, Tuteja N (2005) Cold, salinity and drought stresses: an overview. Arch Biochem Biophys 444: 139–158 Mallona I, Lischewski S, Weiss J, Hause B, Egea-Cortines M (2010) Validation of reference genes for quantitative real-time PCR during leaf and flower development in Petunia hybrida. BMC Plant Biol 10: 4 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 Biol 10: 56 de Oliveira LA, Breton MC, Bastolla FM, Camargo SS, Margis R, Frazzon J, Pasquali G (2012) Reference genes for the normalization of gene expression in Eucalyptus species. Plant Cell Physiol 53: 405–422 Ottow EA, Polle A, Brosche M, Kangasjärvi J, Dibrow P, Zörb C, Teichmann T (2005) Molecular characterization of PeNhaD1: the first member of the NhaD Na+ /H+ antiporter family of plant origin. Plant Mol Biol 58: 75–88 Park SY, Fung P, Nishimura N, Jensen DR, Fujii H, Zhao Y, Lumba S, Santiago J, Rodrigues A, Chow TF, Alfred SE, Bonetta D, Finkelstein R, Provart NJ, Desveaux D, Rodriguez PL, McCourt P, Zhu JK, Schroeder JI, Volkman BF, Cutler SR (2009) Abscisic acid inhibits type 2C protein phosphatases via the PYR/PYL family of START proteins. Science 324: 1068–1071 Pfaffl MW (2001) A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Res 29: e45 Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP (2004) Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper – Excel-based tool using pair-wise correlations. Biotechnol Lett 26: 509–515 Pihur V, Datta S (2009) RankAggreg, an R package for weighted rank aggregation. BMC Bioinformatics 10: 62 Qiu Q, Ma T, Hu Q, Liu B, Wu Y, Zhou H, Wang Q, Wang J, Liu J (2011) Genome-scale transcriptome analysis of the desert poplar, Populus euphratica. Tree Physiol 31: 452–461 Rozen S, Skaletsky H (2000) Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol 132: 365–386 Silveira ED, Alves-Ferreira M, Guimarães LA, da Silva FR, Carneiro VT (2009) Selection of reference genes for quantitative real-time PCR expression studies in the

544

apomictic and sexual grass Brachiaria brizantha. BMC Plant Biol 9: 84 Suzuki T, Higgins P, Crawford D (2000) Control selection for RNA quantitation. Biotechniques 29: 332–337 Tang N, Zhang H, Li X, Xiao J, Xiong L (2012) Constitutive activation of transcription factor OsbZIP46 improves drought tolerance in rice. Plant Physiol 158: 1755–1768 Tang S, Liang H, Yan D, Zhao Y, Han X, Carlson JE, Xia X, Yin W (2013) Populus euphratica: the transcriptomic response to drought stress. Plant Mol Biol 83: 539–557 Taylor S, Wakem M, Dijkman G, Alsarraj M, Nguyen M (2010) A practical approach to RT-qPCR – publishing data that conform to the MIQE guidelines. Methods 50: S1–S5 Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F (2002) Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol 3: 1–12 de Vega-Bartol JJ, Santos RR, Simões M, Miguel C (2011) Evaluation of reference genes for quantitative PCR analysis during somatic embryogenesis in conifers. BMC Proc 5(suppl 7): O44 Wan H, Zhao Z, Qian C, Sui Y, Malik AA, Chen J (2010) Selection of appropriate reference genes for gene expression studies by quantitative real-time polymerase chain reaction in cucumber. Anal Biochem 399: 257–261 Wei L, Miao H, Zhao R, Han X, Zhang T, Zhang H (2013) Identification and testing of reference genes for Sesame gene expression analysis by quantitative real-time PCR. Planta 237: 873–889 Xu M, Zhang B, Su X, Zhang S, Huang M (2011) Reference gene selection for quantitative real-time polymerase chain reaction in Populus. Anal Biochem 408: 337–339 Yang T, Peng H, Whitaker BD, Jurick WM (2013) Differential expression of calcium/calmodulin-regulated SlSRs in response to abiotic and biotic stresses in tomato fruit. Physiol Plant 148: 445–455 Zachová D, Fojtová M, Dvoˇráˇcková M, Mozgová I, Lermontova I, Peška V, Schubert I, Fajkus J, S´ykorová E (2013) Structure–function relationships during transgenic telomerase expression in Arabidopsis. Physiol Plant 149: 114–126 Zhu JK (2002) Salt and drought stress signal transduction in plants. Annu Rev Plant Biol 53: 247–273

Supporting Information Additional Supporting Information may be found in the online version of this article: Table S1. The standard deviations for the CT values under each treatment (Fig. 3A).

Physiol. Plant. 152, 2014

Table S2. CV, SD values and coefficient of determination (r2 ) values calculated using BestKeeper.

Fig. S3. Semi-quantitative RT-PCR data of the three stress inducible genes.

Fig. S1. The detection result of the integrities of 12 RNA samples using microfluidic devices.

Appendix S1. The R script file example of RankAggreg.

Fig. S2. Melting curves of the 16 RGs for qRT-PCR amplification specificity detection.

Edited by V. Hurry

Physiol. Plant. 152, 2014

545