ISRN BiotechnologyVolume 2013 (2013), Article ID 735053, 9 pageshttp://dx.doi.org/10.5402/2013/735053

Research Article

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

Neha Kapila,1,2 Amit Kishore,¹ Monika Sodhi,¹ Ankita Sharma,¹ Pawan Kumar,² A. K. Mohanty,³ Tanushri Jerath,¹ and M. Mukesh¹

¹DNA Fingerprinting Unit, National Bureau of Animal Genetic Resources, Karnal, Haryana 132001, India² Biotechnology Division, Singhania University, Jhunjhnu, Rajasthan 333515, India³Animal Biotechnology Centre, National Dairy Research Institute, Karnal, Haryana 132001, India

Received 3 November 2012; Accepted 23 November 2012

Academic Editors: R. Greiner and J. Jia

Copyright © 2013 Neha Kapila et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. 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–8]. In earlier studies, researchers have relied mostly on GAPDH, ACTB, and RS18 as suitable reference genes [9–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–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–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–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), beta-actin (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).

Table 1: General description, cellular localization, and biological functions for set of candidate genes evaluated as reference genes in the present study.

Gene symbol
Description
Cellular localization
Biological function/component

A2M
Alpha 2-microglobulin
Cytoplasm
A protease inhibitor and cytokine transporter

ACTB
β-actin
Cytoplasm
Cytoskeletal structural protein, nucleotide, and ATP binding

B2M
Beta 2-microglobulin
Golgi membrane, plasma membrane, early endosome membrane, extracellular region
Cytoskeletal protein, immune response, protein binding

EEF1A1
Eukaryotic translation elongation factor 1 alpha 1
Cytoplasm
Translation elongation factor activity

GAPDH
Glyceraldehyde 3-phosphate dehydrogenase
Plasma membrane
Glycolytic enzyme, oxidoreductase in glycolysis and gluconeogenesis

GTP
GTP-binding protein
Cytoplasm and nucleus
Biogenesis of the 60S ribosomal subunit

HMBS
Hydroxymethylbilane synthase
Cytoplasm
Heme synthesis, porphyrin metabolism, transferase activity

HPRT1
Hypoxanthine phosphoribosyltransferase
Cytoplasm
Purine synthesis in salvage pathway

RPL22
Ribosomal protein L22
Cytoplasm
Component of the 60S subunit and encodes a ribosomal protein

RPL4
Ribosomal protein L4
Cytoplasm
Component of the 60S subunit and encodes a ribosomal protein

RPS15
Ribosomal protein S15
Cytoplasm
Protein synthesis/40S subribosome

RPS18
Ribosomal protein S18
Cytoplasm
Component of the 40S ribosome

RPS23
Ribosomal protein S23
Cytoplasm
Protein synthesis/40S subribosome

RPS9
Ribosomal protein S9
Cytoplasm
Protein synthesis/40S subribosome

UBC
Ubiquitin C
Cytoplasm and nucleus
Protein degradation

UXT
Ubiquitously expressed transcript
Cytoplasm and nucleus
Transcriptional activation, ATP binding, microtubule binding, unfolded protein binding

2. Materials and Methods

2.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²⁺-, Mg²⁺-free) (Hyclone) for several times until the solution was pellucid and without milk. The tissue sample was cut into 1 mm³ 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₂. 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.

2.2. 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₂ 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₂ 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.

2.3. Isolation of Total RNA and cDNA Synthesis

Total RNA was extracted from MECs harvested at 30 m, 2 h, 4 h, 8 h, 12 h, 16 h, 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 dT12–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.

2.4. 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.

Table 2: Gene name, GenBank accession numbers, primer sequences, annealing temperature (), amplicon length, and PCR efficiency for the studied genes.

Genes
Accession number
Primers 5′-3′ (forward, reverse)

Amplicon size (bp)
PCR efficiency (%)⁴

A2M 1
CR452243
CACCCAGGACACAGTGGTAGCCCCTGAAGACTGGATGGTCAC
60°C
103
97.20

ACTB 1
AY141970
GCGTGGCTACAGCTTCACCTTGATGTCACGGACGATTTC
60°C
56
90.70

B2M
NM_173893
CTGCTATGTGTATGGGTTCCGGAGTGAACTCAGCGTG
60°C
101
131.37

EEF1A1
BC105315
CATCCCAGGCTGACTGTGCTGTAAGCCAAAAGGGCATGC
60°C
101
115.50

GAPDH 1
BC102589
TGGAAAGGCCATCACCATCTCCCACTTGATGTTGGCAG
60°C
60
96.64

