Identification of Appropriate Reference Genes for qRT-PCR Analysis of Heat-Stressed Mammary Epithelial Cells in Riverine Buffaloes (Bubalus bubalis)

Gene expression studies require appropriate normalization methods for proper evaluation of reference genes. To date, not many studies have been reported on the identification of suitable reference genes in buffaloes. The present study was undertaken to determine the panel of suitable reference genes in heat-stressed buffalo mammary epithelial cells (MECs). Briefly, MEC culture from buffalo mammary gland was exposed to 42 °C for one hour and subsequently allowed to recover at 37 °C for different time intervals (from 30 m to 48 h). Three different algorithms, geNorm, NormFinder, and BestKeeper softwares, were used to evaluate the stability of 16 potential reference genes from different functional classes. Our data identified RPL4, EEF1A1, and RPS23 genes to be the most appropriate reference genes that could be utilized for normalization of qPCR data in heat-stressed buffalo MECs.


Introduction
The riverine buffaloes (bubalus bubalis) exhibit signs of great distress when exposed to direct solar radiations. This is generally attributed to their specific morphological, anatomical, and behavioural characteristics [1]. The effect of heat stress on mammary epithelial cell (MEC), the major cell type in lactating mammary gland, could be one of the prime factors responsible for lower milk production in animals. Understanding of expression profile of these cells from different livestock species during different physiological stages would provide molecular basis of heat stress specific transcriptomic response of mammary gland. Recently, few efforts [2,3] have been made to unravel the transcriptional response of mammary gland to heat stress condition; still, the molecular mechanism of such responses are thought to be too complex to understand. qPCR is a common tool to determine the expression level of any target gene; for accurate quantification of expression level, there is a need to identify the appropriate reference genes under the particular experimental setup. Such approaches are helping a great deal to normalize the real time data for reliable interpretation of expression studies in different species [4][5][6][7][8]. In earlier studies, researchers have relied mostly on GAPDH, ACTB, and RS18 as suitable reference genes [9][10][11][12][13][14][15]. However, Vandesompele et al. and Bustin et al. had shown that use of single reference gene can lower the reliability of expression data and strongly advocated use of multiple reference genes for each experimental setup [10,12].
A number of studies have been conducted to identify the stably expressed candidate genes in different tissues of various livestock species such as pig, sheep, and bovines [5,[16][17][18][19]. These reports suggested that expression levels of commonly used reference genes vary considerably between different tissues and experimental conditions. Such variations reported in different studies [10,17,[19][20][21][22] indicated that there cannot be set universal reference genes for all tissues or experimental conditions, and hence there is a need for identification of tissue specific stably expressed genes [23][24][25][26][27][28][29][30][31]. Environmental heat stress affects the mammary gland that results in low milk production or truncated milk production. Mammary epithelial cells (MECs) are responsible for converting most precursors into milk constituents and transporting them to the mammary lumen, the first line that gets affected by heat stress. As MEC ares the predominant cell types in lactating mammary gland, changes in their genes expression could provide an insight of the mammary gland mechanism. The present study was therefore undertaken to identify a panel of appropriate reference genes for normalization of transcriptional data of heat-stressed buffalo MECs. A total of 16 known reference genes, namely, glyceraldehyde 3-phosphate dehydrogenase (GAPDH), betaactin (ACTB), ubiquitously expressed transcript (UXT), ribosomal protein S15A (RPS15A), beta 2-microglobulin (B2M), alpha 2-microglobulin (A2M), ribosomal protein L-4 (RPL4), ribosomal proteinS18 (RS18), ribosomal protein L-22 (RPL22), ribosomal protein S9 (RPS9), ribosomal protein S23 (RPS23), hydroxymethylbilane synthase (HMBS), hypoxanthine-guanine phosphoribosyltransferase1 (HPRT1), GTP-binding protein (GTP), eukaryotic translation elongation factor 1 alpha 1 (EEF1A1), and ubiquitin C (UBC) belonging to different functional categories, were evaluated (Table 1).

Sampling and Culturing of Mammary Epithelial Cells.
Mammary tissue was obtained from an adult riverine buffalo and was immediately transported to laboratory in Dulbecco's Modified Eagle's Medium/Ham's F12 media (DMEM/F12, 1 : 1 mix) (Hyclone, Logan, UT, USA) containing antibiotics 100 U/mL penicillin-streptomycin (Hyclone). Five grams of tissue was washed with PBS (Ca 2+ -, Mg 2+ -free) (Hyclone) for several times until the solution was pellucid and without milk. The tissue sample was cut into 1 mm 3 cubes and washed again. The smaller pieces of tissue were transferred to collagen-coated cell culture dishes (Corning, USA), containing DMEM/F12 supplemented with 10% fetal bovine serum (PAA), 100 U/mL penicillin-streptomycin (Hyclone), 5 L/mL insulin, 50 M hydrocortisone, 1 g/mL -estradiol, 5 g/mL holotransferrin, and 1 g/mL Progesterone. (Sigma-Aldrich) and incubated at 37 ∘ C, 5% CO 2 . Initially, the basal media was replaced after every 12 h and then after every 48 h with fresh media until cells were visibly spread across the bottom of the culture dish. Cells were detached with 0.25% trypsin containing 0.02% EDTA (Sigma-Aldrich) and transferred to T75 culture flasks (Corning, USA). The process was repeated up to 10th passage, and pure mammary epithelial cells were obtained for further experimentation.

Heat Stress Treatment of Mammary Epithelial Cells. The
MECs were transported to collagen-treated 12-well plates in two sets with one plate assigned as control (kept at 37 ∘ C throughout the time course) and another plate as treated (exposed to 42 ∘ C). Initially, all the cells were incubated at 37 ∘ C with 5% CO 2 for 30 m to stabilize the culture. Subsequently, the plate marked as treated was exposed to 42 ∘ C for 1 h to simulate the heat stress condition. After completion of 1 h, the cells were allowed to recover at 37 ∘ C with 5% CO 2 and harvested at different time intervals (30 m, 2 h, 4 h, 8 h, 12 h, 16 h, 24 h, and 48 h). The samples from control plates were also harvested at same time points corresponding to the treated plates. After assessing the viability, cells were transferred to chilled Trizol reagent (Invitrogen Corp., CA) and stored at −80 ∘ C until RNA extraction.  24 h, and 48 h after heat stress using ice-cold Trizol (Invitrogen Corp., CA). RNeasy Mini Kit columns providing on with column digestion by RNAse-free DNase enzyme (Qiagen, Germany) were used to remove the traces of genomic DNA. Total RNA concentration and purity were measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies). The quality check for all samples was performed using experion automated electrophoresis System (Biorad). All the extracted RNA samples were stored at −80 ∘ C and utilized within one month. cDNA was synthesized using 100 ng RNA, 1 L dT 12−18 (Invitrogen Corp. CA), 1 L 10 mmol/L dNTP mix (Invitrogen Corp., CA), 1 L random primers (Invitrogen Corp., CA), and 10 L DNase-/RNase-free water. The mixture was incubated at 65 ∘ C for 5 min and kept on ice for 3 min. A total of 6 L of master mix composed of 4.5 L 5X First-Strand Buffer, 1 L 0.1 M DTT, 0.25 L (50 U) of SuperScript III RT (Invitrogen Corp., CA), and 0.25 L of RNase inhibitor (10 U, Promega, WI) was added. The reaction was performed in an Eppendorf Gradient cycler using the program: 25 ∘ C for 5 min, 50 ∘ C for 60 min and 70 ∘ C for 15 min. cDNA was then diluted 1 : 4 (v : v) with DNase-/RNase-free water.

Primer Designing and Validation.
To facilitate the real time PCR analysis, primers were either selected from the literature or designed using Primer Express 3.0 software (Applied Biosystem) with minimum amplicon size ranging between 50 and 115 bp and limited 3 G+C content. Primer details for all genes are given in Table 2. To check the sequence specificity, primers were aligned against publicly available databases at NCBI and UCSC's Cow (Bos taurus) genome browser gateway using BLASTN. Prior to qPCR, primer specificity was further confirmed in a 20 L PCR reaction using the same protocol described for qPCR except for the final dissociation protocol. Five L of the PCR product was evaluated in 2% agarose gel stained with ethidium bromide. The accuracy of primer pairs was also ensured by the presence of a unique peak during the dissociation step at the end of qPCR.

Real Time Quantitative PCR (RT-qPCR)
. qPCR reaction was performed using LightCycler 480 instrument (Roche, Germany) in a 96-well white plate (Roche, Germany). Each reaction was comprised of 4 L diluted cDNA combined with 6 L of a mixture composed of 5 L 2X LightCycler 480 SYBR Green I master mix (Roche, Germany), 0.4 L each of 10 pmole forward and reverse primers, and 0.2 L DNase-/RNase-free water. For each gene, samples were analyzed in duplicate (technical replicates) along with 6-point relative standard curve and the nontemplate control. Following amplification conditions were used: 2 min at 50 ∘ C, 10 min at 95 ∘ C, 40 cycles of 15 s at 95 ∘ C (denaturation), and 1 min at 60 ∘ C (annealing + extension). A dissociation protocol with incremental temperatures of 95 ∘ C for 15 s plus 65 ∘ C for 15 s was used to investigate the specificity of qPCR reaction and presence of primer dimers. The qPCR expression data for each reference gene was extracted in the form of crossing points. The data was acquired using the "second derivative maximum" method as computed by the LightCycler Software 3.5 (Roche Diagnostics) and subjected for subsequent analysis.

Evaluation of Expression Stability.
The expression stability of each of the studied 16 genes was evaluated using three independent statistical applications: geneNorm [10], NormFinder [34], and Bestkeeper [35]. The geNorm, a Microsoft Excel-based application, was used to measure the expression stability as value which is based on overall pairwise comparison among the reference genes. The calculated value is inversely correlated to gene expression stability and ranks the reference genes accordingly. NormFinder, also a model-based approach, was also used to determine the optimal reference genes and the combination of two genes for a two-gene normalization factor with its corresponding stability value. Bestkeeper, another software used in the study, is based on pairwise comparisons of raw cycle threshold (Ct) values of each gene. The analysis assumes that the genes which are stably expressed should be highly correlated to each other. The 10 most stable genes as identified by geNorm analysis were considered for Bestkeeper analysis.

Results
The total RNA extracted from individual MECs samples exhibited high purity as determined by mean (±SEM) 260/280 ratio of 2.06 ± 0.014. The bioanalyzer-based RQ value of >8 also indicated sufficiently good quality of each extracted RNA. The qPCR performance for each gene in  2 Pérez et al. [32], and 3 Hernandez et al. [33]. 4 qPCR efficiencies for each primer pair were calculated from six-point standard curves using fivefold dilution series of pooled cDNA from control and heatstressed samples.    terms of coefficient of determination, ( 2 ), and efficiency of amplification ( = 10 −1/slope ) on the basis of slope of six-point standard curve are summarized in Table 2. The efficiency of PCR reactions ranged from 90.70% for ACTB to 131.37% for B2M. The characteristics of individual 16 genes based on their cycle threshold (Ct) values are shown as box whisker plot (Figure 1).

Analyses of Gene Expression Stability by GeNorm.
The value obtained under geNorm analysis ranged from 0.089 (EEF1A1 and RPL4) to 0.415 (A2M) ( Table 3). All candidate genes performed well displaying values below the accepted limit of 1.5. The genes were ranked from the most stable (lowest value) to the least stable(highest value): RPL4, (Figure 2). Additionally, pairwise variation termed as "V value" was also calculated to determine the optimal Determination of the optimal number of control genes for normalization Figure 3: Determination of optimal number of reference genes for normalization by calculation of pairwise variation ( ) of normalization factor ratios for different number of genes. number of genes necessary to calculate normalization factor. As suggested by Vandesompele et al., genes showing V value below the cutoff limit of 0.15 were selected as optimal number of genes for normalization [10]. In this analysis, we started with the two most stably expressed genes and then sequentially included less stably expressed genes. The lower the pairwise variation the better the combination of genes is. From this perspective, the V value was also calculated by adding the third and fourth less stable genes, that is, 3/ 4 and 4/ 5 combinations. The contribution of each gene to the variance of normalization factor ratio was calculated to illustrate the effect of adding or removing a particular gene from the final set of reference genes. Our result showed that a combination of the two most stable genes (RPL4 and EEF1A1) gave V value of 0.037 which is well within the acceptable limit ( Figure 3). Similar approach has been used in a number of other studies to find out the suitable housekeeping genes [16,19,20,36]. However, as the aim in such studies is to achieve an overall variation of less than 0.15 with minimum number of genes, EEF1A1 and RPL4 gene-combination could be appropriate to normalize the target genes expression in heat-stressed buffalo MECs.

Analyses of Gene Expression Stability by NormFinder.
In addition to geNorm, we also utilized NormFinder software to find out expression stability of 16 genes. The analysis identified the same set of genes (RPL4, EEF1A1, GTP, UXT, RPS23) being most stable as revealed by geNorm only with slight change in their ranking order (RPL4 > EEF1A1 > GTP > UXT > RPS23) ( Table 4, Figure 4). A2M and RPL22 genes were found to be the least stable while RPL4 as the most stable gene.

Analyses of Gene Expression Stability Using BestKeeper.
BestKeeper analysis revealed the stability order as RPL4, B2M > UXT > RPS15A, EEF1A1 >, and RPS23 with the crossing point standard deviation (SD [±CP]) value of 0.11, 0.11, 0.12, 0.13, 0.13, and 0.15, respectively (Table 5). These values were well within the acceptable range of fold change expression (<2). On the other hand, RPS9 was the least stable with SD value of 2.06. In addition, intergene relation for the 10 most stable reference gene pairs was also estimated. The highly correlated reference genes were combined into BestKeeper index, and the correlation between each reference gene and BestKeeper was analyzed. The Pearson correlation coefficient ( ), coefficient of determination ( 2 ), and the values were estimated to describe the correlation between reference genes and BestKeeper index ( Table 6). The best correlation was observed for RPS9 ( = 0.983) and RPS23 ( = 0.840) followed by EEF1A1, RPL4, HMBS, and UXT with values of 0.679, 0.567, and 0.565, respectively. The three different algorithms (geNorm, NormFinder and BestKeeper) employed in the present study to identify best suitable reference genes lead to cumulative result, and stability ranking showed almost similar trend.

Discussion
For accurate interpretation of transcriptional studies, there is a widespread realization about the importance of suitable panel of reference genes for every species/tissue under study. In this study, our focus was to identify suitable reference genes in heat-stressed buffalo MECs as no such information is available in buffaloes. The findings of the present study will be a step forward to initiate transcriptional studies in this important dairy species of Indian subcontinent. The present analysis using different statistical applications clearly revealed a panel of the most stable reference genes (RPL4, EEF1A1, and RPS23) at different time points under heat-stressed condition of buffalo MECs. Our results for ranking of the most stable reference genes utilizing these different algorithms were quite comparable, albeit not identical. Similar to our findings, RPL4 was also reported to be one of the most suitable genes in multiple tissues, namely, udder, muscle, liver, and kidney of water buffalo [37]. On the other hand, our results showed unstable expression with respect to the most commonly used reference genes, for example, ACTB and GAPDH. In the past, investigators have relied mostly on GAPDH, ACTB, and RPS18 as reference genes [9][10][11][12][13][14][15]. Historically, GAPDH and ACTB genes have been used quite frequently as single control gene in more than 90% of the studies [21] and were considered good references for many years. But in the recent past, the expression of these reference genes has been shown to be affected by experimental condition [30,38]. Several other studies have also shown variation in their transcription level [10,17,[19][20][21][22] making them unsuitable for transcriptional studies. Proper evaluation of these two reference genes in any cell type or tissue of interest is therefore mandatory for correct interpretation of qPCR results.
In the present analysis, several genes were found to fulfill the criteria of suitable normalizer gene based on values in geNorm (Figure 2), stability index in NormFinder (Figure 4), and Ct values in BestKeeper (Table 6). Overall, RPL4, EEF1A1, and RPS23 were selected to be the three most stable reference genes for normalization of gene expression data in heat-stressed buffalo MECs under in-vitro conditions. Vandesompele et al. [10] and Bower and Johnston [39] have also recommended the use of geometric average of the most stable reference genes for accurate normalization of qPCR data. As this is the first paper on testing the reference genes in heat-stressed buffalo MECs, our analyses were restricted to some of the well-known housekeeping genes reported in other livestock species. In conclusion, RPL4, EEF1A1, and RPS23 are the best reference genes for heat-stressed studies in buffalo MECs, and their geometric means would provide accurate normalization factor. For other experimental settings involving buffalo MECs, the use of these reference genes may be carefully evaluated as their expression may change in other specific experimental conditions. However, the panel reference genes identified in the present study would certainly be useful for accurate normalization of buffalo MECs expression data during heat challenge experiments.