Selection of a reference gene for studies on lipid‐related aquatic adaptations of toothed whales (Grampus griseus)

Abstract Toothed whales are one group of marine mammals that has developed special adaptations, such as echolocation for predation, to successfully live in a dynamic aquatic environment. Their fat metabolism may differ from that of other mammals because toothed whales have acoustic fats. Gene expression in the metabolic pathways of animals can change with respect to their evolution and environment. A real‐time quantitative polymerase chain reaction (RT‐qPCR) is a reliable technique for studying the relative expressions of genes. However, since the accuracy of RT‐qPCR data is totally dependent on the reference gene, the selection of the reference gene is an essential step. In this study, 10 candidate reference genes (ZC3H10, FTL, LGALS1, RPL27, GAPDH, FTH1, DCN, TCTP, NDUS5, and UBIM) were initially tested for amplification efficiency using RT‐qPCR. After excluding DCN, the remaining nine genes, which are nearly 100% efficient, were selected for the gene stability analysis. Stable reference genes across eight different fat tissue, liver, and muscle samples from Grampus griseus were identified by four algorithms, which were provided in Genorm, NormFinder, BestKeeper, and Delta CT. Finally, a RefFinder comprehensive ranking was performed based on the stability values, and the nine genes were ranked as follows: LGALS1 > FTL > GAPDH > ZC3H10 > FTH1 > NDUS5 > TCTP > RPL27 > UBIM. The LGALS1 and FTL genes were identified as the most stable novel reference genes. The third‐ranked gene, GAPDH, is a well‐known housekeeping gene for mammals. Ultimately, we suggest the use of LGALS1 as a reliable novel reference gene for genomics studies on the lipid‐related aquatic adaptations of toothed whales.


