Transcriptome Analysis of Circulating PBMCs to Understand Mechanism of High Altitude Adaptation in Native Cattle of Ladakh Region

Ladakhi cattle is native population of Leh and Ladakh region and constantly exposed to hypobaric hypoxia over many generations. In present study, transcriptome signatures of cattle from Ladakh region (~5500 m) and Sahiwal cattle from tropical regions were evaluated using Agilent 44 K microarray chip. The top up-regulated genes in Ladakhi cows were INHBC, ITPRI, HECA, ABI3, GPR171, and HIF-1α involved in hypoxia and stress response. In Sahiwal cows, the top up-regulated genes eEF1A1, GRO1, CXCL2, DEFB3 and BOLA-DQA3 were associated with immune function and inflammatory response indicating their strong immune potential to combat the pathogens prevalent in the tropical conditions. The molecular pathways highly impacted were MAPK signaling, ETC, apoptosis, TLR signaling and NF- kB signaling pathway indicating signatures of adaptive evolution of these two cattle types in response to diverse environments. Further, qPCR analysis revealed increased expression of DEGs such as HIF-1, EPAS-1, VEGFA, NOS2, and GLUT-1/SLC2A1 in cattle types from high altitude suggesting their pivotal role in association with high altitude adaptation. Based on data generated, native cattle of Ladakh region was found to be genetically distinct from native cattle adapted to the tropical region of India.

The adaptation of livestock to high altitude hypoxic environment is important as they play an important role in sustaining the socio-economic condition of the local residents. It is a well-known fact that high-altitude adaptation is an evolutionary process in mammals with considerable physiological changes so as to make them survive and perform optimally at the extreme environmental condition. Conversely, 'acclimatization' is an immediate physiological response to changing environments. The term "high-altitude adaptation" can thus be described as irreversible, long-term physiological response associated with heritable, behavioural and genetic changes. One of the characteristic environmental feature of high altitude region is sustained hypobaric hypoxia as reflected by lower oxygen (O 2 ) pressure that results in the insufficient supply of O 2 to the tissues 1 . On the basis of O 2 concentration level, hypoxia can be defined as moderate (5-8% O 2 ) or anoxic (less than 1% O 2 ) 2 , whereas; on the basis of time scale, hypoxia can be defined as acute (lasting for seconds to minutes) or chronic (lasting for hours to days) 3 . Low O 2 pressure along with low temperature presents numerous physiological challenges to the animals living at high altitude.
Animals living under harsh conditions of high altitude have the unique ability to adapt themselves to chronic hypoxia and low ambient temperature in comparison to the animals living at sea level 4 , by undergoing modifications at physiological, molecular and cellular levels. They respond to hypoxia by maintaining O 2 delivery through increased respiration rate, red blood cell mass and blood volume exhibiting increased blood O 2 carrying capacity. They also possess greater ability to enhance O 2 uptake and its delivery to tissues by suppressing their metabolic activity 3 , and/or enhancing flux capacity of the O 2 transport 5 . Further, animals offset their O 2 transport system in order to maintain the tissue O 2 level for their growth, development and reproduction 6 . The high altitude animals also initiate several oxidative mitochondrial metabolism and signal transduction pathways which cascade into a series of biochemical and physiological alterations that enable the animal to survive in hypoxic condition. Amongst all the mammals, yak represents the best example of high altitude adaptation. They are genetically adapted to hypoxia, low ambient temperature, strong solar radiation and alpine environment 7 through natural selection over millennia. Hypoxia in animals stimulates the expression of hypoxia inducible factor-1α (HIF-1α), a heterodimeric transcription factor which further regulates the transcription of several hypoxia related genes associated with several processes like erythropoiesis, angiogenesis, and glycolysis 8 . Prolonged hypoxia also leads to enhanced generation of reactive oxygen species (ROS), and lipid peroxidation. In such a situation, antioxidants such as superoxide dismutase (SOD), glutathione-peroxidase (GPX) and catalase (CAT) are activated. These antioxidants reduce the oxidative stress by preventing damage to important cellular components caused by ROS such as free radicals, peroxides, lipid peroxides and heavy metals 9 . Thus, high altitude environment presents numerous physiological challenges to the animal.
In India, one such region Leh and Ladakh also known as "land of high passes" is situated at an altitude of 3,500-5,500 m. It harbors difficult terrain, barren lands and little vegetation with cold arid and hypoxic conditions. This region in the state of Jammu and Kashmir lies between the Kunlun mountain range in the north and the main Great Himalayas to the south. The temperature here drops below −20 °C in winter season which lasts for about six months. Leh and Ladakh region is blessed with several native domesticated livestock species such as yak, yak and cattle cross (dzo, dzomo), mountain goats, sheep, cattle, horse, donkeys and double hump camel. Each of these species has developed an effective mechanism to survive at high altitude, particularly to combat the low temperature and low O 2 condition. These animal genetic resources are very important in the life of the local people of Ladakh region as land resources are meagre and their livelihood is mainly dependent on these animals. Under such harsh climatic conditions, the survival and performance of exotic breeds is not a viable option as the lack of O 2 and fodder/feeds at higher altitudes, only allow the well adapted animal genetic resources to thrive and perform.
The local cattle from Leh and Ladakh region is one such native population that has evolved naturally over the years and got adapted to the harsh high altitude hypoxia environment. This cattle is the lifeline of the local people and plays a pivotal role in fulfilling their nutritional needs. The subsistence on poor quality feed/fodder and better resistance capabilities to withstand environmental stress make this cattle distinct and preferred population by local people. Further, naturally evolved Ladakhi cattle could be a rich source of gene pool, as it might have inherited beneficial traits which allow them to combat the extreme environmental conditions.
Keeping in view the remarkable adaptive traits and importance, the Ladakhi population is an interesting resource to mine the genes and pathways associated with high altitude adaptation. In the present study, the transcriptome signature of peripheral blood mononuclear cells (PBMCs) of Ladakhi cows adapted to high altitude vis a vis Sahiwal cows adapted to the arid/semi-arid region at mean sea level was established using bovine expression microarray chips.

