Identification of candidate reference genes for qRT-PCR normalization studies of salinity stress and injury in Onchidium reevesii

Animals and treatments

Shanghai Ocean University of Leicester granted Ethical approval to carry out the study within its facilities (Approval number: Shou-DW-2019-010). Sexually mature O. reevesii adults were obtained from the coast of Shanghai, China, and housed in the laboratory according to Shen’s method (Shen et al., 2004) in aquaculture tanks (each tank containing no more than 50 individuals), with enough fine silt and seawater (10 ppt) to simulate their natural living environment and shelters for the animals to hide. Feeding (corn flour and diatoms), removal of faeces and fresh seawater replacement were performed in a timely manner. The animals were allowed to acclimate to the new environment for 7 days before the experiments. For salinity stress, 150 O. reevesii individuals of the same size and weight were divided into five salinity groups (0 ppt, 5 ppt, 15 ppt, 25 ppt, 35 ppt). The water used in the experiment was saline, and the volume of water used ensured that all samples were retained. For injury stress, a surgical blade was used to inflict wounds of ∼5–7 mm long and 2–3 mm deep on the dorsal skin of 30 O. reevesii individuals; the experiment was divided into 5 groups according to the weight of the individuals: 8 g, 11 g, 15 g, 18 g and 22 g. To eliminate the influence of temperature change, the experimental temperature was maintained at approximately 25 °C, the temperature at which O. reevesii exhibits the best activity in the field.

Sampling

Dorsal skin tissues were used as samples in salinity stress and injury experiments; sampling was performed at 2 h, 4 h, 12 h, 24 h, 48 h and 7 days. Three individuals from each group were sampled at every time point (Bai et al., 2014). Dorsal skin tissues were also utilized as samples for the groups of different weights. Six tissues were used for assessing gene expression, including the dorsal skin muscle, intestine, lip, pleopod, liver pancreas and gonad. All of the samples were placed in freezer tubes, frozen in liquid nitrogen and stored at −80 °C.

Total RNA extraction and first-strand cDNA synthesis

Samples were rapidly ground in liquid nitrogen, and total RNA was extracted using TRIZOL (TaKaRa, Otsu, Japan). All RNA samples were resuspended in RNase and DNase-free ddH₂O (TaKaRa, Otsu, Japan). The integrity of total RNA was determined by 1% agarose gel electrophoresis, and a NanoDrop 2000c spectrophotometer (Thermo Scientific, Wilmington, DE, USA) was employed to evaluate the RNA quality and concentration. RNA with a 260/280 ratio between 1.8 and 2.2 and 260/230 ratio >1 and <3 was considered satisfactory for use in experiments.

Selection of internal control genes

Seven commonly used reference genes were identified from the transcriptome of O. reevesii. The following are some of the criteria we applied for selecting candidate reference genes: to minimize the effect of co-regulation, the nucleic acid sequence should encode a protein that plays various roles in cellular metabolism with different molecular functions; sequences that have previously been examined for stability in, to a certain extent, similar biological contexts; for reliability, the sequence should be tested specifically to verify the accuracy of the data (Die et al., 2017). The previously published candidate reference genes (Purohit et al., 2015; Altmann et al., 2015; Etich et al., 2016) are cyclophilin (CYC), beta actin (ACTB), elongation factor-1 alpha (EF1a), ribosomal protein L28 (RPL28S), β-tubulin (TUBB), ubiquitin (Ubiq) and 18S ribosomal RNA (18S RNA). To select corresponding sequences, all candidate genes selected from the transcriptome were analysed in the NCBI database using BLAST, and sequences were uploaded to the NCBI database to obtain the GenBank accession number (Table 1).

Primer design and real-time qPCR assays

The software Primer 5.0 was used to design gene-specific primers. Optimized primer pairs were selected based on their amplification efficiencies and specificities (D’Haene, Vandesompele & Hellemans, 2010). The specificity of the PCR primers utilized was evaluated using the melting curve produced by the Applied Biosystems™ QuantStudio™ 6 Flex Real-Time PCR System (Thermo Fisher, Waltham, MA, USA). The fragment size of the PCR product was determined by 2.0% agarose gel electrophoresis. Primer standard curves were created using a cDNA dilution series, and amplification efficiencies were calculated using the following equation: E = (10^−1∕slope − 1)∗100% (Radonić et al., 2004). A correction for different amplification efficiencies was introduced in the sample quantification process (Marino & Cook PMiller, 2003).

NovoStart^® Reverse Transcriptase (NOVOPROTEIN, Shanghai, China) was used to reverse transcribe 2000 ng of total RNA from each sample into cDNA in a 20 µl volume. To allow for increasing the volume of template cDNA for subsequent experiments, stock cDNA samples were diluted 5-fold for real-time PCR. Applied Biosystems™ QuantStudio™ 6 Flex Real-Time PCR System was employed for real-time PCR using a 20 µl reaction system containing the following components: 2 µl cDNA sample, 10 µl 2 ×NovoStart^® SYBR qPCR SuperMix (NOVOPROTEIN, Shanghai, China), 0.8 µl each primer, 0.4 µl ROX Reference Dye II and 6.0 µl ddH₂O. Following denaturation at 94 °C for 5 min, 40 cycles of melting at 95 °C for 15 s, annealing at 57 °C for 20 s and extension at 72 °C for 30 s were carried out. A melting curve analysis was also performed.

