Identification of Internal Reference Genes in Peripheral Blood Mononuclear Cells of Cattle Populations Adapted to Hot Arid Normoxia and Cold Arid Hypoxia Environments

To estimate gene expression in a reliable manner, quantitative real-time polymerase chain reaction data require normalisation using a panel of stably expressed reference genes (RGs). To date, information on an appropriate panel of RGs in cattle populations reared at cold arid high-altitude hypoxia and hot arid tropical normoxia environments is not available. Therefore, the present study was carried out to identify a panel of stably expressed RGs from 10 candidate genes (GAPDH, RPL4, EEF1A1, RPS9, HPRT1, UXT, HMBS, B2M, RPS15, and ACTB) in peripheral blood mononuclear cells (PBMCs) of cattle populations reared at cold arid high-altitude hypoxia and hot arid normoxia environments. Four different statistical algorithms: geNorm, NormFinder, BestKeeper, and RefFinder were used to assess the stability of these genes. A total of 30 blood samples were collected: six adult heifers each of Ladakhi (LAC) and Holstein Frisian crosses (HFX) and 4 Jersey (JYC) cows from cold arid high-altitude hypoxia environments (group I) and five adult heifers each of Sahiwal (SAC), Karan Fries (KFC), and Holstein Friesian (HFC) cows from hot arid normoxia environments (group II). Combined analysis of group I and group II resulted in identification of a panel of RGs like RPS9, RPS15, and GAPDH that could act as a useful resource to unravel the accurate transcriptional profile of PBMCs from diverse cattle populations adapted to distinct altitudes.


INTRODUCTION
India has been blessed with several cattle breeds of Bos indicus lineage adapted to various agroclimatic zones from high land to hot tropical regions. Leh-Ladakh, also known as a "COLD DESERT" is a part of the western Himalayan agro-climatic and high-altitude temperate sub-agro-climatic zone in India. Ladakh is situated at an altitude of 3,500-5,500 m with difficult terrain and harsh climate conditions such as extreme temperature (−40°C in winter and 35°C in summer), low humidity (25-40%), low precipitation (80-300 mm), and low oxygen level (nearly 60-70% of the oxygen concentration at sea level). In spite of harsh weather, this region is blessed with several unique livestock populations such as yak, cattle, horses, sheep, goat, donkey, and double-hump camel. Over thousands of years of the evolutionary process, these animals have developed the special ability to survive in cold and hypoxia environments prevalent in Ladakh. Amongst all the livestock species, the native cattle of Ladakh known as "Ladakhi cattle" are the major livestock species that play an important role in the agriculture economy of the region. The cattle from the Trans-Himalayan region of Ladakh are short in stature and well-adapted to the high-altitude environment. Based on morphometric data on 275 animals and genetic characterization using microsatellite markers (unpublished data), this cattle population was observed to be highly distinct from native cattle breeds adapted to other agroclimatic zones of India. Recently, our group was successful in delineating the distinct transcriptome signatures of peripheral blood mononuclear cells (PBMCs) of high-altitude-adapted Ladakhi cattle and tropically adapted Sahiwal cattle (Verma et al., 2018a;Verma et al., 2018b). In the last few years, the purity of Ladakhi cattle is believed to have declined due to widespread intermixing with Jersey cattle. However, considering the unique hypoxiatolerant characteristics of Ladakhi cattle, preserving its purity will be a key for long-term conservation and sustainable utilization.
On the other hand, India has also been blessed with a huge native cattle genetic resource base which has adapted to hot and tropical conditions. For example, native cattle breeds like Sahiwal, Tharparkar, Rathi, Gir, Ongole etc. are known for their superior thermotolerance as compared to their exotic counterparts of the Bos taurus lineage. The superior cellular tolerance ability of PBMCs from native cattle has been shown in a few studies published by our group Kishore et al., 2014;Sharma et al., 2019) using PBMCs as the cellular model. However, native cattle breeds are also facing genetic dilution due to crossbreeding with exotic germplasm in order to enhance milk production. These native cattle populations that are adapted to distinct altitudes might have acquired a distinct gene pool during the course of the evolutionary process. Such genetic resources with remarkable adaptive traits could be an interesting resource to mine gene expression and the mechanism underlying changes associated with adaptation to cold arid and hot arid environments. It would be interesting to define the importance of various genes in conferring adaptation to cattle populations adapted to diverse altitudes through targeted gene expression analysis. Quantitative real-time polymerase chain reaction (qRT-PCR) has been widely employed to quantify the expression of target genes of interest in different tissues/cells exposed to a variety of experimental conditions. In spite of many of its advantages, this tool is prone to analytical variations arising due to differing amounts of starting material, pipetting errors, and differing efficiencies of RNA extraction and reverse transcription (Vandesompele et al., 2002;Huggett et al., 2005;Bustin et al., 2009). To overcome the limitations of experimental variation, the use of appropriate internal control genes (ICGs) or reference genes (RGs) to successfully normalise the RT-qPCR data has been reported in several studies (Bustin 2010;Castigliego et al., 2010;Crookenden et al., 2017;Die et al., 2017;Sang et al., 2018). The approach to identify a suitable panel of RGs during various experimental/physiological conditions has also been reported in different livestock species (Aggarwal et al., 2013;Kapila et al., 2013;Zhu et al., 2015;Jatav et al., 2016;Li et al., 2016;Kaur et al., 2018). However, to the best of our knowledge, no comparative data on suitable RGs are available for cattle populations reared at cold arid high-altitude and hot arid tropical regions. The present study was planned to identify a panel of stably expressed RGs in PBMCs of six cattle populations from the high-altitude cold arid region of Leh-Ladakh and hot arid tropical climate of India.

