Journal of Animal Reproduction and Biotechnology 2019; 34(4): 280-291
Published online December 31, 2019
Copyright © The Korean Society of Animal Reproduction and Biotechnology.
Sang Yoon Lee1, Choul-Ji Park2 and Yoon Kwon Nam1,*
1Department of Marine Bio-Materials and Aquaculture, Pukyong National University, Busan 48513, Korea
2Genetics and Breeding Research Center, National Institute of Fisheries Science, Geoje 53334, Korea
Potential utility of 14 candidate housekeeping genes as normalization reference for RT-qPCR analysis with developmental samples (fertilized eggs to late veliger larvae) in Pacific abalone
Keywords: abalone Haliotis discus hannai, developmental samples, reference genes, RT-qPCR assay
Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) is a versatile and routine tool for gene expression analysis (Udvardi et al., 2008; Bustin et al., 2009). However, the accuracy of RT-qPCR results may be affected by both biological and technical variations, and an effective normalization using appropriate reference gene(s) is essential for reliable interpretation of the expression of target genes (Guénin et al., 2009). Although several traditional housekeeping genes such as cytoskeletal β-actin (ACTB) and glyceraldehyde 3’-phosphate dehydrogenase (GAPDH) have been widely used for RT-qPCR assays, a considerable number of reports has also claimed that those traditional reference genes would not always show a stable expression under certain conditions (Li and Shen, 2013; Taylor et al., 2013; Yuan et al., 2014; Lee and Nam, 2016a). Ideally, reference genes must be invariantly expressed and non-regulated by experimental or biological conditions. However, it is widely agreed that no given single reference gene can be universally applied to all experimental conditions and the utility of a reference gene could be largely variable under different conditions (Guénin et al., 2009; Nakayama et al., 2018). Hence, in most situations, it is often unavoidable that the most appropriate reference gene(s) may be empirically determined in a case-specific manner.
Developmental samples (embryos and early larvae) of animals have been often reported to display dynamic modulations of various gene sets related with morphogenesis and organ development (Deschamps and Duboule, 2017; Praggastis and Thummel, 2017). Zygotes of many aquatic animals representing external fertilization should undergo a ‘maternal-zygotic transition’ of gene expression, which makes the interpretation of gene expression profile complex and complicated (Blaxter, 2013; Liu et al., 2014; Romney and Podrabsky, 2017). Newly hatched larvae would be protected from embryonic membrane no longer and they should mount significant changes of gene expression to enter new environment. During early ontogeny, larvae often undergo a process involving conspicuous and relatively abrupt change in body structure, which is called metamorphosis. Metamorphosis is known to be accompanied with expression changes of diverse genes related with cell growth and differentiation (Das et al., 2006; Alves et al., 2016; Zhao et al., 2016). Taken together, careful efforts are needed to precisely interpret the expression pattern of target genes with developmental samples, and undoubtedly the use of appropriate reference genes is a critical component for the accurate quantification of the target gene of interest. Selection of unstable housekeeping genes that are involved in developmental and ontogenic processes may lead to serious misinterpretations.
The objective of this study was to assess suitable reference genes for RT-qPCR normalization with developmental and larval samples in Pacific abalone,
Mature female (n = 5-8) and male (n = 5) abalones were induced to release eggs and sperm based on the conventional method of air exposure followed by ultraviolet (UV)-irradiated seawater treatments. Eggs were washed three times with 1 μm-filtered seawater at 18-19°C and inseminated with sperm using wet methods. Fertilized eggs were placed on static incubator containing 1 μm-filtered seawater at 20°C until hatch. Fertilization rate was estimated with at least 330 randomly chosen embryos as percentage of embryos showing successful progress initial cleavages (at 2 hours post insemination: HPI) out of initial number of eggs inseminated. Hatching success was also calculated as percentage of hatched larvae out of initial number of eggs with 200-300 eggs. Both fertilization rates and hatching success were estimated in triplicates per egg batch. After hatch, about 300,000-500,000 hatchlings were incubated in rectangular tank containing 15 tons of 1 μm-filtered at 19 ± 1°C with a daily water exchange rate of 120%. Rearing of swimming larvae were continued until late veliger stage. Approximately 10,000 embryos were sampled at just fertilized (0 HPI), 2-4 cell stage (2 HPI), 8-16 cell stage (4 HPI), morula stage (5 HPI), gastrula stage (7.5 HPI) and hatching (hatched trocophore; 15 HPI). In addition, about 10,000 swimming larvae at early veliger stage (20 HPH) and late veliger stage (45 HPH) were netted for sampling (Lee and Nam, 2016b). Three independent spawning trials were prepared using different abalone broods in order to prepare three biological replications for each developmental stage. Within each biological replication, triplicate samplings were prepared at each developmental stage as triplicate technical replications.
Candidate reference genes were selected based on previous reports on abalone housekeeping genes (Wan et al., 2011; López-Landavery, 2014; Lee and Nam, 2016a). With cDNA templates from abalone developmental samples, specific amplification of each gene segment was confirmed using end-point RT-PCR and ethidium bromide staining of the RT-PCR band (data not shown). Amplified RT-PCR product were sequenced and correct annotation of each gene was validated with BLAST search against NCBI GenBank (http://www.ncbi.nlm.nih.gov). A total of 14 candidate reference genes were selected: cytoskeletal β-actin (ACTB), β-tubulin (B-TU), elongation factor 1α (EF-1A), glyceraldehyde 3-phosphate dehydrogenase (GAPDH), peptidyl-prolyl cis-trans isomerase B (PPIB), ribosomal protein isoforms L3 (RPL3), L4 (RPL4), L5 (RPL5), L7 (RPL7), L7A (RPL7A), L8 (RPL8), L36 (RPL36), S18 (RPS18), and ubiquitin-conjugating enzyme E2 (UBE2). For GenBank accession numbers, refer to Lee and Nam (2016a). qPCR primers used in this study have 20 bp in length, 50-60% of G + C content, 58-63°C of melting temperature, 156-168 bp in amplicon size.
Upon being sampled, embryos and larvae were immediately frozen on dry ice and stored at -85°C until used. Total RNA was extracted using Trizol Reagent (Invitrogen; Thermo Scientific, Wilmington, DE, USA) and further purified using RNeasy Plus Micro Kit (Qiagen, Hilden, Germany) according to the manufacturers’ protocols. The quality and quantity of the purified RNA were examined by measuring the absorbance ratios at 260/230 nm and 260/280 nm using a Nano Drop ND1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA).
A sample of total RNA was reverse transcribed into cDNA using Omniscript RT Kit (Qiagen) according to the manufacturer’s instructions. Synthesized cDNA was 8-fold diluted with sterile distilled water and 2 μL of the diluted cDNA was used as a template for qPCR amplification. qPCR amplification reactions were performed using Light Cycler 480 (Roche Applied Science, Penzberg, Germany) and 2X SYBR Mater Mix (Roche) under the following thermal cycling conditions: 45 cycles of denaturation at 95°C for 5 s, 58°C for 10 s and 72°C for 15 s with an initial denaturation step at 95°C for 2 min. At the end of each run, melting curve analysis was performed from 65°C to 95°C to confirm the specificity of amplification reaction. To validate PCR efficiency (
Because triplicate samplings were made at each of eight developmental stages from every biological replication, a total of 72 cDNA templates were finally prepared (24 cDNA samples per biological replication). The three samples at each stage from each spawning batch were treated as technical replications. Consequently, within each spawning batch, 336 amplification reactions were carried out to generate 112 median Cq values (8 stage-specific Cq values for each of 14 reference genes). Within this context, a total of 1,008 amplification reactions were made to create 336 Cq values for analysis of expression stability in this study.
The stability of gene expression was assessed using geNorm (Vandesompele et al., 2002; Hellemans et al., 2007), BestKeeper (Pfaffl et al., 2004), NormFinder (Andersen et al., 2004), and comparative ΔCT method (Silver et al., 2006). Then, an overall final ranking of expression stability for each reference gene was determined by calculating geometric mean value of rankings from the four different algorithms. The geNorm analysis was performed using qbaseplus (Biogazelle). Raw Cq values were converted into relative quantity values (RQ) via the formula RQ = E-ΔCq, where E is the validated amplification efficiency of each gene, and ΔCq is the difference between the Cq value and the minimum Cq of each gene among the samples. Expression stability (geNorm M value) of each candidate reference gene was estimated based on the average pairwise variation for the reference gene with all of the other genes. A low M value represented stable gene expression with the cut-off value of 1.5. The pairwise variation (Vn/n + 1) between the two sequential normalization factors (NFn and NFn + 1) was also estimated to define the optimal number of reference genes needed for developmental and larval samples. geNorm results were compared with three different software algorithms. Microsoft Excel-based BestKeeper calculated stability of reference genes based on the standard deviation (SD) of untransformed Cq values, where the candidate reference gene with the lowest SD is considered as the most stable gene. The NormFinder add-in for Microsoft Excel was used with the same RQ values for geNorm in order to give a stability value with the lowest value as most stable expression. The ΔCT method compares the relative expression of pairs of genes within each sample and calculated the average SD against the other reference genes to confidently identify useful reference genes.
Effects of choice of different reference genes on relative expression levels of target gene were tested. Based on our previous study, metallothionein (MT) gene was selected as a target gene because it has been reported to represent a significant difference in its mRNA levels between early embryos and swimming larvae (Lee and Nam, 2016b). A segment of MT transcripts (amplicon = 189 bp) was amplified with a primer pair (forward: 5´-GGTACCGACTGCAAGTGTAA-3´ and reverse: 5´-TCATCGGAAGTCATGTGAGC-3´) (Lee and Nam, 2016b). The cDNA templates were chosen from two developmental stages, just fertilized and late veliger. Expression level of MT transcripts in late veliger larvae relative to that in fertilized embryos was normalized with the three most stable and three most unstable reference genes determined in the overall stability ranking. Relative expression levels of MT genes determined by different references were tested using ANOVA (followed by Duncan’s multiple range tests) with a significance level at
Fertilization rates of all of three spawning batches were higher than 88% and averaged to be 92.0 ± 2.4% [mean ± SD] (p > 0.05). Hatching success was slightly more variable among batches ranging to 66% to 81% [75.5 ± 6.4%] (p > 0.05) (Fig. 1). No notable sign for the outbreak of mass mortality was observed during the rearing of swimming larvae, although we didn’t quantitatively estimate the larval viability.
Mean fertilization rate and hatching success of artificially propagated Pacific abalone
The PCR efficiency of each gene was at least higher than 0.94. The slopes of the standard curves ranged from -3.213 to -3.098, and the correlation coefficients (
Raw Cq profile of 14 candidate reference genes in all developmental samples in Pacific abalone. In the box plot, the lower and upper boundaries of each box indicate the 25th and 75th percentiles, respectively. A line within the box marks the median value. Whiskers above and below the box indicate the 90th and 10th percentiles, respectively. Outliers are indicated by dots above and/or below the box. Value in parenthesis is the percent coefficient of variation (% CV) for each reference gene.
Expression modulations in developing embryos and swimming larvae were variable among reference genes. In a broad sense, the patterns of 14 genes were divided into four categories. For Pattern-I, expression levels of reference genes were progressively increased with developments, consequently giving rise to significant differences in expression levels between early embryonic stages and swimming larval stages. Reference genes such as B-TU, ACTB, EF1A and PPIB belonged to this pattern. The Pattern-II could be characterized by gradual decrease of expression levels after early cleavage stages, but after hatch, increase of expression could be observable either at early veliger or late veliger larval stage. This pattern is typically exemplified by RPL7A and GAPDH genes. The Pattern-III is similar, in overall, with the Pattern-II. However, unlike the Pattern-II, rebound or re-increase of expression at swimming larval stages was hardly seen. Typical examples were RPL7 and RPL8. Finally, RPL3 showed a rapid drop of expression level immediately after fertilization that was followed by a quick rebound at morula stage. Afterward, its expression was gradually declined with the progress of development (Pattern-IV) (Fig. 3).
Expression patterns of 14 candidate reference genes with developmental progress. Time course expressions of reference genes were manually categorized into 4 modulatory patterns based on normalized relative quantity values determined with qbaseplus. Abbreviations of developmental stages are fertilized eggs (FE), early cleavage-1 (EC-1; 2-4 cells), early cleavage-2 (EC-2; 8-16 cells), morula (MO), gastrula (GA), hatched trocophore (HT), early veliger (EV), and late veliger (LV).
With geNorm-based analysis, RPL5 and RPS18 had the lowest M values (0.154 and 0.155, respectively), indicating these two genes are most stably expressed with developmental samples examined in this study. These lowest expression stability values were followed by RPL8 (M = 0.185), RPL7 (0.212) and UBE2 (0.237). Several ribosomal protein genes (RPL36, RPL4, and PRL7A) and GAPDH showed moderate expression stability. Conversely, B-TU (1.041), ACTB (0.847), EF1A (0.596), PPIB (0.500) and RPL3 (0.411) were the five most unstable reference genes in geNorm analysis (Fig. 4A). The optimal number of reference genes determined using geNorm showed a range of V values between 0.038 to 0.155, where the V2/3 value was only 0.081 (Fig. 4B).
Expression stability of candidate reference genes determined with geNorm algorithm. (A) geNorm M value, (B) geNorm V value. In (B), the proposed cut-off V value (0.15) is indicated with dashed line.
BestKeeper identified RPL7 as the most stable gene, and also RPL8, UBE2, and RPL5 as fairly stable references. On the other hand, B-TU, ACTB, EF1A were recognized as the three least stable reference genes. According to NormFinder algorithm, RPL4, RPL7A, PPIB, RPS18 and RPL5 were the five most stable reference genes. Stability ranking of other genes by NormFinder was as follow: RPL3, GAPDH, EF1A, UBE2, RPL8, RPL36, RPL7, ACTB and B-TU. From the comparative ΔCT method, the three most stable genes were RPL4, RPL5 and RPS18, whereas the three most unstable genes were B-TU, ACTB and EF1A (Fig. 5).
Analysis of expression stability of candidate reference genes determined with BestKeeper, NormFinder and comparative ΔCT method.
Expression stability ranking of reference genes determined by individual software programs were shown in Fig. 6A. Different algorithms represented different orders of reference genes according to stability ranking. However, ACTB and B-TU were invariantly identified as the least stable references by all the algorithms (two lowest rankers). Besides these two genes, RPS18, GAPDH and RPL5 showed less variations in rankings across algorithms than did other references, whereas RPL4, PPIB and RPL7 exhibited relatively large fluctuations in ranking depending on algorithms used.
Stability ranking of 14 candidate reference genes. (A) Ranking determined by individual algorithms (geNorm, NormFinder, BestKeeper and ΔCT). (B) Overall ranking calculated with geometric means of rankings from individual programs. In (B), values in parentheses are geometric means for the most stable (RPL5) and least stable (B-TU), respectively. (C) Calculated geometric standard deviation for the final ranking of each reference gene.
An overall final ranking was calculated by integrating rankings from the four different algorithms (geNorm, BestKeeper, NormFinder, and ΔCT) (Fig. 6B). Based on the geometric mean, most and least recommended references were RPL5 and B-TU, respectively. Ranking from the most to the least stable expression was as follow: RPL5 (geometric mean = 2.51), RPL4 (2.74), RPS18 (3.31), RPL8 (4.16), RPL7 (4.56), UBE2 (5.54), RPL7A (5.58), GAPDH (6.29), RPL36 (7.50), PPIB (7.95), EF1A (10.84), ACTB (13.00) and B-TU (14.00). Calculation of geometric standard deviations for the ranking of each reference showed that ACTB, B-TU and RPS18 received relatively consistent or similar ranks by different algorithms, whereas rankings of RPL7, PPIB and RPL4 were more variable according to algorithms used (Fig. 6C).
MT expression of veliger larvae relative to fertilized embryos were measured to be 10.9-, 7.9- and 10.7-fold increase, when normalized with the three most recommended references RPL5, RPL4 and RPS18, respectively. Average fold increase was 9.84 1.26 relative to fertilized embryos. However, normalizations with the three least stable references (EF1A, ACTB and B-TU) showed in significantly different results, in which relative expression level in veliger larvae were 1.6, 0.3 and 0.1 normalized with the reference genes EF1A, ACTB and B-TU, respectively. Calculated average expression level with these three references was 0.70 ± 0.60, indicating 1.4-fold downregulation at veliger larvae relative to fertilized embryos (Fig. 7).
Relative expression levels of metallothionein transcripts normalized with the three most stable and three least stable reference genes. Statistical differences in expression levels of MT with different normalization references were tested with ANOVA followed by Duncan’s multiple range tests at
Biological and technical errors introduced throughout the RT-qPCR process should be normalized using appropriate reference genes by internally controlling for such errors as well as input RNA amounts. High RNA quality and purity are mandatory requirements for reliability and reproducibility of gene expression data (Schmittgen and Livak, 2008; Udvardi et al., 2008). We used only the RNA samples showing the 1.91-2.15 ratios for both 260/280 nm and 260/230. Reasonable PCR efficiencies (higher than 0.94 in this study) with fairly uniform
In the present study, we were unable to identify the universal reference gene that were consistently ranked as the best reference across the four software programs. Although the RPL4 was redundantly identified as the most stable gene in NormFinder and ΔCT method, it was only ranked 7th and 8th recommended gene in geNorm (top ranker = RPL5) and BestKeeper (top ranker = RPL7), respectively. This discrepancy may come from different algorithms used in these programs such as ‘pairwise comparison’ or ‘model-based approach’ (Nakayama et al., 2018). To manage this situation, we evaluated overall ranking of each gene based on geometric mean value of individual rankings from the four algorithms (De Santis et al., 2011; Lee and Nam, 2016a; Nakayama et al., 2018). Geometric mean-based ranking identified RPL5 as the best candidate, which was already recommended as the most stable gene by geNorm (lowest M value = 0.154). The RPL5 was ranked 2nd to 5th in other programs. The second gene is RPL4 owing to strong recommendations by NormFinder and ΔCT method as mentioned above. The 3rd recommended gene is the RPS18 with a relative uniformity in rankings across the four different algorithms, although none of program identified it as the best one. Calculation of geometric standard deviation of ranking indicates that the RPS18 (3.64) shows a lower value than all others (5.65-138.78) except ACT and B-TU (both 1.00 corresponding to zero in arithmetic standard deviation) showing invariant ranking irrespective of algorithms. Selection of RPS18 as an internal control for RT-qPCR normalization with developmental samples is congruent with the previous report on other mollusk species (Crassostrea gigas) (Du et al., 2013).
In contrast to variations for top rankers, genes identified to be very unstable represented relatively consistent rankings across software programs. In particular, both B-TU (14th) and ACTB (13th) were always recognized as the two least stable references according to all the four algorithms. Their expression levels are significantly elevated with the progress of development as expressed by normalized RQ values. Fold changes of normalized RQ of B-TU and ACTB between early cleavage and veliger larval stages were more than 15 times, which was apparently dissimilar from other references (maximum 2.3-fold upregulation and 3.4-fold downregulation). Upregulation of these two genes during developmental progress may reflect the increasing demand of cytoskeletons for massive cell growth and morphogenesis particularly in larval period, because swimming abalone larvae undergo dynamic changes of shape and body architecture (Searcy-Bernal et al., 2007; Bunnell et al., 2011; Breuss et al., 2017). Meanwhile, EF1A is an another unstable gene that is consistently ranked 12th by all the algorithms except NormFinder (8th ranked), and changing pattern of normalized RQ during development is similar with those of ACTB and B-TU. GAPDH is one of widely used traditional references (Barber et al., 2005; Zainuddin et al., 2010). However, it was not robustly recommended as a suitable reference gene by any of algorithms (ranked between 7th or 8th in three programs and 4th ranked only in ΔCT method) in this study, which is similar with previous argument on the unsuitability of GAPDH as a reliable reference (Glare et al., 2002; Cho et al., 2008; Sikand et al., 2012). The RPL3 is suggested as the least stable gene out of all ribosomal protein genes tested in this study. Expression pattern of RPL3 during development could be characterized as the downregulation in cleavage stages (2-4 cell stage and further decrease at 8-16 cell stage) relative to fertilized eggs, followed by a quick rebound at morula stage. This change could be explained by a ‘maternal-zygotic transition’ pattern where maternally provisioned transcripts in eggs were replaced with de novo synthesized transcripts by zygotes after early cleavage phase (Harvey et al., 2013; Liu et al., 2014). Such a phenomenon is not clearly observed with other ribosomal protein genes examined in this study, and RPL3-specific fluctuation of expression levels may attribute to its relatively lower rankings than other ribosomal proteins.
Based on the pairwise variation V between two sequential normalization factors containing an increasing number of genes (
In conclusion, potential suitability of 14 candidate housekeeping genes as the normalization reference for RT-qPCR analysis with developmental samples in Pacific abalone was evaluated using four different statistical algorithms. Our data suggest that RPL5 and RPS18, either singly or in combination, are appropriate for normalizing gene expression in developmental samples of this abalone species, whereas ACTB, B-TU and EF1A are less stable and not recommended as references. In addition, our finding also suggests that not only geometric mean-based ranking but also standard deviation in the ranking should be taken into account together for final selection of the best combination of reference gene set. This study could be a useful basis to facilitate the generation of accurate and reliable RT-qPCR data with developmental samples in Pacific abalone.
No potential conflict of interest relevant to this article was reported.
This research was supported by the grant from the Golden Seed Project (GSP), Ministry of Oceans and Fisheries to YKN (#213008-05-2-CG600).