GTP 1
AK074976
CTTGGAATCCGAGGAGCCACCTGGGATCACCAGAGCTGT
60°C
101
102.91

HMBS 2
BC112573.1
CTTTGGAGAGGAATGAAGTGGAATGGTGAAGCCAGGAGGAA
60°C
80
101.4

HPRTI 3
BC103248
GAGAAGTCCGAGTTGAGTTTGGAAGGCTCGTAGTGCAAATGAAGAGT
64°C
190
99.70

RPL22 1
BC114880
AAGATGGCGCCGAAGAAAG TTTCCCGAATCAAAAATTCCA
60°C
101
102.73

RPL4
NM_001014894
TTGGAAACATGTGTCGTGGGGCAGATGGCGTATCGCTTCT
60°C
101
93.72

RPS15A
BC108231
GAATGGTGCGCATGAATGTCGACTTTGGAGCACGGCCTAA
60°C
101
119.08

RS18
DQ222453.1
TTGCCTTTGCCATCACTGCTTGTATTGGCGTGGATT
60°C
158
129.97

RPS23 1
BC102049
CCCAATGATGGTTGCTTGAACGGACTCCAGGAATGTCACC
60°C
101
115.5

RPS9 1
DT860044
CCTCGACCAAGAGCTGAAGCCTCCAGACCTCACGTTTGTTC
60°C
54
96.45

UBC 2
BE668033
TCCCTACCTGCATCATGTGCGGAATTTGGGCCAGTGCTC
59°C
71
104.71

UXT 1
CR452243
TGTGGCCCTTGGATATGGTTGGTTGTCGCTGAGCTCTGTG
60°C
101
99.74

1Bionaz and Loor [19],²Pérez et al. [32], and ³Hernandez et al. [33].

4qPCR efficiencies for each primer pair were calculated from six-point standard curves using fivefold dilution series of pooled cDNA from control and heat-stressed samples.

2.5. 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.

2.6. 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.

3. Results

The total RNA extracted from individual MECs samples exhibited high purity as determined by mean (±SEM) ratio of . The bioanalyzer-based RQ value of >8 also indicated sufficiently good quality of each extracted RNA. The qPCR performance for each gene in terms of coefficient of determination, , and efficiency of amplification 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).

Figure 1: Overall expression pattern of 16 genes evaluated in heat-stressed MECs. The data is represented as qPCR cycle threshold (Ct) values of each gene in the box and whisker diagram. The median is shown as a dashed line across the box. The boxes represent median and 1st and 3rd quartiles ranges, while whiskers indicate maximum and minimum values.

3.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, EEF1A1 > RPS23 > GTP > UXT > RPS9 > RPS15A > HMBS > B2M > HPRT1 > UBC > RPS18 > GAPDH > ACTB > RPL22 > A2M (Figure 2). Additionally, pairwise variation termed as “V value” was also calculated to determine the optimal 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, and 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.

Table 3: Ranking of reference genes based on their expression stability.

Ranking order
Genes
value

1
EEF1A1
0.089

2
RPL4
0.089

3
RPS23
0.109

4
GTP
0.130

5
UXT
0.142

6
RPS9
0.151

7
RPS15A
0.157

8
HMBS
0.164

9
B2M
0.174

10
HPRT1
0.182

11
UBC
0.201

12
RS18
0.219

13
GAPDH
0.236

14
ACTB
0.251

15
RPL22
0.305

16
A2M
0.415

Figure 2: Average expression stability measures ( value) for 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.

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.

3.2. 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.

Table 4: Standard deviation (SD) and accumulated standard deviation (Acc. SD) for each gene analyzed through NormFinder.

Genes
SD
Acc. SD

RPL4
0.0115
0.0115

EEF1A1
0.0446
0.023

GTP
0.1029
0.0376

UXT
0.1056
0.0386

RPS23
0.1397
0.0416

RPS9
0.1457
0.0424

B2M
0.1578
0.0427

RPS15A
0.176
0.0434

UBC
0.1762
0.0432

HMBS
0.2081
0.0441

HPRT1
0.2147
0.0446

GAPDH
0.2591
0.0463

ACTB
0.2735
0.0476

RS18
0.3278
0.05

RPL22
0.6443
0.0634

A2M
1.1554
0.0936

Figure 4: Bar plot showing gene variability in 16 candidate reference genes by NormFinder.

3.3. 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 , and the values were estimated to describe the correlation between reference genes and BestKeeper index (Table 6). The best correlation was observed for RPS9 and RPS23 followed by EEF1A1, RPL4, HMBS, and UXT with values of 0.679, 0.567, and 0.565, respectively.

