Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Selection and evaluation of new reference genes for RT-qPCR analysis in Epinephelus akaara based on transcriptome data

  • Huan Wang,

    Affiliations State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China, State-Province Joint Engineering Laboratory of Marine Bioproducts and Technology, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China

  • Xiang Zhang,

    Affiliations State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China, State-Province Joint Engineering Laboratory of Marine Bioproducts and Technology, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China

  • Qiaohong Liu,

    Affiliations State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China, State-Province Joint Engineering Laboratory of Marine Bioproducts and Technology, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China

  • Xiaochun Liu,

    Affiliation State Key Laboratory of Biocontrol, Institute of Aquatic Economic Animals and Guangdong Provincial Key Laboratory for Aquatic Economic Animals, School of Life Sciences, Sun Yat-Sen University, Guangzhou, Guangdong, China

  • Shaoxiong Ding

    sxding@xmu.edu.cn

    Affiliations State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China, State-Province Joint Engineering Laboratory of Marine Bioproducts and Technology, College of Ocean and Earth Sciences, Xiamen University, Xiamen, Fujian, China

Abstract

Groupers are an economically important fish species in world fishery markets. Because many studies using RT-qPCR have addressed gene expression in groupers, appropriate reference genes are required to obtain reliable and accurate results. In this study, the most suitable reference genes were identified from eleven candidate genes of one of the most valuable species, Epinephelus akaara, in a range of different experimental conditions. Using the software packages geNorm, NormFinder, BestKeeper and refFinder, three traditionally used reference genes, β-actin (β-ACT), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and beta-2-microglobulin (B2M), were identified as not suitable for E. akaara gene expression studies, whereas two newly identified reference genes, conserved oligomeric Golgi complex subunit 5 (Cog5) and brefeldin a-inhibited guanine nucleotide-exchange protein 1 (ARFGEF1), could be universally applied under all the tested conditions. These data provide the foundation for more precise results in RT-qPCR studies of gene expression in E. akaara and other Epinephelus species.

Introduction

Currently, real-time quantitative reverse transcription polymerase chain reaction (RT-qPCR) is a ubiquitous mainstay technology for quantifying gene expression because of its advantages, which include accuracy, specificity, sensitivity, reproducibility and convenience [1, 2]. Despite the numerous advantages, normalization to a constitutively expressed gene is required to avoid experimental errors caused by different sample amounts, variations in RNA, enzymatic efficiency for cDNA transcription and PCR efficiency [3, 4]. Traditionally, the reference genes that have been used are primarily housekeeping genes, such as β-actin (β-ACT), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), tubulin (TUB), elongation factor 1-α (EF1α), polyubiquitin (UBQ) and ribosomal RNAs (18S or 28S rRNA). However, many studies show that the transcript levels of these genes vary considerably across cellular conditions [59], which indicate that they should not be commonly used in all gene expression analysis situations. Therefore, the selection of a suitable reference gene to study specific objects or experimental conditions is an important factor for gene expression analysis by RT-qPCR. Recently, with the rapid development of new biotechnology, transcriptome sequencing has been applied to identify new reference genes in various species, including monkey [10], buckwheat [11], crazyweed [12], deer [13], rice [14], and oilseed rape [15], and the results show that the new reference genes are more stably expressed and that the orthologs of these genes are suitable reference genes in different species that have a close relationship [11]. Therefore, examination of transcriptome-sequencing data may be a useful method to identify more reliable reference genes for RT-qPCR studies on specific fauna or under different experimental conditions.

Grouper (Perciformes: Epinephelidae) have high economic value in Asian fishery markets, and grouper mariculture has rapidly developed due to the insufficient supply of wild caught fish, which does not satisfy the high demand of markets [16]. However, similar to other aquaculture species, the grouper aquaculture industry is also hindered by several problems that remain unresolved. One constraint of grouper culture is the production of fertilized eggs due to the lack of a standardized method for controlling sex change and the unavailability of mature male brood stock caused by the protogynous hermaphrodite mode of reproduction [17]. Larval production is another challenge for grouper farming because of the high level of mortality of embryos and larvae, which may be associated with biological aspects, nutritional deficiencies, genetic makeup of the parent population, variations in ecological factors and outbreaks of infectious disease [18, 19]. The variation in ecological factors and the diseases caused by pathogens are also large problems for farming grouper to commercial size. Salinity is an important ecological factor and not only influences the egg fertilization and incubation, early embryogenesis and larval growth but is also a key factor in controlling growth [20, 21]. Vibrio alginolyticus, a causative agent of outbreaks of vibriosis, is one of the most serious diseases and a common problem during various stages of grouper culture that causes fish mortality and serious economic losses [22, 23]. The natural disease caused by V. alginolyticus in fish includes symptoms of septicaemia, haemorrhage, dark skin, and ulcers on the skin surface in some cases [24]. In recent years, to contribute to the development of artificial grouper breeding and the protection of wild resources, an increasing number of investigations use RT-qPCR to examine the functional genes involved in the aspects discussed above. However, in the gonadal transcriptome data of an important grouper species, Epinephelus akaara, we found variable expression levels of some commonly used reference genes, such as β-ACT. Therefore, these genes may not be suitable as reference genes for analyses of gene expression in grouper and could lead to the erroneous normalization of RT-qPCR data. Thus, the identification of suitable reference genes in grouper studies is crucial for the accurate profiling of gene expression.

The objective of this study was to evaluate the appropriateness of traditional reference genes and new candidate reference genes derived from our previous transcriptome sequencing data of E. akaara and identify the most suitable reference genes for different experimental conditions. To the best of our knowledge, this report is the first assessment of valid reference genes for RT-qPCR studies on groupers. Because the Epinephelus species has the characteristics of a relatively short evolutionary time and closely related phylogeny [25], the identification of these new reference genes for the normalization of RT-qPCR data not only will be an essential tool for future gene expression studies in E. akaara but also may be applicable to other commercially important Epinephelus species, such as E. coioides, E. lanceolatus, Cromileptes altivelis and Anyperodon leucogrammicus.

Materials and methods

Sample preparation and experimental procedures

