Identification and evaluation of new reference genes in Gossypium hirsutumfor accurate normalization of real-time quantitative RT-PCR data

Gene expression analysis is increasingly important in many fields of biological research. Understanding patterns of expressed genes is crucial to provide insights into complex regulatory networks and will lead to the identification of genes relevant to new biological processes [1].

Reverse transcription real-time quantitative polymerase chain reaction (qPCR) is a robust method to study gene expression changes [2]. The main advantages of qPCR when compared to other experimental techniques used to evaluate gene expression levels, such as Northern blot hybridization and reverse transcription-polymerase chain reaction (RT-PCR), are its higher sensitivity, specificity, and broad quantification range of up to seven orders of magnitude [3]. Therefore, qPCR analysis has become the most common method for validating the whole-genome microarray data or a smaller set of genes and molecular diagnostics [4]. Although being extremely powerful technique, qPCR suffers from certain pitfalls, noteworthy the use of unreliable reference genes for the normalization step [5]. Normalization is necessary for the correction of non-specific variations, such as inaccurate quantification of RNA and problems in the quality of RNA that can trigger variable reverse transcription and PCR reactions. A number of strategies have been proposed to normalize qPCR data but normalization remains one of the most important challenges concerning this technique [5].

The expression of reference genes used for normalization in qPCR analysis should remain constant between the cells of different tissues and under different experimental conditions; otherwise, it can lead to erroneous results. Recent reports have demonstrated that some of the most well-known and frequently used reference genes are inappropriate for normalization in qPCR analysis due to expression variability [6–8]. The importance of reference genes for plant qPCR analysis has been recently emphasized even though the identification of these genes is quite laborious [9, 10]. Microarray datasets can also be a rich source of information for selecting qPCR reference genes [6], but unfortunately, this tool is still not available for most of plant species, including cotton.

The classical housekeeping genes involved in basic cellular processes such as 18 S rRNA, ubiquitin, actin, β-tubulin, and glyceraldehyde-3-phosphate dehydrogenase have been recurrently used as internal controls for gene expression analysis in plant as they are supposed to have a uniform expression all samples and experimental conditions tested. However, several reports demonstrated that the transcript levels of these genes also vary considerably under different experimental conditions and are consequently unsuitable for gene expression studies [6, 11]. Statistical algorithms such as geNorm [1], NormFinder [12] and BestKeeper [13] have been developed for the evaluation of best suited reference gene(s) for normalization of qPCR data in a given set of biological samples. Recognizing the importance of reference genes in normalization of RT-qPCR data, various housekeeping genes have been evaluated for stable expression under specific conditions in various organisms. Many works have been carried on animal and human health [3, 14] field that describe the identification of multiple reference genes for normalisation of qPCR data, but similar reports are scarce in plant research [4, 15, 16]. Czechowski et al. (2005) employed a new strategy for the identification of reference genes in Arabidopsis thaliana. Based on the microarray data of Affymetrix ATH1, several new reference genes were revealed in Arabidopsis [6]. Some of these genes have no previous information about function in Arabidopsis or any other organism. The list of new Arabidopsis reference genes revealed by Czechowski and collaborators was successfully employed to search reference genes in unrelated species such as Vitis vinifera by sequence homology [9]. Recently, our group was also successful in providing new reference genes for qPCR in Coffea arabica and Brachiaria brizantia using the same strategy employed in V. vinifera [17, 18].

Cotton (Gossypium spp.) is the world's most important source of natural fiber and also an important source of edible oil [19]. Because of its unique reproductive developmental aspects and speciation history, G. hirsutum has attracted considerable scientific interest, not only among plant breeders and agricultural scientists, but also among taxonomists, developmental geneticists, and evolutionary biologists [20–24]. In spite of this, qPCR analyses in cotton are still hampered by the use of inappropriate references genes.

In this study, we report the validation of housekeeping genes to identify the most suitable internal control gene(s) for normalization of qPCR data obtained in different plant organs and floral verticils and also during flower and fruit development. In addition, to illustrate the usefulness of the new reference genes, we provided a detailed expression analysis of two MADS-box transcription factors in cotton, putative homologues of Arabidopsis AGAMOUS and SEPALLATA3 genes.