Table 5: Parameters based cycle point (CP) values for the 10 most stable housekeeping genes.

EEF1A1
RPL4
RPS23
GTP
UXT
RPS9
RPS15A
HMBS
B2M
HPRT1

9
9
9
9
9
9
9
9
9
9

GM [CP]
23.64
22.02
19.81
24.94
25.90
19.15
19.15
25.81
22.07
22.57

AM [CP]
23.64
22.02
19.81
24.94
25.90
19.54
19.15
25.81
22.07
22.57

min [CP]
23.31
21.79
19.42
24.54
25.70
10.27
18.80
25.31
21.86
22.14

max [CP]
23.91
22.30
20.02
25.27
26.30
20.87
19.43
26.06
22.37
23.00

SD [±CP]
0.13
0.11
0.15
0.16
0.12
2.06
0.13
0.18
0.11
0.20

CV [%CP]
0.57
0.52
0.76
0.64
0.47
10.54
0.67
0.71
0.49
0.89

min [x-fold]
−1.26
−1.18
−1.31
−1.32
−1.15
−470.28
−1.28
−1.42
−1.16
−1.34

max [x-fold]
1.21
1.21
1.16
1.25
1.31
3.30
1.21
1.19
1.23
1.35

SD [±x-fold]
1.10
1.08
1.11
1.12
1.09
4.17
1.09
1.13
1.08
1.15

: number of samples; GM [CP]: geometric mean of cycling point; AM [CP]: arithmetic mean of CP; min [CP] and max [CP]: extreme values of CP; SD [±CP]: standard deviation of the CP; CV [%CP]: coefficient of variation expressed as a percentage on the CP values; min [x-fold] and max [x-fold]: extreme values of expression levels expressed as absolute x-fold over or under coefficient; SD [±x-fold]: standard deviation of the absolute regulation coefficients.

Table 6: Repeated pairwise correlation analysis among genes and with BestKeeper index (BI).

Genes
versus
EEF1A1
RPL4
RPS23
GTP
UXT
RPS9
RPS15A
HMBS
B2M
HPRT1

RPL4
R
0.893
—
—
—
—
—
—
—
—
—

P value
0.001
—
—
—
—
—
—
—
—
—

RPS23
R
0.841
0.800
—
—
—
—
—
—
—
—

P value
0.005
0.010
—
—
—
—
—
—
—
—

GTP
R
0.790
0.592
0.718
—
—
—
—
—
—
—

P value
0.011
0.094
0.029
—
—
—
—
—
—
—

UXT
R
0.776
0.426
0.519
0.793
—
—
—
—
—
—

P value
0.014
0.251
0.151
0.011
—
—
—
—
—
—

RPS9
R
0.712
0.608
0.766
0.351
0.464
—
—
—
—
—

P value
0.032
0.082
0.016
0.355
0.207
—
—
—
—
—

RPS15A
R
0.507
0.250
0.592
0.788
0.634
0.363
—
—
—
—

P value
0.163
0.518
0.094
0.012
0.067
0.337
—
—
—
—

HMBS
R
0.658
0.630
0.775
0.791
0.413
0.430
0.812
—
—
—

P value
0.054
0.069
0.014
0.011
0.269
0.248
0.008
—
—
—

B2M
R
0.206
0.056
0.072
0.179
0.407
−0.042
0.384
0.168
—
—

P value
0.593
0.885
0.855
0.646
0.277
0.915
0.308
0.666
—
—

HPRT1
R
0.382
0.305
0.273
0.639
0.368
0.146
0.676
0.743
−0.013
—

P value
0.312
0.423
0.478
0.064
0.329
0.708
0.045
0.022
0.977
—

BestKeeper Index
R
0.804
0.679
0.840
0.503
0.565
0.983
0.497
0.567
0.035
0.273

P value
0.009
0.044
0.005
0.167
0.113
0.001
0.172
0.112
0.931
0.478

0.646
0.461
0.706
0.253
0.319
0.966
0.247
0.321
0.001
0.075

: Pearson correlation coefficient; : coefficient of determination.

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.

4. 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–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–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.

Acknowledgments

The work was supported by the Indian Council of Agricultural Research, New Delhi, under the National Fellow Scheme. The authors duly acknowledge the Director of NBAGR for providing the research facilities for this work.

References

I. F. M. Marai and A. A. M. Haeeb, “Buffalo's biological functions as affected by heat stress—a review,” Livestock Science, vol. 127, no. 2-3, pp. 89–109, 2010.

R. J. Collier, G. E. Dahl, and M. J. Vanbaale, “Major advances associated with environmental effects on dairy cattle,” Journal of Dairy Science, vol. 89, no. 4, pp. 1244–1253, 2006.