All E. akaara used in this study were obtained from a local commercial aquatic farm in Fujian Province, China. This study was conducted in strict accordance with the guidelines for the Care and Use of Laboratory Animals. The Institutional Animal Care and Use Committee of Xiamen University approved the protocol. All fish surgery procedures were conducted under MS-222 (Tricaine Methanesulfonate) to induce sedation and anaesthesia. All eggs, larvae, and tissue samples were incubated in an RNA fixer (Aidlab, Beijing, China) and stored at -80°C for RNA isolation.

All experiments were performed according to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines [26]. A detailed description of the sample preparation and experimental procedures, including information on sample collection under different experimental conditions, RNA isolation, cDNA synthesis, primer design, stand curve analysis, and RT-qPCR amplification conditions, is found in the S1 File.

Selection of the candidate reference genes

For the selection of new candidate reference genes, we analysed our previous gonadal transcriptome sequencing data (PRJNA277894) from the following five different gonadal development phases of E. akaara: undifferentiated-phase (UN), developing female (DF), developing male (DM), mature female (MF) and mature male (MM) [27]. To estimate the expression stability of each gene, we analysed the raw data for all genes as follows: (1) sorting the genes expressed in all of the tested tissues; (2) calculating the mean expression values (MVs), standard deviations (SDs) and coefficients of variation (CVs); (3) screening the genes that satisfied the conditions described by de Jonge et al [28]; and (4) ordering the gene list according to the lowest CV value.

Based on these selection procedures for the transcriptome sequencing data, 10 genes that had a minor variation in expression were selected. Of these genes, 8 were newly identified candidate reference genes, including brefeldin a-inhibited guanine nucleotide-exchange protein 1 (ARFGEF1), conserved oligomeric Golgi complex subunit 5 (Cog5), putative ATP-dependent RNA helicase DHX30 (DHX30), neuron navigator 3 (Nav3), homeodomain-interacting protein kinase 3 (Hipk3), probable E3 ubiquitin-protein ligase MYCBP2 (Mycbp2), E3 ubiquitin-protein ligase MGRN1 (Mgrn1), and 60S ribosomal protein L17 (RPL17), and two were traditionally used reference genes, GAPDH and beta-2-microglobulin (B2M). The most commonly used reference gene, β-ACT, in grouper was also evaluated in our study.

Primer specificity and RT-qPCR efficiency analysis

We used full-length unigene sequences from gonadal transcriptome data to design specific primers for RT-qPCR. The RT-qPCR products ranged from 86 to 243 bp. The specificities of all the primers were demonstrated by the single bands of expected size in agarose gel electrophoresis and by the single-peak melting curves of the RT-qPCR products (S1 Fig). The Sanger sequencing of amplicons with alignment confirmed the specificity and accuracy of the designed primers. The Cq values of all the primers were greater than 35 in the non-template controls (NTC). Therefore, minute amounts of primer dimers did not affect the fluorescence level of the amplified target gene. The amplifications detected in non-reverse transcription controls (NRT) were similar to the NTCs, which indicated that there was no genomic DNA contamination. The E-values of the eleven genes in this study were between 91% and 99%, which were within the acceptable range. The primer sequences and relevant amplification information are presented in Table 1.

thumbnail
Table 1. Eleven selected candidate reference genes, gene function, primers and different parameters derived from RT-qPCR data analysis.

https://doi.org/10.1371/journal.pone.0171646.t001

Gene stability analysis

The PCR data were analysed using three approaches, geNorm [4], NormFinder [29] and BestKeeper [30]. For all three programs, Ct values were transformed into Cq values by the standard formula Cq = Log (2)/Log (E), where E is the efficiency of the amplification of each primer pair. Then, for both geNorm and NormFinder, the Cq values were transformed to relative quantities using the delta-Cq formula Q = 2 delta-Cq, where delta-Cq is the difference between the sample with the lowest Cq and the Cq value of the sample in question. The geNorm program calculates gene expression stability (M-value), which is based on the pairwise variation of a single reference candidate gene relative to all other tested control genes. Genes with an M-value below the threshold of 1.5 were those stably expressed and lower M-values indicated greater stability of a gene. The geNorm program also estimates the pairwise variation (Vn/Vn+1) between two sequential normalization factors, NFn and NFn+1, to determine the optimal number of reference genes required. The optimal number of reference genes is necessary to satisfy Vn/Vn+1 below the cut-off value of 0.15. NormFinder uses an ANOVA-based model to identify the most stable reference genes. NormFinder calculates intra- and inter-group variations, and genes with the least variations are those that the most stable. BestKeeper calculates the SD and CV based on the Cq values of each reference gene. Genes with an SD value of < 1.0 are stable, and the gene with the lowest SD and CV values was identified as the most stable. Additionally, refFinder, a user-friendly comprehensive tool, was used to integrate the results from the three programs and select the best reference gene.

Validation of the recommended reference genes

To validate the stability of the identified reference genes, the expression levels of two target genes were analysed in different experimental conditions. A non-mammalian myostatin, growth and differentiation factor-8 (MSTN), was used to validate the expression levels of the most stable, the universal, and the least stable reference genes for E. akaara during gonad development and early ontogenetic development phases. An important immune responsive cytokine, interleukin (IL)-1β, was used to validate the expression levels of the most stable, the universal, and the least stable reference genes for salinity treatment and V. alginolyticus challenged E. akaara. Relative quantification of these two genes in different samples was performed using the 2-ΔΔCt method [31].

Results

Expression stability of candidate reference genes during different gonad development phases of E. akaara

As shown in S1 Table and Fig 1A, the Cq values of the eleven candidate reference genes during E. akaara gonad development ranged between 16.0 and 30.5. RPL17 was the most highly expressed gene, with Cq values ranging between 16.0 and 17.7, followed by β-ACT and GAPDH, which exhibited Cq values from 17.9 to 21.7 and from 17.1 to 21.8, respectively. The remaining genes had Cq values greater than 22, with the highest Cq value (30.5) exhibited by Hipk3. For expression stability analysis, the results from the geNorm, NormFinder, and BestKeeper programs were combined, and seven genes, Hipk3, ARFGEF1, B2M, Mycbp2, RPL17, Nav3, and Cog5, were stably expressed during different gonad development phases of E. akaara (Fig 2A, 2B, and 2C). The overall ranking order of stability of candidate reference genes determined by refFinder was Cog5>Nav3>RPL17>Mycbp2>B2M>ARFGEF1>Hipk3>DHX30>GAPDH>Mgrn1>β-ACT (Fig 2D). The results of the pairwise variation calculation showed that the V2/3 value was less than 0.15, which indicated that two reference genes were required for normalization (Fig 3A). Thus, we recommend Cog5/Nav3 as the best reference gene pair for the normalization of E. akaara gonad development.