In order to compare the expression levels of target genes in different tissues at the same time, it is crucial to normalize all the samples by the same set of reference genes. For the evaluation of potential reference, a well known housekeeping gene, poly-ubiquitin (GhUBQ14), was included in the qPCR experiments [26]. We selected eight new candidates to housekeeping genes (GhACT4, GhEF1α5, GhFBX6, GhPP2A1, GhMZA, GhPTB, GhGAPC2, GhβTUB3) in G. hirsutum. These genes are putative homologues of eight Arabidopsis genes included in the list of 27 best reference genes for qPCR analysis (Table 1) [6]. For the selection of the putative cotton housekeeping genes, we searched in the Cotton DB for homologues to the Arabidopsis referenced genes, only eight candidates that showed very high similarities (E-value > 1e-75) were included in the final list. The eight genes found in the cotton databanks belong to different functional classes based on Arabidopsis sequence information, which reduce the chances of co-regulated expression occurrence among these genes (Table 1). The gene name, accession number, A. thaliana homologue locus, A. thaliana annotation, similarity end identity, gene size, and primer sequence, are provided in Table 1. The nine cotton candidate reference genes were evaluated for gene expression stability by qPCR in a set of 23 cotton samples grouped into five different experimental sets. The first experimental set was composed of plant organs: leaves, stem, branch, root, flower buds (RNA pools of stages 2 to 12 mm) and fruits (RNA pools of stages 10 to 15 to fruits larger than 30 mm). The second set included floral meristem and size selected flower buds, based on their diameter of 2, 4, 6, 7, 8, 10 and 12 mm. The third experimental set was composed of the floral verticils: sepal, petal, stamen, carpel and pedicel. The fourth experimental set consists of four stages of fruit development based on it diameter: 10 to 15 (1), 16 to 20 (2), 21 to 30 (3) and larger than 30 mm (4). Finally, in the fifth set, we included all the tissues samples used in this study (23 distinct biological samples).

Total RNA was isolated from different tissue samples and reverse transcribed. The RNA quality for all samples was checked by gel eletrophoresis analisys and spectrophotometric assays (data not shown). Within a biological replicate, for a tissue sample, the same cDNA pool was used for qPCR analysis of each of the nine genes using gene-specific primers. qPCRs were performed in triplicate for each of the 23 cDNA pools along with a no template control in parallel for each gene. The melting-curve analysis performed by the PCR machine after 40 cycles of amplification and agarose gel electrophoresis showed that all the 9 primer pairs amplified a single PCR product of desired size from various cDNA (results not shown). Primer efficiencies for all primer combinations were higher than 0.90 (90%) in all experimental sets. Although, two primers pairs presented efficiencies below 90% in four samples: GhGAPC2 in flower buds and floral and plant organs and GhPP2A1 in floral organs (Table 2). The mean Cq value (average of 6 values from the two biological replicates) in a tissue sample for each gene is shown in Table 2. Cq values were in the range of 15.30 and 29.17. GhACT4, GhUBQ14 and GhEF1α5 are the top three most expressed genes in all sets followed by GhMZA, GhβTUB3, GhPP2A1 and GhPTB. GhFBX6 and GhGAPC2 genes present the lowest expression levels in all samples.

We used geNorm v3.5 software, to analyze the expression stability of the tested genes in all samples, and ranked them accordingly to gene stability measure (M). The M value is obtained by the use of relative expression values for each cDNA sample as input for the geNorm algorithm based on the geometric averaging of multiple control genes and mean pairwise variation of a gene from all other control genes in a given set of samples. Therefore, genes with the lowest M values have the most stable expression. The results obtained with geNorm algorithm are presented in the Figure 1 and summarized in Table 3. The geNorm algorithm also determines the pairwise variation Vn/n + 1, which measures the effect of adding further reference genes on the normalisation factor (that is calculated as the geometric mean of the expression values of the selected reference genes). It is advisable to add additional reference genes to the normalisation factor until the added gene has no significant effect. Vandesompele et al. (2002) used 0.15 as a cut-off value, below which the inclusion of an additional reference gene is not required. Pairwise variation analysis (Figure 2) showed that the ideal number of reference genes may be different for distinct set of samples. For instance, for the normalization of the floral organ set, only two genes are necessary. On the other hand, five genes are required for the normalization of the plant organ set. When evaluating all the pairwise variation, the least stable housekeeping gene was GhGAPC2 followed by GhβTUB3 since they significantly increased the pairwise variation during the whole assay by increasing the V value as shown in Figure 2. However, Vandesompele and collaborators recommend the use of at least three reference genes whenever this result obtained in our analysis is observed [1].