S. Tao, J. W. Bubolz, B. C. do Amaral et al., “Effect of heat stress during the dry period on mammary gland development,” Journal of Dairy Science, vol. 94, pp. 5976–5986, 2011.

A. de Ketelaere, K. Goossens, L. Peelman, and C. Burvenich, “Technical note: validation of internal control genes for gene expression analysis in bovine polymorphonuclear leukocytes,” Journal of Dairy Science, vol. 89, no. 10, pp. 4066–4069, 2006.

S. Tramontana, M. Bionaz, A. Sharma et al., “Internal controls for quantitative polymerase chain reaction of swine mammary glands during pregnancy and lactation,” Journal of Dairy Science, vol. 91, no. 8, pp. 3057–3066, 2008.

L. Castigliegoa, A. Armania, X. Lia, G. Grifonib, D. Gianfaldonia, and A. Guidi, “Selecting reference genes in the white blood cells of buffalos treated with recombinant growth hormone,” Analytical Biochemistry, vol. 403, pp. 120–122, 2010.

C. R. Galiveti, T. S. Rozhdestvensky, J. Brosius, H. Lehrach, and Z. Konthur, “Application of housekeeping npcRNAs for quantitative expression analysis of human transcriptome by real-time PCR,” RNA, vol. 16, no. 2, pp. 450–461, 2010.

M. Sodhi, A. Kishore, K. Khate et al., “Evaluating suitable internal control genes for transcriptional studies in heat-stressed mammary explants of buffaloes,” Journal of Animal Breeding and Genetics. In press

J. Huggett, K. Dheda, S. Bustin, and A. Zumla, “Real-time RT-PCR normalisation; strategies and considerations,” Genes and Immunity, vol. 6, no. 4, pp. 279–284, 2005.

J. Vandesompele, K. de Preter, F. Pattyn et al., “Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes,” Genome Biology, vol. 3, no. 7, p. RESEARCH0034, 2002.

E. Deindl, K. Boengler, N. van Royen, and W. Schaper, “Differential expression of GAPDH and β-actin in growing collateral arteries,” Molecular and Cellular Biochemistry, vol. 236, no. 1-2, pp. 139–146, 2002.

S. A. Bustin, “Absolute quantification of mrna using real-time reverse transcription polymerase chain reaction assays,” Journal of Molecular Endocrinology, vol. 25, no. 2, pp. 169–193, 2000.

S. Musters, K. Coughlan, T. McFadden, R. Maple, T. Mulvey, and K. Plaut, “Exogenous TGF-β1 promotes stromal development in the heifer mammary gland,” Journal of Dairy Science, vol. 87, no. 4, pp. 896–904, 2004.

S. R. Thorn, S. Purup, W. S. Cohick, M. Vestergaard, K. Sejrsen, and Y. R. Boisclair, “Leptin does not act directly on mammary epithelial cells in prepubertal dairy heifers,” Journal of Dairy Science, vol. 89, no. 5, pp. 1467–1477, 2006.

S. R. Thorn, M. J. Meyer, M. E. van Amburgh, and Y. R. Boisclair, “Effect of estrogen on leptin and expression of leptin receptor transcripts in prepubertal dairy heifers,” Journal of Dairy Science, vol. 90, no. 8, pp. 3742–3750, 2007.

D. Garcia-Crespo, R. A. Juste, and A. Hurtado, “Selection of ovine housekeeping genes for normalisation by real-time RT-PCR; analysis of PrP gene expression and genetic susceptibility to scrapie,” BMC Veterinary Research, vol. 1, article 3, 2005.

N. A. Janovick-Guretzky, H. M. Dann, D. B. Carlson, M. R. Murphy, J. J. Loor, and J. K. Drackley, “Housekeeping gene expression in bovine liver is affected by physiological state, feed intake, and dietary treatment,” Journal of Dairy Science, vol. 90, no. 5, pp. 2246–2252, 2007.

K. Svobodová, K. Bílek, and A. Knoll, “Verification of reference genes for relative quantification of gene expression by real-time reverse transcription PCR in the pig,” Journal of Applied Genetics, vol. 49, no. 3, pp. 263–265, 2008.

M. Bionaz and J. J. Loor, “Identification of reference genes for quantitative real-time PCR in the bovine mammary gland during the lactation cycle,” Physiological Genomics, vol. 29, no. 3, pp. 312–319, 2007.

B. S. Stamova, M. Apperson, W. L. Walker et al., “Identification and validation of suitable endogenous reference genes for gene expression studies in human peripheral blood,” BMC Medical Genomics, vol. 2, article 49, 2009.