Data analysis

Three software tools, geNorm, BestKeeper, and NormFinder, were used to evaluate reference gene stability according to their respective (Vandesompele et al., 2002; Andersen, Jensen & TF, 2004; Pfaffl et al., 2004). The geNorm program ranks the most stable reference genes based on the average pairwise variation of a reference gene with other selected housekeeping genes and sorts reference genes using their expression stability value (M). In general, the lower the M value, the higher the expression stability. BestKeeper predicts “ideal” reference genes according to pair-wise correlation analysis among all pairs of candidate reference genes. NormFinder calculates the entire variation of candidate reference genes in all samples and also performs intragroup and intergroup comparisons. RefFinder (https://omictools.com/reffinder-tool) is a user-friendly web-based comprehensive tool developed for evaluating and screening reference genes from extensive experimental datasets. It integrates the currently available major computational programs (geNorm, NormFinder, BestKeeper) to compare and rank candidate reference genes. Based on the rankings from each program, RefFinder assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall final ranking (Kim et al., 2010).

Real-time PCR amplification efficiencies

To ensure comparability among the seven housekeeping genes evaluated, PCR amplification efficiencies were calculated based on the slopes resulting from the measurement of cDNA serial dilutions. All tested reference gene amplification efficiencies were in the range of 90–110% (RPL28S (94.3%), ACTB (97.03%), TUBB (97.32%), CYC (96.61%), EF1α (92.03%), 18S RNA (95.54%) and Ubiq (100.89%). The melting curves for all amplification products demonstrated a single peak, verifying primer-specific amplification (Fig. S1).

Expression stability of the tested genes

Based on the statistical results of expression analyses for each housekeeping gene in the seven different tissue types according to BestKeeper (Table 2), genes showing the lowest expression level were CYC and Ubiq, with Ct mean (AM[Ct]) values of 28.96 and 31.09 cycles, respectively. The genes exhibiting the highest expression levels were ACTB, TUBB and EF1α, with AM[Ct] values of 22.16, 23.06 and 19.75 cycles, respectively. RPL28S displayed the least variation in expression among the analysed tissues [SD (±Ct) = 0.48], whereas 18S RNA [SD(±Ct) =1.72] and ACTB [SD(±Ct) = 2.67] showed the greatest variation. Candidate reference genes with an SD value higher than 1 must be considered unsuitable (Mehta et al., 2010), and the M value of ACTB according to the geNorm program was clearly higher than 1.5. Overall, the ranking of these candidate housekeeping genes by BestKeeper, NormFinder and geNorm programs was consistent. In addition, ACTB and 18S RNA appeared to be less stable, whereas RPL28S, EF1 α, and TUBB were the most stable genes (Table 2; Fig. 1 Group1, Fig. 2A). According to geNorm, a V2/3 value of 0.15 is the proposed cut off value, under which the inclusion of an additional reference gene is not necessary (Fig. S2A). In this study, the V2/3 value was 0.265, and the values of V3/4, V4/5, V5/6, and V6/7 were all greater than 0.15, indicating no need to include another gene in as a normalization factor. Thus, the optimal number of reference genes for normalization in this example is one. The comprehensive ranking of the seven candidate genes given by RefFinder indicated RPL28S as the best reference gene. EF1a was also found to be an ideal reference gene when abundant reference gene expression is required.

Expression stability in muscle of the seven genes under different salinity treatments

Considering that the salinity of O. reevesii’s habitat often fluctuates greatly due to the influence of tides, river flow, rainfall and other factors, multiple salinity gradients were established to observe changes in housekeeping gene expression to explore the response of O. reevesii to changes in salinity (osmotic pressure regulation ability).

0 ppt salinity

Under conditions of varying salinity, BestKeeper, NormFinder, and geNorm were employed to examine the expression level of each housekeeping gene. According to BestKeeper (Table 3), the expression stability ranking of the seven reference genes for the seven sampling time points was as follows: 18S RNA >RPL28S >ACTB >TUBB >Ubiq >EF1a >CYC. Additionally, geNorm analysis showed that TUBB, EF1a and ACTB had the highest stabilities (Fig. 1 Group 2). NormFinder analysis revealed ACTB to be the gene with the greatest stability (Fig. 2B). Because the Vn/n+1 values provided by geNorm were all greater than 0.15, multiple housekeeping genes were not required as internal references (Fig. S2B). Based on the results of a comprehensive analysis by RefFinder, we recommend that TUBB be used as an internal reference for qRT-PCR analyses of O. reevesii under this condition. However, when the reference gene should exhibit a high level of expression, ACTB is the most appropriate choice.

5 ppt salinity

According to BestKeeper (Table 4), the expression stability ranking of the seven candidate genes was 18S RNA >RPL28S >EF1a >ACTB >TUBB >Ubiq >CYC, an order that was slightly different from the result for the 0 ppt salinity condition. However, except for 18S RNA (0.68), the SD[±Ct] values for the housekeeping genes were greater than 1. NormFinder (Fig. 1 Group3) and geNorm (Fig. 2C) indicated EF1a, RPL28S and TUBB to be the most stable genes, with Vn/n+1 values all greater than 0.15 (Fig. S2). Although the results from the three software programs are inconsistent and not ideal, according to RefFinder, RPL28S should be used as a reference gene under this condition. Due to the insufficient level of 18S RNA, RPL28S and EF1a expression, ACTB should be used as a reference gene under the stress of this level of salinity.

15 ppt salinity

BestKeeper results (Table 5) showed a stability ranking of expression for the seven candidate genes of CYC>ACTB >EF1a>RPL28S >TUBB >18S RNA >Ubiq, which was significantly different from that of the 0 ppt and 5 ppt conditions, even though SD[ ±Ct] for all housekeeping genes was lower than 1. According to NormFinder (Fig. 1 Group4) and geNorm (Fig. 2D), EF1a, ACTB, CYC and TUBB were the most stable reference genes, and the V2/3 values of 0.127 suggested that the normalization factor should preferably consist of two housekeeping genes (Fig. S2D). The best combination of two genes was TUBB and EF1a (Fig. 2D). RefFinder indicated that TUBB was the most stable internal reference gene. Ubiq exhibited slightly lower expression, but EF1a, ACTB and CYC can all be used as internal reference genes for this treatment.

25 ppt salinity

For this condition, BestKeeper results (Table 6) revealed an expression stability ranking of Ubiq >ACTB >EF1a >RPL28S >CYC >TUBB >18S RNA; however, SD[ ±Ct] for 18S RNA was lower than 1. Different from the above conditions, the NormFinder (Fig. 1 Group5) and geNorm (Fig. 2E) results were similar to those of BestKeeper. V2/3 values were 0.156, indicating no need for multiple housekeeping genes as internal references (Fig. S2E). Although the Ubiq gene was found to be the most stable, its expression level was far lower than that of the ACTB gene. Therefore, the best reference gene is ACTB.

35 ppt salinity

TUBB >CYC >RPL28S >Ubiq >ACTB >EF1a >18S RNA was the expression stability ranking according to BestKeeper (Table 7). As with the condition of 25 ppt salinity, SD[±Ct] values were lower than 1, except for 18S RNA, and the NormFinder (Fig. 1 Group6) and geNorm (Fig. 2F) results were similar to those of BestKeeper. As V2/3 values were over 0.15, multiple housekeeping genes would not be needed (Fig. S2F). The expression abundance of all candidate genes was within the acceptable range, after considering both expression stability and abundance, the best reference gene was found to be TUBB.

Expression stability analysis of the seven genes in muscle after injury

For each of the seven housekeeping genes, transcript abundance was assessed in three independent muscle pools collected at time points ranging from 2 h to 7 days after skin damage. According to BestKeeper (Table 8), the lowest level of expression gene was displayed by Ubiq, with an AM[Ct] value of 29.85 cycles. In contrast, ACTB and EF1α showed the highest levels, with AM[Ct] values of 18.75 and 18.92 cycles, respectively. The BestKeeper program indicated SD(±Ct) values lower than 1 for the reference genes other than CYC, Ubiq and 18S RNA; the geNorm M values of the seven reference genes were also lower than 1.5. BestKeeper indicated that the gene exhibiting the least variation in gene expression was ACTB [SD(±Ct) = 0.46]. However, NormFinder (Fig. 1 Group7) and geNorm (Fig. 2G) showed TUBB to be the most stable. geNorm analysis V2/3 values were 0.137, suggesting that the normalization factor should comprise additional housekeeping genes (Fig. S2G). The most suitable two-gene combination was TUBB plus EF1a (Fig. 2G). If a single gene is used as an internal reference, RefFinder indicated EF1a as the most appropriate choice.

Expression stability analysis of the seven genes in individuals of different weights

Regarding groups of O. reevesii of different weights, BestKeeper analysis (Table 9) revealed an SD[±Ct] lower than 1 only for CYC gene. In contrast, NormFinder (Fig. 1 Group8) and geNorm (Fig. 2H) showed RPL28S, TUBB and EF1a to be significantly more stable than the other candidate genes, with Vn/n+1 values greater than the threshold of 0.15 (Fig. S2H). Based on RefFinder comprehensive analysis, RPL28S is more suitable as an internal reference gene when compared to the other candidates, similar to the results for the skin injury condition. Because a reference gene requires a high level of expression, ACTB is the most appropriate reference gene for this experimental condition.