Selection of housekeeping genes as internal controls for quantitative RT-PCR analysis of the veined rapa whelk (Rapana venosa)

Larvae culture and sample collection

Egg capsules of R. venosa were collected naturally from Laizhou Bay, Laizhou, China. Following published methods, larvae were incubated in appropriately sized tanks at Blue Ocean Co. Limited (Laizhou, China) (Pan et al., 2013). Newly hatched pelagic larvae were transferred to 2.5 m × 2.5 m × 1.5 m tanks with a density range of 0.3–0.05 ind/mL, determined by developmental stage. Larvae were fed a mixture of microalgae containing Platymonas subcordiformis, Isochrysis galbana, and Chlorella vulgaris (13.0 × 10⁴ cells/mL daily). Seawater was treated by sand filtration and UV irradiation before samples were cultured. Seawater temperature was below 25 ± 1 °C. Larvae samples were examined by microscope to ensure synchronous growth in developmental stages including blastula, juvenile, and adult stages. Samples were collected and washed with distilled water, frozen in liquid nitrogen, and stored at −80 °C until use. We selected five biological replicates from 12 larval stage ((blastula (L), gastrula (M), trochophore (N), early intra-membrane veliger (R), mid intra-membrane veliger (S), late intra-membrane veliger (T), one-spiral whorl larvae (C), two-spiral larvae (D), early three-spiral whorl larvae (F), late three-spiral whorl larvae (G), four-spiral whorl larvae (J), and juvenile stage (Y)) and all tested tissues (gill, hemocyte, intestine, Leiblein’s gland, liver, kidney, mantle, and gonad) were aseptically dissected from five adult specimens. Hemolymph was extracted from the pericardial cavity using a 1 mL medical injector. Hemocytes were obtained by centrifugation at 4 °C and 1,000× g for 10 min.

Total RNA extraction and cDNA synthesis

Total RNA was extracted from tissue samples of gills, intestines, Leiblein’s glands, livers, kidneys, and gonads, and from larvae of different developmental stages using MiniBEST Universal RNA Extraction Kit (TaKaRa, Tokyo, Japan), and from hemocyte and mantle using RNAiso Plus (TaKaRa) according to manufacturer’s instructions. RNA integrity was confirmed by gel electrophoresis based on the predicted product size. RNA from each sample was diluted with nuclease-free water and 0.1 µg RNA of each sample was used as the template for cDNA synthesis using a PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa). Prior to qRT-PCR, cDNA was diluted 10-fold.

Selection of candidate internal controls

According to RNA-seq transcriptome data of developmental samples, which were derived from stages C to Y and performed in triplicate for each stage (Song et al., 2016), genes with similar expression patterns were identified and classified to different clusters. Candidate housekeeping genes were selected from clusters exhibiting expression stability based the RPKM value in different developmental stages. In total, 13 genes were selected using a previously published RNA-seq data (Song et al., 2016).

Primer design and qRT-PCR

The primers for qRT-PCR were designed using Primer Premier 5 (PREMIER Biosoft, USA) and are listed in Table 1. qRT-PCR was performed using a SYBR Green^® real-time PCR assay consisting of a SYBRPrimeScript™ RT-PCR Kit II (TaKaRa) with a Mastercycler^® ep realplex S (Eppendorf; Hamburg, Germany). Amplifications were carried out in a total volume of 20 μL (10 μL of SYBR Green Master Mix, 0.4 μL of each forward and reverse primer (10 μmol∕L), 1 μL of diluted cDNA, and 8.2 μL RNase-free water) as follows: 95 °C for 2 min followed by 40 cycles of 95 ° C for 15 s, the respective annealing temperature (Tm) for 15 s, and 68 °C for 20 s. The Tm for GAPDH, EF-1α, ACT, COX1, NDUFA7, and RL5 is 58.5 °C, the Tm for RL28, TUBB, RS25, RS8, UBE2, and PPIA is 57.5 °C, and the Tm for HH3 is 56.5 °C. Melting-curve analysis of the amplification products was performed and following electrophoresis each gel picture was analyzed to confirm the product by the predicted size. Each assay was performed in triplicate and Cq values were recorded for further analysis.

Analysis of gene expression stability

