Impact of Heat Stress on Cellular and Transcriptional Adaptation of Mammary Epithelial Cells in Riverine Buffalo (Bubalus Bubalis)

The present study aims to identify the heat responsive genes and biological pathways in heat stressed buffalo mammary epithelial cells (MECs). The primary mammary epithelial cells of riverine buffalo were exposed to thermal stress at 42°C for one hour. The cells were subsequently allowed to recover at 37°C and harvested at different time intervals (30 min to 48 h) along with control samples (un-stressed). In order to assess the impact of heat stress in buffalo MECs, several in-vitro cellular parameters (lactate dehydrogenase activity, cell proliferation assay, cellular viability, cell death and apoptosis) and transcriptional studies were conducted. The heat stress resulted in overall decrease in cell viability and cell proliferation of MECs while induction of cellular apoptosis and necrosis. The transcriptomic profile of heat stressed MECs was generated using Agilent 44 K bovine oligonucleotide array and at cutoff criteria of ≥3-or ≤3 fold change, a total of 153 genes were observed to be upregulated while 8 genes were down regulated across all time points post heat stress. The genes that were specifically up-regulated or down-regulated were identified as heat responsive genes. The upregulated genes in heat stressed MECs belonged to heat shock family viz., HSPA6, HSPB8, DNAJB2, HSPA1A. Along with HSPs, genes like BOLA, MRPL55, PFKFB3, PSMC2, ENDODD1, ARID5A, and SENP3 were also upregulated. Microarray data revealed that the heat responsive genes belonged to different functional classes viz., chaperons; immune responsive; cell proliferation and metabolism related. Gene ontology analysis revealed enrichment of several biological processes like; cellular process, metabolic process, response to stimulus, biological regulation, immune system processes and signaling. The transcriptome analysis data was further validated by RT-qPCR studies. Several HSP (HSP40, HSP60, HSP70, HSP90, and HSPB1), apoptotic (Bax and Bcl2), immune (IL6, TNFα and NF-kβ) and oxidative stress (GPX1 and DUSP1) related genes showed differential expression profile at different time points post heat stress. The transcriptional data strongly indicated the induction of survival/apoptotic mechanism in heat stressed buffalo MECs. The overrepresented pathways across all time points were; electron transport chain, cytochrome P450, apoptosis, MAPK, FAS and stress induction of HSP regulation, delta Notch signaling, apoptosis modulation by HSP70, EGFR1 signaling, cytokines and inflammatory response, oxidative stress, TNF-alpha and NF- kB signaling pathway. The study thus identified several genes from different functional classes and biological pathways that could be termed as heat responsive in buffalo MEC. The responsiveness of buffalo MECs to heat stress in the present study clearly suggested its suitability as a model to understand the modulation of buffalo mammary gland expression signature in response to environmental heat load.


