With advantages of relatively accurate quantification, high sensitivity and high throughput, quantitative reverse transcription polymerase chain reaction (qRT-PCR) has become one of the most widely used techniques to detect changes in gene expression (Johnson et al., 2014; Phillips et al., 2009). Indeed, the use of qRT-PCR has increased tremendously in nearly all branches of biology (Jian et al., 2008; Yang et al., 2018; Leidinger et al., 2016; Chervoneva et al., 2017). However, there are inevitably a number of influencing factors that affect the efficiency of the reaction (Bustin et al., 2009), such as discrepancy in pipetting, the RNA quality and concentration, the efficiency of reverse transcription and amplification among different samples, the PCR procedures and primer amplification efficiency. Therefore, to avoid variations or errors, it is fundamental to standardize the level of target gene expression by utilizing in parallel a reference gene as an internal control (Huggett et al., 2005). In general, an ideal reference gene should demonstrate a consistent level of expression across all tested tissues or conditions (Bustin & Bustin, 2002). Nonetheless, there is increasing evidence that expression of assumed reference genes can vary observably with experimental conditions such as developmental stage and chemical treatment, significantly affecting relative quantification and qPCR result interpretation.
Certain housekeeping genes, such as 18S ribosomal RNA (18SrRNA), 28S ribosomal RNA (28SrRNA), β-actin (ACTB), cyclophilin (CYP2), elongation factor 1-alpha (EF1a), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), translation elongation factor (TEF), tubulin (TUBB) and polyubiquitin (UBQ) (Huan, Wang & Liu, 2017; Song et al., 2017; Gao et al., 2017) are commonly used as reference genes. As these housekeeping genes are related to basic metabolic processes and are essential for normal cell growth, their expression levels are thought to be stable. However, many recent studies have revealed that fluctuations in the expression levels of housekeeping genes may be largely influenced by the tissue, individual developmental stage or experimental conditions (Hu, Xie & Yao, 2016). In fact, there is no clear evidence to show that a single universal reference gene is suitable for different tissues and varying experimental conditions (Vandesompele et al., 2002). Thus, it is crucial to determine one or two stably expressed reference genes before they are applied for normalizing the expression levels of target genes based on qRT-PCR.
Onchidium reevesii (Mollusca, Gastropoda, Pulmonata, Systellommatophora, Onchidiidae) is a brackish water amphibious sea slug that resides mainly in river ports of the intertidal zone and coastal tidal flats, reed and mangroves ecosystems (Wang et al., 2018; Sun et al., 2014). The species is rich in natural collagen, which has high nutritional and medicinal value; thus, it may constitute an excellent specialty aquatic product for humans if issues pertaining to their breeding in captivity can be resolved (Guan et al., 2013). Moreover, O. reevesii is an important part of the biodiversity of wetlands and is considered to be a model organisms for studying the evolution of marine invertebrates from sea to land (Sun et al., 2016). O. reevesii has been investigated with regard to morphological, physiological and active substance aspects, and further molecular biology research is urgent (Wang, 2014; Shen et al., 2010; Cheng et al., 2015). Additionally, suitable reference genes are important for accurately interpreting expression analyses of functional genes and for evaluating RNAi efficiency. O. reevesii is a euryhaline organism, and the salinity of its habitat varies greatly under the influence of tides, rainfall, river flow and other factors. Due to its complex geographical environment and biodiversity, these organisms are inevitably attacked and injured by natural enemies. In this study, an extreme environment (salinity) and stress conditions were imposed on O. reevesii in an effort to determine stably expressed internal reference genes under such treatment and to provide basic data for future studies on its ability to regulate osmotic pressure, to adapt to extreme environments and to repair damage.
Material and Methods
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.
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 ddH2O (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).
|Gene abbreviation||Gene name||Biological function||Primer sequence (5′–3′)||Product size (bp)||Accession no.|
|CYC||cyclophilin||Immunosuppressant, protein folding||F: GTGGGATGTTCCTCTTTACC||269||KY593122|
|RPL28||Ribosomal protein l 28||60S ribosomal subunit||F: CTGGCACAGGCAAAGTGTCC||196||KY593120|
|ACTB||β-actin||Cytoskeletal structural protein||F: GTCCACCGCAAGTGCTTCT||214||KY593121|
|TUBB||β-tubulin||Structural protein||F: GTGCTGTTGCCGATGAAAG||157||KY593125|
|EF1α||elongation factor 1-alpha||Essential component of the eukaryotic translational apparatus||F: GGAGATGCCAGCCTCAAAC||165||KY593123|
|Ubiq||Ubiquitin||Protein degradation||F: GCCGAGGCTACATTCCAGT||203||MF680835|
|18S RNA||18S rRNA||rRNA in the ribosome||F: TCCGCAGGAGTTGCTTCGAT||142||FJ843070|
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 ddH2O. 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.
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.
In this study, seven genes, CYC, RPL28S, ACTB, TUBB, EF1a, Ubiq and 18S RNA, were selected as candidates for reference gene screening for use in O. reevesii. When using a relative quantification technique to analyse mRNA expression of a target gene, the final results are more accurate and scientific when appropriate internal reference genes are employed. In general, screening of internal reference genes should to meet the following criteria: it should be widely involved in all aspects of the organism’s cellular metabolism; its expression should maintain a high level and stability in a range and with low sensitivity to changes in the external environment; candidate genes should be expressed stably under different experimental conditions.
Because of advantages of high timeliness, sensitivity and convenience, qRT-PCR is often used to screen reference genes. At present, there are three programs commonly used: BestKeeper, geNorm and NormFinder. However, due to the different modes in which data are processed by the respective software, there are some inconsistencies regarding the output. RefFinder comprehensively evaluates the results of the above three programs and provides a relatively reasonable internal control gene rank to help the experimenter ultimately determine the optimal choice.
We examined different salinity stress levels, skin muscle injury, different tissues of O. reevesii adults and individuals with different weights to identify reference genes. The most suitable internal reference genes in different tissues of Haliotis rufescens are reportedly RPL5S and CYC; the most stable in O. reevesii was found to be the RPL28S (Sun et al., 2012), but EF1a can serve as an alternative. When analysing differences in target gene expression in various tissues of O. reevesii, RPL28S should be used as the reference gene if the abundance of target gene expression is low, whereas EF1a is recommended if target gene expression is high. Previous studies have compared the most stable internal reference genes of flatworms under salinity stress, with GM2-activator protein (GM2AP) and ACTB indentified; these results are similar to those for O. reevesii, in which the most suitable internal control genes are TUBB and ACTB (Plusquin et al., 2012). However, our candidate gene expression analysis showed that at low salinity (0 ppt, 5 ppt) (low osmotic pressure environment) under laboratory conditions, expression of these genes was significantly lower than that at 15 ppt and 25 ppt. One reason for this may be that a low osmotic pressure environment leads to excess moisture entering tissue, decreasing cellular metabolism and eventually housekeeping gene expression (Mizuno & Ogawa, 2011; Orskov, 2010). Nonetheless, at high salinity (35 ppt), the housekeeping genes examined essentially maintained a normal level of expression, which may be due to the strong ability of O. reevesii to obtain water from the external environment, ensuring the stability of its internal environment. It can be inferred from the above that O. reevesii is not strongly tolerant to low salinity and osmotic pressure and that its optimal living salinity is approximately 15-25 ppt. Sun found the reference genes RPL13S and RPL32S to be the most stably expressed in contused rat skeletal muscle; however, EF1a plus TUBB constitutes the most suitable internal reference combination for O. reevesii after tissue damage (Sun et al., 2012). Heavy metal ion stress in organisms is also a major focus of current research, and as tidal flat inhabitants, O. reevesii feeds on the surface soil of these flats and may very likely serve as an indicator of heavy metal ion pollution. This is a future research direction of our laboratory (Jáuregui et al., 2015; Aydın-Önen, 2016; Authman et al., 2015). Although we compared data for O. reevesii of different sizes, because there is no consensus regarding the relationship between its growth stage and body weight, these experimental data need further confirmation. This is the first study to screen internal reference genes for O. reevesii under different conditions, and the results will be useful for relative quantitation of gene expression in the future.
In this study, we first ascertained and evaluated the expression stability of seven housekeeping genes for qRT-PCR normalization in O. reevesii tissues and under conditions of salinity stress and tissue injury. (1) In our assessment of different tissues, RPL28S and EF1a were found to be the most suitable and stable internal genes among the six tissue samples tested as well as among individuals of different weights. (2) The results suggest that ACTB and TUBB are the most stable genes, with high expression levels when assessing salinity stress. (3) Regarding muscle injury, EF1a is the most stable candidate gene. (4) Among all experimental groups, data analysis of two groups (15 ppt, injury) revealed TUBB plus EF1a to constitute a suitable reference combination. Based on our results, we propose that the three housekeeping genes ACTB, TUBB and EF1a be the first choice of reference genes for qRT-PCR. Our experimental data indicate that O. reevesii has low tolerance to low osmotic pressure and that a salinity range of approximately 15-25 ppt is the most suitable living environment for this organism. To our knowledge, this study is the first to select and evaluate optimal reference genes for O. reevesii, and the findings are expected to provide theoretical data support for future experiments involving qPCR. Although the optimal internal reference gene differs among treatments, such as during salinity stress and tissue injury, it is important to understand the importance of the selection of these genes. Overall, for different experimental studies of O. reevesii, the selection of appropriate reference genes should be taken into consideration, and our results provide basic experimental data for this purpose.
Melting curves of seven candidate reference genes
(a) to (g) represent seven candidate reference genes, CYC, RPL28S, ACTB, TUBB, EF1a, Ubiq, 18S RNA, respectively.
Pairwise variation (V) of candidate reference genes, as calculated by geNorm software
Vn/n+1 values were used to determine the optimal number of reference genes. A cutoff of 0.15 (Vn value) is usually applied. A different tissues, B dorsal muscle tissue under 0 ppt salinity; C dorsal muscle tissue under 5 ppt salinity; D dorsal muscle tissue under 15 ppt salinity; E dorsal muscle tissue under 25 ppt salinity; F dorsal muscle tissue under 35 ppt salinity; G dorsal muscle tissue after injury; H dorsal muscle tissue from animals of different weights.