| 17143
SENEVIRATHNA ET Al. phylogeny of complete mitochondrial genomes (Senevirathna et al., 2021). Evolutionary forces and environmental factors may affect the shaping of rare dolphin populations. These major evolutionary forces consist of genetic mutations, natural selection, genetic drift, and gene flow (Saeb & Al-Naqeb, 2016). Previous research has examined the evolution of fitness rates due to the interaction of social and genetic factors in a bottlenose dolphin population (Frère et al., 2010). In addition, an analysis of dolphin genomes has observed an adaptive evolution of nervous system genes and a slow metabolic rate (McGowen et al., 2012). Therefore, evolutionary studies on metabolic genes are vital to identify special aquatic adaptations of these cetaceans.
We theorize that genes involved in lipid metabolism may have a stronger evolutionary influence on cetaceans. Furthermore, evolutionary or selective pressure can affect gene evolution for cetacean aquatic adaptations. The intake of resources from the living environment impacts the adaptive evolution of metabolism and diversification and the synthesis of fatty acids in marine animals (Twining et al., 2021). One important ecological feature of the Risso's dolphins is their carnivorous feeding behavior (Baird, 2009). Feeding mainly on squids, such as cephalopods, is a common habit of toothed whales, which may relate to their lipid metabolism.
Generally, toothed whales exhibit unique metabolic adaptations, such as heat regulation to survive in the cold marine environment (Yuan et al., 2021), and metabolic changes can be caused by anthropogenic and climatic factors. Toothed whales consist of several types of fat deposits, including the acoustic fats in the head region, which are involved in the special aquatic adaptation of echolocation.
These acoustic fats predominantly contain unusual wax esters and triglycerides (Koopman et al., 2006;Norris, 1974). However, their fatty acid composition varies based on the species, age, and type of tissue (Koopman, 2018). There is an increasing interest in studying toothed whales at the genomic and transcriptomic levels to reveal potential genes and distinctive metabolic pathways that are implicated in this specialized process of fat-metabolism-related aquatic adaptations. A recent supportive study has identified positively selected genes for lipid metabolism in Cetacea as well as unique features for examining functional modifications of multiple genes (Endo et al., 2018). However, environmental stressors can alter lipid metabolism pathways and directly affect lipid oxidation in animals (Koelmel et al., 2020). Lipidomics is also relevant in environmental toxicology research; it can be applied to lower organisms and higher organisms, such as mammals, and is an ideal molecule for analyzing animal health (Aristizabal-Henao et al., 2020). Another study has confirmed the importance of identifying a stable reference gene to clarify biomarkers for exploring environmental stresses on the cellular adaptations of animals using Nacella spp. under different environmental conditions (Koenigstein et al., 2013). Therefore, environmental lipidomics represents the next level of molecular application in wild animal research, and investigations of lipid-related reference genes in toothed whales will be key for various future ecological studies concerning toxicology and lipidomics.
A reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) has numerous biological applications for targeted gene expression and is especially advantageous for non-model species (Fassbinder-Orth, 2014). This technique requires an accurate and reliable reference or housekeeping gene (HKG) to normalize specific gene expression data for relative quantification (Almeida-Oliveira et al., 2017). In certain tissues of a particular species, the expression level of genes may differ significantly compared with a selected reference gene (Radonić et al., 2004). Therefore, it is crucial to identify a reliable reference gene before conducting a comparison with the genes of interest. Specifically, in toothed whales, the adipose tissue metabolic pathway is still not clear, and supportive lipid metabolism genes have not been investigated at the transcriptomic level. It is possible that there are unique de novo biosynthesis pathways of lipids in toothed whales for various ecological adaptations (Koopman, 2018). Therefore, identifying a stable reference gene is necessary for future RT-qPCR assays of the functional genomics of lipid metabolic pathways to understand molecular evolutionary implications.
The present study used tissue samples from Risso's dolphins ( Figure 1). The experimental design emphasized the identification of a novel reference gene for adipose tissues (Figure 2). Ten candidate genes were initially selected based on transcriptomic fragments per kilobase of transcript per million mapped reads (FPKM) values, the coefficient of variance (CV), gene function, and relevance to metabolism. Then, the identified genes were evaluated for amplification efficiency through RT-qPCR. Only the most efficient genes were retained for further stability analysis using eight types of fat tissues, liver tissue, and muscle tissue with the aim of identifying the most stable reference gene. The null hypothesis, which predicted no significant difference in selected reference gene stability, was tested by four statistical algorithms. This work ultimately validates stable novel genes in the selected tissue types. Future research can utilize these findings for gene expression analysis in marine mammal species and particularly to study fat metabolism in toothed whales. Ten types of tissue from three male Risso's dolphins (sample ID: 19TK409,19TK410,and 19TK411) were received from the Taiji Fisheries Association with the cooperation of the biological surveys of the National Research Institute of Far Seas Fisheries and the Japan Fisheries Research and Education Agency. The tissues included four parts of the melon, two types of jaw fat, two parts of blubber, liver, and muscle. All of these tissue types are non-reproductive; therefore, the use of only male samples for this analysis does not risk any sex-related bias. These tissues were preserved in an RNAlater solution at the site before being transported to the laboratory, where they were stored at −80°C until RNA extraction.