T. Suzuki, P. J. Higgins, and D. R. Crawford, “Control selection for RNA quantitation,” BioTechniques, vol. 29, no. 2, pp. 332–337, 2000.

Z. Tong, Z. Gao, F. Wang, J. Zhou, and Z. Zhang, “Selection of reliable reference genes for gene expression studies in peach using real-time PCR,” BMC Molecular Biology, vol. 10, article 71, 2009.

P. D. Lee, R. Sladek, C. M. T. Greenwood, and T. J. Hudson, “Control genes and variability: absence of ubiquitous reference transcripts in diverse mammalian expression studies,” Genome Research, vol. 12, no. 2, pp. 292–297, 2002.

A. Pachot, J. L. Blond, B. Mougin, and P. Miossec, “Peptidylpropyl isomerase B (PPIB): a suitable reference gene for mRNA quantification in peripheral whole blood,” Journal of Biotechnology, vol. 114, no. 1-2, pp. 121–124, 2004.

J. Biederman, J. Yee, and P. Cortes, “Validation of internal control genes for gene expression analysis in diabetic glomerulosclerosis,” Kidney International, vol. 66, no. 6, pp. 2308–2314, 2004.

B. Etschmann, B. Wilcken, K. Stoevesand, A. von der Schulenburg, and A. Sterner-Kock, “Selection of reference genes for quantitative real-time PCR analysis in canine mammary tumors using the GeNorm algorithm,” Veterinary Pathology, vol. 43, no. 6, pp. 934–942, 2006.

F. J. Hoerndli, M. Toigo, A. Schild, J. Götz, and P. J. Day, “Reference genes identified in SH-SY5Y cells using custom-made gene arrays with validation by quantitative polymerase chain reaction,” Analytical Biochemistry, vol. 335, no. 1, pp. 30–41, 2004.

M. Jung, A. Ramankulov, J. Roigas et al., “In search of suitable reference genes for gene expression studies of human renal cell carcinoma by real-time PCR,” BMC Molecular Biology, vol. 8, article 47, 2007.

F. Ohl, M. Jung, A. Radonić, M. Sachs, S. A. Loening, and K. Jung, “Identification and validation of suitable endogenous reference genes for gene expression studies of human bladder cancer,” Journal of Urology, vol. 175, no. 5, pp. 1915–1920, 2006.

A. Radonić, S. Thulke, H. G. Bae, M. A. Müller, W. Siegert, and A. Nitsche, “Reference gene selection for quantitative real-time PCR analysis in virus infected cells: SARS corona virus, Yellow fever virus, Human Herpesvirus-6, Camelpox virus and Cytomegalovirus infections,” Virology Journal, vol. 2, article 7, 2005.

N. Silver, S. Best, J. Jiang, and S. L. Thein, “Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR,” BMC Molecular Biology, vol. 7, article 33, 2006.

R. Pérez, I. Tupac-Yupanqui, and S. Dunner, “Evaluation of suitable reference genes for gene expression studies in bovine muscular tissue,” BMC Molecular Biology, vol. 9, article 79, 2008.

L. L. Hernandez, S. W. Limesand, J. L. Collier, N. D. Horseman, and R. J. Collier, “The bovine mammary gland expresses multiple functional isoforms of serotonin receptors,” Journal of Endocrinology, vol. 203, no. 1, pp. 123–131, 2009.

C. L. Andersen, J. L. Jensen, and T. F. Ørntoft, “Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets,” Cancer Research, vol. 64, no. 15, pp. 5245–5250, 2004.

M. W. Pfaffl, A. Tichopad, C. Prgomet, and T. P. Neuvians, “Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper—excel-based tool using pair-wise correlations,” Biotechnology Letters, vol. 26, no. 6, pp. 509–515, 2004.

A. M. Brunner, I. A. Yakovlev, and S. H. Strauss, “Validating internal controls for quantitative plant gene expression studies,” BMC Plant Biology, vol. 4, article 14, 2004.

V. Terzi, C. Morcia, M. Spini et al., “Identification and validation of reference genes for gene expression studies in water buffalo,” Animal, vol. 4, no. 6, pp. 853–860, 2010.

E. M. Glare, M. Divjak, M. J. Bailey, and E. H. Walters, “β-actin and GAPDH housekeeping gene expression in asthmatic airways is variable and not suitable for normalising mRNA levels,” Thorax, vol. 57, no. 9, pp. 765–770, 2002.

N. I. Bower and I. A. Johnston, “Selection of reference genes for expression studies with fish myogenic cell cultures,” BMC Molecular Biology, vol. 10, article 80, 2009.