Introduction
Heat stress is a significant issue for many livestock enterprises, particularly for dairy. The negative impact of heat stress on dairy animals have been well documented [1][2][3][4][5][6] that includes reduced feed intake, milk production, milk quality, fertility and increased respiratory, heart rates, panting activity, peripheral blood flow, sweating etc. These impacts can result in significant loss of income and increase in management costs. Estimated average losses resulting from heat stress for the US dairy industry were about US$900 million per annum [7]. Such apprehensions have increased even more in recent years with the increased global warming and strive of dairies to maximize milk production.
In India, dairy industry is mainly dependent on cattle and buffaloes with buffaloes contributing more than 50% of milk production in addition to draught power and meat. Depending upon the habitat and chromosome number, the domesticated water buffaloes has been classified into two main categories, namely riverine and swamp. The riverine buffaloes (2n = 50) primarily developed for milk, meat and draught purposes, are mainly found in India and countries to the west of India like Pakistan, Iran, Iraq, Egypt, Bulgaria, Italy etc. The swamp buffaloes (2n = 48) on the other hand have been primarily developed for draught purposes and are mainly found in South East Asia and China. The riverine buffalo is established as an economically important species not only in India but in many Asian and Mediterranean countries and hence its genetic improvement ranks high amongst agricultural research needs of these countries.
Of the total world buffalo population, 97.1% buffaloes are distributed predominantly in Asian countries, while only 2.9% are found in rest of the world. The largest buffalo populated countries are India, Pakistan and China, of which India and Pakistan account for almost 65% of the total world buffalo population. India is bestowed with immense richness of buffalo diversity with over 96.9 million buffaloes constituting more than half (56.8%) of the total world (170.4 million) and 58.9% of Asian (165.4 million) buffalo population, respectively [8].
The role of buffaloes as a major milk producing species is well recognized in the Indian subcontinent, especially in India and Pakistan where 92% of the world buffalo's milk is produced. With higher total solids (protein, fat, minerals) of 18-23% as compared to 13-16% in cow milk, buffalo milk confers advantage in the preparation of cheese, curd and other dairy products [8].
Although buffaloes are well suited to hot and humid climates and muddy terrain, but they exhibit signs of great distress when exposed to direct solar radiation or while working in the sun during hot weather. It is reported that milk yield, growth rate and fertility of buffaloes get reduced during periods of high ambient temperature [9]. This is due to the fact that buffaloes absorb a great deal of solar radiations because of their dark skin and sparse coat or hair. In addition to that they possess less efficient evaporative cooling system due to their rather poor sweating ability as buffalo skin has one-sixth of the density of sweat glands than that of cattle skin. Under heat stress, buffalo's body temperature, pulse rate, respiration rate and general discomfort increase quickly [9]. It further evokes a series of drastic changes in biological functions including decreased intake efficiency and utilization of feed; disturbances in metabolism of water, protein; and changes in energy and mineral balances, enzymatic reactions, hormonal secretions and blood metabolites. All this result in impairment of growth, reproductive as well as productive performance and hence overall well-being [10].
Even though buffaloes play a major role in sustaining the economy of India's agriculture, still the full genetic potential of this species has not yet been well exploited, primarily due to lack of basic data on buffalo genome. The progress in the development of cattle, pig, sheep and chicken genome database continues worldwide, but similar efforts for buffalo are relatively of lesser magnitude. Previously, efforts have been made understand response of buffalo to heat stress on the basis of anatomical differences and physiological parameters; however, genetic components, gene regulatory pathways, alterations in gene expression and molecular mechanisms underlying changes in heat stress response are not well established in this species. The present scenario calls for mining of genomic data on this important genetic resource to understand about the genes, pathways and molecular mechanism of heat stress response. Microarray technology proven to be a powerful tool for the simultaneous expression analysis of thousands of genes in a tissue could be the most appropriate tool for this purpose [11].
Heat stress is an important environmental stimulus that can affect the expression of mRNA in different animal tissues or cells. Mammary gland, an important economical organ of dairy animals has always drawn attention of the scientists for over a century because of special function for milk synthesis and secretion. To understand the physiological function of mammary gland, mammary tissue cells or explants have been widely-used as in vitro models. However, when using tissue explants, it is inherently difficult to distinguish between primary mitogens and secondary regulators of mammary gland function and development. To circumvent most of these difficulties, emphasis has been placed on cell culture methodologies to study growth regulation, hormonal responsiveness, or biochemical properties of mammary epithelial cells [12]. 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 are the predominant cell types in lactating mammary gland, changes in their genes expression could provide an insight of the mammary gland mechanism. Hence buffalo MECs could be an interesting in-vitro model to delineate the genes whose expression is significantly modulated due to heat stress challenge. To the best of our knowledge, no systematic initiative has been attempted so far to highlight the molecular mechanism or identify gene networks and molecular pathways associated with heat stress response in buffaloes using MEC. Transcriptomic adaptations of buffalo mammary epithelial cells during heat stress can be easily and efficiently identified utilizing bovine arrays. Considering the above issues, the present study was planned to generate global expression profile of buffalo MECs during normal and heat stressed state and identify molecular pathways significantly regulated in heat stressed MECs

Buffalo MECs primary culture and heat treatment
The buffalo mammary gland tissues of approximately 5gm were obtained from a healthy adult buffalo from Gazipur abattoir 28.734190N and 77.272830E, New Delhi, India. The primary MECs were cultured using DMEM/F12, supplements and growth conditions as described earlier [13]. After several passages, 80% confluent buffalo MECs on 10th passage were distributed in collagen treated 12-wellplates (Corning, USA) in two sets with one plate assigned as control (kept at 37°C all the time) and other plate as treated (exposed to 42°C). Initially, cells were incubated at 37°C with 5% CO 2 to stabilize the culture. Subsequently, the plate marked as treated was exposed to 42°C for one hour to simulate heat stress (HS) condition. After 1h, the cells were allowed to recover at 37°C, 5% CO 2 and harvested by trypsinization at different time points (30m, 2h, 4h, 8h, 12h, and 24h). The samples from control (CTR) plates were also trypsinized and harvested at the same time points corresponding to the treated plates. Followed by exposure to heat stress, cell viability and growth characteristics of buffalo MECs in normal and heat treated samples were determined using commonly used trypan blue dye exclusion method.

Estimation of cellular proliferation towards heat stress to MECs
The induction and inhibition of proliferation of buffalo MECs under normal and heat stress condition in in vitro model was evaluated using MTT assay kit (Cayman, Ann Arbor). Cells were seeded in triplicate with a density of 5x10 3 cells/well in 100 μl of culture medium in 96 well plates (Corning, USA) and cultured for 24-48 h at 37°C, 5% CO 2 . Cells in control plates were maintained at 37°C, 5% CO 2 throughout the time-course, while those in treatment plates were exposed at 42°C, 5% CO 2 for 1 h and then shifted to 37°C, 5% CO 2 . The post heat treated cells were harvested at different time points (0h, 30 m, 2h, 4h, 6h, 8h, 12h, 24h and 48h) and cell proliferation assay was performed following manufacturer's instructions.

