Skip to main content

Identification and validation of quantitative real-time reverse transcription PCR reference genes for gene expression analysis in teak (Tectona grandis L.f.)

Abstract

Background

Teak (Tectona grandis L.f.) is currently the preferred choice of the timber trade for fabrication of woody products due to its extraordinary qualities and is widely grown around the world. Gene expression studies are essential to explore wood formation of vascular plants, and quantitative real-time reverse transcription PCR (qRT-PCR) is a sensitive technique employed for quantifying gene expression levels. One or more appropriate reference genes are crucial to accurately compare mRNA transcripts through different tissues/organs and experimental conditions. Despite being the focus of some genetic studies, a lack of molecular information has hindered genetic exploration of teak. To date, qRT-PCR reference genes have not been identified and validated for teak.

Results

Identification and cloning of nine commonly used qRT-PCR reference genes from teak, including ribosomal protein 60s (rp60s), clathrin adaptor complexes medium subunit family (Cac), actin (Act), histone 3 (His3), sand family (Sand), β-Tubulin (Β-Tub), ubiquitin (Ubq), elongation factor 1-α (Ef-1α), and glyceraldehyde-3-phosphate dehydrogenase (GAPDH). Expression profiles of these genes were evaluated by qRT-PCR in six tissue and organ samples (leaf, flower, seedling, root, stem and branch secondary xylem) of teak. Appropriate gene cloning and sequencing, primer specificity and amplification efficiency was verified for each gene. Their stability as reference genes was validated by NormFinder, BestKeeper, geNorm and Delta Ct programs. Results obtained from all programs showed that TgUbq and TgEf-1α are the most stable genes to use as qRT-PCR reference genes and TgAct is the most unstable gene in teak. The relative expression of the teak cinnamyl alcohol dehydrogenase (TgCAD) gene in lignified tissues at different ages was assessed by qRT-PCR, using TgUbq and TgEf-1α as internal controls. These analyses exposed a consistent expression pattern with both reference genes.

Conclusion

This study proposes a first broad collection of teak tissue and organ mRNA expression data for nine selected candidate qRT-PCR reference genes. NormFinder, Bestkeeper, geNorm and Delta Ct analyses suggested that TgUbq and TgEf-1α have the highest expression stability and provided similar results when evaluating TgCAD gene expression, while the commonly used Act should be avoided.

Background

The flux of information from DNA to protein is connected by mRNA, and the level of mRNA transcription is one of the factors determining the degree of gene expression [1]. Changes in gene expression are critical for cell development [2], integration of metabolism [3] and resistance to biotic and abiotic stresses [4, 5], and as such are a research area of great interest to the fields of medicine, pharmacy, life sciences and agronomy.

Methods currently available for gene expression assessment include microarray analysis, Northern blotting, in situ hybridization, RNase protection assay, RNA sequencing (RNA-seq), qualitative RT-PCR, competitive RT-PCR, and quantitative real-time reverse transcription RT-PCR (qRT-PCR). The qRT-PCR is considered an efficient, safe (free of radioactive reagents), fast, affordable, reproducible, reliable and specific for quantifying levels of transcripts [6]. However, some variables such as the integrity, amount and purity of the RNA used as well as enzyme efficiency during cDNA synthesis and PCR amplification make an additional step to normalize the data necessary [7]. Normalization requires the use of one or more reference genes (also called internal control genes) for which expression is constant and stable at different developmental stages, nutritional conditions or experimental conditions [8]. Unfortunately, a gene has not been found for which expression is absolutely stable under all circumstances or across species that can be used indiscriminately for qRT-PCR analysis.

Bioinformatics tools have been developed to assess and identify the most suitable reference genes for qRT-PCR data normalization. geNorm shows expression stability throughout a set of housekeeping candidates [9], the Normfinder algorithm chooses the best candidate reference genes according to its calculations [10], while the Excel-based tool called BestKeeper determines the best candidate of pair-wise correlations [11]. Other statistical approaches used include Delta Ct [12] and “Stability index” methods [13].

Reference genes commonly used that present sufficiently stable expression are those related to cell maintenance such as actin, tubulin, glyceraldehyde-3-phosphate dehydrogenase, elongation factor 1-α and 18S ribosomal RNA[8, 14]. New genes have been studied as internal controls in model or commercial plants, such as Arabidopsis[15], Populus[16] and Brachypodium[17]. Tests for the selection of reference genes for qRT-PCR in teak have not been published yet.

Teak is a deciduous tree, native to countries of southeast Asia such as Myanmar, Thailand, India, Laos and Java [18]. Its wood is known internationally for its beauty, weightlessness, durability and weather resistance and it is used in the building of ships, furniture, house floors and walls, and general carpentry [19, 20]. Currently, the wood market has a great interest in teak extractives such as naphthoquinones and anthraquinones, which have shown remarkable antifungal and antitermitic effect [21, 22]. Additionally, teak populations serve significant environmental roles, as they can be used in agroforestry systems and forest recovery [22]. These characteristics make teak one of the most widely grown and economically profitable trees around the world [23]. Despite the great economic importance of teak, there are no studies of gene expression, the genome sequence is not available and sequenced genes are limited.

To select suitable qRT-PCR internal control genes for teak, this study analyzed the expression levels of candidate reference genes in different tissues and organs such as leaves, flowers, seedlings, roots, stem and branch secondary xylem of trees. Eight candidate reference genes were identified by their orthologous genes in model plants. These candidates were cloned, sequenced and tested. The selected genes are involved in different biological functions such as the formation of cellular cytoskeleton (Actin and β-Tubulin), elongation phase of translation (Elongation factor 1-α), DNA packaging (Histone 3) protein modification (Ubiquitin), intracellular transport (Clathrin adaptor complexes medium subunit family), vesicular transport (SAND family), protein biosynthesis (Ribosomal protein 60s) and carbohydrate metabolism (Glyceraldehyde-3-phosphate dehydrogenase).

Finally, in order to validate our results, the most stable reference genes were used to assess the TgCAD gene expression levels in different tissues and organs.