thumbnail
Fig 1. Expression levels of the 11 candidate reference genes.

Graphical representation of absolute Cq values for each gene analysed during different gonad development phases (A), during early ontogenetic development stages (B), for the salinity treatment of E.akaara (C), and for V.alginolytics challenged E.akaara (D). Values are Cq values. A line across the box depicts the median. The box indicates the 25th and 75th percentiles. Whiskers represent the maximum and minimum values.

https://doi.org/10.1371/journal.pone.0171646.g001

thumbnail
Fig 2. Determination of the expression stability during gonad development phases according to different programs.

(A) geNorm, the red line indicates the geNorm cut-off value of 1.5. (B) NormFinder. (C) BestKeeper, the red line indicates the BestKeeper cut-off value of 1.0. (D) refFinder, the red indicates the gene expression was unstable.

https://doi.org/10.1371/journal.pone.0171646.g002

thumbnail
Fig 3. Determination of the optimal number of reference genes.

Pairwise variations (V) of the 11 candidate reference genes were calculated by geNorm to determine the optimal number of reference genes for the accurate normalization during the gonad development phases, during the entire early ontogenetic development stage and in the separate embryonic ontogenesis stage and larval ontogenesis stage, and for salinity treatment samples (A) and in different tissues after V. alginolyticus challenge (B). The black line indicates the geNorm cut-off value of 0.15.

https://doi.org/10.1371/journal.pone.0171646.g003

Expression stability of candidate reference genes during E. akaara early ontogenetic development stage

When the entire early ontogenetic development stage (from two-cell to 32 days after hatching (dah)) was considered, the Cq values of all eleven genes ranged from 18.2 to 36.4 (S1 Table and Fig 1B). Unlike highly expressed RPL17 (Cq ranging from 18.2 to 29.2) and weakly expressed Hipk3 (Cq ranging from 29.4 to 35.4), the other candidate reference genes were expressed at moderate levels. Combining the results from geNorm, NormFinder, and BestKeeper programs, Nav3 and Mgrn1 were stably expressed genes during the entire early ontogenetic development stage of E. akaara (Fig 4A, 4B and 4C). The comprehensive ranking of stability according to the integrative refFinder was Mycbp2>ARFGEF1>Cog5>Mgrn1>DHX30>Nav3>β-ACT>Hipk3>B2M>RPL17>GAPDH (Fig 4D). However, the pairwise variation showed that all Vn/Vn+1 were larger than 0.15, which indicated that no optimal number of reference genes was determined (Fig 3A). Therefore, no ideal choice for normalization during the E. akaara early ontogenetic development stage could be determined, although Nav3/Mgrn1 might be the best option.

thumbnail
Fig 4. Determination of the expression stability during the entire early ontogenetic development stage and in the separate embryonic ontogenesis stage and larval ontogenesis stage according to different programs.

(A) geNorm, the red line indicates the geNorm cut-off value of 1.5. (B) NormFinder. (C) BestKeeper, the red line indicates the BestKeeper cut-off value of 1.0. (D) refFinder, the red indicates the gene expression was unstable. G: GAPDH, R: RPL17, H: Hipk3, D: DHX30, A: ARFGEF1, β: β-ACT, N: Nav3, C: Cog5, Mg: Mgrn1, My: Mycbp2, B: B2M.

https://doi.org/10.1371/journal.pone.0171646.g004

When the early ontogenetic development stage was divided into the embryonic ontogenesis stage (from two-cell to newly hatched larvae) and larval ontogenesis stage (from first feeding larvae to 32 dah), consistent with other reports [32, 33], better results were obtained. According to the analysis of geNorm, NormFinder and BestKeeper, Nav3, Mgrn1, Cog5, Mycbp2, β-ACT and ARFGEF1 were stably expressed during the embryonic ontogenesis stage, whereas Cog5, RPL17, Nav3, B2M, GAPDH, β-ACT, ARFGEF1, DHX30 and Mgrn1 were stably expressed during the larval ontogenesis stage (Fig 4A, 4B, and 4C). According to refFinder, the recommended comprehensive ranking of stability was Nav3>Cog5>Mgrn1>Mycbp2> ARFGEF1>β-ACT>DHX30>B2M>Hipk3>GAPDH>RPL17 during the embryonic ontogenesis stage and Cog5>RPL17>Nav3>B2M>ARFGEF1>GAPDH>Mgrn1> β-ACT>Mycbp2>DHX30>Hipk3 during the larval ontogenesis stage (Fig 4D). For the number of genes required, V3/4 and V2/3 were less than 0.15 during the embryonic ontogenesis stage and larval ontogenesis stage, respectively, suggesting that three and two reference genes were sufficient for normalization, respectively (Fig 3A). Therefore, we recommend Nav3/Cog5/Mgrn1 and RPL17/Cog5 as the optimum suitable reference genes for the embryonic ontogenesis stage and the larval ontogenesis stage, respectively.

Expression stability of candidate reference genes for salinity treatment of E. akaara

For the salinity treatment group, the Cq values of the eleven genes ranged between 17.3 and 33.4. RPL17 and β-ACT were the most highly expressed genes, with Cq values ranging from 17.3 to 18.6 and 18.3 to 19.1, respectively, whereas DHX30 and Hipk3 were the least expressed, with Cq values ranging from 30.8 to 33.2 and 29.1 to 33.4, respectively (S1 Table and Fig 1C). Because the Cq values of the eleven candidate reference genes between the head kidney and the gills did not differ significantly under salinity treatment, we combined the two tissues to evaluate the gene expression stability. The integrative analysis of the geNorm, NormFinder, and BestKeeper programs showed that most of the candidate reference genes were stable, except for Hipk3 and GAPDH (Fig 5A, 5B and 5C). The refFinder recommended ranking of stability was Nav3>Mgrn1>RPL17>Mycbp2>β-ACT>ARFGEF1>B2M>DHX30>Cog5>Hipk3>GAPDH (Fig 5D). The geNorm analysis showed that the V3/4 value was less than 0.15; therefore, three reference genes were required (Fig 3A). In conclusion, we recommend Nav3/Mgrn1/ RPL17 as the most suitable reference genes for E. akaara under salinity treatment.

