Selection of reference genes for quantitative gene expression normalization in flax (Linum usitatissimum L.)

Flax is one of the earliest cultivated plants known to man and the remains of seeds found in different archeological sites suggest that it was also used before domestication, most probably for textiles [1]. Cultivated flax (Linum usitatissimum L.) is an annual, diploid, self-pollinated crop used as a source of both high-quality, cellulose-rich bast fibers and oil. The main areas of fiber production are currently found in France (77,000 ha) with an annual turnover of over 200 million euros followed by Russia and China (80,000 ha). In contrast, oil varieties (linseed) are mainly grown in India (930,000 ha), Canada (811,000 ha) and China (570,000 ha) [2].

Although flax fibers have now been largely replaced in textiles by cotton or synthetic fibers, they still are used in the fabrication of high-quality linen material and are also incorporated in the polymeric matrix during the fabrication of biocomposites in order to improve their mechanical properties [3].

In the flax stem, fibers are derived from the procambial cells of the protophloem [4] and are referred to as bast fibers. Individual bast fibers are extremely-long (≤ 77 mm) single cells [5] characterized by the presence of a very thick secondary cell-wall containing high amounts of cellulose and very low amounts of lignins as compared to the shorter fibers found in the xylem [6]. Bast fibers provide mechanical support to the plant and are organized in bundles that occupy the great majority of the external tissues between the epidermis and the vascular cambium. Flax fiber extraction is initiated during the retting process when stems are left on the ground in order to promote microbiological-mediated separation of the fiber bundles from the surrounding tissues [7].

Secondary cell-wall formation is most often studied in woody species such as poplar, pine and eucalyptus, or else in model species such as Arabidopsis and tobacco. More recently, growing interest in the unusual structure of the secondary cell wall of bast fibers (high crystalline cellulose, low lignin) is stimulating transcriptomic and genomic studies in flax and other bast fiber species [8–10]. A better understanding of the molecular regulation controlling fiber cell-wall metabolism should increase our knowledge of the very complex coordination necessary for the production and assembly of different polymers within the cell wall.

Flax is a suitable model for genetic and functional genomics as it is autogamous, easy to grow, has a relatively short vegetative stage, and can be genetically transformed [11–13]. The first studies on flax DNA structure were undertaken almost thirty years ago [14, 15] and the first cDNA sequence was deposited in Genbank in 1993 [16]. Previous analyses of the flax transcriptome have mainly involved studies of single gene expression by Northern blot [17–23] or multiple gene expression by classical RT-PCR [24]. Only relatively few studies [10, 25] have reported the use of qRT-PCR in flax. Quantitative real-time RT-PCR (qRT-PCR) is currently the most accurate method for detecting low abundant mRNAs and is capable of detecting slight variations in gene expression in different tissues of the same plant [26]. When compared to traditional methods used to evaluate transcript accumulation, the main advantages of qRT-PCR are its very high sensitivity and specificity. However, in order to obtain reliable and reproducible results, it is necessary to include suitable internal controls. Indeed, the high sensitivity of qRT-PCR can lead to misinterpretation of expression data [27], especially when biological samples show similar gene expression levels. Slight variations in the RNA integrity or quantity, in the reverse transcription yield, or more generally in the transcriptional activity of the studied organ can all have a major impact on quantification. It is therefore absolutely necessary to develop a normalisation strategy before undertaking gene expression studies. The identification of suitable control genes is thus an important challenge in transcriptome analyses and the major aim is to identify genes ubiquitously expressed in every tissue independently of the experimental context [28].

In the past a number of different genes were commonly used for normalising expression, especially in classical RT-PCR approaches. Most of these genes proved to be suitable for such studies at the time, mainly because of the low sensibility of RT-PCR. More recently, the same reference genes have also been used for qRT-PCR and it has been assumed that their expression is perfectly stable even though this has not necessarily been experimentally verified. Examples of such commonly-used genes include GAPDH (glyceraldehyde-3-phosphate dehydrogenase) [29], cyclophilin and 18S rRNA [30], actin [31], EF1alpha [32], ubiquitin [33] and beta tubulin [34].