The expression stability of 13 candidate housekeeping genes among the different RNA samples was calculated using the Excel-based tool GeNorm v3.4 (https://genorm.cmgg.be/) and overall stability of these candidate genes was determined using RefFinder (http://fulxie.0fees.us/?type=referencei=1). We used this suite of tools to ensure a statistically thorough analysis and robust identification of housekeeping genes for use in qRT-PCR of R. venosa.

GeNorm evaluates gene stability (M) of inputted genes using a statistical algorithm according to geometric averaging of multiple control genes and means pairwise variation of a gene from those remaining control genes in all provided samples. Vandesompele et al. (2002) proposed a value of 1.5 as a cut-off for suitability as an endogenous control, especially with heterogeneous samples such as different cell types or tissues. Based on this approach, genes with the lowest M value have the highest expression stability. The best combinations of two internal control genes with a constant level were selected by stepwise exclusion of the gene with the highest M value followed by a recalculation of new M values for all of the remaining genes (Vandesompele et al., 2002). In addition, GeNorm is used to determine the optimal number of housekeeping genes by pairwise number variation analysis. It computes the geometric mean of the selected genes that are expressed steadily for accurate normalization. A pairwise variation below 0.15, which is determined by Vn∕n + 1, means that an added control gene (n + 1) would not further improve the normalization factor (Vandesompele et al., 2002).

RefFinder is an online analysis tool, which includes GeNorm v3.4, NormFinder v20, BestKeeper v1, and delta CT, that is used to avoid one-sidedness and the potential limitations of relying on a single tool or algorithm for stability analysis to identify reference genes of interest (Pfaffl et al., 2004). Andersen, Jensen & Ørntoft (2004) developed NormFinder to estimate both overall and subgroup variation of the sample set of candidate genes for normalization factors (NFs) in a gene expression study. Three candidate housekeeping genes and two tested samples per group are required to set the minimum input data. Raw Cq values are first log-transformed and used as input. Random co-regulated genes would not bias the results of the software. NormFinder ranks the best candidate reference genes according to the lowest expression stability values and the lowest variation values by combining intra- and inter-group variability.

BestKeeper calculates descriptive statistics of the Cq values and Pearson correlation coefficient (Pfaffl et al., 2004). Internal control genes with stable expression for use as a housekeeping gene are identified based on highly correlated expression levels. The correlation between each candidate reference gene and the BestKeeper index, which is determined by calculating the geometric mean of the Cq values of the candidate genes, is estimated by the Pearson correlation coefficient (r), the coefficient of determination (r²), and P value. BestKeeper ranks the candidate housekeeping genes according to Cq variation, which is displayed as standard deviation (SD), r, and as r² with the BestKeeper index value. An SD threshold value less than 1.0 is recommended by Pfaffl et al. (2004) and the closer r² is to 1, the better.

Selection of housekeeping genes to be used as internal controls

Eight different subclusters that exhibited various gene expression patterns were identified. Genes that have similar gene expression patterns were evaluated based on their correlation and classified as a single subcluster (Fig. 1). We found that each subcluster had between 179 and 8,222 genes. The gene expression patterns of subclusters 1 and 8 were highly variable in six developmental stages, whereas genes in subclusters 2, 3, and 4 indicated the greatest expression stability. Therefore, genes from these latter subclusters were selected as eligible candidate housekeeping genes. TUBB, which is a common housekeeping gene already used as an internal control, was selected from subcluster 2. ACT, UBE2, and GAPDH were selected from subcluster 3, while the other nine candidate housekeeping genes were from subcluster 4.

Real-time PCR amplification of housekeeping genes

Single peaks of the melting curves in different samples confirmed primer accuracy and gene-specific amplification (Fig. S1). In addition, agarose gel electrophoresis exhibited a single band for each amplified gene and PCR product were confirmed based on the expected size.

Gene expression stability analysis in tissues

In our study, we identified 13 candidate reference genes and tested them for expression stability in eight different tissues, specifically, gill, hemocyte, intestine, Leiblein’s gland (Leiblein), liver, kidney, mantle, and gonad from five adult individuals.

We analyzed the raw quantification cycle data (Cq values) obtained from qRT-PCR and the determined variation among these candidate housekeeping genes in trial samples. In all tissues, we found that the average Cq values of the 13 genes ranged from 14.71 to 33.45 (Table 2). We found that COX1 had the lowest mean Cq values, which represent the highest expression levels, both in all tissues and in individual tissues, whereas HH3 had the highest mean Cq values, which indicates the lowest expression levels, in various tissues with the exception of the mantle. We found that each housekeeping gene displayed minor variability in its expression level in the various tissues under the same conditions. According to computed values of standard error (SE), we found that RL28 (SE = 0.37), EF-1α (SE = 0.38), NDUFA7 (SE = 0.39), and RS25 (SE = 0.39) have the least varying transcript abundance values when all of the tissues were analyzed together; however, we found that different housekeeping genes displayed variable levels of transcript abundance in different tissues. The genes with the lowest SE in the eight tissues examined was RS25 (SE = 0.26) in gill, TUBB (SE = 0.35 and SE = 0.72) in hemocyte and Leiblein respectively, HH3 (SE = 0.44) in intestine, RL5 (SE = 0.44) and RL28 (SE = 0.42) in liver, PPIA (SE = 0.56) in kidney, ACT (SE = 0.44) and RS25 (SE = 0.45) in mantle, and EF-1α (SE = 0.46) and NDUFA7 (SE = 0.47) in gonad. These findings demonstrate that no single candidate housekeeping gene is expressed at a stable level on the basis of Cq value only in these eight tissues from different adult R. venosa samples, and therefore, it is necessary to select better-suited housekeeping genes using additional statistical analyses.

The GeNorm-derived M values of candidate housekeeping genes in all tissues and for each tissue are shown in Fig. 2. We found that EF-1α and RL5 (both M = 0.52) showed the highest stability in all tissues as well as in gonads, kidneys, and livers. However, we found that EF-1α and RL28 were the best combination of two internal control genes for gill, intestine, and mantle, whereas the best control gene pairs for hemocyte and Leiblein were PPIA and HH3, and RL28 and UBE2, respectively. These results indicate that EF-1α , RL5, and RL28 are the most appropriate internal control genes for most tissues.

Figure 2 also shows a ranking of the stability values calculated by NormFinder for tissue-specific housekeeping genes. The best housekeeping genes on the basis of all tissues together was EF-1α with a stability value of 0.63 as well as in Leiblein and mantle (0.16 and 0.10, respectively), whereas in the gill and intestine it was COX1 (0.08 and 0.20, respectively), in hemocyte and gonad it was RL28 (0.68 and 0.84, respectively), and in liver and kidney it was both RL5 and EF-1α (both 0.17 and 0.05, respectively). Based on these findings, EF-1α is an appropriate stable internal control gene for all tissues analyzed either together or separately. In addition, RL5, COX1, and RL28 can be used for their respective tissue-specific analyses.

Using SD values generated from BestKeeper for each of the various housekeeping genes, we found major differences in each tissue (shown in Tables 3 and 4). RL28 was identified as the best gene for all tissues together and in livers, whereas RS25 was identified for use in gills and mantles, TUBB in hemocyte and Leiblein, HH3 in intestines, PPIA in kidneys, and NDUFA7 in gonads. However, as ranked by r, we found that EF-1α was the most stably expressed gene in all tissues together and in most separate tissues, namely, gills, Leiblein, mantle, gonads, and kidneys, although EF-1α and RL5 had the same rank position in kidneys. In addition, RL5 was the best gene in hemocyte, whereas COX1 and PPIA were ideal for intestines and livers, respectively. Therefore, based on our findings with r, EF-1α was the most stable housekeeping gene in most tissue samples.

Gene expression stability analysis in developmental stages

In this study, using five specimens, we identified 13 candidate reference genes and tested them for expression stability in 12 different developmental larval stages. For data description and display in tables and figures, we merged these 12 developmental larval stages into five groups/stages according to their developmental characteristics: stage 1 (L, M, N); stage 2 (R, S, T); stage 3 (C, D, F); stage 4 (G, J); and stage 5 (Y).

We analyzed the Cq values obtained from qRT-PCR and calculated variation among the candidate housekeeping genes evaluated in the samples. In all stages combined, the average Cq of the 13 genes ranged from 19.26 to 30.36 (Table 5). We found that EF-1α and COX1 had the lowest mean Cq values, which represented the highest expression level, in both the total for all stages and separately for each of the five stages, whereas we found that HH3 had the highest mean Cq values in all of the different stages, which indicates it had the lowest expression levels. Interestingly, stages 1 and 2 showed more variation in gene expression compared to that found in the other stages, and that integral Cq values decreased gradually from stage 1 to stage 3 and were then constant (Table S1). According to computed SE values, we found that HH3 (SE = 0.32), EF-1α and NDUFA7 (both SE = 0.34), and GAPDH (SE = 0.35) had the least varying transcript abundance when all stages were analyzed together. However, the tested housekeeping genes showed different expression states in each stage under the same condition. The gene with the lowest SE value for each of the five stages is RL28 (SE = 0.31) and COX1 (SE = 0.32) in stage 1, PPIA (SE = 0.62) in stage 2, TUBB and RS25 (both SE = 0.15) in stage 3, PPIA and TUBB (both SE = 0.10) in stage 4, and RS25 (SE = 0.19) in stage 5. These findings illustrate that no single candidate housekeeping gene was expressed at a stable level on the basis of Cq value only across five stages from different R. venosa larvae samples, and therefore, further statistical analyses were required to identify the best housekeeping genes.

GeNorm-derived M values of the candidate housekeeping genes for all stages together and stage-specific are shown in Fig. 3. COX1 and RL28 (both M = 0.34) have the highest stability in all the stages together. However, we found that UBE2 and PPIA were the most suitable combination of two internal controls for stage 1, RL5 and RL28 were optimal for stage 2, NDUFA7 and RL28 in stage 3, RS25 and PPIA in stage 4, and COX1 and RL5 in stage 5. These results indicate that RL28 is the most appropriate internal control gene for most of the different stages examined. In addition, Figure 3 shows the ranked stability value calculated by NormFinder of the tested candidate housekeeping genes, and we found that the best housekeeping gene for all combined stages was RL28 with a stability value of 0.188. It was also the best gene for stages 1 and 3 (0.277 and 0.228, respectively), whereas RL5 was the best gene for stages 2 and 5 with respective stability values of 0.142 and 0.112, and TUBB (0.022) for stage 4. Based on these findings, RL28 and RL5 are the appropriate stable internal control genes for most developmental stages.

Based on SD values generated from BestKeeper, we found that various housekeeping genes manifested major differences in expression in each stage (Tables 4 and 6). We identified HH3 as the best gene for all stages together, whereas the best stage-specific genes were RL28 for stage 1, NDUFA7 for stage 2, TUBB for stage 3, and RS25 for stages 4 and 5. However, when ranking by r, we found that RL28 was the most stably expressed gene when all stages were combined and for stage 3 with PPIA, RL5, TUBB, and GAPDH being the ideal genes for stages 1, 2, 4, and 5, respectively.

Determination of the optimal number of internal controls for normalization

In regards to the tissue-specific pairwise variation calculated by GeNorm (Fig. 4), we found that in most tissues (gill, gonad, intestine, Leiblein, kidney, and mantle) the V4/5 value of 0.141 indicated EF-1α and RL5 are insufficient for normalization and that RL28 and NDUFA7 should be included; however, the V2/3 value was under 0.15, which suggests two internal control genes were sufficient. In liver, the V2/3 value (0.148) was comparable to the cut-off value although the V3/4 value (0.092) was sufficiently low, which indicates the inclusion of a third reference gene is needed to improve stability of normalization (Fig. S2). We found that none of the pairwise variations determined before V10/11 was less than 0.15 in hemocyte, and therefore, over 10 genes are suitable as reference genes based on the determined conditions (Fig. S4). In terms of pairwise variation for all developmental stages (Fig. 4), we found that the V2/3 value was below threshold (0.15), which suggests RL28 and COX1 (on the basis of the M value) are sufficient for normalization, and that inclusion of an additional reference gene is not required in most stages with the exception of stage 1. We found in stage 1 that the V2/3 value exceeds threshold, whereas V3/4 is below threshold, which indicates the necessity of adding a third reference gene (RS25 based on M value) to improve the robustness of normalization (Fig. S3).