To assess the validity of the procedure for the selection of control genes detailed above, the relative expression level of two cotton genes that belong to MADS-box family were inspected. After the search in Cotton db using the MADS-box consensus sequence, 18 ESTs were found with high similarity to MIKC MADS box family (E-value ≤ 7 × 10^-13) (Data not shown). The reduced number of cotton MIKC type genes is expected since the ESTs sequencing efforts in cotton are very limited when compared to other species such as Arabidopsis and rice. In spite of the low number of MADS-box genes, the phylogenetic analysis identified good candidates to homologous genes of Arabidopsis AGAMOUS (AG) and SEPALLATA3 (SEP3) (data not shown). The homologue of AG, was previously characterized by RT-PCR and named GhMADS3 [36]. RT-PCR analysis suggests that GhMADS3 expression is restricted to stamens and carpels. Ectopic expression in Nicotiana tabacum L. indicates that it is the cotton orthologous gene to AG [36]. The Arabidopsis thaliana SEP3 is expressed in the three inner whorls of organs throughout flower development, but there is no information of the putative homologue of cotton (GhSEP-like1), identified by our phylogenetic analysis [37]. The expression of GhMADS3 and GhSEP-like1 was estimated in different plant organs, during flower development and in the floral organs of 6 mm flower buds. The qPCR analysis employed the control genes recommended by NormFinder program for the normalization of gene expression. The analysis revealed that G. hirsutum GhMADS3 and GhSEP-like1 genes very similar expression profiles of AG and SEP3 genes from Arabidopsis (Figure 3). However, we also observed unexpected expression patterns: GhSEP-like1 is expressed in cotton fruits and the GhMADS3 in pedicels of 6 mm flower buds.

The qPCR is broadly accepted as the method of choice for accurate and sensitive quantification of gene transcript levels, even for those genes whose transcript levels are low. For valid qPCR analysis, accurate normalization of gene expression against an appropriate internal control is required. The ideal control gene should have similar expression regardless of experimental conditions, including different cell types, developmental stages, and/or sample treatment. However, no one gene has a stable expression under every experimental condition, as numerous studies reported that expression of housekeeping genes can also vary considerably with experimental conditions. Consequently, normalization of gene expression with a single reference gene can trigger erroneous data and, consequently, misinterpretation of experiment results. Therefore, it is necessary to validate the expression stability of a control gene under specific experimental conditions prior to its use in qPCR normalization.

Normalisation with multiple reference genes is becoming the golden standard, but reports that identify such genes in plant research are limited [3, 4, 17, 18, 38, 39], even though algorithms are available to test the expression stability of candidates [1, 12, 13] and a number of candidate reference genes for Arabidopsis have been proposed [6]. To obtain a solid basis for normalization of our gene expression data when studying the flower development in cotton, we evaluated the expression stability of nine candidate reference genes, including one traditional "housekeeping" gene in five different experimental sets. Candidate genes were selected according to the level of DNA sequence similarity to genes previously identified as reference genes in Arabidopsis and cotton. This strategy has been successful in finding good reference genes in other species such as grape [39] and it was already employed by our group in coffee and B. brizantha [17, 18]. Another strategy used to identify bona fide qPCR reference genes is to check housekeeping genes previously used in Northern and RT-PCR studies [40, 41]. However, it has be shown that the expression of traditional reference genes may vary enormously depending on the test condition [6]. In cotton, Tu and collaborators tested six putative constitutive genes (Histone3, UBQ7, Actin, Cyclophilin, Gbpolyubiquitin-1 and Gbpolyubiquitin-2), two of them (Gbpolyubiquitin-1 and Gbpolyubiquitin-2) from previously published data [42]. In contrast to the present work, roots, floral stages and verticils samples were not included in the final set of samples [41]. The reference genes evaluation was performed using exclusively geNorm and the value obtained for the pairwise variation with the best control genes was above the cut-off value of 0.15 suggested by Vandesompele et al. [1]. Moreover, the expression in the fiber developmental series of the all six putative reference genes varied greatly, hampering their use for qPCR [41].

We elected the NormFinder as the preferential method for the selection of the best references genes since it considers intra- and inter-group variations for the normalization factor (NF). However, geNorm was also important to compose the final set of references genes for the experimental conditions tested in this work. Our analysis has shown that each experimental condition tested demands a specific set of reference genes (Table 3 and 4). This result emphasizes the importance of reference genes validation for each experimental condition, especially when samples belong to very different groups, e.g. different organs.

When plant organs and all samples were tested, GhUBQ14 and GhPP2a1 were considered the most appropriate reference genes. GhUBQ14 and GhPP2a1 should avoid error transferences since NormFinder chose them as the best combination of genes. NormFinder chose GhACT4 and GhUBQ14 as the best combination of two genes in flower buds. Both programs ranked GhACT4 as the most stable gene, conferring higher robustness to the NF. Our analyses of different floral organs revealed that GhACT4 and GhFBX6 are the most appropriated genes for qPCR normalization, since they represent the best combination of genes considered by NormFinder to improve NF. GhFBX6 was ranked by both algorithms as the most stable gene in the floral organs set. Finally, fruit development GhMZA was considered as the most stable gene in both the NormFinder and geNorm programs, and NormFinder chose GhMZA and GhPTB as the best combination of genes.