Ethics Statement and Animal Selection
The blood sampling of animals was performed in accordance with the relevant guidelines and regulations as approved by the Institutional Animal Ethics Committee (IAEC) of ICAR-National Bureau of Animal Genetics Resources (NBAGR), Karnal. The study has included three cattle populations from the cold arid high-altitude region of Leh-Ladakh, viz., Ladakhi cattle (native), Jersey cattle (exotic), and Frieswal cattle {Sahiwal x Holstein Frisian cross}, and three cattle populations from the hot dry and semi-arid condition of Haryana state, viz., Sahiwal cattle (native), Holstein Frisians cattle (exotic), and Karan Fries cattle {Tharparkar x Holstein Frisian}. Ladakhi cattle (LAC) are the native breed of the Bos indicus lineage and are naturally adapted to the high-altitude region of Ladakh. Jersey cattle (JYC) are of the Bos taurus lineage and originated from temperate regions, while Frieswal cattle (HFX) are a cross-bred cattle population. Both Jersey and Frieswal cattle are non-native to Leh-Ladakh and have been reared in limited organized farms in the region since the last 2 decades. Amongst the cattle populations selected from the hot semi-arid region, Sahiwal cattle (SAC) are a very popular native cattle breed of Bos indicus lineage and are known for their adaptation potential to hot dry and tropical conditions. Holstein Frisian cattle (HFC) are of Bos taurus lineage and are non-native to India. However, Holstein Frisian cattle have been widely used in India in various cross-breeding programmes to enhance milk production of local breeds. Karan Fries cattle (KFC) are a popular cross-breed that were developed several decades ago by crossing Tharparkar cattle (native) with Holstein Frisian. Therefore the study has included two populations of native cattle: one from high-altitude (Ladakhi cattle; LAC) and the other from hot arid tropical regions (Sahiwal cattle; SAC); and two populations of cross-breeds and two populations of exotic cattle from two extreme altitudes. The geographical coordinates of the sampling site representing the hot arid climate were latitude-29°3′ 56.7828″ N, and longitude 76°2′ 25.7892″ E. The geographical coordinates of the sampling site from the high-altitude region of Ladakh were latitude-34°9′ 9.3168″ N, and longitude 77°34′ 37.3764″ E. About 7-8 ml of whole blood samples was collected aseptically from the external jugular vein of the animals in sterile EDTA-coated vacutainer tubes.
In total, 30 blood samples were collected from adult heifers; five each of Sahiwal (SAC), Karan Fries (KFC), and Holstein Friesian (HFC) cows from the hot arid normoxia environment and 5 heifers each of Ladakhi (LAC), Holstein Frisian crosses (HFX), and Jersey (JYC) cows from the cold arid high-altitude hypoxia environment. The entire workflow of the qPCR experiment is depicted in Figure 1.

