Suitable reference gene for quantitative real-time PCR analysis of gene expression in gonadal tissues of minnow Puntius sophore under high-temperature stress

Background High ambient temperature is known to affect fish gonadal development and physiology in a variety of ways depending on the severity and duration of exposure; however, the underlying molecular mechanisms are poorly understood. Gonadal gene expression influence the gonadal development, physiology and the quality of egg/sperm produced in teleosts and the mechanistic understanding of spatio-temporal changes in the gonadal gene expression could be instrumental in controlling the fate of egg/sperm and the quality of seed produced. Real time-quantititative polymerase chain reaction (RT-qCR), is a high throughput, sensitive and reproducible methodology used for understanding gene expression patterns by measuring the relative abundance of mRNA transcripts. However, its accuracy relies upon a suitable reference gene whose expression levels remain stable across various experimental conditions. In the present study, we evaluated the suitability of ten potential reference genes to be used as internal controls in RT-qPCR analysis in gonadal tissues (ovary and testis) of minnow Puntius sophore exposed to high temperature stress for different time periods (7 days, 60 days). Expression analysis of ten different constitutively expressed genes viz. 18S ribosomal RNA (18S rRNA), beta actin (βactin), β-2 microglobulin (b2mg), eukaryotic elongation factor-1 (eef1), glyceraldehyde-3phosphate dehydrogenase (gapdh), glucose-6-phosphate dehydrogenase (g6pd), ribosomal binding protein L13 (rpl13), tubulin (tub), tata box binding protein (tbp), ubiquitin (ubi) was carried out by using RT-qPCR and the stability in their expressions were evaluated by using four different algorithms; namely, delta Ct, BestKeeper, geNorm and NormFinder. Results In ovary, eef1 was found to be the most suitable reference gene in all the algorithms used. In testis, b2mg was found to be the most suitable reference gene in delta Ct, BestKeeper, NormFinder analysis while tbp and eef1 were found to be the most suitable reference genes in geNorm analysis. Conclusions In conclusion, eef1 and b2mg were found to be the most suitable reference genes in ovary and testis, respectively, of Puntius sophore exposed to high temperature stress, and could be used as internal controls for gene expression analysis in gonadal tissues of Puntius sophore.


Background
Gonadal gene expression influence sex determination, growth, and survivability of gametes in teleosts [1]. The ability to fully control sexual maturation, spawning and production of high quality seeds are some of the primary requirements for the successful development of aquaculture and in depth understanding of gene expression mechanism and the ability to manipulate it with high precision are important aspects of the reproductive physiology and biochemistry [2]. The ability to regulate sexual differentiation, maturation and reproduction provides fish breeders control over breeding processes thereby affecting largescale industrial production [3]. The quality of fish gametes are determined by expression of an array of genes influenced by several factors (age, management, feeding, chemical and physical factors, water quality, etc.) that have an impact on the survivability of embryos, larvae and/or fry [4]. In recent times, fish gonadal gene expression studies are being used to assess the sperm/egg quality. Thus, gonadal gene expression studies have profound implications on the fisheries and aquaculture.
Gonadal gene expression in fish is affected by environmental factors such as water temperature, pH, dissolved oxygen content etc. [3]. Temperature is an important environmental factor that influences the sex ratio, quality of eggs and can have long lasting effects [5,6] and increase in global temperature has been a cause of concern in the climate change regime [7,8]. As fish exhibit enormous diversity in their morphology, habitat occupancy, and biology, it is necessary to extensively study the impacts of external temperature on the overall physiology and in the reproductive physiology in particular to assess species that will adapt to the changing environmental temperature and studying the gene expression pattern of fish species would be helpful in carrying out such assessments.
Gene expression study using quantitative real time polymerase chain reaction (RT-qCR), has become a high throughput, sensitive and reproducible methodology for measuring the relative abundance of mRNA transcripts and has been one of the most widely used method of gene expression analysis due to its accuracy, fastness, efficiency, sensitivity, reproducibility and broad dynamic range [9,10]. However, its accuracy relies upon a suitable reference gene whose expression levels remain stable across various experimental conditions. In recent times, extensive effort are being made to find suitable reference genes that maintain constant expression levels under different experimental conditions [10][11][12]. However, no single gene has been found to have constant level of expression in all fish species, different tissues and experimental conditions. For instance, 18S rRNA, b2m and elfa have been found to be suitable reference genes in gonadal tissues of zebrafish (Danio rerio) at different developmental stages whereas rpl7 has been found to be the most suitable reference gene in gonadal tissues of Oryzias latipes [13,14]. Therefore, it is necessary to identify species and organ specific reference genes that could be used for studies on transcriptomic responses of organisms under various conditions. In this context, the objective of the present study was to identify suitable reference gene(s) that can be used for RT-qPCR analysis in gonadal tissues (both ovary and testis) of minnow Puntius sophore exposed to high temperature stress for different time periods, both short term (7 days) and for long term (60 days). Puntius sophore is a close relative of Danio rerio, a widely used animal model and is a nutrient-dense small indigenous fish that is especially rich in micronutrients. It is available in most of the freshwater aquatic ecosystems in many tropical and subtropical countries like India, China, Bangladesh, Thailand, and Myanmar [15,16] and owing to its high nutritive value, attempts are being made to bring it into aquaculture fold [17,18]. The study identified suitable candidate reference genes in gonadal tissues of Puntius sophore which could be useful in evaluating different aspects of its reproductive physiology.