| RNA extraction, sequencing, and cDNA synthesis
The RNA was extracted by the RNAiso Plus (total RNA extraction reagent) method according to the manufacturer's instructions (Code 9108, Takara Bio Inc., Japan). In summary, 20 mg of each sample was homogenized for 2 min at 2616 g by a tissue homogenizer (Precellys ® 24) with RNAiso plus and beads. The homogenized samples were transferred to a centrifuge tube and incubated for 5 min at room temperature (RT). Then, the samples were centrifuged at 12,000 g for 5 min at 4°C, and supernatants were collected in new tubes.
Chloroform was added to each sample in the amount of one-fifth of the amount of RNAiso plus, and the solution was mixed until it became milky. The tubes were kept at RT for 5 min and subsequently centrifuged at 12,000 g for 15 min at 4°C. The top layer, which contained RNA, was transferred to new tubes. The RNA pellets were precipitated by adding one-half of the amount of isopropanol and centrifuging at 12,000 g for 10 min at 4°C. The precipitates were washed with 75% cold ethanol and centrifuged at 7500 g for 5 min at Paired-end raw reads were generated by RTA2 and bcl2fastq2-v2-20.0 (Macrogen NGS service).
The RNA clean-up was conducted with the NucleoSpin RNA Clean-up XS kit (MACHEREY-NAGEL, Düren, Germany) following the protocol (View, 2014). First, DNA in the crude RNA extracts was digested by adding a one-tenth mixture of rDNase (Macherey-Nagel, 740963) and a reaction buffer and incubating for 10 min at 37°C.
An equal volume of buffer RCU was then added to each sample and mixed for 2 × 5 s. The solutions were added to NucleoSpin RNA XS columns with collection tubes and centrifuged for 30 s at 11,000 g.
Next, 400 µl of buffer RA3 was added to each tube and centrifuged for 30 s at 11,000 g. The flow-through was discarded, and 200 µl of buffer RA3 was introduced to the columns and centrifuged for 2 min at 11,000 g to dry the membrane containing purified RNA.
Subsequently, each column was placed in a nuclease-free collection tube (1.5 ml, supplied), and RNA was eluted in 30 µl of RNase-free water by centrifuging at 11,000 g for 30 s. The RNA concentrations were checked with the Qubit ® RNA Assay Kit in Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and the RNA quality was verified by means of the Agilent 2200 Tape Station (Agilent Technologies, Santa Clara, CA, USA). The RNA samples were stored at −80°C for further use. F I G U R E 1 Study species, Risso's dolphin (Grampus griseus) F I G U R E 2 Diagram of experimental design. CV%, coefficient of variance calculated by FPKM values; AE%, amplification efficiency based on RT-qPCR; SA%, stability analysis by statistical algorithms First-strand cDNA was synthesized for each RNA sample using the PrimerScript™ RT master mix kit (Perfect Real Time) (Code RR037A, Takara Bio Inc., Japan) with 2 µl of master mix, up to 500 ng of purified RNA, and RNase-free water in a 10 μl solution on ice followed by incubation on a thermal cycler at 37°C for 15 min and 85°C for 5 s. The Qubit was employed to check the cDNA concentrations for all samples, and the presence of cDNA was checked by normal PCR using a selected primer before RT-qPCR. The cDNA samples were stored at −20°C until further use.