Results

Normalization of gene expression experiments, especially of qRT-PCR using a set of reference genes is currently a critical procedure when analyzing expression levels of target genes in different tissues or under different conditions. In the present study, nine potential reference genes for qRT-PCR of teak were assessed. A total of 36 cDNA samples including several organ types (leaf, root, and flower) and secondary xylem tissues from stems and branches of different ages were analyzed (Figure 1).

Figure 1
figure 1

Teak tissue and organ sample set. A = leaf, B = flower, C = seedling, D = root, E = stem secondary xylem, F = branch secondary xylem.

Identification and cloning of references genes in teak

As teak does not have the relevant genetic sequence information available in databases, it was necessary to design degenerate primers to amplify, clone and sequence the reference genes according to the most common genes used for qRT-PCR analysis in trees such as Platycladus orientalis[24], Vernicia fordii[25], Quercus suber[26], Populus euphratica[27] and Pyrus pyrifolia[28]. GAPDH (FN431982.1) was the only teak sequence available in GenBank (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/genbank). Therefore, we performed multiple nucleotide sequence alignment of the reference genes of different species for the remaining genes (Additional file 1). For each gene, at least four sequences were used in the alignments. Degenerate primers were designed to amplify the most conserved domains and at least 250 bp of the teak cDNA (Table 1, Additional file 2).

Table 1 Candidate reference genes, primers used to amplify in teak and their PCR parameters

The degenerate primers were able to produce specific amplicons ranging from 271 to 1440 bp using cDNA of teak leaves as template (Table 1). After gel purification, PCR fragments were cloned into the pJET1.2/Blunt vector (Thermo Scientific, USA) and transformed into DH5α™ competent cells (Life Technologies, USA). Recombinant colonies were selected to extract plasmid DNA for sequencing.