Peripheral Blood Mononuclear Cells
Isolation, RNA Extraction, and cDNA Synthesis Immediately after collection, the blood samples were processed to isolate peripheral blood mononuclear cells (PBMCs) using the density gradient centrifugation method as described by Verma et al. (2018b). In brief, blood was diluted in a 1:1 ratio with 1 × PBS (Ca 2+ and Mg 2+ free; Hyclone, Utah) and was gently over laid on a Histopaque-1077 (Sigma-Aldrich Inc., United States) followed by centrifugation at 4000 RPM for 30 min at RT. After removing the buffy coat in a separate 15 ml tube, cells were treated with 2 ml of chilled RBC lysis buffer for 5 min at RT and washed twice with 1 × PBS (Ca 2+ and Mg 2+ free; Hyclone, Utah). Finally, the cells were suspended in 1.0 ml of ice cold Trizol reagent (Invitrogen, Carlsbad, California), homogenized, and stored at −80°C. Total RNA was extracted using Trizol reagent according to the manufacturer's instructions (Invitrogen, Corp., CA, United States). The extracted RNA was further purified using RNeasy mini kit columns (Qiagen, Germany). To remove the traces of genomic DNA, rnase free DNase enzyme was used according to the manufacturer's instructions (Qiagen, Germany). RNA concentration and purity were checked using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, United States) and Experion Bio-analyzer (Bio-Rad, United States). The OD 260 /OD 280 absorption ratio for different samples varied from 1.92 to 2.10. The RNA integrity number (RIN) values for all purified RNA samples were in the satisfactory range (6.5-8.5). The RNA integrity of each sample was also confirmed by visualizing 28S and 18S ribosomal bands on 1.5% agarose gel. The cDNA was synthesized using a Revertaid First strand cDNA synthesis kit (Fermentas, United States) as described in our previous studies Verma et al., 2018b). Briefly, first strand cDNA was synthesized with 200 ng of purified RNA, oligo-dT (18) primer, dNTP mix, random primers, RiboLockRNase inhibitor, and M-MuLV reverse transcriptase supplied with RevertAid First Strand cDNA Synthesis (Thermo Scientific, CA, United States), using the program: 25°C for 5 min, 50°C for 60 min, and 70°C for 15 min. Before using them as templates for qPCR, each of the cDNA samples was diluted 1:4 (v:v) with DNase/RNase free water. The quality of 30 cDNAs was confirmed by amplifying the GAPDH gene using a similar protocol to that described for qPCR except for the final dissociation protocol. A small aliquot of amplified products for all the samples was run on 2.5% agarose gel to check the primer specificity and amplification quality. FIGURE 1 | The entire workflow of the present study conducted using PBMCs of cattle breeds from cold arid high-altitude hypoxia and hot arid normoxia groups. The map of India in Figure 1 was created with the help of Smartdraw software-https://www.smartdraw.com).