Surprisingly, there are relatively few reports concerning the study of stably-expressed genes, referred to as housekeeping genes (HKGs), in different plant species, even though qRT-PCR is routinely used in many research projects. Suitable HKGs have been identified in model species with sequenced genomes such as rice [35], Arabidopsis [36], grapevine [37] and poplar [38], as well as in economically-important crops such as wheat [39], barley [40], Coffee [41], tomato [42], potato [43] and soybean [44]. In most cases, gene expression stability was determined by using statistical approaches such as geNorm [45] and NormFinder [46]. The geNorm algorithm allows the identification of the most suitable reference gene(s) and the optimal number of genes that should be used. It relies on a pairwise comparison and evaluates the variation of relative quantity ratios for each gene pair in a set of expression data. The NormFinder algorithm depends on a statistical and mathematical model that not only estimates the overall expression variation of a candidate gene, but also considers the variation between the chosen subgroups.

In the present study, we have evaluated the expression profiles of 13 putative HKGs during the development of flax plants with a special focus on the stem which contains two major cell-wall tissue types: 1) inner tissues (lignin-rich secondary cell walls of xylem) and 2) outer tissues (lignin-poor, cellulose-rich secondary cell walls of bast fibers). The two algorithms geNorm and NormFinder were used in order to determine the best reference genes needed for normalisation. The selected HKGs were then used to normalise gene expression in flax for an investigation of the expression of LuMYB1 - a flax MYB transcription factor highly similar to AtMYB59 [47].

MYB transcription factors are of particular interest in studies on secondary cell wall development since they have been shown to regulate the expression of genes coding enzymes of the phenylpropanoid biosynthetic pathway responsible for the production of lignin monomers [48, 49]. We decided to investigate the expression profiles of LuMYB1 since initial transcriptome studies in our laboratory (data not shown) had indicated that the gene was most highly expressed in inner stem tissues and therefore potentially associated with secondary cell wall formation in flax xylem tissues. In addition, AtMYB59 has been shown to undergo alternative splicing leading to different transcripts in rice and Arabidopsis, and therefore represents an interesting model for gene regulation studies.

Over the last years several very powerful techniques have been developed to detect differences in gene expression levels between different cell types, tissues and organs. Among these, qRT-PCR is now commonly used in many laboratories to undertake accurate expression profiling of different candidate genes that have been previously identified, either because they are known to be involved in a specific biological process, or else because they have shown an interesting expression profile in global expression analyses (microarrays).

However, reliable relative quantification can only be performed if an accurate normalization is performed following the choice and verification of suitable (i.e. stable) HKGs. Although microarray data have been analyzed to determine the most stably-expressed genes in Arabidopsis [36] and rice [51], it is not always possible to directly transpose suitable HKGs identified in one species to another species. For example UBQ10 is very stable in Arabidopsis [36] but not in rice [52], nor in soybean [44] nor in Brachypodium [53]. Therefore, since there are no universally-suitable reference genes, it is necessary to verify the expression levels of the candidate reference genes under the same experimental conditions as those used for the gene of interest.

In this study, we have measured the expression levels and stability of candidate genes in different tissues at three developmental stages in flax plants. Our major research interest in this species concerns the cell wall formation and development of the phloem (bast) fibers located in the external tissues of the stem. The polymer composition of the secondary cell walls in these long fiber cells (as well as that of the shorter xylem cells located in the inner tissues) has been previously analyzed [6] and we have shown that the abundance and structure of the phenolic polymer lignin varies according to the age and physiological state of the plant. In order to further investigate cell wall formation (including lignification) in flax, we are currently undertaking gene expression profiling by qRT-PCR on different tissues. However, due to important differences between the structure and the metabolic state of cells in flax inner and outer tissues, it is necessary to identify the most relevant combination of reference genes. Normalization of gene expression in flax has previously been done with either a homolog to an Arabidopsis NADH Ubiquinone oxidoreductase determined by cDNA-AFLP [25], or with an Elongation Factor EF1A homolog [10].

We originally selected 20 candidate HKGs from the literature and after determination of amplification efficiencies for the selected primers we retained 13 candidates. We then evaluated their normalization potential using 2 commonly-used algorithms (geNorm [45] and NormFinder [46]) in a systematic study of their expression stability in different flax tissues including stems, flowers, roots and leaves. The NormFinder algorithm [46] generally ranked GAPDH as the most stably expressed gene in all the flax samples irrespective of whether the samples were assembled into 1 main group or divided into several sub-groups. In only 1 case did NormFinder rank GAPDH as the second most stable gene. These results are consistent with those observed recently in Brachypodium distachyon [53] during cold/heat stress and in Coffea arabica [41] in different organs and tissues. In contrast to NormFinder, the geNorm algorithm indicated that we should use 3 HKGs (EFIA, ETIF5A, Ubiquitin) for studies involving a mix of root, leaf, stem and flower tissues. However, the same program indicated that only 2 genes (ETIF1, ETIF4F) should be used when just stem tissues were considered. While EF1A is often described as a stable gene and has been used as a reference in many species [37, 43, 52, 54], there are only relatively few reports of Eucaryotic Translation Initiation Factor genes being used as reference genes. For example in soybean, ETIF1A and ETIF1B genes were ranked as the most stable genes tested when 21 samples were pooled, or when photoperiodic treatments were modified [44]. In contrast, expression patterns of several translation initiation factors were shown to be unstable in wheat [55], and have been excluded as good candidate reference genes in poplar [38] and Darnel ryegrass [54].