thumbnail
Fig 5. Determination of the expression stability for salinity treatment samples according to different programs.

(A) geNorm, the red line indicates the geNorm cut-off value of 1.5. (B) NormFinder. (C) BestKeeper, the red line indicates the BestKeeper cut-off value of 1.0. (D) refFinder, the red indicates the gene expression was unstable.

https://doi.org/10.1371/journal.pone.0171646.g005

Expression stability of candidate reference genes for Vibrio alginolyticus challenged E. akaara

As shown in S1 Table and Fig 1D, the Cq values of the eleven genes ranged between 17.2 and 36.5. The most highly expressed gene was RPL17 in the liver, intestine, stomach and spleen; β-ACT in the head kidney, heart and gills; and GAPDH in white muscle. The least expressed gene was Mgrn1 in the liver, intestine, spleen, heart, head kidney, gills and white muscle and DHX30 in the stomach. For the optimal number of genes required for data normalization, the V2/3 values were lower than 0.15 in the intestine, heart, and gills, whereas the V3/4 values were lower than 0.15 in the liver, stomach, spleen, head kidney, and white muscle; therefore, two and three reference genes were required for normalization, respectively (Fig 3B). The detailed stability ranking according to the integrative refFinder is shown in Fig 6. Based on these results, the most suitable reference genes for E. akaara were RPL17/Mgrn1/B2M in liver, ARFGEF1/Hipk3 in intestine, GAPDH/ARFGEF1/RPL17 in stomach, ARFGEF1/Nav3/Cog5 in spleen, Mycbp2/DHX30 in heart, Mgrn1/ARFGEF1/β-ACT in head kidney, ARFGEF1/Nav3 in gills, and DHX30/Mycbp2/RPL17 in white muscle,. Although more genes were stably expressed in each tissue, the genes that could be universally used in all the tested tissues were ARFGEF1, β-ACT and Cog5. Considering the complexity and susceptibility of the experimental studies and the V values for each tissue, three reference genes were most appropriate for the normalization of different tissues following exposure to a pathogen. Therefore, we recommend ARFGEF1/β-ACT/Cog5 as the most suitable reference genes for all the tested tissues of E. akaara challenged by V. alginolyticus.

thumbnail
Fig 6. Geomean of ranking values of the 11 candidate reference genes calculated by refFinder for V. alginolyticus challenged E. akaara.

The red diamond indicates the gene expression was unstable as determined by geNorm, NormFinder, or BestKeeper. Ranking of the gene expression stability was performed for the (A) liver, (B) intestine, (C) stomach, (D) spleen, (E) heart, (F) head kidney, (G) gills, and (H) white muscle.

https://doi.org/10.1371/journal.pone.0171646.g006

The universality of candidate reference genes for E. akaara in all tested conditions

A universally used reference gene would be popular because of its convenience. To find a reference gene that could be universally applied for E. akaara, we listed the stability status and evaluated the consensus ranking performance of eleven candidate reference genes for all tested conditions. In Table 2, the candidate reference genes are divided into two groups: one group with stable gene expression, and one group with unstable expression. The stably expressed genes were screened according to geNorm, NormFinder and BestKeeper, which must simultaneously satisfy the stability thresholds of these programs. The consensus ranking of each gene was calculated from the geometric mean of the refFinder ranking under all tested conditions. As shown in Fig 7, ARFGEF1, Nav3, Mgrn1, and Cog5 were the top four genes, whereas in Table 2, only Cog5 and ARFGEF1 were listed in the column of stably expressed genes among all the experimental conditions of E. akaara. Therefore, we recommend Cog5 and ARFGEF1 as the optimal choices for universal reference genes in E. akaara.

thumbnail
Fig 7. Consensus ranking of the 11 candidate reference genes.

https://doi.org/10.1371/journal.pone.0171646.g007

thumbnail
Table 2. Expression stability ranking of the 11 candidate reference genes.

https://doi.org/10.1371/journal.pone.0171646.t002

Validation of the recommended reference genes

To examine the availability of the identified most stable reference genes, recommended universal reference genes, and the least stable reference genes, two extensively studied genes, MSTN and IL-1β, were used as the target genes in different experimental conditions. The relative quantification of these two genes varied according to the different reference genes used during normalization (Fig 8).

thumbnail
Fig 8. Validation of the recommended reference genes.

Expression profiles of MSTN and IL-1β were analysed using the most stable, the recommended universal and the least stable reference genes during different gonad development phases (A), during the early ontogenetic development stage (B), for the salinity treatment of E. akaara (C), and for V. alginolyticus challenged E. akaara (D). Bars represent the means and standard error of three biological replications. HS: high salinity, LS: low salinity, HK: head kidney, GS: gills.

https://doi.org/10.1371/journal.pone.0171646.g008

Comparing the relative expression levels of MSTN and IL-1β revealed expression differences of the target genes when normalized with the most stable, the recommended universal and the least stable reference genes. Compared with the least stable genes, the most stable genes and the recommended universal genes showed clearly similar levels of gene expression. Therefore, the recommended universal genes could be used for normalization of target genes under all experimental conditions.

Discussion

Normalization is essential to obtain accurate and reliable quantitative gene expression data from RT-qPCR analyses. Therefore, the systematic validation of reference genes for different experimental conditions is necessary [34]. During the past few years, the number of articles reporting the validation of reference genes in fishes has increased [3540], but most studies explore the stability of traditional reference genes such as β-ACT, GAPDH, B2M and EF1α. However, transcriptome-wide surveys of gene expression stability in other species, such as monkey [10], buckwheat [11], and oilseed rape [15], show that many newly identified reference genes are more stably expressed. In this study, large-scale gonadal transcriptome data from E. akaara were used as the source for reference gene selection. The results of RT-qPCR analysis of the different gonad development phases were consistent with the transcriptome sequencing, and the newly identified candidate reference genes were more stably expressed than the traditional genes in some situations, which demonstrated that transcriptome data are an accurate and useful source for candidate reference gene screening and represented an important strategy for large-scale reference gene selection for non-model organisms.