The GhACT4, GhEF1α5, GhFBX6, GhPP2A1, GhMZA, GhPTB, GhGAPC2, GhβTUB genes were identified as novel reference genes in A. thaliana through microarray experiments and were validated by qPCR [7]. Among them, GhGAPC2 gave poor results in our analysis in cotton. GhUBQ14, a traditional reference gene in cotton [26] was well evaluated by NormFinder ranking in the best combination in three of the five experimental sets. Although, evaluations of a traditional reference genes by the same procedures used in this work not always give support to their frequent use. For instance, UBQ10 gene shows highly stable expression in Arabidopsis [6] whereas its putative homologue has been shown unsuitable for normalization of different tissues at different developmental stages in rice and soybean [4, 43].

Other commonly used housekeeping gene, GhβTUB, displayed inappropriate expression variability limiting its use as internal control in cotton. A similar result was also observed for the β-tubulin of B. brizantha when male and female reproductive tissues, spikelets, roots and leaves were evaluated [17]. On the other hand, β-TUB is one of most stably expressed genes in poplar (Populus ssp) tissue samples among the 10 reference genes tested [10]. GAPDH, another traditional reference gene, was considered the most appropriate reference gene when coffee leaves drought-stressed vs. control plants and different coffee cultivar leaves were analyzed [18]. Taken together, these results suggest that the housekeeping genes are regulated differently in different plant species and may exhibit differential expression patterns. This may partly be explained by the fact that housekeeping genes are not only implicated in the basal cell metabolism but also may participate in other cellular functions [11].

The programs employed to evaluate reference genes in our study (geNorm and NormFinder) use the same input data, i.e. non-normalized relative quantities, and Cqs need to be transformed considering primer pair efficiencies. In our experience, it is crucial to evaluate primer pair efficiencies for each sample tested since primer efficiency varies depend on the according to biological sample. The importance of this step can be well illustrated by the primer efficiency variation of GhGACP2 in flower buds compared to fruits (Table 2).

The values of Cq presented here should not be considered alone, but they may help in the selection of best combination of reference genes when there is previous data about target gene expression levels. Similar expression levels of the reference and target genes are considered an important issue regarding qPCR normalization [1]. Indeed, references genes with excessively high/low expression levels compared to target genes can trigger problems for data analysis [44, 45].

As suggested by Remans and collaborators [7], biological replicates were submitted to geNorm and NormFinder as independent samples. This procedure increased the credibility of the most suitable cotton reference genes because it takes into account possible variations in reference gene expression that are not due to different treatments, but intrinsic to the gene itself.

To illustrate the suitability of the reference genes revealed in the present study, two putative cotton homologues to AG and SEP3 (GhMADS3 and GhSEP-like1) had their expression profile evaluated in different plant organs, during flower development and in floral organs at flower buds of 6 mm (Figure 3). As it is observed to AG and SEP3, the GhMADS3 and GhSEP-like1 genes are highly expressed in flower buds, but GhSEP-like1 also shows a high expression in fruits. GhMADS3 also is expressed in higher levels after stage of 2 mm and throughout cotton flower development. The low expression of GhMADS3 in floral meristem is expected as well a high expression level in stamen and carpels of 6 mm flower bud. The AG gene is expressed in few cells during the initial flower development to establish organ identity and is also important at later stages of stamens and carpels development [46, 47]. The GhMADS3 expression observed in pedicels may be the result of contamination of material derived from carpels. These two organs are merged, which hamper a perfect separation during flower bud dissection to collect the samples. Our analysis of GhMADS3 expression by RT-qPCR refined the information of the previous study adding accuracy, spatial and temporal information to the expression during floral development [36]. In addition, it revealed that this GhMADS3 is also expressed in cotton fruits (Figure 3).

The high expression in fruits of GhSEP-like1 contrast to the spatial and developmental expression pattern of SEP3 in Arabidopsis, former AGL9 (Figure 3) [37]. SEP3 is expressed in all Arabidopsis flower verticils throughout development but no transcripts are found in siliques. However, PPERSEP3, a putative Prunus persica homologue of SEP3, is expressed during fruit development [48]. In addition, GhMADS4 and GhMADS7, genes from AGAMOUS subclass in cotton, are also expressed during fiber development [23]. Taking together, these results suggest that besides flower identity SEP3 and AG-like genes in cotton may be involved in an additional developmental process during fruit development.

This work constitutes the first in-depth study to validate the optimal control genes for the quantification of transcript levels in different cotton organs and during flower and fruit development. The use of the new cotton reference genes combined with size collected flower buds and floral organ dissection allowed a precise spatial and temporal characterization of two MADS-box genes in cotton plants. In summary, the new cotton reference genes will enable more accurate and reliable normalization of qPCR results for gene expression studies in this important crop plant.