Identification of DEGs in PBMCs of Ladakhi and Sahiwal cows.
In the present study, the transcriptome profile of PBMCs of Ladakhi cows native to Leh and Ladakh region (high altitude) was generated. To understand the transcriptomic variation caused by difference in altitude, transcriptome profile of PBMCs of Sahiwal cows adapted to tropical condition (mean sea level) was also generated. The transcriptome data was generated using bovine gene expression microarray from Agilent Technologies (Product number: G2519F). The Agilent bovine oligonucleotide array used in this study was in 4 × 44 K format with design ID V2:023647 and 43,803 probes. The extracted RNA was intact and sufficiently good in quality as indicated by electropherogram with two characteristics peaks of 18S rRNA & 28S rRNA. (Supplementary Fig. S1). To ensure successful microarray hybridization, the yield and specific activity of Cyanine3 (Cy3)-labeled complementary RNA (cRNA) was measured in terms of Cy3 dye concentration (pmol/μl), RNA absorbance ratio (260 nm/280 nm) and cRNA concentration (ng/μl). The concentration of cRNA (ng/μl) was used to determine the μg cRNA yield and concentrations of cRNA (ng/ μl) and Cy3 (pmol/μl) was used to determine the specific activity. In all the samples of Ladakhi cows and Sahiwal cows, the yield of cRNA was >1.65 μg and specific activity was >9.0 pmol Cy3 per μg cRNA, indicating the good quality of cRNA obtained. The scanned images analysed using Agilent feature extraction software showed successful hybridization of the cRNA samples. The microarray raw data files were subjected for quantile normalization to remove unwanted technical variation. Based on intensity distribution, the box whisker plot showed a consistent distribution of normalized intensity values in samples of Ladakhi and Sahiwal cows (Fig. 1A). In addition principle component analysis plot (Fig. 1B) revealed distinct grouping of Ladakhi and Sahiwal cows based on transcriptome data. The normalized transcriptome data was filtered based on multiple testing corrected P values, P < 0.05 (FDR < 0.05) and different fold change criteria. Line plot of transcripts showing differential genes expression of between Ladakhi and Sahiwal cows PBMCs is shown in Fig. 1C.
Overall microarray data analysis revealed a total of 8417 differentially expressed genes (DEGs) between Ladakhi and Sahiwal PBMCs with multiple testing (FDR < 0.05; Fig. 1D). With additional cut off criteria i.e., fold change of 2 or more, a total of 3910 genes were found to be differentially expressed between the two cattle types. Venn analysis showed that out of these 3910 genes, 1207 were up-regulated while 2703 genes were down-regulated in Ladakhi cows in comparison to Sahiwal cows (Fig. 1E). Overall analysis revealed substantial differences in transcriptome signature between the two cattle types.
Clustering of transcriptome data. The normalized transcriptome data was partitioned using hierarchical and k-means clustering tools. For creating a hierarchical clustering, a total of 3910 DEGs were used (Fig. 1F). The k-mean clustering yielded a total of 6 major clusters (Fig. 1G). These clusters helped to differentiate the genes on the basis of their expression pattern between Ladakhi and Sahiwal cows. The heat maps thus generated helped to judge the similarities/patterns between genes and between samples. The clustering tools and heat map collectively revealed distinct transcriptome signature for high altitude adapted Ladakhi cows and tropically adapted Sahiwal cows (Fig. 1H).
Top up-regulated DEGs in PBMCs of Ladakhi cows. Out of 3910 DEGs, a total of 1207 genes were found to be up-regulated in Ladakhi cows. List of top 50 DEGs with their absolute fold change values (Fc), log fold change values (Log Fc), gene symbol, gene description are presented in Table 1.
The transforming growth factor TGF-beta (TGFBI) gene that is a part of transforming growth factor signaling pathway and is involved in many cellular processes including cell growth, differentiation and apoptosis, cellular homeostasis was also up-regulated (26. Another category of genes that showed upregulation in Sahiwal was PHD finger protein genes including, PHF2 (2.83 fold), PHF8 (2.03 fold), PHF17 (4.45 fold) and PHF20 (6.75 fold). These genes are related to methyl lysine-binding protein, a component of the MOF histone acetyltransferase protein complex. The higher expression of these genes might be attributed to the presence of higher atmospheric O 2 that is must for their activity.  Supplementary Table S1. Under GO term 'response to stimulus' the major subcategories included were response to stress (GO: 00006950; p value: 0.5850) and immune response (GO: 0006955; p value: 0.8423). Under term 'response to stress' (Fig. 2C) several other terms like cellular response to stress; defense response; response to wounding, oxidative stress, starvation and hypoxia were found to be enriched. The genes represented under each subcategory of response to stress in Ladakhi cows are listed in Supplementary Table S2. The top molecular functions identified in Ladakhi cows are mentioned in Supplementary Table S3 included Supplementary Table S4.
In tropically adapted Sahiwal cows, the major GO terms identified under biological processes are shown in Supplementary Fig. S2A  stimulus, (GO: 0050896; p value: 0.0788), cellular process (GO: 0009987; p value: 0.0634), biological regulation (GO: 0065007; p value: 0.0525) were the most common GO terms observed in the data set. Under response to stimulus, the important sub-ontologies identified were: response to stress, biotic stimulus, abiotic stimulus and immune response. Under response to stress (Fig. S2B) two sub-ontologies identified were cellular response to stress and oxidative stress. Under immune response (Fig. S2C), the major sub-ontologies enriched were innate immune response, adaptive immune response and cell activation involved in immune response. List of genes identified under sub-categories of response to stress and immune response in Sahiwal cows is mentioned in Supplementary Tables S5 and S6. Additionally, the major GO terms identified are shown in Supplementary  Fig. S2D

Validation of gene expression data by qPCR.
To validate the microarray results quantitative-real time PCR (qPCR) was performed for some hypoxia related target genes such as HIF-1α, EPAS-1, VEGF-A, ECE-1, NOS2, GLUT-1, HK2, TNFα, GRα and heat shock proteins; HSP27, HSP70 and HSP90. For data normalization GAPDH, RPS9, HMBS and RPS15 were identified as the most stably expressed internal control genes (ICGs) in all the 30 PBMC samples. The identified ICGs were subsequently utilized to normalize the qPCR data generated for hypoxia and heat stress related genes in PBMCs of different cattle types. The qPCR performance for each ICG and target genes in terms of slope of six-point standard curve, coefficient of determination of standard curve (R 2 ) and efficiency of amplification (E = 10 (−1/slope) ) is mentioned in Supplementary Tables S9 and S10.
The qPCR results showed an increased expression of hypoxia related target genes viz; HIF-1α, EPAS-1, VEGF-A, ECE-1, NOS2, GLUT-1, HK2, TNFα and GRα in cattle types from high altitude region (LAC, HFX, JYC) than cattle types from low altitude (SAC, HFC, KFC). Among high altitude cattle types, Ladakhi cattle showed least expression of all the target genes in comparison to Holstein Friesian crosses and Jersey cattle residing at high altitude (Fig. 3). In comparison to high altitude cattle, the qPCR data showed a significant (p < 0.05) increase in mRNA expression of HSP27, HSP70 and HSP90 in cattle from low altitude than cattle from high altitude region.

Discussion
Ladakhi cattle population is hypoxia tolerant, living in an extremely inhospitable high-altitude environment with low temperature, compared to Sahiwal cattle living in extreme hot tropical environment. Out of 3910 DEGs, the number of down-regulated DEGs was quite high in PBMCs of Ladakhi cows (2703) in comparison to Sahiwal cows (1207). This could probably be explained by the fact that high altitude hypoxia environment has a down-regulatory effect on overall transcriptional machinery of Ladakhi cows living at high altitude in comparison to Sahiwal cows living at tropical/normoxic conditions. Variation in the number of DEGs in Ladakhi and Sahiwal cows clearly show the difference in transcriptome profile of both the cattle from two different altitudes. Further, hierarchical clustering, k-means clustering and distinct heat maps also revealed that the transcriptome signature of Ladakhi cows is highly distinct from that of the Sahiwal cows.
The genes that were up-regulated in Ladakhi cattle might be responsible for high altitude adaptation. Amongst these, Inhibin beta c (INHBC) that was the top most up-regulated gene (197.63 Fold) in Ladakhi cows is involved in regulating cell growth, proliferation, differentiation 11 and a number of diverse functions such as hypothalamic and pituitary hormone secretion, gonadal hormone secretion, germ cell development and maturation, erythroid differentiation, insulin secretion, nerve cell survival and bone growth. Inositol 1, 4, 5-triphosphate receptor type 1 (ITPR1), the second most up-regulated (75.37 Fold) gene in Ladakhi cows is a target of HIF-2α 12 and plays an important role in ER stress-induced apoptosis. Some of the important hypoxia related genes such as HIF-1 HIF-3α, VEGFA, VEGFB, VEGFC, NOS2, GLUT-1, GLUT-3 and HK2 identified in present data set have been marked as important candidate genes in previous studies. Hypoxia inducing factor-alpha (HIF-1α), up-regulated  15.66 fold in Ladakhi PBMCs act as a key regulator of O 2 homeostasis at the cellular and systemic level 13 . The activation of HIF-1α is mainly dependent on the hypoxia induced stabilization of the α-subunit, whereas its β-subunit can act independent of the concentration of O 2 . Under hypoxia, activation of the anaerobic glycolysis mainly occurs to compensate for the energy deficit in the cell. HIF-1α is known to up-regulate the expression of enzymes involved in glucose uptake (glucose transporters) and glycolysis 14 . In parallel, HIF-1-dependent factors  modify the vascular tonus resulting in an improved blood circulation. Till now, about 100 target genes have been shown to be regulated by this gene 15 . Our data showing significant induction of HIF-1α gene strongly support the role of HIF gene in maintaining the O 2 homeostasis and other processes necessary for survival of high altitude adapted Ladakhi cows. Simultaneously, hypoxia inducing factor 3A (HIF3A) that acts as a transcriptional regulator in adaptive response to low O 2 tension was also up-regulated 18.53 fold in PBMCs of Ladakhi cow. HIF is also known to activate vascular endothelial growth factor (VEGF), which encodes for protein that induces angiogenesis 16 . VEGFs are essential for growth of new blood vessels promoting angiogenesis and play an important role in hypoxia response and high altitude adaptation 17 . All the three vascular endothelial growth factors, VEGFA, VEFGB and VEFGC were up-regulated in Ladakhi PBMCs by 6.35, 4.75 and 6.93 fold, respectively. Up-regulated expression of VEGFs in Ladakhi cows suggest them to be candidate genes in high altitude adaptation. Similarly, Endothelial PAS domain-1 (EPAS-1), another major transcription factor was up-regulated (19.72 fold) in Ladakhi cows. The association of EPAS1 gene in high altitude adaptation is well documented in humans and dogs [18][19][20] and has been related to the haemoglobin concentration in Andeans and Tibetans [21][22][23] . The nitric oxide synthase (NOS2) was also up-regulated (17.46 fold) in Ladakhi PBMCs, further substantiating the previous studies that high altitude hypoxia increases the expression of inducible nitric oxide synthase 24 . The nitric oxide (NO) generated from NOS is associated with several physiological processes at higher altitude. In the presence of NO, blood vessels inner lining (endothelium) signal the surrounding smooth muscle to relax which results in vasodilation and increase in blood flow to enhance the O 2 supply under hypoxic condition. The enhanced expression of NOS2 gene may be responsible for higher production of NO in Ladakhi cattle so as to help in its pulmonary vasculature vasodilation under hypobaric hypoxia condition aiding in better survivability. Calcium channel, voltage-dependent, T type, alpha 1G subunit (CACNA1G) that encodes one of the voltage-sensitive calcium channels mediating the entry of calcium into cell was upregulated by 26.51 fold in Ladakhi cattle. This gene is an integral part of calcium signaling pathway and reported to play an important role in high altitude adaptations to hypoxia 25 . Forkhead box M1 (FOXM1) and regulator of telomere elongation helicase (RTEL1) genes upregulated in Ladakhi cows have been reported to play an essential role in repairing oxidative DNA damage 26,27 . This gene was also marked as an important candidate gene in Snub nosed monkeys 28 and Tibetan pigs 29 .
Further, up-regulation of Myosin heavy chain 2 (MYH2) gene in Ladakhi cows may be related to skeletal muscle contraction and regulation of processes like glycolysis and glucose metabolism. This gene has also been reported to have higher expression in yak in comparison to cattle of low altitude 30 . Similarly, higher expression of hexokinase 2 (HK2) in Ladakhi cows suggested increased glucose metabolism in hypoxia environment as HK2 is known to be an important glycolytic enzyme involved in the metabolism of glucose to glucose-6-phosphate 31 . Several studies have reported a close association of HK2 expression with glucose transport 32,33 . Other genes involved in anaerobic metabolism like glucose transporters (GLUT1; GLUT3), lactate dehydrogenase A (LDHA), pyruvate dehydrogenase kinase 1 (PDK) were also present in higher abundance in Ladakhi cows. Several hypoxia-responsive genes namely, MMP1, MUC1, APLN, FZD3 were also upregulated in Ladakhi cows. Amongst these, APLN and FZD3 encode for factors with proangiogenic functions as reported under hypoxic conditions 34 . PRKAA1, another important upregulated gene in Ladakhi cows has been reported to be important in achieving genetic adaptation and maintaining metabolic homeostasis at high altitude 35 . Furthermore, the transcript for SENP2 an enzyme implicated in erythropoiesis 36 ; MICU1, VAV3, ARNT2 linked to hemoglobin concentration; CXCL17 involved in vascular physiology and PAFAH1B3 linked to hypoxia 37,38 were also observed to be upregulated.
Interestingly, most of the genes from immune functional classes were down-regulated in Ladakhi cows PBMCs. The lower expression of these genes in Ladakhi cows could be explained by the fact that at high altitude environment, there is lesser load of pathogens/microbes. This scenario may reduce the risk of opportunistic infection in Ladakhi animals, hence these immune genes were relatively lowly expressed in comparison to Sahiwal cows that belongs to tropical conditions. In comparison to high altitude adapted Ladakhi cows, genes like eEF1A1 responsible in stabilizing the formation of functional ribosome and in translation initiation; chemokine genes such as CXCL3, CCL4, CCL16 related to several immune functions and host defence were present several fold higher in Sahiwal PBMCs. The chemokines are known to induce strong immune response in cattle breeds and play fundamental role in the homeostasis, and functioning of immune system 39 . Another major immune related gene CSF3R, a member of cytokine receptors family that helps in cell surface adhesion or in migration of leukocytes was observed to be up regulated in Sahiwal PBMCs. The protein encoded by this gene is the receptor for colony stimulating factor 3, a cytokine that controls the production, differentiation, and function of granulocytes. The series of beta defensin genes viz., DEFB3, DEFB1, DEFB4A and DEFB10 involved in host defence system was observed to be upregulated in Sahiwal PBMCs. These antimicrobial peptides are thought to be evolved under the pressure of natural selection to maintain a host-pathogen balance in cattle breeds adapted to tropical region 40 . The higher abundance of these antimicrobial peptides can possibly be explained with the fact that cattle types from tropical region are exposed to diversity of pathogens. Additionally, the enrichment of genes like BOLA-DQA3, BoLA, BOLA DQA1, BOLA DQB and BOLA DRB3 in PBMCs of Sahiwal cows in comparison of Ladakhi cows signifies the immune potential and active host defence mechanism in Sahiwal cattle to combat the pathogens normally present in the tropical conditions. Overall, these genes would enable the immunological defences of Sahiwal cows under harsh tropical environmental condition. The differences in expression of BOLA genes in Sahiwal and Ladakhi cows might be due to the differences in immune system activation at two distinct altitudes and types of pathogens prevailing in high and low altitude environment. In addition, increased expression of molecular chaperons such as DNAJC2, DNAJC27, HSPA4, HSPB11 in Sahiwal than Ladakhi, could be due to environmental heat load present in the tropical conditions in comparison to environmental condition (cold arid) prevalent in Leh and Ladakh. The overexpression of cluster of genes related to DEAD Box polypeptides (DDX), DDX47, DDX55, DDX6, DDX18, DDX21, DDX23, DDX54, and DDX59 was also significant in Sahiwal than Ladakhi cows. These genes are known to be highly conserved in nature and are generally involved in RNA metabolism. The results clearly suggest that there could be substantial differences in the processes or expression of proteins involved in RNA metabolism between tropically adapted Sahiwal and high altitude adapted Ladakhi cows. The lower abundance of DDX genes in Ladakhi cows from high altitude is similar to few other studies that have shown that under hypoxic condition, there is down-regulation of DDX expression 41 . Also significant increase in expression of several ribosomal proteins viz; RPL34, RPL36, RPS12, RPS24, and mitochondrial ribosomal proteins MRPS15 MRPL44, MRPL46, MRPS17 could be due to enhanced need of ribosomal biogenesis to maintain normal protein translation and essential cellular functions in Sahiwal cows. Overall, the transcriptome data generated for Sahiwal cows PBMCs indicated the comparative enrichment of genes associated with immune function, inflammatory response, heat stress, biological process regulation and ribosomal biogenesis related genes.
Gene ontology analysis revealed that genes involved in biological process like cellular response to stimulus  expected. This could be attributed to the typical arid or semiarid environmental condition exist in the region of sampling for Sahiwal cows and constantly exposed to tropical pathogens.
The molecular pathways impacted in both the cattle types adapted to extreme agro-climatic regions included MAPK signaling pathway, Electron transport chain, apoptosis, IL2 signaling, IL3 signaling, IL6 signaling, TLR signaling pathway, Delta Notch signaling pathway, Myometrial relaxation and contraction and TNF-alpha and NF-kB signaling pathway. Among these MAPK signaling pathway known to be involved in triggering the cellular response to environmental stress was the most impacted pathway. A total of 39 genes (14 up and 25 down regulated) were matched in the pathway entity list (Fig. 4). The upregulated genes TAB2, TGFB2, MAP4K1, MINK1, ATF2, MAPT, ARRB2, MAPK3, KRAS, CDC25B, PTPN7, JUND, MAP4K3 of this particular pathway in Ladakhi cows may be associated with the high-altitude adaptation. Some of the recent studies have also shown indication of MAPK signaling pathway in high altitude adaptation especially in Tibetans 42 . Further, in recent past, relationship of this pathway with hypoxic response has also been documented 43 . Taken together, the series of pathways impacted in the present data set along with distinct transcriptome profile under hypoxic (high altitude) and normoxic (mean sea level) conditions indicated signatures of adaptive evolution of these two cattle types in response to diverse environments.
The qPCR analysis also substantiated higher expression of hypoxia related target genes such as HIF-1α, EPAS-1, VEGFA, NOS2, ECE-1, GLUT1, HK2, TNFα, GRα in PBMCs of Ladakhi cattle in comparison to breed from low altitude region suggesting pivotal role of these genes in high altitude adaptation. The enrichment of genes reported for high altitude adaptation in few other species like, Tibetan antelopes, Tibetan wild boars, yaks and ground tits seems to be corroborated with the present findings, indicating existence of common evolutionary mechanism for high altitude adaptation [44][45][46][47] . Based on the generated data, local native Ladakhi cattle are better adapted to withstand the hypoxic condition than exotic or cross-bred populations. The data indicated that hypoxia related genes get accumulated under hypoxic conditions and probably are essential adaptive component for the animals surviving at high altitude. On the other hand, significantly higher expression of molecular chaperons HSP27, HSP70 and HSP90 genes in tropically adapted Sahiwal cattle than Ladakhi and other exotic breeds from high altitude could be due to induction of molecular chaperons at high ambient temperature prevalent under tropical habitat. Overall, the expression pattern of several target genes generated by qPCR corroborated well with the microarray data set.
In conclusion, the comparative data analysis indicated distinct transcriptome signature of two cattle types native to extreme climatic conditions. Based on data generated in the present study, native cattle of Ladakh region was found to be genetically distinct from native cattle adapted to tropical region of India

Methods
Ethics statement and animal selection. The blood sampling of animals were performed in accordance with the relevant guidelines and regulations as approved by Institutional Animal Ethics Committee (IAEC) of National Bureau of Animal Genetics Resources (NBAGR), Karnal. In the present study, transcriptome profile was generated in a total of 11 PBMCs samples; 5 heifers of Ladakhi cattle adapted to high altitude hypoxic condition; and 6 heifers of Sahiwal cattle from low altitude normoxic condition. The sampling of both the cattle types was accomplished in the month of August. About 25-30 ml of whole blood was collected aseptically from external jugular vein in sterile EDTA coated vacutainer tubes. The blood samples were then taken to laboratory in an icebox and were subjected for PBMC isolation.
PBMC isolation and RNA extraction. The PBMCs were isolated from the whole blood using density gradient centrifugation method as described in our previous reports 48 . In brief, whole blood was diluted with 1 × PBS (Ca 2+ and Mg 2+ free; Hyclone, Utah) in 1:1 ratio and gently over laid on Histopaque-1077 (Sigma-Aldrich Inc., USA) and centrifuged at 400 RPM for 30 mins at room temperature. After removing the buffy coat, lymphocyte pellet was treated with 2 ml of chilled RBC lysis buffer. The reaction was stopped by adding 8.0 ml of 1 × PBS (Ca 2+ and Mg 2+ free; Hyclone, Utah) followed by centrifugation at 260 g for 10 min. After getting a clear white pellet, the cells were washed twice with 1 × PBS and were suspended in 1.0 ml of ice cold Trizol reagent (Invitrogen, Carlsbad, California). Total RNA was extracted and purified from 11 PBMC samples using ice cold Trizol reagent according to the manufacturer's instructions (Invitrogen, Corp., CA). The traces of genomic DNA were removed by RNase free DNase treatment (Qiagen, Germany) using RNeasy Mini kit columns (Qiagen, Germany). The quality and concentration of extracted RNA was measured using NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies) and bioanalyzer (Bio-Rad, USA). cRNA labeling and hybridization with bovine 44 K oligonucleotide chip. The transcriptome signature was generated using Agilent whole genome bovine 44 K chip harbouring 60 mer oligos and 21520 entities. The entire microarray related procedures viz; labeling, hybridization, washing and scanning were followed as described in our previous study 49 . The entire procedure is shown in Fig. 5. Labeling and microarray processing was performed as per "One color protocol using Agilent's low input quick amp labeling kit" (Agilent Technologies, Santa Clara, CA). 2 µl of total RNA (100 ng/µl) samples were labeled with T7 promoter primer at 65 °C for 10 min and incubated in ice for 5 min. cDNA was constructed from labeled RNA samples after adding cDNA master mix (5× First Strand Buffer, 0.1 M DTT, 10 mM dNTP mix and Affinity Script RNase Block Mix) followed by incubation at 40 °C for 2 h, 70 °C for 15 min and final incubation on ice for 5 min. The cDNA was reverse transcribed to synthesize cRNA (complimentary RNA) and was amplified by adding transcription master mix (5× Transcription Buffer, 0.1 M DTT, NTP mix, T7 RNA Polymerase Blend and Cyanine 3-CTP) followed by incubation at 40 °C for 2 h. The amplified labeled cRNA was purified (RNeasy mini column kit, Qiagen, Germany), and quantified to obtain the yield of cRNA {μg cRNA yield = (Concentration of cRNA) ×30 μl (elution volume)/1000} and specific activity of Cy3 {pmol Cy3 per μg of cRNA = (Concentration of Cy3/Concentration of cRNA) ×1000}. For hybridization, 1.65 μg of linearly amplified and Cy3-labeled cRNA were fragmented using fragmentation mix (10× Blocking Agent and 25× Fragmentation Buffer) and incubated at 60 °C for exactly 30 min. After incubation the fragmented RNA samples were immediately transferred on ice for one min and added 2× GEx hybridization buffer to stop the reaction. The fragmented samples were carefully dispensed to the array on gasket slide placed in hybridization chamber base without introducing air bubble. Slowly placed an array onto the gasket slide with "Agilent labeled barcode facing down and the numeric barcode facing up. Placed the assembled slide chamber in rotisserie in a hybridization oven allowing it to rotate at 10 rpm for 17 hrs at 65 °C. Post hybridization, the microarray slides were disassembled in gene expression (GE) wash buffer 1 (pre warmed overnight at 37 °C) containing 0.005% Triton-X to reduce the possibility of array wash artifacts. This was followed by second wash with GE Wash Buffer 2 for 5 min at room temperature.
Scanning, data acquisition and analysis. The slides were scanned immediately after washing on Agilent DNA Microarray Scanner (G2505B) using one colour scan setting. The signal intensities were checked from the digit images using Agilent Feature Extractor Software Version 9.5 (Agilent Technologies). The microarray raw data files in.txt format were imported to GeneSpring GX 13.1 (Agilent Technologies) and subjected for quantile normalization which is highly effective in reducing variation between arrays. Normalized data was analysed for statistically significant gene expression differences between Ladakhi cows and Sahiwal cows to the fold change statistic that uses overall gene expression variation to calculate a gene-specific variance. To keep the number of false positives to alpha = 0.05, a false discovery rate adjustment was used 50 . Genes that were statistically significant in expression between Ladakhi cows and Sahiwal cows were grouped according to fold expression (up-vs. down-regulation). The data was further subjected to hierarchical clustering and k-means clustering based on entities (genes), gene ontology enrichment analysis and pathway analysis.  Holstein Frisian crosses (HFX), Jersey cows (JYC) from high altitude region; and Sahiwal cows (SAC), Karan Fries cows (KFC), Holstein Friesian cows (HFC) from tropical region using qPCR. For optimal gene expression analysis, qPCR data of target genes was normalized using a total of 10 previously known candidate reference genes from different functional categories viz., GAPDH, RPL4, EEF1A1, RPS9, HPRT, UXT, HMBS, B2M, RPS15 and ACTB. The primers for each of the ICGs and target genes were either selected from literature or designed using Primer express 3.0 software (Applied Biosystem). Primer details for all the ICGs and target genes are provided as supplementary information (Tables S7 and S8). The accuracy of primer pairs was also ensured by the presence of a unique peak during the dissociation step at the end of qPCR cycle. The qPCR was performed using Light Cycler 480 instrument (Roche, Germany) as described in our previous reports 51 . The data was acquired using the 'second derivative maximum' method as computed by the Light Cycler Software 3.5 (Roche Diagnostics) and subjected for subsequent analysis.
Data normalization. qPCR data of target genes in high altitude and tropically adapted cattle was normalized by selecting best suited ICGs utilizing geNorm, NormFinder and Best keeper softwares [52][53][54] .The mean value was calculated using relative quantification 2 −ΔΔCT method 55 . The differences between groups were tested by Tukey's Multiple Comparison Test. P values less than 0.05 and 0.01 were considered significant. Statistical test was performed using GraphPad PRISM version 5.0 (La Jolla, CA, USA).