Many statistical methods are used to determine the stability of gene expression and assess the most suitable reference genes [11, 41, 42]; however, no consensus agreement has emerged on which approach is optimum. To improve the reliability and accuracy of the assessment, the combination of different statistical algorithms is the most reliable method. In this study, geNorm, NormFinder and BestKeeper, the most frequently used statistical programs, were used to assess the expression variation of candidate reference genes. We found that the results generated from these three software applications were very similar, despite some nuances in the stability ranking order that might be caused by the different algorithms. Because of the subtle deviations among software programs, refFinder, which compares and ranks the tested candidate reference genes based on the rankings from each program, was used to integrate the three methods.

To avoid erroneous data that may be obtained using a single reference gene, normalization with multiple reference genes is becoming the most common technique [39, 43]. However, because available reference genes are limited, most previous studies typically apply a single reference gene. For the optimal number of reference genes required, the V value calculated by geNorm can provide a reference. Based on our results, the V2/3 or V3/4 value was lower than the threshold 0.15 of geNorm under all the different experimental conditions, which indicated that 2 or 3 reference genes were required for normalization. Considering the complexity, susceptibility, economic cost and workload of an experimental study, the number of reference genes for normalization is appropriate and practical, and such a strategy has been widely used in many studies [32].

In previous studies on gene expression using the RT-qPCR method, including investigations on grouper functional genes, β-ACT, GAPDH, and B2M are commonly used as internal controls. However, the stability of expression of these genes in grouper remains unclear. In our study, although β-ACT is the most frequently used reference gene in grouper [4447], this gene was not the most suitable choice in most situations, and different results from other studies have made this gene a controversial reference gene. β-ACT was a stable gene in 12 tissues of the half-smooth tongue sole following LPS or bacterial challenge [48], but the expression varied to different extents in turbot tissues before and after bacterial challenge [49], during Kuruma shrimp growth and developmental stages [50], and under environmental stresses in Ciona savignyi [32]. The different performances of β-ACT may be due to its multiple functions not only as a ubiquitous cytoskeleton protein involved in cell structure and motility [1] but also because it participates in diverse biological processes such as transcription and stress response [51]. Although GAPDH and B2M were identified as stable reference genes by the gonadal transcriptome data, these genes were unusable under certain situations in the present study. For example, GAPDH was identified but could not be used in the gonad development phase, embryonic ontogenesis stage, for the salinity treatment, and in many tissues after V. alginolyticus challenge, and B2M could not be used in the embryonic ontogenesis stage and most tissues after V. alginolyticus challenge. The instability of GAPDH has been reported in numerous other systems, including different developmental stages of Atlantic halibut [35], embryonic development in zebra fish [38], virus-infected salmon [52], and different tissue types after bacterial infection in Japanese flounder [39]. The performance of GAPDH may be due to its diverse range of functions in the glycolysis cycle, nuclear RNA export, DNA replication and repair, protein phosphotransferase/kinase reactions, membrane transport and fusion, translational regulation and phosphotransferase activity [6], making the expression sensitive to many perturbations in cellular homeostasis [53]. B2M was identified as the least stable gene across all tested tissue types after V. alginolyticus challenge. Consistent with our results, several studies report considerable variation of B2M across tissues, including in turbot and flounder when infected by megalocytivirus and Edwardsiella tarda [39, 43, 49]. The instability of B2M was induced by infection, which may reflect its role in antigen binding and presentation by the major histocompatibility complex (MHC), hence its limitation as a reference gene in such studies [54]. In the current study, none of these three commonly used housekeeping genes showed high stability, indicating that their expression levels were influenced by both the developmental process and the immune response; thus, these housekeeping genes unsuitable as internal controls.

Clearly, a universal reference gene would be more acceptable to researchers because of its convenience. Of the eleven candidate reference genes, based on the integration of the stability results across all tested conditions, only the stability parameters of Cog5 and ARFGEF1 satisfied the stability threshold for all the experimental conditions. Although Cog5 was a newly identified reference gene in fish, previous studies of human non-small cell lung cancer and thyroid cancer used this gene as an internal control [5557]. Furthermore, in a Pacific oyster microarray, Cog5 was also identified as a candidate reference gene because its mRNA level appeared stable across all analysed tissues [58]. The stability of Cog5 was closely tied to its function, which is based on the conserved oligomeric Golgi (COG) complex subunit that is involved in the Golgi-associated membrane trafficking and glycoconjugate synthesis [59]. The other gene, ADP-ribosylation factor, ARFGEF1, is ubiquitous in eukaryotic cells and functions as a regulator of vesicular traffic and actin remodelling [60]. Previous studies reveal that ARFGEF1 is one of the most stable reference genes in plants and animals, including wheat, barley, rye, citrus, and beetle [6164]. Moreover, the combination of Cog5 and ARFGEF1 showed similar normalization results with the most stable genes when the relative expression of MSTN and IL-1β was analysed. Therefore, we recommend that Cog5 and ARFGEF1 be used as reference genes for E. akaara under all tested conditions. However, for a specific experimental condition, the most suitable reference genes are strongly recommended because of the better stability and applicability.

Conclusions

In summary, the primary conclusions of this study were as follows: (1) the three traditionally used reference genes, β-ACT, GAPDH and B2M, were not applicable for E. akaara gene expression studies in most situations; (2) the two newly identified reference genes, Cog5 and ARFGEF1, could be universally applied in E. akaara under all tested conditions; and (3) the most suitable reference genes were identified for each specific experimental conditions of E. akaara. Additionally, our results might be applicable to other groupers.

Supporting information

S1 Fig. Melting curve analyses of 11 candidate reference genes using RT-qPCR.

https://doi.org/10.1371/journal.pone.0171646.s001

(TIF)

S2 Fig. Expression stability values (M) of the 11 candidate reference genes calculated by geNorm for V. alginolyticus challenged E. akaara.

Low M values indicate more stable expression. The least stable genes are on the left, and the most stable genes are on the right. The red line indicates the geNorm cut-off value of 1.5. Ranking of the gene expression stability was performed for the (A) liver, (B) intestine, (C) stomach, (D) spleen, (E) heart, (F) head kidney, (G) gills, and (H) white muscle.

https://doi.org/10.1371/journal.pone.0171646.s002

(TIF)

S3 Fig. Expression stability V values of the 11 candidate reference genes calculated by NormFinder for V. alginolyticus challenged E. akaara.