Selection of Candidate Reference Genes and Real-Time Quantitative PCR Primers
In the present study, 10 candidate genes belonging to different functional groups were selected for evaluation as suitable RGs ( Table 1). The candidate genes included in the study were; glyceraldehyde 3-phosphate (GAPDH), ribosomal protein L-4 (RPL4), eukaryotic elongation factor 1 alpha (EEF1A1), ribosomal protein S9 (RPS9), hypoxanthine guanine phosphoribosyl transferase 1 (HPRT1), ubiquitin expressed transcript (UXT), hydroxyl methylbilane (HMBS), beta 2-microglobulin (B2M), ribosomal protein S15 (RPS15), and beta actin (ACTB). The primers specific for these 10 RGs were available in the laboratory and have been utilized successfully in several of our previous studies Jatav et al., 2016;Kaur et al., 2018). The information about sequences, amplicon length, and annealing temperature for each primer pair are summarized in Table 1.

Real-Time Quantitative PCR Reference Genes Transcripts
Quantitative PCR was performed in a 10 μL reaction volume containing 4 μl of diluted cDNA and 6 μl of master mix composed of 5 μL of 2X LightCycler 480 SYBR Green (Roche Life Science, Germany), 0.4 μL each of 10 μM forward and reverse primers, and 0.2 μL of DNase/RNase free water. Each qRT-PCR reaction was performed in duplicate to check the quality by assessing intra-assay variation. The amplification was carried out in a 96-well block using a LightCycler 480-II real-time PCR instrument (Roche Life Science, Germany) with the following conditions; 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). In order to evaluate the quality of qPCR reactions in terms of nonspecific amplification and primer-dimer formation, a dissociation curve for each gene was obtained by increasing the temperature from 60°C to 95°C. A six-point relative standard curve was prepared for each gene by using five-fold serial dilutions of pooled cDNA samples in duplicate. The amplification specificity for each primer was checked by the presence of a single band of expected size on 2.5% agarose gel (Supplementary Figure S1), and also by observing the single melt curve peak after completion of qPCR (Supplementary Figure S2). The qPCR data for each gene were extracted using the "second derivative maximum" method (Rasmussen 2001) as computed by Light Cycler software 3.5 for subsequent analysis.