We also evaluated the stability of tubulin and actin genes that are very often used to normalize expression data. Although actin has been shown to be a suitable normalization gene in developmental studies [44], it also appears to be unstable in many biological processes [27]. In the same way, tubulin has been shown to be stable during development in orobranche [56], but is apparently unstable (geNorm analyses) during development [53] and abiotic stress [54]. Our results would also suggest that these 2 genes are among the most instable for expression profiling during development in flax.

Altogether, our results showed that the geNorm and NormFinder algorithms came to different conclusions concerning the best candidate HKG(s) to use for expression normalization in flax. In order to see whether the use of these different HKGs modified the determined expression profiles (and hence possible interpretations of biological role) of a gene of interest, we analyzed the expression of LuMYB1 in different tissues and organs of Flax. LuMYB1 is a homolog of AtMYB59, a gene that has been shown to undergo alternative splicing in Arabidopsis and rice [47]. We decided to investigate the expression profile of LuMYB1 since our team is interested in gaining a better understanding of the regulation of cell wall formation in flax and different MYB transcription factors have been previously shown to play an important role in controlling this process [48]. In addition, preliminary transcriptome data in our lab had indicated that LuMYB1 was highly expressed in stem inner tissues. Such an observation could suggest that this gene is involved in the transcriptional regulation of cell-wall genes during the formation of xylem.

Our results showed that the 2 different flax spliceforms (LuMYB1-1 and LuMYB1-2) are strongly expressed in the inner stem tissues, but only expressed at low levels in external stem tissues, leaves and roots. In addition, only LuMYB1-2 is expressed in flowers. Such a result is interesting since previous work [47] has shown that the corresponding Arabidopsis gene is mainly expressed in roots, leaves and seedlings, but only poorly expressed in stem tissues. More recently [57], AtMYB59 has been shown to be specifically expressed during the S phase of the cell cycle in Arabidopsis cell suspensions, suggesting that this gene plays a key role in cell cycle regulation. Our observation that the corresponding flax gene (both spliceforms) are highly expressed in inner stem tissues could also suggest a role for LuMYB1 in cell cycle events since these flax tissues are the site of intense mitotic activity associated with the division of vascular cambium initials to form new xylem cells. However, further work is obviously necessary to confirm this hypothesis.

In addition, our results also revealed that both the choice of HKGs and the decision or not to analyze samples together or in sub-groups modified the determined expression profiles for the gene of interest LuMYB1. Although such choices did not affect the identity of the 4 (stem) tissues showing the highest expression levels, they did affect the ranking of these different tissues. Such an observation suggests that care should be taken when interpreting the biological significance of small differences in gene expression, and that in this case, the researchers should probably consider using different normalizations to verify their data. In addition, it is clear that functional approaches are also necessary to fully understand the role of a given gene in a particular biological process. Nevertheless, our results showing important differences in LuMYB1 expression levels (regardless of the normalization type) would suggest that this gene is potentially associated with the formation of xylem tissue in flax stems.

A considerable quantity of transcriptional data is currently available for major model plant species and in silico analyses can therefore be used to identify suitable HKGs for gene expression normalization in these species. However, for most other plant species, the suitability of potential HKGs identified in the literature or in heterologous databases must be verified by qRT-PCR. In this study, we have identified several suitable reference genes for studying developmental gene expression in flax, with a special focus on stem tissues. The use of HKGs identified by both geNorm and NormFinder algorithms allowed us to determine the expression profiles of a gene of interest (LuMYB1) in a large range of different flax tissues. Our results also showed that certain classically-used reference genes such as actin and tubulin are not necessarily the most suitable for studies of quantitative gene expression in flax. In conclusion, our study has identified suitable HKGs for future developmental transcriptome studies in flax. These genes also represent potentially interesting targets for normalizing gene expression in flax under stress conditions, although it would obviously be necessary to verify their stabilities in this case.