Low V values indicate more stable expression. The least stable genes are on the left, and the most stable genes are on the right. Ranking of the gene expression stability was performed for the (A) liver, (B) intestine, (C) stomach, (D) spleen, (E) heart, (F) head kidney, (G) gills, and (H) white muscle.

https://doi.org/10.1371/journal.pone.0171646.s003

(TIF)

S4 Fig. Expression stability values (SD) of the 11 candidate reference genes calculated by BestKeeper for V. alginolyticus challenged E. akaara.

Low SD values indicate more stable expression. The least stable genes are on the left, and the most stable genes are on the right. The red line indicates the BestKeeper cut-off value of 1.0. Ranking of the gene expression stability was performed for the (A) liver, (B) intestine, (C) stomach, (D) spleen, (E) heart, (F) head kidney, (G) gills, and (H) white muscle.

https://doi.org/10.1371/journal.pone.0171646.s004

(TIF)

S1 File. Sample preparation and experimental procedures.

https://doi.org/10.1371/journal.pone.0171646.s005

(DOCX)

S1 Table. Average Cq values ± SD of the 11 candidate reference genes for all the different experimental samples.

https://doi.org/10.1371/journal.pone.0171646.s006

(DOCX)

Acknowledgments

We wish to thank Prof. Michael Masao Miyamoto for his valuable comments and advice on the manuscript.

Author Contributions

  1. Conceptualization: HW XZ SD.
  2. Formal analysis: HW XZ.
  3. Funding acquisition: XL SD.
  4. Investigation: HW QL.
  5. Resources: HW QL.
  6. Validation: HW XZ QL.
  7. Visualization: HW SD.
  8. Writing – original draft: HW.
  9. Writing – review & editing: HW SD.