Analysis of Expression Stability of Candidate Reference Genes
The qPCR data recorded for each gene were subsequently analysed to evaluate the expression stability. Four statistical approaches, viz., geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper (Xu et al.,. qPCR, efficiencies for each primer calculated pair-wise from a six-point standard curve using a five-fold dilution series of pooled DNA of Ladakhi and Sahiwal cow PBMCs R 2: correlation coefficient of the slope of the standard curve. Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 730599 2016), and the RefFinder web tool (Xie et al., 2012) were used to calculate the expression stability and ranking of RGs in the highaltitude hypoxia cold arid group (LAC, JYC, and HFX) and normoxia hot arid group (SAC, KFC, and HFC). The C q values of each RG were exported to an Excel work sheet and modified as per the requirement of the software. Like for creating an input file for geNorm and NormFinder analysis, the C q values were first transformed into relative quantities by using formula 2− ΔCT , in which ΔC q corresponding C q value -minimum C q value. geNorm calculated the expression stability of individual genes on the basis of the M value which indicates the stability in expression of a gene. The genes with smaller M values (<1.5) are considered to have higher expression stability . In addition, geNorm was also used to conduct pair-wise variation analysis (Vn/Vn+1) in order to select the optimal number of RGs to normalise the target gene expression data. The cut off value of Vn/Vn+1 < 0.15 was used to decide the optimal number of RGs to be employed for calculating the normalisation factor (Vandesompele et al., 2002). This analysis is based on the principal that the expression ratio of the two best RGs will always remain similar across samples. The NormFinder software calculated the stability values of each RG based on inter-and intra-group variations. Similar to geNorm, for Normfinder analysis, the C T values were first transformed into relative quantities. However, for Bestkeeper analysis, C T values were not transformed into relative quantities. The BestKeeper algorithm was used to calculate gene expression variation based on cycle threshold values (C q ), crossing point standard deviation [{SD (±CP)} <1], and coefficient of variance (CV [%CP]). In BestKeeper analysis, genes with low SD (<1), low CV, and high coefficient of correlation (r) are generally considered stably expressed and vice versa. Finally, the RefFinder tool (https://www.heartcure.com.au/ reffinder/) was also employed to estimate the overall ranking of the 10 RGs by assigning an appropriate weight to each gene (Xie et al., 2012). The RefFinder analysis integrates the outcome of geNorm, NormFinder, and BestKeeper tools to provide an overall final ranking of RGs.

Primer Specificity, Amplification Efficiency, and Descriptive Statistics
This study evaluated the expression stability of 10 RGs (GAPDH, RPL4, EEF1A1, RPS9, HPRT1, UXT, HMBS, B2M, RPS15, and ACTB) in PBMCs of cattle types from cold arid high-altitude hypoxia and hot arid normoxia environments. The specificity of each primer pair was ascertained by the presence of the single specific amplicon of expected size on 2.5% agarose gel (Supplementary Figure S1). The melt curve plot of each RG showed a single peak suggesting the highly specific nature of primer pairs used (Supplementary Figure S2). The amplification efficiencies estimated from the six-point standard curve (generated from five-fold serial dilution of pooled cDNA) ranged from 88 to 110%. The slope values of the standard curve for different RGs ranged between −3.05 and −3.5 which was within the acceptable limit (−2.96 to −3.6). Based on the overall evaluation of the melt curve, amplification efficiencies, and slope values, it could be safely assumed that RT-qPCR data for each primer pair were of high quality ( Table 1).
The average raw Cq values of individual RGs across all PBMCs in cold arid and hot arid groups (combined analysis) are summarized in Table 2. The average Cq values of individual RGs were quite variable and ranged from 17.66 (EEF1A1) to 27.86 (HMBS). On the basis of their average C t scores, the 10 RGs were classified into group-I (abundantly expressed), group-II (moderately expressed), and group-III (least expressed). Group-1 included EEF1A1, B2M, RPS15, and RPL4 genes with a high expression level and average Cq scores of 17. 66, 18.18, 19.35, and 19.68 Table 2). Considering the distribution of average raw Cq scores and interquartile range, EEF1A1 exhibited the lowest coefficient of variations (lowest variability across samples). On the other hand, the HMBS gene showed the highest coefficient of variation. Based on this parameter, EEF1A1, B2M, RPS15, and RPL4 RGs were most stable while UXT, HPRT1, and HMBS genes were the least stable across the combined dataset ( Figure 2).
Similar to the combined dataset, in cold arid as well as hot arid normoxia groups, EEF1A1 showed the highest maximum and HMBS showed the least expression level. The average Cq values of individual RGs in cold arid hypoxia and hot arid normoxia groups are summarized in Supplementary Tables S1,S2 and Supplementary Figures S3,S4 respectively.