| Selection of candidate reference genes using RNA-seq data and designation of primers
Raw RNA sequencing data were obtained from the company Macrogen Japan Corp. (order no. 1908JNHX-0031). To obtain highquality clean reads, raw reads were used to trim adapters and filtered by removing the poly-A sequence that contained more than five bases at the 3′ end and had a Q Phred quality score of 20 or less. (https://github.com/c2997 108/OpenP ortab lePip eline). Cufflink was used for mapped read splicing, and reads were then compared with genomic annotation information to explore new genes (Trapnell et al., 2010). Furthermore, the number of reads mapped to a particular gene was calculated with a ballgown algorithm (Frazee et al., 2015) in the pipeline. Finally, the gene count data and normalized FPKM values were obtained.
Ten candidate reference genes were selected, including novel genes, that displayed a low CV value, gene function in metabolism, and high FPKM values in all tissues. These candidate genes were chosen based on unpublished RNAseq data, which are available from the corresponding author upon request. Table 1 specifies the gene IDs, gene names, and functions of the candidate reference genes.
The expression stability of the selected genes was calculated and summarized as follows (Wang et al., 2014). First, using the FPKM values of gene expression in 30 samples, mean expression values and standard deviations (SDs) were calculated for each gene. Second, the CV of each gene was ranked according to the CV value ( Table 2).
Standard curves were acquired for each primer pair by amplification in serial dilutions, such as 1:1, 1:10, 1:100, 1:1000, and 1:10,000, for all samples (Figure 4). The correlation coefficient (R 2 ) and the PCR amplification efficiency (E) for each gene was calculated from the slope of a standard curve, E = (10 (−1/slope) − 1) × 100% (Bustin et al., 2009). Based on the closer to 100% of E and closer to 1 of R 2 values, nine reference genes were selected for further analysis of gene stability to test for all biological samples in three replicates.

| Gene stability analysis
A normality test for Ct data was conducted in the SPSS 17 statistical package (SPSS Inc, 2008) using the Kolmogorov-Smirnov test, and the normality of different genes was evaluated based on the p-value.
The statistical algorithms of Genorm (Vandesompele et al., 2002), the comparative Delta CT method (Silver et al., 2006), and RefFinder The BestKeeper algorithm was applied to calculate the SD; genes were rated on the basis of variance, which also ascribes low SD and CV values to stable genes. The Delta CT values were obtained from pairwise differences between Ct values of reference genes to determine the stability of genes by repeatability among all samples.
RefFinder considered the results of all of the above algorithms and generated a complete final ranking of the reference genes. The most stable reference genes were selected and subsequently analyzed for gene stability in 10 different tissues by using data from three replicates. Comprehensive graphs were prepared for the discussion.

| RE SULTS
In this study, a total of 30 samples (10 tissues from three individual Risso's dolphins) were employed to identify potential reference genes for the normalization of RT-qPCR data across fat tissues and to validate transcriptomic data. Based on the RNA-seq data (see Table A1), the study identified 10 highly expressed (FPKM > 1000) genes with CVs ranging from 46.14% to 175.75% (Table 2).  Table A2.

| Primer specificity, amplification efficiency, and expression profiles
The specificity of each primer pair was confirmed by 2% agarose gel electrophoresis with a single band (Figure 3a). All genes presented a single melting peak in the qPCR, which indicates specific amplification ( Figure 3b). The cDNA-free NTC (negative control) samples did not display any melting curve product. These results suggest an absence of errors in the RT-qPCR amplification. RT-qPCR performance was good according to the coefficient (R 2 ) and efficiency (E) results (see Table 3 and Figure 4). As the standard curves illustrate, nine genes showed good slopes (nearly −3.4) and efficiencies (nearly 80%-120%), including six genes (90%-110% efficiency). Therefore, these genes were selected for gene stability analysis and to finally identify a potential reference gene for toothed whales. The DCN gene was not considered for further stability analysis because of its significantly low efficiency.
The expression profiles of the nine selected genes were tested by the mean quantification cycle (Ct) values for each sample with BestKeeper analysis (

| Validation of reference genes and identification of the most stable reference genes
The Genorm statistical algorithm evaluated the stability of the selected reference genes by using the M value. The lowest value was assigned to the most stable gene, LGALS1, and FTL. The RPL27A gene had the highest M value, which suggests that it is the least stable (see Table 5 and Figure 5). NormFinder and Delta CT indicated the high-

| Stability of top three reference genes in different tissues
Based on the comprehensive gene stability values, the FTL gene was the most stable in the inner jaw fat sample (JF1), as it had the lowest stability value (1.414), and was the least stable in the outer jaw fat sample (JF2). Similar results were obtained for the GAPDH gene.
However, interestingly, the LGALS1 gene was more stable in muscle samples than in other fat tissues and was least stable in BL1 and JF2 ( Figure 6).

| D ISCUSS I ON
In this research, all of the samples were collected at the same time.
The analysis used 30 biological replicates of 10 types of tissue from three animals. All samples were preserved in RNAlater solution.
Therefore, we believe that the total RNA of the various types of tissue are mainly involved with aquatic adaptations. Liver and muscle tissues were used as controls, as they usually express many types of genes related to metabolism.
From comparative genomics analysis, it is evident that many genes are expressed in dolphins for aquatic adaptations, blubber, and fat storage (Mancia, 2018). The fat tissues of dolphins contain numerous fatty acids and fatty alcohol combinations; thus, the metabolism of these fats is still ambiguous. Considering all of these aspects, we determined candidate genes with the lowest CV percentage. Reference genes for qPCR can be identified through highthroughput transcriptomic data by next-generation sequencing (Bao et al., 2020;Gao et al., 2018;Pombo et al., 2018). All of our selected candidate reference genes exhibited a high level of expression by more than 500 FPKM values, and, according to multiple sources, their functions include regulatory and maintenance functions, metabolic processes, and catalysis in cells. Finally, this experiment was designed to detect 10 candidate reference genes, including eight uncommon reference genes and two well-known HKGs. None of the tissues or genes is dependent on sex, which presumably precludes bias in the selected candidate genes due to the selection of only males for this study.
A range of ecological and evolutionary research has applied qPCR technologies during the last decade. One study has concluded that the qPCR approach is more sensitive and reproducible than the conventional polymerase chain reaction (cPCR) method for detecting target DNA in the prey-predation relationship (Yang et al., 2020). The qPCR assay has also been used in a study of the ecology and evolution of cryptic nematodes in a marine environment (Derycke et al., 2012). Furthermore, researchers have employed RT-qPCR to study age-and sex-related opsin gene expression in guppies compared with reference genes such as COI, β-actin, and Myosin-HC (Laver & Taylor, 2011). This study indicated that the variation in gene expression depends on various factors, including age and sex, as the adults of a guppy species expressed higher LWS compared with juveniles, which relates to reproductive fitness and male coloration. Therefore, the identification of a suitable reference gene for the normalization of qPCR-based gene expression is important for ecological and evolutionary molecular studies of toothed whales. A stable and suitable reference gene is critical for correctly analyzing qPCR data (Linardic & Braybrook, 2021). quantitative or complex traits in organisms (e.g., genetic variance and heritability; Farries et al., 2019). Consequently, reference genes are also crucial for toothed whales' quantitative and population genetics.
RT-qPCR has become a widely used, accurate, and sensitive method to determine gene expression under various conditions and for many cell types in a range of organisms. Reference genes are essential for determining the level of expression of target genes, such as up-regulation or down-regulation by normalization (Freitas et al., 2019;Kozera & Rapacz, 2013). Several reference genes can be diagnosed by statistical algorithms, including NormFinder, BestKeeper, Genorm, and Delta CT (Xie et al., 2012). Therefore, the use of multiple HKGs in the normalization of genes of interest is currently a common practice. The LGALS1 and FTL genes were identified as suitable novel genes for future functional genomics analysis of toothed whales from myriad ecological, evolutionary, and metabolic perspectives. The use of two reference genes for qPCR data analysis has been recommended in other research (Linardic & Braybrook, 2021).
For the selection of stable genes, it is advised to employ more than one algorithm, which allows diverse analyses of data. Additionally, housekeeping genes also maintain normal cellular functions, which are confirmed in our experiment by the stable expression of all selected genes in muscle and liver tissues.
In our study, the LGALS1 gene, one of the galectins, was ranked first by all algorithms in the stability tests. Functions of this gene include the modulation of immune responses (Nishi et al., 2008) and RAS protein signaling (Ruvolo et al., 2020). According to the recent updates, the LGALS1 gene is highly expressed in fat tissues (NCBI,

| CON CLUS IONS
This study has contributed to research on toothed whales by systematically selecting and evaluating stable reference genes for RT-qPCR using several types of tissue, including melon fat, jaw fat, blubber fat, liver, and muscle, from Risso's dolphins. Despite many limitations, the study has effectively checked 10 highly expressed candidate reference genes for efficiency. Two novel genes, LGALS1 and FTL, have been identified as stable reference genes in toothed whales. The LGALS1 gene was the most optimal reference gene and is therefore important for future ecological and evolutionary genomics studies on toothed whales. It also seems plausible that this study provides an accurate normalization factor for expression data on the genes of interest in various tissues of toothed whales.

ACK N OWLED G M ENTS
We would like to thank the Taiji

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Selected transcriptomic data will be deposited in the DDBJ (DNA Data Bank of Japan) under the BioProject Accession PRJDB11720 and available on request from the corresponding author. All other data will also be publicly available with this article.