References

  1. 1. Bustin SA. Absolute quantification of mRNA using real-time reverse transcription polymerase chain reaction assays. Journal of molecular endocrinology. 2000;25(2):169–193. pmid:11013345
  2. 2. Bustin SA, Benes V, Nolan T, Pfaffl MW. Quantitative real-time RT-PCR-a perspective. Journal of molecular endocrinology. 2005;34(3):597–601. pmid:15956331
  3. 3. Huggett J, Dheda K, Bustin S, Zumla A. Real-time RT-PCR normalisation; strategies and considerations. Genes and Immunity. 2005;6(4):279–284. pmid:15815687
  4. 4. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome biology. 2002;3(7):1–12.
  5. 5. Yperman J, De Visscher G, Holvoet P, Flameng W. Beta-actin cannot be used as a control for gene expression in ovine interstitial cells derived from heart valves. Journal of Heart Valve Disease. 2004;13(5):848–853. pmid:15473489
  6. 6. Bustin SA. Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems. Journal of molecular endocrinology. 2002;29:23–39. pmid:12200227
  7. 7. Ruan W, Lai M. Actin, a reliable marker of internal control? Clinica Chimica Acta. 2007;385(1–2):1–5.
  8. 8. Thorrez L, Van Deun K, Tranchevent L C, Van Lommel L, Engelen K, Marchal K, et al. Using ribosomal protein genes as reference: a tale of caution. PloS one. 2008;3(3):e1854. pmid:18365009
  9. 9. Glare E M, Divjak M, Bailey M J, Walters E H. Beta-actin and GAPDH housekeeping gene expression in asthmatic airways is variable and not suitable for normalising mRNA levels. Thorax. 2002;57(9):765–770. pmid:12200519
  10. 10. Park SJ, Kim YH, Huh JW, Lee SR, Kim SH, Kim SU, et al. Selection of new appropriate reference genes for RT-qPCR analysis via transcriptome sequencing of cynomolgus monkeys (Macaca fascicularis). PloS one. 2013;8(4):e60758. PubMed Central PMCID: PMC3626658. pmid:23613744
  11. 11. Demidenko NV, Logacheva MD, Penin AA. Selection and validation of reference genes for quantitative real-time PCR in buckwheat (Fagopyrum esculentum) based on transcriptome sequence data. PloS one. 2011;6(5):e19434. PubMed Central PMCID: PMC3093374. pmid:21589908
  12. 12. Zhuang H, Fu Y, He W, Wang L, Wei Y. Selection of appropriate reference genes for quantitative real-time PCR in Oxytropis ochrocephala Bunge using transcriptome datasets under abiotic stress treatments. Frontiers in plant science. 2015;6. PubMed Central PMCID: PMC4484982.
  13. 13. Liu M, Yao B, Zhang H, Guo H, Hu D, Wang Q, et al. Identification of novel reference genes using sika deer antler transcriptome expression data and their validation for quantitative gene expression analysis. Genes & Genomics. 2014;36(5):573–582.
  14. 14. Narsai R, Ivanova A, Ng S, Whelan J. Defining reference genes in Oryza sativa using organ, development, biotic and abiotic transcriptome datasets. BMC plant biology. 2010;10:56. PubMed Central PMCID: PMC2923530. pmid:20353606
  15. 15. Yang H, Liu J, Huang S, Guo T, Deng L, Hua W. Selection and evaluation of novel reference genes for quantitative reverse transcription PCR (qRT-PCR) based on genome and transcriptome data in Brassica napus L. Gene. 2014;538(1):113–122. pmid:24406618
  16. 16. Pierre S, Gaillard S, Prévot-D'Alvise N, Aubert J, Rostaing-Capaillon O, Leung-Tack D, et al. Grouper aquaculture: Asian success and Mediterranean trials. Aquatic Conservation: Marine and Freshwater Ecosystems. 2008;18(3):297–308.
  17. 17. Huntsman GR, Schaaf WE. Simulation of the Impact of Fishing on Reproduction of a Protogynous Grouper, the Graysby. North American Journal of Fisheries Management. 1994;14(1):41–52.
  18. 18. Lim LC. Larviculture of the Greasy Grouper Epinephelus tauvina F. and the Brown‐Marbled Grouper E. fuscoguttatus F. in Singapore. Journal of the World Aquaculture Society. 1993;24(2):262–274.
  19. 19. Kohno H, Ordonio-Aguilar R S, Ohno A, Taki Y. Why is grouper larval rearing difficult?: an approach from the development of the feeding apparatus in early stage larvae of the grouper, Epinephelus coioides. Ichthyological Research. 1997;44(2–3):267–274.
  20. 20. Boeuf G, Payan P. How should salinity influence fish growth? Comparative Biochemistry and Physiology Part C:Toxicology & Pharmacology. 2001;130(4):411–423.
  21. 21. Gracia-López V, Kiewek-Martínez M, Maldonado-García M. Effects of temperature and salinity on artificially reproduced eggs and larvae of the leopard grouper Mycteroperca rosacea. Aquaculture. 2004;237(1–4):485–498.
  22. 22. Yin Z-X, He W, Chen W-J, Yan J-H, Yang J-N, Chan S-M, et al. Cloning, expression and antimicrobial activity of an antimicrobial peptide, epinecidin-1, from the orange-spotted grouper, Epinephelus coioides. Aquaculture. 2006;253(1–4):204–211.
  23. 23. Rasheed VM. Diseases of Cultured Brown-Spotted Grouper Epinephelus tauvina and Silvery Black Porgy Acanthopagrus cuvieriin Kuwait. Journal of Aquatic Animal Health. 1989;1(2):102–107.
  24. 24. Egidius E. Vibriosis: pathogenicity and pathology. A review. Aquaculture. 1987;67(1):15–28.
  25. 25. Ma KY, Craig MT, Choat JH, van Herwerden L. The historical biogeography of groupers: Clade diversification patterns and processes. Molecular phylogenetics and evolution. 2016;100:21–30. pmid:26908372
  26. 26. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clinical chemistry. 2009;55(4):611–622. pmid:19246619
  27. 27. Liu M, Wang YY, Shan XJ, Kang B, Ding SX. Primary male development of two sequentially hermaphroditic groupers, Epinephelus akaara and Epinephelus awoara (Perciformes: Epinephelidae). Journal of fish biology. 2016;88:1598–1613. pmid:26935897
  28. 28. De Jonge H J, Fehrmann R S, de Bont E S, Hofstra R M, Gerbens F, Kamps W A, et al. Evidence based selection of housekeeping genes. PloS one. 2007;2(9):e898. pmid:17878933
  29. 29. Andersen C L, Jensen J L, Ørntoft T F. Normalization of Real-Time Quantitative Reverse Transcription-PCR Data: A Model-Based Variance Estimation Approach to Identify Genes Suited for Normalization, Applied to Bladder and Colon Cancer Data Sets. Cancer research. 2004;64(15):5245–5250. pmid:15289330
  30. 30. Pfaffl M W, Tichopad A, Prgomet C, Neuvians T P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper-Excel-based tool using pair-wise correlations. Biotechnology letters. 2004;26(6):509–515. pmid:15127793
  31. 31. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–408. pmid:11846609
  32. 32. Huang X, Gao Y, Jiang B, Zhou Z, Zhan A. Reference gene selection for quantitative gene expression studies during biological invasions: A test on multiple genes and tissues in a model ascidian Ciona savignyi. Gene. 2016;576(1 Pt 1):79–87.
  33. 33. Sun M, Lu MX, Tang XT, Du YZ. Exploring valid reference genes for quantitative real-time PCR analysis in Sesamia inferens (Lepidoptera: Noctuidae). PloS one. 2015;10(1):e0115979. PubMed Central PMCID: PMC4293147. pmid:25585250
  34. 34. Nolan T, Hands RE, Bustin SA. Quantification of mRNA using real-time RT-PCR. Nature protocols. 2006;1(3):1559–1582. pmid:17406449
  35. 35. Fernandes JM, Mommens M, Hagen O, Babiak I, Solberg C. Selection of suitable reference genes for real-time PCR studies of Atlantic halibut development. Comparative biochemistry and physiology Part B, Biochemistry & molecular biology. 2008;150(1):23–32.
  36. 36. Olsvik PA, Lie KK, Jordal AE, Nilsen TO, Hordvik I. Evaluation of potential reference genes in real-time RT-PCR studies of Atlantic salmon. BMC molecular biology. 2005;6:21. PubMed Central PMCID: PMC1314898. pmid:16293192
  37. 37. Overgard AC, Nerland AH, Patel S. Evaluation of potential reference genes for real time RT-PCR studies in Atlantic halibut (Hippoglossus Hippoglossus L.); during development, in tissues of healthy and NNV-injected fish, and in anterior kidney leucocytes. BMC molecular biology. 2010;11:36. PubMed Central PMCID: PMC2882370. pmid:20459764
  38. 38. Tang R, Dodd A, Lai D, McNabb WC, Love DR. Validation of Zebrafish (Danio rerio) Reference Genes for Quantitative Real-time RT-PCR Normalization. Acta Biochimica et Biophysica Sinica. 2007;39(5):384–390. pmid:17492136
  39. 39. Zheng WJ, Sun L. Evaluation of housekeeping genes as references for quantitative real time RT-PCR analysis of gene expression in Japanese flounder (Paralichthys olivaceus). Fish & shellfish immunology. 2011;30(2):638–645.
  40. 40. Purohit GK, Mahanty A, Mohanty BP, Mohanty S. Evaluation of housekeeping genes as references for quantitative real-time PCR analysis of gene expression in the murrel Channa striatus under high-temperature stress. Fish physiology and biochemistry. 2015.
  41. 41. Garcia-Vallejo JJ, Hof BV, Robben J, Van Wijk JA, Van Die I, Joziasse DH, et al. Approach for defining endogenous reference genes in gene expression experiments. Analytical biochemistry. 2004;329(2):293–299. pmid:15158490
  42. 42. Szabo A, Perou CM, Karaca M, Perreard L, Quackenbush JF, Bernard PS. Statistical modeling for selecting housekeeper genes. Genome biology. 2004;5(8):R59. pmid:15287981
  43. 43. Zhang J, Hu YH, Sun BG, Xiao ZZ, Sun L. Selection of normalization factors for quantitative real time RT-PCR studies in Japanese flounder (Paralichthys olivaceus) and turbot (Scophthalmus maximus) under conditions of viral infection. Veterinary immunology and immunopathology. 2013;152(3–4):303–316. pmid:23332581
  44. 44. Luo Y-S, Hu W, Liu X-C, Lin H-R, Zhu Z-Y. Molecular cloning and mRNA expression pattern of Sox9 during sex reversal in orange-spotted grouper (Epinephelus coioides). Aquaculture. 2010;306(1):322–328.
  45. 45. Huang W, Zhou L, Li Z, Gui J-F. Expression pattern, cellular localization and promoter activity analysis of ovarian aromatase (Cyp19a1a) in protogynous hermaphrodite red-spotted grouper. Mol Cell Endocrinol. 2009;307(1):224–236.
  46. 46. Liu Q, Lu H, Zhang L, Xie J, Shen W, Zhang W. Homologues of sox8 and sox10 in the orange-spotted grouper Epinephelus coioides: sequences, expression patterns, and their effects on cyp19a1a promoter activities in vitro. Comparative biochemistry and physiology Part B, Biochemistry & molecular biology. 2012;163(1):86–95.
  47. 47. Tang X, Liu X, Zhang Y, Zhu P, Lin H. Molecular cloning, tissue distribution and expression profiles of thyroid hormone receptors during embryogenesis in orange-spotted grouper (Epinephelus coioides). Gen Comp Endocrinol. 2008;159(2–3):117–124. pmid:18805420
  48. 48. Li Z, Yang L, Wang J, Shi W, Pawar RA, Liu Y, et al. Beta-Actin is a useful internal control for tissue-specific gene expression studies using quantitative real-time PCR in the half-smooth tongue sole Cynoglossus semilaevis challenged with LPS or Vibrio anguillarum. Fish & shellfish immunology. 2010;29(1):89–93.
  49. 49. Dang W, Sun L. Determination of internal controls for quantitative real time RT-PCR analysis of the effect of Edwardsiella tarda infection on gene expression in turbot (Scophthalmus maximus). Fish & shellfish immunology. 2011;30(2):720–728.
  50. 50. Sellars MJ, Vuocolo T, Leeton LA, Coman GJ, Degnan BM, Preston NP. Real-time RT-PCR quantification of Kuruma shrimp transcripts: a comparison of relative and absolute quantification procedures. Journal of biotechnology. 2007;129(3):391–399. pmid:17350129
  51. 51. Percipalle P, Visa N. Molecular functions of nuclear actin in transcription. The Journal of cell biology. 2006;172(7):967–971. PubMed Central PMCID: PMC2063754. pmid:16549500
  52. 52. Jorgensen SM, Kleveland EJ, Grimholt U, Gjoen T. Validation of reference genes for real-time polymerase chain reaction studies in Atlantic salmon. Marine biotechnology. 2006;8(4):398–408. pmid:16676145
  53. 53. Small BC, Murdock CA, Bilodeau-Bourgeois AL, Peterson BC, Waldbieser GC. Stability of reference genes for real-time PCR analyses in channel catfish (Ictalurus punctatus) tissues under varying physiological conditions. Comparative biochemistry and physiology Part B, Biochemistry & molecular biology. 2008;151(3):296–304.
  54. 54. Rodrigues PN, Dixon B, Roelofs J, Rombout JH, Egberts E, Pohajdak B, et al. Expression and temperature-dependent regulation of the beta2-microglobulin (Cyca-B2m) gene in a cold-blooded vertebrate, the common carp (Cyprinus carpio L.). Journal of Immunology Research. 1998;5(4): 263–275.
  55. 55. Beau-Faller M, Ruppert A M, Voegeli A C, Neuville A, Meyer N, Guerin E, et al. MET gene copy number in non-small cell lung cancer: molecular analysis in a targeted tyrosine kinase inhibitor naive cohort. Journal of Thoracic Oncology. 2008;3(4):331–339. pmid:18379349
  56. 56. Bell DW, Lynch TJ, Haserlat SM, Harris PL, Okimoto RA, Brannigan BW, et al. Epidermal growth factor receptor mutations and gene amplification in non-small-cell lung cancer: molecular analysis of the IDEAL/INTACT gefitinib trials. Journal of clinical oncology: official journal of the American Society of Clinical Oncology. 2005;23(31):8081–8092.
  57. 57. Bentzien F, Zuzow M, Heald N, Gibson A, Shi Y, Goon L, et al. In vitro and in vivo activity of cabozantinib (XL184), an inhibitor of RET, MET, and VEGFR2, in a model of medullary thyroid cancer. Thyroid: official journal of the American Thyroid Association. 2013;23(12):1569–1577. PubMed Central PMCID: PMC3868259.
  58. 58. Dheilly N M, Lelong C, Huvet A, Favrel P. Development of a Pacific oyster (Crassostrea gigas) 31,918-feature microarray: identification of reference genes and tissue-enriched expression patterns. BMC genomics. 2011;12(1):1.
  59. 59. Oka T, Vasile E, Penman M, Novina CD, Dykxhoorn DM, Ungar D, et al. Genetic analysis of the subunit organization and function of the conserved oligomeric golgi (COG) complex: studies of COG5- and COG7-deficient mammalian cells. The Journal of biological chemistry. 2005;280(38):32736–32745. pmid:16051600
  60. 60. Kim JH, Kim TW, Kim SJ. Downregulation of ARFGEF1 and CAMK2B by promoter hypermethylation in breast cancer cells. BMB Reports. 2011;44(8):523–528. pmid:21871176
  61. 61. Paolacci AR, Tanzarella OA, Porceddu E, Ciaffi M. Identification and validation of reference genes for quantitative RT-PCR normalization in wheat. BMC molecular biology. 2009;10:11. PubMed Central PMCID: PMC2667184. pmid:19232096
  62. 62. Gimenez MJ, Piston F, Atienza SG. Identification of suitable reference genes for normalization of qPCR data in comparative transcriptomics analyses in the Triticeae. Planta. 2011;233(1):163–173. pmid:20960006
  63. 63. Carvalho K, de Campos MKF, Pereira LFP, Vieira LGE. Reference gene selection for real-time quantitative polymerase chain reaction normalization in "Swingle" citrumelo under drought stress. Analytical biochemistry. 2010;402(2):197–199. pmid:20363209
  64. 64. Shi XQ, Guo WC, Wan PJ, Zhou LT, Ren XL, Ahmat T, et al. Validation of reference genes for expression analysis by quantitative real-time PCR in Leptinotarsa decemlineata (Say). BMC Res Notes. 2013;6:93. PubMed Central PMCID: PMC3600673. pmid:23497596