Selection of Reference Genes by geNorm Analysis
The first analysis to determine the expression stability of 10 candidate RGs in PBMCs of heifer cows from the hypoxia cold arid group (LAC, JYC, and HFX) and normoxia hot arid group (SAC, KFC, and HFC) was based on the geNorm algorithm. The expression stability of individual RGs was evaluated first by combining RT-qPCR data for all 30 PBMCs samples (cold arid hypoxia and hot arid normoxia groups). All the RGs in combined group analysis showed stable expression with a stability index within the acceptable range (<M value < 1.5). In the combined analysis, RPS9 and RPS15 showed the highest expression stability (M 0.464), followed by RPL4 (M 0.527) and GAPDH (M 0.539), whereas HPRT1 was the least stable gene (M 1.228). On the basis of average expression stability measure, RGs were arranged from most stable (lowest M value) to the least stable (highest M value): RPS9 RPS15 > RPL4> GAPDH > HMBS > EEF1A1> B2M > UXT > ACTB > HPRT1 ( Figure 3A). Since the samples were from two distinct altitudes, we tried to determine the optimal number of RGs on the basis of pair-wise variation (Vn/n+1 were either > or < the threshold of 0.15 ( Figure 3B). It is generally considered that when Vn/n+1 is < 0.15 (threshold value), inclusion of an additional gene is not required for  Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 730599 6 calculating the normalisation factor. In the combined dataset, combination V3/V4 resulted in a Vn/n+1 value < 0.116 which was less than the cut off value of 0.15 suggesting that the first three most stable RGs will be sufficient for accurate normalisation of RT-qPCR data. Therefore, based on M (stability) and V values (pair-wise variation) derived from geNorm analysis, the most stable RGs for the combined dataset were RPS9, RPS15, and RPL4 ( Figure 4 and Table 3).
In addition to combined analysis, the geNorm tool was used to evaluate the expression stability of individual RGs separately in the cold arid hypoxia group (LAC, JYC, and HFX) and hot arid normoxia group (SAC, KFC, and HFC). The M values for all the 10 FIGURE 3 | GeNorm analysis for ranking of genes based on average expression stability measure (M value) and pair-wise variation (Vn/Vn+ 1) between the normalisation factors NFn and NFn + 1 to determine the optimal number of reference genes. Analysis in the combined group (cold arid and hot arid) (A and B, respectively), cold arid hypoxia group (C and D, respectively), and hot arid normoxia group (E and F, respectively). RGs in the hypoxia cold arid group were <1.5. The ranking of RGs in order of expression stability within this group were: RPL4 EEF1A1 >RPS9> HMBS > RPS15 > B2M > GAPDH > HPRT11 > UXT > ACTB ( Figure 3C; Table 3). RPL4 and EEF1A1 showed the highest expression stability (M 0.255) followed by RPS9 (M 0.285) and HMBS (M 0.312), whereas ACTB was least stable (M 0.940). The pair-wise variation analysis for different combinations of RGs was well below the threshold value of 0.15. The results obtained for V2/3 (0.092), V3/4 (0.075), V4/5 (0.097), and V5/6 (0.077) combinations indicated that RPL4 and EEF1A1 would provide the most accurate normalisation ( Figure 3D). The expression stability of the three best RGs in combined analysis is shown in Figure 4B.

Selection of Reference Genes by NormFinder Analysis
In the NormFinder-based intergroup (combined) analysis covering both the conditions and all samples (cold arid hypoxia and hot arid normoxia), RPL4, RPS9, and RPS15 were found to be most stable with stability values of 0.282, 0.289, and 0.292 respectively. On the other hand, UXT, ACTB, and HPRT1 were least stable with stability values of 0.605, 0.737, and 1.288, respectively ( Table 3). The graph showing intragroup variation of RGs in the combined dataset is shown in Figure 5A. Based on stability values, the RGs were ranked as RPL4> RPS9> RPS15 > GAPDH > HMBS > B2M > EEF1A1> UXT > ACTB > HPRT1. Within the cold arid hypoxia group, RPS15, EEF1A1, and RPS9 were the most stable RGs with stability values of 0.140, 0.140, and 0.155, respectively, whereas, ACTB was the least stable gene with the highest variability value of 0.808 ( Table 3). Figure 5B. The RGs were ranked as RPS15> EEF1A1> RPS9> RPL4> B2M > GAPDH > HMBS > HPRT1>UXT > ACTB. Within the hot arid normoxia group, NormFinder analysis identified RPS15, RPS9, and HPRT1 as the three most stable RGs with stability values of 0.107, 0.134, and 0.136, respectively ( Table 3). The graph showing intragroup variation analysis of RGs in the combined dataset is shown in Figure 5C. Based on stability values, the RGs were ranked as RPS15 > RPS9> HPRT1> HMBS > GAPDH > RPL4> EEF1A1> UXT > ACTB > B2M. Overall, there was a good agreement in Frontiers in Genetics | www.frontiersin.org February 2022 | Volume 12 | Article 730599 9 geNorm and NormFinder outcomes for all the datasets ( Table 3).

Selection of Reference Genes by RefFinder Analysis
Additionally, the RefFinder algorithm was used to evaluate the comprehensive ranking of individual RGs in combined, cold arid

Validation of Selected Reference Genes
To evaluate the reliability of the best suitable and worst panel of RGs, a validation qPCR experiment was performed using some of the known candidate target genes associated with high-altitude hypoxia and heat stress response. The qPCR data for two target genes associated with hypoxia and high altitude such as HIF-alpha and EPAS1 and two target genes associated with heat stress response such as HSP70 and HSP27 were generated in PBMC samples of highaltitude-adapted and tropically adapted cattle populations. As shown in Figure 6A, it is quite evident that the panel of best reference genes (RPS19, RPS15, and GAPDH) normalised the target gene data more accurately. We expected higher expression of high-altitudeassociated genes in PBMCs of high-altitude cattle populations. On the other hand, genes related to heat stress response such as HSPs should have higher expression in PBMCs of cattle populations from the tropical region. As shown in Figure 6A, the expression of target genes such as HIF-alpha and EPAS1 was higher in highaltitude (HA)-adapted cattle while genes like HSP70 and HSP27 were more expressed in low-altitude (LA) cattle breeds. Further, the relative expression, standard deviation (SD), and standard error (SE) of the four target genes also supported the high quality of the panel of reference genes used (Supplementary Table S3).
On the other hand, normalisation of same target genes with least stable genes (HPRT and ACTB) resulted in an unexpected   pattern of expression. As shown in Figure 6B, HSPs showed higher expression in PBMCs of high-altitude cattle populations. Further, the use of least stable reference genes resulted in higher SD and SE values (Supplementary Table S4). Considering the above facts, it could be stated that the panel of RGs identified in the present study will be applicable for accurate normalisation of target genes in studies involving high and low-altitude cattle populations.

DISCUSSION
In an era of high-throughput platforms, quantitative PCR (qRT-PCR) is being employed as a most preferred tool to validate gene expression data. Even though qRT-PCR is the most sensitive technique, it suffers from several analytical variations like differences in the amount of starting material, RNA extraction, and efficiency of the reverse transcription process. To a great extent the effect of these nonbiological variables can be nullified by normalising gene expression data by a panel of stable reference genes (RGs). As suggested in the "Minimum information for publication of Quantitative Real-time PCR Experiments" (MIQE) guideline (Bustin et al., 2009), accuracy in gene expression is largely governed by the availability of reliable RGs. At present there is no consensus for the set of reference gene(s) that can be used universally for the normalisation purposes. In the past, several studies have utilized RGs without proper validation or arbitrarily selected commonly used reference genes such as ACTB and GAPDH. Unfortunately, the use of such empirical RGs might not provide accurate normalisation and create doubt on the reliable estimation of gene expression. Moreover, many reports showed the variable expression of commonly used reference genes in different cells, tissue, and conditions (Kim and Yun 2011;Zhao et al., 2012;Thomas et al., 2014). The major challenge in any biological experiment involving different tissues, cell types, disease state, and physiological and or developmental stages is the knowledge about appropriate RGs whose expression remains constant without any observable variations across samples (Bustin 2002;Radonić et al., 2004). It has been seen that a particular RG appropriate in one condition might have variable expression in another set of conditions. Thus, identification and validation of proper RGs is the prerequisite for any specific experimental condition. Our group has successfully identified a panel of appropriate RGs for various cellular types and experimental conditions involving zebu cattle and riverine buffaloes Sood et al., 2017;Lagah et al., 2019). The present study was also conducted on similar lines to identify suitable reference genes for cattle populations from diverse environmental conditions. The 10 RGs (RGs selected for the present work were part of our earlier studies, wherein these genes were evaluated for their suitability as normalisers in different cell types, cattle breeds, and experimental conditions Jatav et al., 2016;Kaur et al., 2018)).