Cell apoptosis detection by flow cytometry
The effect of heat stress on cell apoptosis of buffalo MECs was determined using ANNEXIN V-FITC / 7-AAD Kit (Beckman Coulter, France). Buffalo MECs were cultured, heat stressed and harvested at different time points as mentioned previously. The cell apoptosis assay was conducted as per manufacturer's instructions. For Annexin V-FITC, the maximum fluorescence absorption and emission were recorded at wavelength of 490nm and 525nm, respectively, whereas, for 7-AAD viability dye, the maximum fluorescence absorption and emission for DNA/7-AAD complex were recorded at 543 nm and 655 nm, respectively.

Real-time quantitative PCR (qPCR)
Total RNA extraction from MECs harvested at CTR, 30m, 2h, 4h, 8h, 12h, 24h after heat stress and cDNA synthesis was performed as described in our previous studies [13]. Primer details for all the heat shock protein and apoptotic gene family are provided as supplementary information (S1 Table). 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. qPCR was performed using Light Cycler 480 instrument (Roche, Germany) as described in our previous reports [13]. 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. In our previous study, RPL4, EEF1A1and RPS23 were observed to be most stably expressed genes in heat stressed buffalo MEC [13] and identified as appropriate panel of reference genes for normalization of qPCR data utilizing geNorm, NormFinder and Best keeper softwares [14][15][16] The geometric averages of these genes were f utilized for normalization of qPCR data generated in the present study. To measure the relative changes in gene expression data, 2 -ΔΔCT method was used [17].