The teak nucleotide identities were checked by BLAST [29] against NCBI non redundant sequences (http://blast.ncbi.nlm.nih.gov) and the results showed that all the clones contained the expected fragments. The most conserved genes were TgAct and TgEf-1α with 92% of similarity, followed by TgUbq and Tgβ-Tub with 91% and 87%, respectively (data not shown). All genes showed at least 79% of similarity. Translated amino acid sequences were obtained by Expasy Translation Tool (http://web.expasy.org/translate/) and used to check for the presence of the expected domains in Pfam Database (http://pfam.sanger.ac.uk). Thereafter, teak amino acid sequences were compared against NCBI protein sequences with the algorithm tBLASTn (http://blast.ncbi.nlm.nih.gov). Results of the in silico analysis showed that all teak putative protein sequences possess the predicted domains, presenting high similarity with the selected reference genes (Additional file 3). At protein level, the most conserved genes were TgAct, TgEf-1α, TgHis3, Tgβ-Tub and TgUbq presenting 99% of similarity (data not shown).

Primer specificity and PCR efficiency

Real-time PCR primers (Table 2, Additional file 4) were designed to amplify the teak sequences of the eight clones (Tgrp60s, TgCac, TgAct, TgHis3, TgSand, Tgβ-Tub, TgUbq and TgEf-1α) and TgGAPDH (Table 2), and were used to detect transcript levels. Primer specificity was evaluated with a single peak in all ten melting curves (Figure 2) and as a single band in the agarose gel analysis (Additional file 5). qRT-PCR efficiency (E) varied from 91.4% for Tgrp60s to 108.5% for TgSand and correlation coefficients (R2) oscillated from 98.2% for TgΒ-Tub to 99.9% for TgEf-1α (Table 2). The acceptable range for PCR efficiencies calculated using standard curve serial dilution experiments is 90–110% (i.e. a slope between 3.1 and 3.58) [30]. The annealing temperature of 65°C was effective for all primers; nevertheless, its choice can impact on the efficiency of the reaction. Altogether, the results showed that the chosen primers accurately amplified the candidate reference genes.

Table 2 Candidate reference genes, TgCAD target gene, specific qRT-PCR primers and different parameters derived from qRT-PCR analysis
Figure 2
figure 2

Specificity of qRT-PCR amplification. Melting curves (dissociation curves) of the 10 amplicons (rp60s, Act, Cad, Cac, Ef-1α, GAPDH, His3, Sand, B-Tub, Ubq genes) after the qRT-PCR reactions, all showing one peak.

To compare the differences in transcript levels between reference genes, the Cq range was determined and the coefficient of variance was calculated for each gene across all samples based on the interquartile range (25-75% percentiles). The average Cq values of the different genes ranged from 22 to 34 cycles (Figure 3, Table 3). TgUbq, Tgβ-Tub and Tgrp60s showed the narrowest variance (lowest Cq dispersion), while TgAct and TgCac exhibited widest variance (highest dispersion). The gene with the most abundant transcript level was TgEf-1α while TgCac was the least abundant, reaching mean threshold fluorescence with 23 and 31 amplification cycles, respectively.

Figure 3
figure 3

Expression levels of candidate reference genes in different plant samples. Expression data displayed as Cq values for each reference gene in all T. grandis samples. The horizontal lines of the box indicate the 25th and 75th quartiles. The central horizontal line across the box is depicted as the median. The whisker caps represent the maximum and minimum values. Dots represent outliers. Genes are in order from the most (lower Cq, on the left) to the least abundantly expressed (higher Cq, on the right).

Table 3 Descriptive statistics and expression level obtained by BestKeeper

Expression stability of the nine candidate reference genes

To evaluate the reference genes’ expression stability, four different methodologies were used: geNorm, NormFinder, BestKeeper and Delta Ct.

geNorm

geNorm was used to rank the reference genes by calculating the gene expression stability value M, which corresponds to the average pairwise variation (V) of a particular gene with all other control genes [9]. The most stable reference gene has the lowest M value, while the least stable has the highest M value. To identify reference genes with stable expression, geNorm indicates genes with M values below the threshold of 1.5, however Vandesompele et al. [9] suggests M values lower than 1.0 to ensure the selection of the most stable genes. When all 36 samples were analyzed together with geNorm (Figure 4), eight genes had M < 1.0, with TgUbq and TgEF1α showing the highest expression stability (M = 0.295) in different tissues. Act was the only gene with M > 1.0, with the lowest expression stability of 1.035.

Figure 4
figure 4

Gene expression stability of the candidate reference genes calculated by different statistical methods. Ranking of each candidate reference gene (rp60s, Act, Cac, Ef-1α, GAPDH, His3, Sand, B-Tub, Ubq) calculated by NormFinder, BestKeeper, geNorm and Delta Ct methods, for all tested cDNA samples (leaf, flower, seedling, root, stem secondary xylem, branch secondary xylem).

To obtain reliable results from qRT-PCR studies, two or more reference genes should be used for data normalization. The optimal number of reference genes can be determined by calculating the pairwise variation (Vn/n+1) using the geNorm algorithm [9]. It is calculated between the two sequential normalization factors (NF), NFn and NFn+1, for all the samples analyzed. Slight variations mean addition of another gene has a low effect on the normalization. Vandesompele et al. [9] proposed 0.15 as the cut-off value for V, below which the inclusion of an additional control gene is not required. This means that if Vn/n+1 < 0.15, it is not necessary to use ≥ n + 1 reference genes for normalization. In this study, the paired variable coefficients indicated that the inclusion of the third reference gene (i.e. TgUbq, TgEF1α and TgGAPDH) would be useful for normalization when considering total samples and only lignified samples, whereas two stable reference genes (TgGAPDH and Tgrp60s) can be employed when analyzing non-lignified tissues (Figure 5). However, if the samples of stem secondary xylem (the most lignified tissue) are excluded from the analysis, two reference genes (TgUbq and TgEF1α) would be optimal for normalizing gene expression (Figure 6).

Figure 5
figure 5

Gene expression differences among the candidate reference genes analyzed by NormFinder. Black circles represent the log-transformed gene expression levels. Vertical bars give a confidence interval for the inter-tissue variation. Top and bottom lines from the graphic represent the maximum standard deviation of the candidate reference genes, with the difference log expression levels between 0.2 and -0.2.

Figure 6
figure 6

Pairwise variations (V) calculated by geNorm to determine the optimal number of reference genes. The average pairwise variations Vn/n+1 was analyzed between the normalization factors NFn and NFn+1 to indicate the optimal number of reference genes required for qRT-PCR data normalization in all the samples, lignified tissues (root, stem secondary xylem, branch secondary xylem), non-lignified tissues (leaf, flower, seedling) and in all samples minus stem secondary xylem.

NormFinder

NormFinder is a Microsoft Excel-based Visual Basic application that allows estimation of stability values of single candidate reference genes. The algorithm is based on intra- and inter-group variations and combines both results into a stability value for each candidate reference gene [10]. The results of the NormFinder analysis were somewhat similar to those of geNorm. Both methods ranked TgUbq, TgEF1α and Tgrp60s as among the four most stable reference genes (Figure 5) and TgAct and TgHis3 as the least stable (Figure 4). However, Tgβ-Tub emerged as the third most stably expressed using NormFinder, whereas it was ranked seventh by geNorm. These discrepant results could be explained due to inter-tissue expression variations detected by NormFinder analysis, which is not take account for gene stability calculations in the geNorm algorithm. When considering only intra-tissue variations, Tgβ–Tub was the most stable gene in lignified tissues (i.e. roots, branches and stems) (Table 4). However, in non-lignified tissues (leaves, flowers and seedlings) Tgβ–Tub was ranked eighth of nine genes, corroborating results obtained by geNorm.

Table 4 NormFinder intragroup expression stability for teak candidate reference genes

BestKeeper

The Bestkeeper software was adopted for descriptive analysis. The program is an Excel-based software tool that estimates gene expression stability based on the coefficient of correlation (r) between each reference gene and an index, defined as the geometric mean of all candidate reference gene Ct (or CP) values [11]. The BestKeeper also calculates the CP standard deviation (SD) and the coefficient of variance (CV) of each candidate gene. Reference genes with SD values >1 are considered not stable and should be avoided. Results of analysis are shown in Table 3. Similarly to geNorm and NormFinder, BestKeeper ranked Act and His3 among the three least stable reference genes with SD values >1.0 (1.42 and 1.02 respectively). In addition, Tgβ–Tub presented unstable expression with SD value of 1.21, being one of the least stable genes as observed in the geNorm analysis (Figure 4). The best reference genes are those that have the lowest coefficient of variance and standard deviation. In this study, Tgrp60s and TgUbq had CV ± SD values of 1.78 ± 0.51 and 2.23 ± 0.60, respectively, displaying a stable expression in all samples. The results of the BestKeeper analysis showed a similar pattern of stability to those obtained from geNorm, which described TgUbq, Tgrp60s, TgEF1α and TgGAPDH as the four best reference genes for the normalization of qRT-PCR data in teak. In the NormFinder analysis, TgGAPDH was replaced by Tgβ–Tub among the top four regarding stability (Figure 4).

Delta Ct

The Delta Ct method is based on the ‘pairs of genes’ comparison using a simple ΔCt approach [12]. The formula used in this method is similar to the standard comparative Ct method (ΔΔCt) [31] except that no endogenous reference gene is incorporated since the purpose is to define stably expressed genes to normalize. In this approach, all pairs of genes are compared to each other and the genes are ranked according to the ΔCt values, from lowest to highest [12]. As observed in geNorm and Bestkeeper analysis, TgEF1α, TgUbq and TgGAPDH were the best reference genes, as they had the highest expression stability (lowest ΔCt values) (Figure 4). As was detected in all used programs, TgHis3 and TgAct were the least stable internal controls for gene expression normalization.

Validation of TgUbq and TgEF-1a as internal controls to assess expression of the teak cinnamyl alcohol dehydrogenase gene in lignified tissues

The use of different reference genes to evaluate relative expression data has an important impact on the final normalized results. As TgUbq and TgEF-1a showed the best stability values by geNorm, Delta Ct, NormFinder and BestKeeper analyses (Figure 7), they were used to evaluate the transcript level of a gene of interest, the teak cinnamyl alcohol dehydrogenase (TgCAD). TgCAD was identified and cloned using the same methodologies described for the reference genes. To validate the selected reference genes, the transcript levels were quantified in leaves from four month-old greenhouse grown teak, and in lignified tissues and organs such as stem secondary xylem from 60 year-old trees, stem from 1 year-old plants and branch secondary xylem from 60 and 12 year-old teak trees. Results showed that no matter which gene is used (TgUbq or TgEF-1a), TgCAD expression decreased in the following order: leaf > stem from 1 year-old plants > branch secondary xylem from 12 year-old trees > branch secondary xylem from 60 year-old trees > stem secondary xylem from 60 year-old trees (Figure 8). On the other hand, between tissues, leaf was the tissue with highest expression and stem and branch secondary xylem from 60 year-old teak trees were the tissues with lowest expression of the TgCAD gene.

Figure 7
figure 7

Venn diagrams. (A) the most stable reference genes present in the first four positions and (B) the least stable genes present in the last five positions identified by the NormFinder, BestKeeper, geNorm and Delta Ct methods. Diagrams were performed with the Smartdraw® program.

Figure 8
figure 8

Expression levels of the TgCAD gene. It was used different tissues and organ ages of teak tree, using the best validated reference genes (TgUbq and TgEF-1a) for normalization and the results are represented as mean fold changes in relative expression compared to stem secondary xylem from 60 years-old trees. Bars are mean standard deviation calculated from the 3 biological replicates.

Discussion

In plant molecular biological research, qRT-PCR has improved the detection and quantification of expression profiles of target genes due to its sensitivity, specificity and accuracy. For correct qRT-PCR measurements, reference genes are used as endogenous controls for gene expression normalization when analyzing the expression of genes of interest [5, 24]. Therefore, a careful choice of reference genes is essential to obtain an accurate quantification of the target gene transcript levels [25].

Currently, teak is one of the most important trees worldwide due to its wood’s properties. Despite the growing importance of teak wood and extractives to the world market, the number of studies adopting techniques of modern biology for teak improvement is still quite limited. As far as has been documented, this is the first study of cloning and expression stability of qRT-PCR reference genes in teak tissues. In total, eight candidate genes (Tgrp60s, TgCac, TgAct, TgHis3, TgSand, TgΒ-Tub, TgUbq, and TgEf-1α) were successfully identified.

One of the challenges of studying gene expression in trees is to ensure good quality of total RNA isolated from stems and branches, which are woody tissues with high lignin contents. In this study, the use of the Salzman protocol [32] for total RNA extraction, followed by a DNAse I (Promega) treatment, provided high quality RNA from all lignified tissues and was chosen as standard method for RNA extraction.

The stability of nine reference genes for qRT-PCR normalization was assessed by four statistical approaches, geNorm, NormFinder, BestKeeper and Delta Ct method, in teak lignified and non-lignified tissues at different developmental stages. In spite of some inconsistencies that are usually observed between these methods [24, 26, 33], our results were quite constant regardless of the algorithm used for analysis. When considering the rank of four most and least stable genes, TgUbq and TgEF1α were selected among the most stable in all methods, while TgAct, TgHis3 and TgCac showed the least expression stability. The only clear discrepancy within the results was the inclusion of Tgβ-Tub in the most stable group by NormFinder and in the least stable group by the other programs (Figure 4). In the NormFinder intra-tissues analysis, Tgβ–Tub was the most stable gene in lignified tissues, whereas it was ranked eighth of the nine genes in non-lignified tissues (Table 4). These results suggest that Tgβ–Tub is a suitable reference gene in lignified tissues and could be used as internal control for quantifying gene expression in them.

Among recent studies in trees searching for suitable reference genes, control genes such as Act, Ubq, Ef1α, α-Tub, Cac, Sand, β-Tub were considered to be stable in various tissues and different conditions [2428]. Act, Ubq and Ef1α were shown to be suitable reference genes for normalization in lignified tissues of Vernicia fordii[25] and Quercus suber[26]. In our analysis, the most stable genes were Ubq and EF1α (Figure 7).

Genes encoding elongation factor-1α and ubiquitin are frequently considered consistent reference genes under different experimental conditions. EF1α has been found to be one of the most stable reference genes in several plants and conditions such as Nicotiana tabacum[34], Lolium perenne[35] and Capsicum annuum[36]. The Ubq gene showed high stability for qRT-PCR normalization in Platycladus orientalis[24] and Brachypodium distachyon[37]. In combination, Ubq and EF1α showed stable expression across different tissues of Vernicia fordii[25] and Dimocarpus longan[33]. However, Ef1α and Ubq were the most variable reference genes in Lycopersicum esculentum[38] and Euphorbia esula[39], respectively, suggesting that these genes might not be suitable for qRT-PCR normalization in some plants and/or conditions. Although Act is one of the most commonly used reference gene in plants, in teak it showed low stability when assessed in different tissue samples and with different statistical methods. Similar results were observed in Nicotiana tabacum plants with viral infections [40] and Glycine max[41].

Studies have shown that the expression of reference genes can vary significantly under different experimental conditions [33, 42]. To mitigate these variations, the use of multiple reference genes to assess target gene expression is appropriate. geNorm analysis of paired variable coefficients suggested the inclusion of a third reference gene (i.e. TgUbq, TgEF1α and TgGAPDH) for normalization when considering the total number of samples (Figure 6). Although the cut off value ≤0.15 is frequently used to confirm the optimal number of reference genes [9], this is not an absolute number because small datasets require fewer reference genes than larger ones and previous studies have reported proper normalization with higher cut-off values [40]. In this study, the combination of the two most stable reference genes (TgUbq and TgEF1α) to evaluate expression stability in all samples provided a coefficient of 0.17 and, thus, can be sufficient for the normalization of qRT-PCR data in teak, especially considering that the use of more than two reference genes in large scale gene expression profiles will significantly increase the costs of analysis.

To validate the utility of TgUbq and TgEF-1a as reference genes, the expression profile of TgCAD was assessed in teak leaves and lignified tissues collected from plants in the field at different development stages (Figure 8). CAD functions in one of the final steps of monolignol biosynthesis in the phenylpropanoid pathway and its study is essential to understand the lignin deposition and cell wall formation in trees. It catalyzes the NAPDH-dependent reduction of cinnamyl aldehydes to cinnamyl alcohols prior to their transport to the secondary cell wall for polymerization into the lignin heteropolymer [43]. In plants, CAD expression may vary according to tissue and development stage [44]. In addition, it has been shown that CAD/CAD-like genes are differentially expressed in plants infected with pests and pathogens [45, 46].

Using TgUbq or TgEF-1a as reference genes, qRT-PCR results showed that TgCAD was strongly expressed in leaves (average 133-fold), followed by stems from 1 year-old plants (40-fold), branch secondary xylem from 12 year-old trees (24-fold), branch secondary xylem from 60 year-old trees (5-fold) and stem secondary xylem from 60 year-old trees (calibrator) (Figure 8). We observed higher expression of TgCAD in younger lignified tissues compared to older ones, probably due to less lignin deposition and secondary wall formation in 60 year-old trees. The TgCAD expression in lignified tissues and leaves presented the same pattern whichever internal control used, indicating that the reference genes identified in this study are suitable for qRT-PCR normalization in different tissues and plant ages.

This is the first attempt to identify qRT-PCR reference genes in several teak tissues. These results suggest the use of TgUbq and TgEF-1a as the best combination of reference genes for gene expression assessment in leaves, flowers, seedlings, roots, and lignified stem and branch secondary xylem of varying ages in teak. In addition, we recommend the researchers to validate the reference genes in their samples of interest before performing any experiment. The different tissues show that TgAct is not a suitable reference gene to normalize gene expression in this tree, highlighting the need to evaluate commonly used reference genes for particular species, conditions, tissues and organs. Finally, they advise that the use of reference genes without validation may reduce precision or produce misleading results.

Conclusions

To the best of our knowledge, this study is the first attempt at cloning, sequencing and evaluating a set of commonly used candidate reference genes for the normalization of gene expression analysis using qRT-PCR in teak. Our data showed that expression stability varied considerably among the nine genes tested in the different samples of teak tissues. Stability analysis using NormFinder, Bestkeeper, geNorm and Delta Ct showed that TgUbq and TgEF-1a are the most stable genes across different tissues and organs, while TgAct was deemed to be unsuitable as a reference gene. TgCAD expression analyses confirmed TgUbq and TgEF-1a stability for correct normalization in teak. Consequently, they can be used in future gene expression studies of target genes in different teak tissues.

Methods

Plant material

Roots, seedlings and leaves were obtained from fifteen four month-old greenhouse grown teak. Flowers and branch and stem secondary xylem (Figure 1) were collected from fifteen twelve year-old teak trees located in Piracicaba, São Paulo State, Brazil. All the harvested tissues were immediately frozen by immersion in liquid nitrogen and stored at -80°C.

Total RNA extraction, purification and quality controls

Frozen tissue samples of 1.0 g were weighed and ground to fine powder in liquid nitrogen using a sterilized mortar and pestle. The fifteen samples from each tissue or organ were divided into three different RNA extractions (five samples for each extraction). Total RNA was extracted following a protocol developed for lignified tissues by Salzman et al. [32]. RNA quality assessment included purity (absence of protein and DNA) and integrity (absence of RNA degradation). 1 μl of each extraction was analyzed spectrophotometrically using a Nanodrop ND-1000 Spectrophotometer (NanoDrop Technologies Inc., USA) and only RNA samples with 260/280 ratio between 1.9 and 2.1 and 260/230 ratio greater than 2.0 were used for subsequent analyses. The concentration of each sample was approximately 2 μg/μl, so they were diluted to a final concentration of 1 μg/μl and 4 μg of total RNA from each sample was treated with DNAse I (Promega). Then, 0.5 μl of each treated sample was analyzed in agarose gels, all displaying clear bands corresponding to rRNA, absence of DNA and no degradation. In addition, PCR control reactions to examine for genomic DNA contamination were performed using total RNA without reverse transcription as template, and negative results (absence of bands) were assessed by electrophoresis on a 1% (w/v) agarose gel with ethidium bromide staining.

cDNA synthesis

Two cDNA samples were synthesized from the three extractions of each tissue or organ from 1.0 μg of the treated RNA using the SuperScript™ III First-Strand Synthesis System for RT-PCR (Invitrogen) according to the manufacturer’s instructions. Each cDNA sample concentration was determined using the Nanodrop ND-1000 Spectrophotometer (NanoDrop Technologies Inc., USA) to be approximately 2000 ng/μl. A concentration of 100 ng/μl (1:20 dilution) and 25 ng/μl (1:80 dilution) was used for PCR amplification and qRT-PCR expression experiments, respectively.

Multiple sequence alignments, PCR and qRT-PCR primer design

Primers (Table 1) were manually designed flanking the conserved domains of rp60s, Cac, Act, His3, Sand, Β-Tub, Ubq, and Ef-1α after doing Clustal alignment (http://www.ebi.ac.uk/Tools/msa/clustalw2) of several orthologous plant sequences obtained from GenBank (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/genbank) to amplify by PCR those genes from teak leaf cDNA (Additional file 2). The eight amplified fragments gel electrophoresis were excised, purified with Fragment CleanUp® (Invisorb, USA) and inserted into the pJET1.2/Blunt vector from the CloneJet™ PCR Cloning Kit (Thermo Scientific, USA) following the manufacturer’s recommendations. Plasmids were cloned in DH5α™ competent cells (Life Technologies, USA) and recombinant colonies were sequenced with the 3100 Genetic Analyzer (Applied Biosystems, USA) using pJET1.2/Blunt vector specific primers. Finally, sequences were blasted using blastx (http://blast.ncbi.nlm.nih.gov) to confirm their percentage amino acid similarity to the conserved domains and were translated (http://web.expasy.org/translate/) to amino acid sequences, which were submitted to PFAM search (http://pfam.sanger.ac.uk/search) (Sanger Institute, England) to confirm the presence of each gene’s canonical protein domains. The primers for qRT-PCR were designed flanking the eight cloned teak sequences and GAPDH (Table 2) with OligoPerfectTM Designer (Life technologies, USA) with default parameters. Teak candidate reference genes, TgCAD target gene, NCBI accession numbers, qRT-PCR primer information and different parameters derived from qRT-PCR analysis are shown in Table 2.

Primer specificity, qRT-PCR Efficiency and R2

Confirmation of primer specificity was based on the dissociation curve at the end of each run (Figure 3). To determine the amplification efficiencies of the candidate genes, it was used cDNA samples from the teak leaf with five dilutions to obtain the standard curve, and then the PCR efficiency for each gene was calculated according to the equation (1 + E) = 10slope. The correlation coefficient (R2) and slope values were obtained from the standard curve (Table 2).

Quantitative real-time reverse transcription PCR

The qRT-PCR mixture contained 5.0 μl of a 1:80 dilution of the six synthesized cDNAs from each tissue or organ, primers to a final concentration of 50 μM each, 12.5 μl of the SYBR Green PCR Master Mix (Applied Biosystems, USA) and PCR-grade water up to a total volume of 25 μl. Each gene reaction was performed in technical replicate. PCR reactions without template were also done as negative controls for each primer pair. The quantitative PCRs were performed employing the StepOnePlus™ System (Applied Biosystems, USA). All PCR reactions were performed under the following conditions: 2 min at 50°C, 2 min at 95°C, and 45 cycles of 15 s at 95°C and 1 min at 65°C in 96-well optical reaction plates (Applied Biosystems, USA). Leaf samples were used as calibrator to normalize the values between different plates.

Analysis of gene expression stability

Gene expression stability was evaluated by applying four different statistical approaches: geNorm [9], NormFinder [10], Bestkeeper [11] and Delta Ct [12]. qRT-PCR data was exported from the StepOnePlus™ System (Applied Biosystems, USA) into an Excel datasheet (Microsoft Excel 2003) as Raw Crossing Point data (Additional file 6) and those values were log transformed by the 2-ΔCt method for further requirements. Each of these approaches generated a measure of reference gene stability, by which each gene was ranked. Venn diagrams were constructed with the Smartdraw® program.

Validation of reference genes

One gene of interest, putatively coding for a cinammyl alcohol dehydrogenase (CAD) (Table 2), an enzyme involved in lignin biosynthesis, one of the terminal steps of the phenylpropanoid pathway, was used to validate the best two reference genes. The relative expression level of the target gene was determined in leaves, and the lignified tissues of stem and branch secondary xylem of 60 year-old trees, stem of 1 year-old trees and branch secondary xylem of 12 year-old trees, expecting a higher expression level in younger tissues with a continuous secondary wall formation. The experimental procedure was the same used for the selection of the reference genes. Stem secondary xylem of 60 year-old tree samples were chosen as calibrator.

Abbreviations

Tg :

Tectona grandis

rp60s :

Ribosomal protein 60s

Cac:

Clathrin adaptor complexes medium subunit family

Act:

Actin

His3:

Histone 3

Sand:

Sand family

Β-Tub:

β-Tubulin

Ubq:

Ubiquitin

Ef-1α:

Elongation factor 1-α

GAPDH:

Glyceraldehyde-3-phosphate dehydrogenase

CAD:

Cinnamyl alcohol dehydrogenase

cDNA:

Complementary DNA

ΔCt:

Delta cycle threshold

mRNA:

Messenger RNA

PCR:

Polymerase chain reaction

qRT-PCR:

Quantitative real-time reverse transcription PCR

CP:

Crossing-point cycle number

GM:

Geometric mean

AM:

Arithmetic mean

SD:

Standard deviation

CV:

Coefficient of variation

bp:

Base pairs

Tm:

Melting temperature

R2:

Correlation coefficient.

References

  1. Hauptman N, Glavac D: MicroRNAs and long non-coding RNAs: prospects in diagnostics and therapy of cancer. Radiol Oncol. 2013, 47 (4): 311-318.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  2. Conte I, Banfi S, Bovolenta P: Non-coding RNAs in the development of sensory organs and related diseases. Cell Mol Life Sci. 2013, 70: 4141-4155.

    Article  PubMed  CAS  Google Scholar 

  3. Ohdan T, Francisco PB, Sawada T, Hirose T, Terao T, Satoh H, Nakamura Y: Expression profiling of genes involved in starch synthesis in sink and source organs of rice. J Exp Bot. 2005, 56 (422): 3229-3244.

    Article  PubMed  CAS  Google Scholar 

  4. Mackay IM, Arden KE, Nitsche A: Real-time PCR in virology. Nucleic Acids Res. 2002, 30 (6): 1292-1305.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  5. Zhu J, Zhang L, Li W, Han S, Yang W, Qi L: Reference gene selection for quantitative real-time PCR normalization in Caragana intermedia under different abiotic stress conditions. PLoS One. 2013, 8 (1): e53196-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  6. Kubista M, Andrade JM, Bengtsson M, Forootan A, Jonak J, Lind K, Sindelka R, Sjoback R, Sjogreen B, Strombom L, Stahlberg A, Zoric N: The real-time polymerase chain reaction. Mol Asp Med. 2006, 27: 95-125.

    Article  CAS  Google Scholar 

  7. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL, Vandesompele J, Wittwer CT: The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009, 55 (4): 611-622.

    Article  PubMed  CAS  Google Scholar 

  8. Dheda K, Huggett JF, Chang JS, Kima LU, Bustin SA, Johnson MA, Rook GAW, Zumla A: The implications of using an inappropriate reference gene for real-time reverse transcription PCR data normalization. Anal Biochem. 2005, 344: 141-143.

    Article  PubMed  CAS  Google Scholar 

  9. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3 (7): 34-

    Article  Google Scholar 

  10. Andersen CL, Jensen JL, Orntoft TF: Normalization of real-time quatitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64: 5245-5250.

    Article  PubMed  CAS  Google Scholar 

  11. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP: Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper-Excel-based tool using pair-wise correlations. Biotechnol Lett. 2004, 26: 509-515.

    Article  PubMed  CAS  Google Scholar 

  12. Silver N, Best S, Jiang J, Thein SL: Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR. BMC Mol Biol. 2006, 7: 33-

    Article  PubMed  PubMed Central  Google Scholar 

  13. Brunner AM, Yakovlev IA, Strauss SH: Validating internal controls for quantitative plant gene expression studies. BMC Plant Biol. 2004, 4: 14-

    Article  PubMed  PubMed Central  Google Scholar 

  14. Radoni A, Thulke S, Mackay IM, Landt O, Siegert W, Nitsche A: Guideline to reference gene selection for quantitative real-time PCR. Biochem Biophys Res Commun. 2004, 313: 856-862.

    Article  Google Scholar 

  15. Czechowski T, Stitt M, Altmann T, Udvardi MK, Scheible W: Genome-wide identification and testing of superior reference genes for transcript normalization in arabidopsis. Plant Physiol. 2005, 139: 5-17.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  16. Xu M, Zhang B, Su X, Zhang S, Huang M: Reference gene selection for quantitative real-time polymerase chain reaction in Populus. Anal Biochem. 2011, 408: 337-339.

    Article  PubMed  CAS  Google Scholar 

  17. Chambers JP, Behpouri A, Bird A, Ng CK: Evaluation of the use of the Polyubiquitin Genes, Ubi4 and Ubi10 as reference genes for expression studies in Brachypodium distachyon. PLoS One. 2012, 7 (11): e49372-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  18. Verhaegen D, Fofana IJ, Logossa ZA, Ofori D: What is the genetic origin of teak (Tectona grandis L.) introduced in Africa and in Indonesia?. Tree Genetics Genomes. 2010, 6: 717-733.

    Article  Google Scholar 

  19. Miranda I, Sousa V, Pereira H: Wood properties of teak (Tectona grandis) from a mature unmanaged stand in East Timor. J Wood Sci. 2011, 57: 171-178.

    Article  CAS  Google Scholar 

  20. Lukmandaru G, Takahashi K: Variation in the natural termite resistance of teak (Tectona grandis Linn. fil.) wood as a function of tree age. Ann For Sci. 2008, 65: 708-

    Article  Google Scholar 

  21. Guerrero-Vásquez GA, Andrade CKZ, Molinillo JMG, Macías FA: Practical first total synthesis of the potent phytotoxic (±)-naphthotectone, isolated from Tectona grandis. Eur J Org Chem. 2013, 2013 (27): 6175-6180.

    Article  Google Scholar 

  22. Healey SP, Gara RI: The effect of a teak (Tectona grandis) plantation on the establishment of native species in an abandoned pasture in Costa Rica. For Ecol Manag. 2003, 176: 497-507.

    Article  Google Scholar 

  23. Hallett JT, Diaz-Calvo J, Villa-Castillo J, Wagner MR: Teak plantations: economic bonanza or environmental disaster?. J For. 2011, 109 (5): 288-292.

    Google Scholar 

  24. Chang E, Shi S, Liu J, Cheng T, Xue L, Yang X, Yang W, Lan Q, Jiang Z: Selection of reference genes for quantitative gene expression studies in Platycladus orientalis (Cupressaceae) using real-time PCR. PLoS One. 2012, 7 (3): e33278-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  25. Han X, Lu M, Chen Y, Zhan Z, Cui Q, Wang Y: Selection of reliable reference genes for gene expression studies using real-time PCR in tung tree during seed development. PLoS One. 2012, 7 (8): e43084-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  26. Marum L, Miguel A, Ricardo CP, Miguel C: Reference gene selection for quantitative real-time PCR normalization in Quercus suber. PLoS One. 2012, 7 (4): e35113-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  27. Wang HL, Chen J, Tian Q, Wang S, Xia X, Yin W: Identification and validation of reference genes for Populus euphratica gene expression analysis during abiotic stresses by quantitative real-time PCR. Physiol Plant. 2014, in press

    Google Scholar 

  28. Imai T, Ubi BE, Saito T, Moriguchi T: Evaluation of reference genes for accurate normalization of gene expression for real time-quantitative PCR in Pyrus pyrifolia using different tissue samples and seasonal conditions. PLoS One. 2014, 9 (1): e86492-

    Article  PubMed  PubMed Central  Google Scholar 

  29. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990, 215: 403-410.

    Article  PubMed  CAS  Google Scholar 

  30. Pfaffl MW: Quantification Strategies in Real-Time PCR. The Real-Time PCR Encyclopedia, A–Z of Quantitative PCR. Edited by: Bustin SA. 2004, La Jolla: International University Line, 87-120.

    Google Scholar 

  31. Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔC T method. Methods. 2001, 25: 402-408.

    Article  PubMed  CAS  Google Scholar 

  32. Salzman RA, Fujita T, Zhu-Salzman K, Hasegawa PM, Bressan RA: An improved RNA isolation method for plant tissues containing high levels of phenolic compounds or carbohydrates. Plant Mol Biol reports/ISPMB. 1999, 17: 11-17.

    Article  CAS  Google Scholar 

  33. Lin YL, Lai ZX: Reference gene selection for qPCR analysis during somatic embryogenesis in longan tree. Plant Sci. 2010, 178: 359-365.

    Article  CAS  Google Scholar 

  34. Schmidt GW, Delaney SK: Stable internal reference genes for normalization of real-time RT-PCR in tobacco (Nicotiana tabacum) during development and abiotic stress. Mol Genet Genomics. 2010, 283: 233-241.

    Article  PubMed  CAS  Google Scholar 

  35. Lee JM, Roche JR, Donaghy DJ, Thrush A, Sathish P: Validation of reference genes for quantitative RT-PCR studies of gene expression in perennial ryegrass (Lolium perenne L.). BMC Mol Biol. 2010, 11: 8-

    Article  PubMed  PubMed Central  Google Scholar 

  36. Bin WS, Wei LK, Ping DW, Li Z, Wei G, Bing LJ, Gui PB, Jian WH, Feng CJ: Evaluation of appropriate reference genes for gene expression studies in pepper by quantitative real-time PCR. Mol Breed. 2012, 30: 1393-1400.

    Article  Google Scholar 

  37. Hong S, Seo PJ, Yang M, Xiang F, Park C: Exploring valid reference genes for gene expression studies in Brachypodium distachyon by real-time PCR. BMC Plant Biol. 2008, 8: 112-

    Article  PubMed  PubMed Central  Google Scholar 

  38. Expósito-Rodríguez M, Borges AA, Borges-Pérez A, Pérez JA: Selection of internal control genes for quantitative real-time RT-PCR studies during tomato development process. BMC Plant Biol. 2008, 8: 131-

    Article  PubMed  PubMed Central  Google Scholar 

  39. Chao WS, Doğramaci M, Foley ME, Horvath DP, Anderson JV: Selection and validation of endogenous reference genes for qRT-PCR analysis in leafy spurge (euphorbia esula). PLoS One. 2012, 7 (8): e42839-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  40. Liu D, Shi L, Han C, Yu J, Li D, Zhang Y: Validation of reference genes for gene expression studies in virus-infected Nicotiana benthamiana using quantitative real-time PCR. PLoS One. 2012, 7 (9): e46451-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  41. Libault M, Thibivilliers S, Bilgin DD, Radwan O, Benitez M, Clough SJ, Stacey G: Identification of four soybean reference genes for gene expression normalization. Plant Genome. 2008, 1 (1): 44-54.

    Article  CAS  Google Scholar 

  42. Barsalobres-Cavallari CF, Severino FE, Maluf MP, Maia IG: Identification of suitable internal control genes for expression studies in Coffea Arabica under different experimental conditions. BMC Mol Biol. 2009, 10: 1-

    Article  PubMed  PubMed Central  Google Scholar 

  43. Trabucco GM, Matos DA, Lee SJ, Saathoff AJ, Priest HD, Mockler TC, Sarath G, Hazen SP: Functional characterization of Cinnamyl alcohol dehydrogenase and caffeic acid O-methyltransferase in Brachypodium distachyon. BMC Biotechnol. 2013, 13: 61-

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  44. Barakat A, Bagniewska-Zadworna A, Frost CJ, Carlson JE: Phylogeny and expression profiling of CAD and CAD-like genes in hybrid Populus (P. Deltoides × P. Nigra): evidence from herbivore damage for subfunctionalization and functional divergence. BMC Plant Biol. 2010, 10: 100-

    Article  PubMed  PubMed Central  Google Scholar 

  45. Bhuiyan NH, Selvaraj G, Wei Y, King J: Gene expression profiling and silencing reveal that monolignol biosynthesis plays a critical role in penetration defense in wheat against powdery mildew invasion. J Exp Bot. 2009, 60 (2): 509-521.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  46. Coelho AC, Horta M, Neves D, Cravador A: Involvement of a cinnamyl alcohol dehydrogenase of Quercus suber in the defence response to infection by Phyutophthora cinnamomi. Physiol Mol Plant Pathol. 2006, 69: 62-72.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge Proteca Biotecnologia Florestal Company for kindly providing teak seeds. EG was recipient of a fellowship from “Coordenação de Aperfeiçoamento de Pessoal de Nível Superior” (CAPES), Brazil. TSV, DAR and HC were recipient of fellowships from “Fundação de Amparo à Pesquisa do Estado de São Paulo” (FAPESP).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Helaine Carrer.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

EG is the primary author of the manuscript. EG, TSV and DAR performed all the sample collection and preparation, experimental procedures, data analysis, qRT-PCR experiments and drafted the manuscript. VF and HC participated in the study design and provided helpful discussions and contributed with reagents, materials and analysis tools. EG, TSV, DAR and HC authors participated in the structuring and editing of the manuscript. HC coordinated the study, supervised the research and provided financial support. All authors read and approved the final manuscript.

Electronic supplementary material

Additional file 1:Information related to the orthologous plant sequences used in this study.(DOCX 24 KB)

13104_2014_2986_MOESM2_ESM.docx

Additional file 2:Clustal alignments used for designing primers to amplify orthologous sequences in teak. Red and green squares mean forward and reverse primers, respectively. (DOCX 200 KB)

Additional file 3:Protein clustal alignments for teak candidate reference genes.(DOCX 38 KB)

13104_2014_2986_MOESM4_ESM.docx

Additional file 4:Teak sequences (with accession numbers) used for designing qRT-PCR primers (Underlined and bolded).(DOCX 17 KB)

13104_2014_2986_MOESM5_ESM.docx

Additional file 5:Agarose gel (2%) electrophoresis showing amplification of a specific PCR product of the expected size for each gene. M represents 50 bp DNA ladder marker (GeneRuler™ 50 bp DNA Ladder, Thermo Scientific, USA) and “-” represents negative control. (DOCX 93 KB)

Additional file 6:Raw CP data used for statistical analysis in this study.(DOCX 21 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Galeano, E., Vasconcelos, T.S., Ramiro, D.A. et al. Identification and validation of quantitative real-time reverse transcription PCR reference genes for gene expression analysis in teak (Tectona grandis L.f.). BMC Res Notes 7, 464 (2014). https://0-doi-org.brum.beds.ac.uk/10.1186/1756-0500-7-464

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/1756-0500-7-464

Keywords