Expression profiles of reference genes
Real-time reverse transcription quantitative polymerase chain reaction (RT-qPCR) is one of the most reliable methods of measuring transcript abundance, representing the spatio-temporal and clinical changes in gene expression. However, the results of RT-qPCR are very much dependent on the normalization of data against stable endogenous controls, often referred as housekeeping genes. Each of the sample was measured in triplicate for each run and nine independent replicates were performed for each set of experimental set up (Control, 7 days and 60 days heat exposed groups).
The expression levels of the candidate reference genes were calculated with respect to mean real-time PCR quantitative cycle (Cq) value. The Cq values of the candidate reference genes in ovary and testis tissues are presented in Table 1. In ovary tissue, 18S rRNA, βactin were found to be the highly expressed genes. gapdh and tbp were found to be least expressed genes in ovary (Cq values ranging between 31 and 36 and 33-35 respectively. Similarly in testis, βactin, rbpl13 and 18S rRNA were found to be the highly expressed genes (Cq values ranging from 10 to 16, 15-23 and 15-24 respectively) whereas g6pd, ubi and gapdh were the least expressed genes.
Gene expression stability of tissue specific reference genes due to thermal stress In order to find out the suitable reference genes for normalization of RT-qPCR data, the stability of expression of ten candidate reference genes, across the control and heat stressed groups in gonadal tissues (ovary and testis) were measured and validated by using four computational methods, such as: comparative ΔCt, geNorm, NormFinder and BestKeeper [19][20][21][22]. These four programs, use different approaches in evaluating stability in gene expression and thus are used in combination to increase confidence over the results obtained.

Comparative ΔCt method
Comparative delta Ct method compares relative expression of a pair of genes; if the ΔCt value between the two genes remains constant when analyzed in different samples, it means both the genes are stably expressed among those samples. However, if the ΔCt fluctuates, then 1 or both genes are variably expressed. Then expression profiles of other genes are taken into consideration to find out the most stably expressed gene [19]. In the present study, eef1 and b2mg were found to be most stably expressed genes in ovary and testis respectively with lowest average standard deviation values. The ranking of other candidate genes for suitability as a reference gene has been provided in Fig. 1.

geNorm analysis
The geNorm software is employed as a mean for determining the expression stability value (M) for each candidate gene on the basis of the average pair wise variation between all genes analysed [10,20]. The gene with the lowest M value is considered to have the most stable expression, while with the highest M value has the least stable expression. All the tested genes were found to be having high stability and stability values were within the geNorm cut off value (1.6) [23] except for βactin, rbl1 and ubi. In testis b2mg, tbp and eef1 were the genes which had stability values within the cutoff limit. geNorm analysis indicated that, the eef1 and tubulin were the most stable genes in ovary while b2mg, tbp were the most stable genes in testis (Fig. 2).

NormFinder analysis
The NormFinder estimates the variance in the Cq value dataset, ranks the genes according to their stability. The more suitable gene expression is indicated by lower average expression stability values [21]. As per the NormFinder  1 Stability in expression of 10 reference genes in a ovary and b testis of Puntius sophore as determined by delta Ct method. eef1 and b2mg were found to be the most suitable reference genes in ovary and testis, respectively analysis, the order of stability in ovary was eef1 > tubulin > b2mg > tbp > g6pd > 18S rRNA > gapdh > βactin > rbpl13 > ubi. Similarly, in testis the order was tbp > βactin > g6pd > 18S rRNA > eef1 > gapdh > rbl13 > tubulin > ubi (Fig. 3).

Recommended comprehensive ranking
Combining the results of the analysis by four different algorithms, a comprehensive ranking was generated by the online tool available at http://leonxie.esy.es/RefFinder/. According to this ranking eef1 was found to be the most suitable reference gene in ovary while b2mg was found to be the most suitable reference gene in testis (Table 3).

Discussion
Gonadal gene expression studies in fish are being employed to assess the quality of egg, spermatid, which in turn has profound implication in fish culture and management systems [5]. The control of gamete quality is of great importance; availability of good quality male and female gametes is necessary to close the lifecycle of a species and obtain subsequent generations and poor egg quality can lead to several types of problems [2]. It is also being used to understand developmental processes, and to assess the environmental effects on reproductive physiology of fish [24]. The RT-qPCR has emerged as a reliable, reproducible and highly sensitive technique for gene expression analysis. However, the reliability of the results largely depend on the selection of reference gene used for normalization of gene expression across different experimental conditions. In this regard, we have evaluated the stability of ten commonly used reference genes using Excel based algorithms deltCt, geNorm, NormFinder and BestKeeper. The results of analysis in all the algorithms used were very similar i. e. in Fig. 2 geNorm stability values of the 10 reference genes in a ovary and b testis of Puntius sophore. In ovary, eef1 and tubulin were the most stable genes whereas in testis b2mg and tbp were the most suitable reference genes. geNorm stability values (M) for each genes was calculated as the average pairwise variation for that gene with all other tested reference genes using the genorm software and genes were ranked in the increasing order of stability values Fig. 3 Normfinder stability values of the 10 reference genes in a ovary and b testis of Puntius sophore. In ovary, eef1 was the most stable gene whereas in testis b2mg and tbp were the most suitable reference genes as determined by the Normfinder algorithm. Normfinder stability values were calculated by combining the inter and intra group variations determined using the excel-Normfinder add-in. The genes were then ranked in the increasing order of stability values ovary eef1 was found to be the most stable gene in all the algorithms used. eef1 performs varieties of functions during cell growth and proliferation, and plays key role in cytoskeleton organization, mitotic apparatus formation, and signal transduction [25][26][27]. It also mediates the recruitment of aminoacyl-tRNA to the A-site of the 80S ribosome during protein synthesis and thus is expressed ubiquitously in all cells [25,26]. Similar to the present study, elf1 has been found to be a suitable reference gene in gonadal tissues of Zebrafish exposed to endocrine disruptors [28] and elfaα has been found to be a suitable reference gene in gonads of rice field eel Monopterus albus in different stages of development [29]. Another translation related protein eif (eukaryotic initiation factor) has been found to be one of the most suitable reference genes in ovarian tissues of Procambarus clarkii [30].
Similarly in testis, b2mg was found to be the most stable gene in deltCt, NormFinder and BestKeeper  algorithms. However, in geNorm analysis tbp and eef1 were found to be the most stable genes. b2mg is a component of MHC class I molecules, which are present on all nucleated cells and it has been found to be a suitable reference gene in mouse hematopoetic stem cells [31], and in K562 cells and leucocytes of normal individuals as well as of malignoma patients [32] and different tissues of Japanese flounder (Paralichthys olivaceus) [33]. However, in some other studies involving human colorectal cancer cells [34] and human reticulocyte cells [19] it has been reported not to be a suitable reference gene. The statistical algorithms used in the present study have distinguished features and thus there could be discrepancies between the results obtained. Therefore, these algorithms are used in combination to increase confidence limits on the results obtained from them. However, in the present study, we obtained same results from all the algorithms used except for the geNorm analysis in testis in which tbp was found to be the most stable gene in contrast to b2mg which was found to be the most stable gene in all the algorithms used.

Conclusion
In conclusion, eef1 and b2mg were found to be the most suitable reference genes in ovary and testis of Puntius sophore exposed to high temperature stress and can be used as internal controls for gene expression analysis in Puntius sophore. Identification of these reference gene would be helpful in carrying out gene expression analysis in gonadal tissues of Puntius sophore.

Ethical statement
The study, including sample collection, experimentation and sacrifice met the ethical guidelines, including adherence to the legal requirements of the study country. The study was approved by the Institute Animal Ethics Committee (IEAC) of ICAR-CIFRI vide approval no. CIFRI/ IEAC-16-17/03.

Sample collection after short and long term thermal exposure
Puntius sophore (weight: 4.22 ± 0.5 g, length: 6.46 ± 0.56 cm) were collected from aquaculture ponds and were acclimatized under laboratory conditions for 30 days. Fishes were kept in glass aquariums of 30 l capacity with temperature control systems. Fish were fed once daily @ 2% of body weight by providing laboratory prepared fish diet.
The fishes were randomly assigned to three experimental groups; control (acclimation temperature 27 ± 0.2°C), fish exposed to heat stress (36°C) for 7 days and fish exposed to heat stress for 60 days. The temperature in the heat stressed groups were initially increased at rate of 2°C/h using the temperature control system attached with the aquarium and maintained at 36 ± 0.4°C for the specified time period (7/60 days). 3/4th of the water was changed daily in each aquarium and were replenished with water of specified temperatures (27 ± 0.2°C in control groups and 36 ± 0.4°C for the heat stressed groups). After completion of the exposure period, fishes were dissected (after euthenization with Tricaine, MS-222, 200 mg/l) and ovary, testis tissues were collected in RNA later (R0901, Sigma) and were kept at −40°C for further use.

RNA extraction and cDNA synthesis
RNA was extracted from tissue samples (9 samples from each experimental group, weighing approximately 50-70 mg of tissues) using RiboZol (Himedia, India) according to the manufacturer's protocol. Following isolation of RNA, the concentration of each RNA sample was measured by Bioanalyzer (Agilent 2100 Bioanalyzer, USA). RNA integrity was confirmed by denaturing agaroseformaldehyde gel (1% w/v agarose, 16% formaldehyde) electrophoresis showing two sharp and intense bands for 18S & 28S ribosomal RNA. RNA samples were treated with the DNase1 (NEB, UK) as per the manufacturer's recommended protocol to remove DNA contamination. 1 μg of DNase treated total RNA was reverse transcribed using M-MLV reverse transcriptase (NEB, UK) according to manufacturer's protocols. The cDNA so obtained were diluted 10 fold with nuclease free water for further amplification in RT-qPCR.

Primer synthesis
Ten candidate reference genes were evaluated in the present study such as βactin, gapdh, rbpl13, tub, 18S rRNA, tbp, b2mg, eef1, g6pd, ubi. Primers for βactin, gapdh, rbpl13, tub, 18S rRNA, and tbp were adapted from Purohit et al. (2016) [10]. Primers for b2mg, eef1, g6pd were adapted from McCurley and Callard [13] and those for ubi were adapted from Olsvik et al. 2008 [35]. Primers with specific nucleotide sequences were obtained from Integrated DNA Technologies, USA. The specificity of the primer sets were confirmed by the presence of a single band of appropriate size obtained after PCR amplification.

RT-qPCR
The real time PCR amplifications were carried out using SYBR Green detection chemistry. cDNA were run in triplicates on a 96 well reaction plates with the CFX Connect real time PCR (Bio-Rad, UK). 20 μl of reaction mixture containing 10 μl of VeriQuest SYBR Green Mix (BioRad, UK), 1.0 μl of each 10 μM of primers and 5 μl of diluted cDNA as template and 3 μl RNase/DNase free sterile water (Thermo Scientific, USA).
The following amplification program were used in all RT-qPCR reactions: 50 0 C for 2 min, 95°C for 10 min, 40 cycles of 15 s 95°C annealing, extension for 45 s at optimized temperatures for specific candidate genes ( Table 4). The specificity of each amplification reaction was verified by a melting curve analysis after 40 cycles. No template controls (NTC) were included for each primer pair to avoid possible contamination of assay reagents and one negative control (without reverse transcription) per sample to identify residual gDNA after DNase digestion. PCR efficiency of the reference genes were determined by a standard curve analysis of cDNA samples according to the MIQE guideline [36]. A series of 10 fold dilution of three replicates of cDNA, were made to determine the gene specific PCR amplification efficiency for each primer pair used in RT-qPCR experiments. Standard curve were derived from the E values by the formula E = 10 -1/slope . The mean efficiency values were obtained for each tissue samples and were used to adjust the quantitative cycle (Ct) values for further analysis.

Data analysis
Stability of mRNA expression of each candidate genes was statistically analyzed with four widely used software programs such as comparative Ct, geNorm, NormFinder and BestKeeper [19][20][21][22]. geNorm stability values (M) for each genes was calculated as the average pairwise variation for that gene with all other tested reference genes using the genorm software and genes were ranked in the increasing order of stability values. The Normfinder stability values were calculated by combining the inter and intra group variations determined using the excel-Normfinder add-in. The genes were then ranked in the increasing order of stability values. Descrptive statistics was carried out by SPSS 16.0 and Bestkeeper excel add-in. The comprehensive ranking of different genes were evaluated using the online resource available at http://leonxie.esy.es/RefFinder/ which assigns an appropriate weight to individual genes and calculates the geometric mean of their weights obtained for the overall final ranking [37].