Generation of microarray based gene expression data
As a first step towards One-Color Microarray-Based Gene Expression Analysis using Low Input Quick Amp Labelling kit (Agilent Technologies, Santa Clara, CA), the RNA samples (CTR, 30 m, 2h, 4h, 8h, 1h, 24 h) were labelled with T7 promoter primer (65°C for 10 min followed by incubation in ice for 5 min). cDNA was constructed from labelled RNA samples after adding cDNA master mix (5X 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 cRNA synthesis and amplification was performed by adding transcription master mix (5X 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 labelled cRNA was purified (RNeasy mini column kit, Qiagen, Germany), quantified {μgcRNA yield = (Concentration of cRNA) x 30 μL (elution volume)/ 1000} and checked for specific activity {pmol Cy3 per μg cRNA = (Concentration of Cy3/ Concentration of cRNA) x 1000}. All the samples exhibited yield and specific activity values higher than the minimum value of 1.65 and 6, respectively using four pack microarray format. For hybridization, 1.65 μg of linearly amplified and cyanine 3-labeledcRNA were fragmented using fragmentation mix (10X Blocking Agent and 25X Fragmentation Buffer) and incubated at 60°C for exactly 30 min. The fragmented RNA samples were immediately transferred on ice for one minute and 55μl of 2x GEx Hybridization Buffer HI-RPM was added to stop the fragmentation reaction. The fragmented samples (110μl volume) were loaded to the array placed in slide chamber. The assembled slide chamber placed in rotisserie was put in a hybridization oven set to rotate at 10 rpm and 65°C. The hybridization was allowed for 17 hours. Followed by hybridization, microarray slide was disassembled in GE wash buffer 1 (pre warm overnight at 37°C) containing Triton X-102 (0.005%),washed with fresh GE wash buffer 1 for 1 min followed by second wash with GE Wash Buffer 2 for 5 min at room temperature. The slide was dried and scanned immediately.
Scanning, data acquisition and data analysis. Immediately after washing, the slides were scanned on GenePix-4000B (Molecular Device) microarray scanner using one colour scan setting with 5μm resolution at wavelength of 532 nm (Cy3). After scanning, the images were collected as 16 bit images and finally stored as tif image files. These files were further subjected to feature extraction using Agilent Feature Extractor Software Version 9.5 (Agilent Technologies) software. GeneSpring GX 12.6 (Agilent technologies) was used to perform the analysis of raw data obtained from feature extraction software. The data was normalized to the 75 th percentile. The microarray analysis workflow employed in the present study consisted of following steps viz., loading of raw data into GeneSpring software, normalization of microarray data, quality check of microarray data, filtration of probe sets, identification of differentially expressed genes based on fold change criteria, classification/clustering of genes, gene ontology enrichment analysis and pathway analysis.

Establishment of primary culture of buffalo MECs
The primary culture of buffalo mammary epithelial cells was achieved following enzymatic digestion of buffalo mammary tissue. In the initial stage heterogeneous population of epithelial and fibroblast-like cells was observed ( Fig 1A) &that was further purified using selective trypsinization procedure. Being more sensitive to trypsin treatment as compared to MECs, fibroblasts cells were removed and homogeneous population of mammary epithelial cells was obtained. On microscopic evaluation the cells showed normal characteristics of MECs (Fig 1B-1E). Upon reaching the confluence stage (5-6 days after seeding), the cells formed a monolayer & aggregated to form typical cobble stone morphology of MECs ( Fig 1C). During the post confluence period, a number of floating dead cells were observed indicating the occurrence of contact inhibition (Fig 1F) The MECs were passaged up to 10 times during continuous sub culturing. Following the similar methodology, primary MEC culture has been established in different livestock species [18][19][20][21][22] and in bovines, established MEC line has also been used to study the effect of heat stress in vitro [12,[23][24][25][26].
Immunofluorescent staining for assessing purity of cultured MECs. Purity of mammary epithelial cells was ensured by evaluating the expression of cytokeratin 18 and vimentin proteins using immune fluorescence staining. The anti-cytokeratin 18 antibody detects the expression of cytokeratins, known to be specific to epithelial cells whereas anti-vimentin antibody detects vimentin specific to fibroblast cells. The staining with anti-cytokeratin 18 antibody revealed strong signals indicating that high percentage of cells are of epithelial lineage (Fig 2A  & 2B). Although, few cells were stained with anti-vimentin, but the signals were week ( [20]. The immuno staining pattern with homogeneity and strong signal for cytokeratin 18, provided the evidence that the buffalo primary culture established in the present study mainly consisted of MECs with no or very little contamination of fibroblast cells. Our findings are in accordance with other reports across species, wherein it is reported that expression of cytoskeleton is specific for epithelial cells [19-20, 22, 27-32].

Effect of heat stress on cell viability, proliferation, cytotoxicity and cellular apoptosis
In order to assess the effect of heat stress on MECs, cell viability & growth parameters were recorded using trypan blue dye, MTT assay & flow cytometric analysis. The cell viability estimation by trypan blue exclusion method showed decreased cell number and viability in heat  (Fig 3). A decrease in cell viability was recorded immediately (30 m) after heat treatment and it remained significantly low up to 8 h of recovery phase. Thereafter, the cell viability increased during later stages of recovery (>12 h). The low percentage of viable cells observed immediately after the heat stress might be attributed to the adverse effect of heat stress on cell membrane structure or/and to cell necrosis or impairment of cell growth, while the recovery in cell viability during the later stages, could be attributed to the restoration of cell survival mechanism by mammary cells. The non-linear response of the cell viability to heat stress observed in the present study is in accordance with other studies where similar response of the cell viability to heat stress in different in vitro cell culture models was examined [23,[33][34][35][36][37].  Further, the effect of heat stress on cell viability & proliferation of MECs was also determined using MTT based assay. The heat stressed MECs(exposed to 42°C for 1 hour) had significantly lower cell viability than unstressed (control) cells. There was significant reduction in cell proliferation immediately after the heat stress & remained low throughout the time course (Fig 4). The lower formazan crystal formation in damaged cells indicated loss of cell proliferation efficiency in post heat stress samples as compared to unstressed cells. Hence decrease in cell proliferation efficiency in heat stressed buffalo MEC might be due to loss of plasma membrane potential. Similar observation has been reported in bovine MECs where thermal stress inhibits the proliferation rate [12,24,36,38].
To check the impact of heat stress on induction of apoptosis or programmed cell death in buffalo MECs, flow cytometric analysis was carried out at various time points of recovery phase post heat stress. The annexinV-FITC /7-AAD dyes were used to detremine the percentage of apoptotic and dead cells. The analysis yielded three different types of cell populations; viable cells with no staining, early apoptotic cells stained positive with annexin-V-FITC dye and late apoptotic/ dead cells stained with both annexin-V-FITC as well as 7-AAD dyes. The cells detected by annexin-V-FITC appearedto be undergoing an early apoptotic process whereas the cells stained with 7-AAD reflected damaged plasma membrane. At 30 min of recovery phase, 5.21% of cells were found to be annexin-V positive suggesting immidiate induction of apoptosis after heat stress. The percentage of MEC undergoing early apoptosis showed progressive increase with increase in time period post heat stress. The percentage of early apoptotic and apoptotic cells was maximum at 24 h post stress with values of 6.79% and 12.95%, respectively. Similarly, the proportion of dead cells also increased with increase in time interval post heat stress (Figs 5 and 6).
Based on present analysis and those reported by other workers [24, 26, 36, and 38] it could be suggested that activation of apoptotic pathway during heat stress could be the major cause of cell death.

mRNA expression profile of heat shock proteins (HSPs)
The expression pattern of heat shock protein genes (HSPs) showed immediate increase in mRNA level of all the analyzed HSPs in response to in vitro heat stress. Each member of the  studied molecular chaperons (HSP27, HSP40, HSP70, HSP60, and HSP90) responded well within 30 min. Most of them showed maximum increase in mRNA expression at 2h-4h, remained elevated till 12h post heat stress and eventually (16h-48h post heat stress) declined to the level of unstressed MECs (S1 Fig). The increase in expression of HSP genes suggested the responsiveness of buffalo MECs to heat stress in vitro. The five major heat shock protein genes were selected for analysis as these are the molecular chaperons that ensures the correct protein folding and apoptosis regulation during physiological stressful conditions [39]. Amongst all HSPs, HSP70 was identified as the most predominant form of transcripts induced in buffalo MECs due to heat stress. The early induction of chaperone activity in buffalo MECs suggested the induction of thermoregulatory mechanism during early stages of cellular response to heat stress. Our findings are in accordance with previous studies conducted across different cell types where it was reported that heat shock remarkably impacts induction of HSP genes; reduction of cellular growth and heat tolerance ability of MECs; and down-regulation of genes associated with cellular metabolism and secretion of milk protein [23,[39][40][41][42][43].
mRNA expression profile of apoptotic genes. The mRNA expression of two genes related to apoptotic pathway, anti-apoptotic (BCL2) and pro-apoptotic (BAX) were measured in heat stressed MECs. The results showed that the amount of BAX mRNA, which is known as apoptotic activator increased immediately (30m) and continued to increase even till48 h of recovery period after heat stress (S2 Fig). Its expression was highest at 48h time point with 3.439 fold greater followed by 24h time point with 2.25 fold greater than unstressed (control) MECs. Conversely, anti-apoptotic gene (BCL2) showed decreased expression up to 24h (0.638 fold) before being recovered close to basal level at 48h (1.037) of heat stress (S2 Fig). The BAX/BCL2expression ratio was highest (3.316 fold) after 48h of heat stress, indicating high rate of apoptosis. The expression kinetics of BAX and BCL2 genes strongly indicate the occurrence of heat stress induced apoptosis in the buffalo MECs. The results were in accordance with flow cytometric based apoptotic data showing highest cell death and apoptosis during late time points recovery stages of post heat stress. Several studies have suggested that pro-apoptotic genes like BAX accelerates programmed cell death by binding to, and antagonizing the apoptosis repressor BCL2 gene. This gene interacts with components of permeability transition pore (PTP) protein complex of mitochondria that consists of the voltage-dependent anion channel in the outer mitochondrial membrane and adenine nucleotide translocase in the inner mitochondrial membrane. Under stress conditions, the interaction of BAX gene with PTP causes the formation of large pore due to conformational changes resulting in the release of cytochrome c and other pro-apoptotic genes that trigger apoptosis. On the other hand, BCL2 gene prevents the opening of PTP and also encodes an integral outer mitochondrial membrane protein that blocks the apoptotic death. The observed expression profile of BAX and BCL2 genes strongly suggests that heat stress induces cell apoptosis by triggering perturbation of mitochondrial function. Based on our data and information available in the literature, it could be stated that mitochondrial functioning is key to regulate apoptosis and cell death in buffalo MECs during heat stress. During heat stress, pro-apoptotic signals trigger a change in mitochondrial permeability which results in release of mitochondrial proteins into the cytoplasm that might be crucial for apoptosis. In line with our study, few others have also shown induction of cellular apoptosis in different cell lines on heat shock treatment [24,44].

Analysis of microarray data for identifying differentially expressed genes
In this study, an attempt was made to obtain a global picture of in vitro heat stress response by investigating transcriptome profile of mammary epithelial cells of buffaloes. The Agilent 44K bovine oligonucleotide array which contain~20,000 probe sets was employed to characterize the gene expression changes in buffalo MEC in response to heat stress. RNA electropherogram profile using bio analyzer (S3 and S4 Figs) revealed that all samples are of good quality and intact. To ensure success in microarray hybridization, the yield and specific activity of Cy3-labeled complementary RNA (cRNA) was measured in terms of cyanine 3 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 μgcRNA yield and concentrations of cRNA (ng/μL) and cyanine 3 (pmol/μL) was used to determine the specific activity In all the MEC samples (CTR, 30 (Fig 8).
In addition, the Venn analysis showed distribution of genes at 3 fold change and identified list of DEG that were commonly upregulated (153) and down regulated (8) across all time points post heat stress (Fig 9). These genes whose transcriptional pattern changed due to heat stress across all time points were termed as heat responsive genes. The Venn diagram analysis helped to identify the genes that are differentially expressed or remained commonly expressed in unstressed (CTR) and heat stressed MEC at individual time points after heat stress. The analysis revealed that there was relatively a greater transcriptional response during the early time points as evident by 894, 2267, 835 DEG at 30 min, 2h, 4h respectively in comparison to 138, 80, 732 and 227 DEG during late hours i.e. 8h, 12h, 16h and 24h (Fig 10). Further, the Venn analysis was also carried out to find out number of genes that were induced specifically during early and late time points post heat stress in comparison to CTR. During early time points viz. 30 min, 2h and 4h, a total of 472, 416 and 118 genes were up regulated respectively in comparison to unstressed (CTR) cells (Fig 11). At later time points viz. 8h, 12h, 16h and 24h, a total of 124, 67, 471 and 172 genes were induced in comparison to unstressed cells.
The top 50 up-regulated and top 50 down-regulated transcripts filtered at 3-fold cut off are listed in Table 1 and Table 2, respectively.
Among most up-regulated genes, BOLA was highly up-regulated gene during heat stress in buffalo MECs, probably because the presence of MHC class I raises the possibility of cells in autoimmune disease state [45]. Jorge et al [46] have also reported the induction of BOLA gene during early logarithmic growth phase of Escherichia coli in response to heat stress.  Furthermore, a significant fold change increase in the expression of BOLA gene after heat stress in biofilm and planktonic stages of growth in Escherichia coli has also been reported [47]. The second most up-regulated gene was MRPL55 (mitochondrial ribosomal protein) that is implicated in protein synthesis within the mitochondrion and cell cycle progression [48]. Similar to our findings, the heat shock to larval stages of Drosophila eye indicated over expression of MRPL55 transcripts [49]. Another most up-regulated gene PFKFB is responsible for maintaining the cellular levels of fructose-2, 6-bisphosphate-a key regulator of glycolysis. Earlier also, significant induction of PFKFB3 mRNA under hypoxic stress was reported in several human and mouse cell lines [50][51]. In addition, PSMC2 (proteasome 26S subunit, ATPase, 2) also showed induction in expression pattern during post heat stress in buffalo MECs. This gene is known to be involved in the ATP-dependent degradation of ubiquitinated proteins. Further, it has been observed that PSMC2 gene also plays an important role in cellular growth and proliferation [52]. The up-regulation of proteosome subunit might indicate the immediate response of PSMC2 during weakening of ubiquitin-proteosome system resulted in accumulation of abnormal proteins that in turn might confer growth and development in buffalo MECs. Additionally, it has been reported for human chronic myelogenous leukemia cell line that HSP70 is involved with the dissociation and reassociation of the 26S proteasome during adaptation to oxidative stress [53]. These findings can be correlated with the present study where HSP70 and PSMC2 were highly up-regulated, explaining the co-relation between both genes during stressful conditions. Another most up regulated gene observed in present study was ARID5Athatplays an important role in development, tissue-specific gene expression, and regulation of cell growth [54]. Along with ARID5A, IL6 was also up regulated conferring the findings that ARID5A has a role in stabilization of IL6 mRNA for promotion of inflammatory responses [55]. Among all down-regulated genes, COL4A1 (collagen, type IV, alpha 1) gene showed high reduction in expression under heat stress. This gene specifically inhibits endothelial cell proliferation and expression of HIF-1alpha and ERK1/2 and plays a major role in p38 MAPK activation. Down regulation of collagen, the main component of the extracellular matrix hasalso been reported for sea anemones during heat stress [56]. The next highly down-regulated gene was IGFBP5 (insulin-like growth factor binding protein 5). It has a role in tissue turnover by reducing the availability of the survival factor IGF-1 as well as increasing extracellular matrix degradation, thereby causing apoptosis and tissue remolding [57]. Similar to our observation, the reduced mRNA expression of IGFBP5has been reported in heat stressed cows as well [58]. Another major affected gene in heat stressed buffalo MECs was CABP2 (calcium binding protein 2), which is a calcium binding protein and is an important component of calcium mediated cellular signal transduction [59]. Along with above described major genes, IREB2 (ironresponsive element binding protein 2) also ranked higher among down-regulated genes. The binding of IREB2 to transferrin receptor mRNA inhibits the degradation of otherwise rapidly degraded mRNA. Its reduced expression was observed in plants during oxidative stress [60] which could be correlated with our study. In addition, FBXO22 gene, a F-box protein, constitute one of the four subunits of the ubiquitin protein ligase complex called SCFs (SKP1-cullin-  F-box), which function in phosphorylation-dependent ubiquitination and are thought to be involved in degradation of specific proteins in response to p53 induction. Similar to our observations, the evolutionarily conserved Arabidopsis thaliana F-box protein showed reduction in transcriptional expression profile during temperature stress [61]. Additionally, NCAM1 gene also showed a reduction pattern in expression level during post heat stress. It is involved in cell-to-cell interactions as well as cell-matrix interactions during development and differentiation and play an important role in immune surveillance. Our findings were in accordance with the reduced expression of NCAM1 gene in adult mice as a consequence of heat stress [62]. Identification of heat stress responsive genes from transcriptome data. Amongst the several hundred genes induced or repressed due to heat stress in vitro, an effort was made to filter out genes associated with; heat shock protein family, apoptosis; immune and oxidative stress ( Table 3). The heat map view for list of genes identified under these categories is presented in Fig 9. As expected, the whole set of genes of heat shock family viz., HSPA6, HSPB8, DNAJB2, HSPA1A etc. were up-regulated in MEC at most of the time points post heat stress. The expression of these genes was more during the early phase of heat stress as compared to late time points post heat stress. Our findings are in accordance with previous studies where heat stress led to induction of HSP genes, and down-regulation of genes associated with cellular metabolism and cellular growth [23,[39][40]. Similar to our study, induction in HSPs expression in different cell/tissue types viz., leukocytes/lymphocytes [41,43,63], bovine endometrial tissue [64,65], bovine concept uses [65], bovine MECs [23], buffalo lymphocytes [66]due to heat stress has also been reported. Heat stress has also been shown to trigger an increase in HSPs in virtually all the vertebrates, including mice [67,68] domestic goats [69], humans [70,71], juvenile Hamadryas baboons [72], common carp [73], domestic chickens [74][75][76][77] and domestic turkey [78]. Our observations thus supported the idea that HSP70 can act as reliable, sensitive biomarker of thermal stress [72,79]. Similarly, several apoptosis related genes were also found to be up-regulated on heat stress. Up-regulation of apoptotic genes could result in disruption of mitochondrion transmembrane potential, thereby causing cytochrome c release leading to the induction in apoptosis. The data on induced expression of apoptotic genes immediately after heat stress suggests that cellular mechanism may not provide protection to the MECs against heat-induced apoptosis in buffalo while during recovery period of heat stress they probably helped the cell to cope with hyperthermia through clearance of damaged proteins. Our findings are in accordance with some previous reports where heat shock showed induction in expression of apoptotic genes in in vitro culture models of buffalo embryos [80] and cat skin fibroblasts [81]. In addition, immune system related genes BOLA, IL1-B, BOLA--DRA, TNF, IL1A, IL10, and CXCL2 and IL6 were also upregulated by heat shock, supporting the RT-qPCR analysis carried out in the present study. Similar to our findings, the induction of immune related gene expression in intestinal mucosa of mice was observed in response to  environmental stress [82]. Furthermore, our findings are supported by increase in expression of IL6 mRNA in mouse macrophages and MEF cells after hyperthermia [83]. It has also been examined that the exposure of heat stress to mice hepatocytes resulted in TNFα induced apoptosis [84]. Along with these findings, effect of heat stress in mammary tissue and peripheral blood mononuclear cells during the dry period in cows revealed altered expression pattern of cytokines exposed to heat stress [85]. Consistent with our study, CXCL2 expression showed an increase in response to high ambient temperatures in bovine [86][87]. Our observations support the idea of induced expression of pro-inflammatory transcripts (IL6, IL8, CXCL2) in porcine intestinal epithelial cells exposed to infections [88]. Abee and Wouters [89] have also cited that stress adaptive genes such as BOLA play a role in controlling cell morphology during heat stress. Our data also showed correlation with increased expression of MHC class genes in porcine intestine [90] and human intestine epithelial cell line [91] during heat stress. In addition, IL1B was found to be most induced during inflammatory response in pigs [92]. Likewise, as reported for heat stressed human skin fibroblasts [93],several genes associated with oxidative stress viz., GSR, DUSP16, GPX7, HMOX1, TXNRD1, GPX4 were specifically up-regulated during the early phase of heat stress response in buffalo MECs. Similar to our results, genes of glutathione peroxidase family were shown to be induced under heat stress condition in Saccharomyces cerevisiae [94] and quail [95]. The family of GPX is known to play an important role in protecting animals and humans against oxidative stress. In addition the elevated expression of HMOX1 and TXNRD1 genes were observed in human melanoma cell culture which confirmed the induction of cellular oxidative stress during harmful insults [96]. Our findings has provided the evidence to suggest the varied expression profile of immune-responsive and oxidative stress related genes in buffalo MECs during heat stress. Thus, in the present study, RT-qPCR analysis validated the transcriptional expression profile of HSPs, apoptotic genes, immune responsive and oxidative stress related genes as observed by microarray gene expression analysis.
Clustering and annotation of early and late transcriptome data. For creating hierarchical clustering, data of differentially expressed genes across all time points was used. This approach has allowed classifying the whole transcriptome data based on variations in gene expression of heat stressed buffalo MEC at different time points. The analysis generated heat maps to judge for the similarities/patterns between genes and between samples. Based on conditions (time points), the analysis revealed 2 major clusters early time points from 30 min to 4 h post heat stress grouped together while late time points covering 8 h to 24 h along with the control (CTR) formed the second cluster. The clustering data reflected the presence of specificity of expression pattern with respect to time points post heat stress. These results provided evidence   16 h and 24 h) at fold change 3 with respect to control. Enrichment score (modified fisher exact p-value) for each gene-set calculated in the Gene Ontology terms using DAVID tool [98]. The p-values were computed by a modified Fisher's exact test for each GO term. We determined the number of genes sharing the same GO terms with a correction for overrepresented (p 0.05) categories based on Gene Ontology data. The top most cluster obtained in early up-regulated time points were generally responsible for regulation of cell cycle (6 genes) with highest EASE score (represents more enrichment of the group), apoptosis (6 genes), chaperon activity (7 genes), transcription factor binding (5 genes). The EASE score for GO terms related to up-and down-regulated genes in early and late time points are given in S2 Table. Whereas the cluster in early down-regulated, identified an enrichment of genes associated with signal transduction (27 genes), oxidation reduction (10 genes), response to stimulus (10 genes). Further, highest scoring GOs in late up-regulated include cell cycle (5 genes), negative regulation of apoptosis (4 genes) and late down-regulated genes were primarily related to immune response (15 genes) and cytoskeleton (10 genes). The gene set level analysis revealed various functional groups of gene and related biological mechanism involved in heat stress. Identification of biological process, molecular functions and pathways affected in buffalo MECs during heat stress. In order to explore the biological significance, detailed annotation of gene function, biological process and cellular distribution of differentially expressed genes (DEG; > = 3 fold change) identified in response to heat stress in vitro in MEC was accomplished by gene ontology (GO) descriptions. Using the data set of all DEG across all time points post heat stress, a total of 32 biological processes were found to be affected across all time points. However, the four GO terms that was most enriched under biological processes included; response to stimulus (638 genes), multicellular organismal process (506 genes), single organism signaling (381) and cellular developmental process (227 genes. The major five molecular functions identified were binding activity (285), molecular transducer activity (156), receptor activity (143), and transporter activity (25). Under molecular transducer activity, signal transducer was the major molecular function. Signal transducer activity basically conveys a signal across a cell to trigger a response in order to change in cell function or state. Under receptor activity, signal receptor activity (129) was found to be major molecular function associated with DEG. Under transport activity, substrate specific transporter activity (25) and transmembrane transporter activity (25) were affected. The three major cellular processes associated with DEG across all time points were extracellular region (440), membrane (367) and cell part (378).
Additionally, REVIGO analysis was performed on DEG (> = 3 fold change) which summarized major GO terms influenced in buffalo MECs under heat stress. REVIGO is a web server that summarizes long, unintelligible lists of GO terms by finding a representative subset of the terms using a simple clustering algorithm that relies on semantic similarity measures. In the present analysis, a higher frequency of several biological process (cell communication, multicellular organismal development, signal transduction & immune response), cellular component (extracellular region, plasma membrane) & molecular functions (protein binding and transporter activity) were obtained after removing redundant GO terms (Table 4). REVIGO analysis corroborated more or less with the GO analysis that was performed using GO module of GenSpring GX software.
Further, 153 genes that were up-regulated in heat stressed MEC across all time points were also assigned with GO terms. A large range of GO categories for biological process were identified including cellular process, metabolic process, single organism process, response to stimulus, biological regulation, immune system processes, signaling etc. Among the GO terms associated with response to stimuli in biological processes, the most significant categories were cellular response to stress, response to chemical stimulus and response to stress (Fig 12). GO terms for  molecular functions were also identified for 153 genes that were commonly over expressed at all time points (Fig 12). They belong to the categories of catalytic activity (18), binding activity (25), transporter activity (3), structural molecular activity (2) and enzyme regulator activity (2). Binding activity was the second major molecular function activity which was enriched. Various sub-categories like ion binding (13), carbohydrate derivative binding (7), protein binding (11), small molecule binding (9) heterocyclic compound binding (14) etc. were affected.
To get better insight into biological function, we further analyzed the differentially expressed genes based on prior knowledge of biological pathways. Several pathways over represented across all time points were; Electron transport chain, Cytochrome P450 pathway, Apoptosis, IL2 signaling, MAPK, FAS and stress induction of HSP regulation, Delta Notch signaling pathway, Apoptosis modulation by HSP70, EGFR1 signaling, Cytokines and Inflammatory response, Nuclear receptors, Oxidative stress, TNF-alpha/ NF-kB signaling pathway and GPCRs pathway. Representative picture of one of the major signaling pathways; Cytokines & Inflammatory response pathways facilitating cell survival and cell death program is shown in Fig 13.

Conclusion
The present work thus presented a suitable strategy to characterize the cellular and transcriptomic adaptation of buffalo mammary epithelial cells to heat stress in-vitro. Use of heterologous bovine Agilent microarray expression chip in the present study proved successful in dissecting the transcriptome of heat stressed and unstressed buffalo MECs. The study thus has identified several heat responsive genes from different functional classes and biological pathways related to chaperons, immune function, cell proliferation and metabolism etc. known to be affected in by heat stress. The present data provided the strong clue about the coordinated transcriptional response of buffalo mammary epithelial cells to heat stress. The responsiveness of buffalo MECs to heat stress in the present study clearly suggested its suitability as a model to understand the modulation of buffalo mammary gland expression signature in response to environmental heat load. In future, such studies could be extended in evaluating the impact of  Table. Candidate target and reference genes evaluated in this study with primer sequences and annealing temperature (Ta). (DOCX) S2 Table. EASE scores for affected GO terms in early and late time points post heat stress. (DOCX)