http://orcid.org/0000-0001-8101-3863
Correction
11 Jan 2018: Kapila N, Sharma A, Kishore A, Sodhi M, Tripathi PK, et al. (2018) Correction: Impact of Heat Stress on Cellular and Transcriptional Adaptation of Mammary Epithelial Cells in Riverine Buffalo (Bubalus Bubalis). PLOS ONE 13(1): e0191380. https://doi.org/10.1371/journal.pone.0191380 View correction
Figures
Abstract
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.
Citation: Kapila N, Sharma A, Kishore A, Sodhi M, Tripathi PK, Mohanty AK, et al. (2016) Impact of Heat Stress on Cellular and Transcriptional Adaptation of Mammary Epithelial Cells in Riverine Buffalo (Bubalus Bubalis). PLoS ONE 11(9): e0157237. https://doi.org/10.1371/journal.pone.0157237
Editor: Yiqun G. Shellman, University of Colorado, School of Medicine, UNITED STATES
Received: November 4, 2015; Accepted: May 26, 2016; Published: September 28, 2016
Copyright: © 2016 Kapila et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Data can be found in GEO database with accession number: GSE84486. Additionally, relevant data can be found in the paper and its Supporting Information files.
Funding: This work was supported under National Fellow Scheme of Indian Council of Agricultural Research, New Delhi. Grant number: 27(3)/2010-HRD, http://www.icar.org.in/files/edu/Award-of-National-Fellow.pdf. MM-Principal Scientist and National fellow received the funding from ICAR.
Competing interests: The authors have declared that no competing interests exist.
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–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 sub-continent, 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
Material and Methods
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% CO2 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% CO2and 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 5x103 cells/well in 100 μl of culture medium in 96 well plates (Corning, USA) and cultured for 24–48 h at 37°C, 5% CO2. Cells in control plates were maintained at 37°C, 5% CO2 throughout the time-course, while those in treatment plates were exposed at 42°C, 5% CO2 for 1 h and then shifted to 37°C, 5% CO2. 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–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 75th 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.
Results and Discussion
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–22] and in bovines, established MEC line has also been used to study the effect of heat stress in vitro [12, 23–26].
A) Mixed population of epithelial and fibroblast cells, B) Formation of islands by purified MECs, at low density seeding, C) Cobble stone morphology shown by confluent MECs, D) Post confluent stage of MECs forming dome structures (arrow), E) MECs from passage 5 forming papillate structures, F) Floating dead cells due to occurrence of contact inhibition during post confluency stage.
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 (Fig 2C & 2D) with staining restricted mainly to the peripheral portion of the cytoplasm. It might be attributed to the presence of networks of intermediate filaments, important for cell to cell communication & polarity [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 stressed cells across different time points post heat stress (30 m, 2 h, 4 h, 6 h, 8 h, 12 h, 24 h & 48 h) whereas the unstressed (control) cells maintained at 37°C did not show reduction in cell viability (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–37].
CTR-unstressed (control) MECs.
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].
CTR-unstressed (control); TRT- heat stress treated MECs.
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).
For each graph, the quadrants display distribution of viable (bottom left), early apoptotic (bottom right), apoptotic (upper right) and dead (upper left) cells.
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–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 min, 2 h, 4 h, 8 h,1 2h, 16 h and 24 h), the yield of cRNA was >1.65 μg and the specific activity was >9.0 pmol Cy3 per μgcRNA, indicating the good quality of cRNA that was obtained for all the samples.
Out of 40,000 probes in Agilent bovine genome array, we obtained a total of 24573 which passed the expression filter as per the criteria mentioned. Box whisker plot showing distribution of normalized intensity values is presented in S5 Fig. In order to identify differentially expressed genes (DEG), fold change values were generated by subtracting the intensity of unstressed cells (CTR) from those at 30 min, 2 h, 4 h, 8 h, 12 h, 16 h, 24 h heat-treated cells and selected the genes showing at least 3-fold up- or down- regulation at all-time points of heat treatment. The line plot showing transcriptional pattern of DEG filtered at 3-fold cutoff criteria is depicted in S6 Fig. In comparison to unstressed (CTR), the line plot showed lot of variation in transcript pattern during the early time points (30 min, 2 h, 4 h) post heat stress.
The fold change data waslog2 transformed and selected for further data analysis. At the cutoff criteria of signed fold change ≥2-or ≤2, a total of 19970 differentially expressed transcripts between unstressed and stressed (all studied time points post stress: 30 min, 2h, 4h, 8h, 12h, 16h, 24h) MECs. At the cutoff criteria of signed fold change ≥3-or ≤3, a total of 15286 DEG was observed across different time points post heat stress. The number of DEG at different fold change criteria (2, 3, 5, 10-fold) is presented in Fig 7. Pair wise comparison with CTR at 3 fold change revealed; 2741 DEG (805/1936 induced/repressed) at 30 min, 3137 DEG (720/2417) at 2h, 2645 DEG at 4h (439/2206), 564 DEG at 8h (399/165), 697 DEG at 12h (377/320), 1769 DEG at 16h (759/1010) and 541 DEG at 24h (385/156) (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 (iron-responsive 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–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 that the transcripts were coordinately regulated in a time dependent manner due to heat stress in vitro. Hierarchical clustering of differentially expressed genes is depicted in S7 Fig.
The genes from hierarchical clustering were further analyzed to establish the functional groups of differentially expressed genes at different time points post heat stress (early and late time points). For making a functional link and understanding about the biological themes that are enriched in the gene set, Expression Analysis Systematic Explorer (EASE) tool available in Database for Annotation Visualization and Integrated Discovery (DAVID) (http://david.niaid.nih.gov/david/ease.htm) was used. EASE [97] is generally used as an application tool for rapid biological interpretation of gene lists that could result from the analysis of microarray, proteomics, SAGE and other high-throughput genomic data. We performed functional enrichment analysis against significantly up (151 genes) & down-regulated (628 genes) from first clusters i.e. early time points (30 min, 2 h and 4 h) and a total of 135 up & 5 down-regulated genes from another cluster of late time points (8 h,12 h, 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 hyperthermia and other physiological stressors in tissue/cell damage and related gene regulation studies to understand buffalo mammary functions.
Supporting Information
S1 Fig. Expression profile of HSP genes (A) HSP27, (B) HSP40, (C) HSP60, (D) HSP70 and (E) HSP90 in buffalo MECs in response to heat stress.
https://doi.org/10.1371/journal.pone.0157237.s001
(TIF)
S2 Fig. mRNA abundance of anti-apoptotic (Bcl2) and pro-apoptotic genes (Bax) in heat stressed buffalo MECs.
https://doi.org/10.1371/journal.pone.0157237.s002
(TIF)
S3 Fig. Virtual gel image of buffalo MEC extracted RNA on bio analyzer (Experion).
L: RNA Ladder; CTR: Control; HS: Heat stress
https://doi.org/10.1371/journal.pone.0157237.s003
(TIF)
S4 Fig. Electropherogram profiles of MEC RNA samples indicating 18S and 28S rRNA peaks.
https://doi.org/10.1371/journal.pone.0157237.s004
(TIF)
S5 Fig. Box whisker plot showing differentially expressed genes (DEGs) after 20-100th percentile normalization showing same median at all time points.
https://doi.org/10.1371/journal.pone.0157237.s005
(TIF)
S6 Fig. Line plot of differentially expressed genes at fold change criteria of ≥3-or ≤3>.
https://doi.org/10.1371/journal.pone.0157237.s006
(TIF)
S7 Fig. Hierarchical clustering across 8 time points with DEGs.
The unstressed (CTR) clusters with later stages of heat stress (8 h to 24 h).
https://doi.org/10.1371/journal.pone.0157237.s007
(TIF)
S1 Table. Candidate target and reference genes evaluated in this study with primer sequences and annealing temperature (Ta).
https://doi.org/10.1371/journal.pone.0157237.s008
(DOCX)
S2 Table. EASE scores for affected GO terms in early and late time points post heat stress.
https://doi.org/10.1371/journal.pone.0157237.s009
(DOCX)
Acknowledgments
The authors duly acknowledge the financial support received for this study under National Fellow project of Indian Council of Agricultural Research, New Delhi. The authors also acknowledge Director, ICAR-NBAGR, Karnal for providing research facilities during the course of this study.
Author Contributions
- Conceived and designed the experiments: MM MS PT.
- Performed the experiments: NK AK AS.
- Analyzed the data: MM AS NK.
- Contributed reagents/materials/analysis tools: MM MS NK AKM.
- Wrote the paper: MM AS NK.
References
- 1.
Chase LE. Climate change impacts on dairy cattle, Climate Change and Agriculture: Promoting Practical and Profitable Responses, Department of Animal Sciences 2006; 8.
- 2. Hansen PJ, Arechiga CF. Strategies for managing reproduction in the heat stressed dairy cows. J Dairy Sci. 1999; 77 (Suppl. 2): 36–50.
- 3. Kadzere CT, Murphy MR, Silanikove N, Malt E. Heat stress in lactating dairy cows: a review. Livestock Prod Sci. 2002; 77: 59–91
- 4. Klinedinst PL, Wilhite DA, Hahn GL, Hubbard KG. The potential effects of climate change on summer season dairy cattle milk production and reproduction. Climate Change. 1993; 23: 21–36.
- 5. Sharma AK, Rodriguez LA, Wilcox CJ, Collier RJ, Bachman KC, Martin FG. Interactions of climatic factors affecting milk yield and composition. J Dairy Sci. 1988; 71:819–825. pmid:3372822
- 6. West JW. Effects of heat stress in dairy cattle. J Dairy Sci. 2003; 86: 2131–2144. pmid:12836950
- 7. St-Pierre NR, Cobanov B, Schnitkey G. Economic losses from heat stress by U.S. livestock industries. J Dairy Sci. 2003; 86: E52–77.
- 8.
FAO. Water Buffalo, An asset undervalued. United Nations Food and Agriculture Organization website: http://www.aphca.org/publications/files/w_buffalo.pdf (2000)
- 9. Marai IFM, Habeeb AAM. Buffalo biological functions as affected by heat stress—a review. Livestock Sci. 2010; 127: 89–109.
- 10.
Basu SK. 1985. Genetic improvement in buffaloes. Kalyani Publishers, Ludhiana.
- 11. Liang T, Spence J, Liu L, Strother WN, Chang HW, Ellison JA, et al. α-Synuclein maps to a quantitative trait locus for alcohol preference and is differentially expressed in alcohol-preferring and non-preferring rats. Proc. Natl. Acad. Sci. 2003; 100: 4690–4695. pmid:12665621
- 12. Hu H, Wang J, Bu D, Wei H, Zhou L, Loor JJ. In vitro culture and characterization of a mammary epithelial cell line from Chinese Holstein dairy cow. PLoS ONE 2009; 4(11): e7636. pmid:19888476
- 13. Kapila N, Kishore A, Sodhi M, Sharma A, Kumar P, Mohanty AK, et al. Identification of appropriate reference genes for qRT-PCR analysis of heat stressed mammary epithelial cells in riverine buffaloes (Bubalus bubalis). ISRN Biotech. 2012; Article ID 735053.
- 14. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002; 3(7): 1–12.
- 15. Andersen CL, Jensen JL, Orntoft TF. Normalization of real time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization applied to bladder and colon cancer data sets. Cancer Res. 2004; 64:5245–5250. pmid:15289330
- 16. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. Biotechnol Lett. 2004; 26:509–515. pmid:15127793
- 17. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using realtime quantitative PCR and the 2ΔΔC (T) Method. Methods. 2001; 25:402–408. pmid:11846609
- 18. Zavizion B, Gorewit RC, Politis I. Subcloning the MAC-T bovine mammary epithelial cell line: morphology, growth properties, and cytogenetic analysis of clonal cells. J Dairy Sci. 1995; 78:515–527. pmid:7540186
- 19. Pantschenko AG, Barber MR, Woodcock-Mitchell J, Bushmich SL, Yang TJ. Establishment and characterization of a caprine mammary myoepithelial cell line (CMMyoEC). In Vitro Cell Dev Biol Anim. 2000; 36(6): 351–356. pmid:10949992
- 20. Pantschenko AG, Woodcock-Mitchell J, Bushmich SL, Yang TJ. Establishment and characterization of a caprine mammary epithelial cell line (CMEC). In Vitro Cell Dev Biol Anim. 2000; 36: 26–37. pmid:10691038
- 21. Kumura H, Tanaka A, Abo Y, Yui S, Shimazaki K. et al. Primary culture of porcine mammary epithelial cells as a model system for evaluation of milk protein expression. Biosci. Biotech Biochem. 2001; 65: 2098–2101.
- 22. Anand V, Dogra N, Singh S, Kumar SN, Jena MK, Malakar D. et al. Establishment and Characterization of a Buffalo (Bubalus bubalis) Mammary Epithelial Cell Line. PLoS ONE. 2012; 7:7.
- 23. Collier RJ, Dahl GE, Van Baale MJ. Major advances associated with environmental effects on dairy cattle. J. Dairy Sci. 2006; 89:1244–1253. pmid:16537957
- 24. Du J, Di HS, Guo L, Li ZH, Wang GL. Hyperthermia causes bovine mammary epithelial cell death by a mitochondrial-induced pathway. J Therm Biol. 2008; 33: 37–47.
- 25. Zhao K, Liu HY, Zhou MM, Liu JX. Establishment and characterization of a lactating bovine mammary epithelial cell model for the study of milk synthesis. Cell Biol Inter. 2010; 34: 717–721.
- 26. Liu HY, Zhao K, Zhou MM, Wang C, Ye JA, Liu JX. Cytoprotection of vitamin E on hyperthermia-induced damage in bovine mammary epithelial cells. J Therm Biol. 2010; 35: 250–253.
- 27. Li P, Wilde CJ, Finch LM, Fernig DG, Rudland PS. Identification of cell types in the developing goat mammary gland. Histochem J. 1999; 31(6): 379–393 pmid:10462224
- 28. Ehmann UK, DeVries JT, Chen MS, Adamos AA, Guzman RC. An in vitro model of epithelial cell growth stimulation in the rodent mammary gland. Cell Prolif. 2003; 36(4): 177–190. pmid:12950387
- 29. Gurzov EN, Nabha SM, Yamamoto H, Meng H, Scharovsky OG, Bonfil RD. Paradoxical antiproliferative effect by a murine mammary tumor-derived epithelial cell line. BMC Cancer. 2007; 7:184. pmid:17908302
- 30. Taylor P, Stampfer M, Bartek J, Lewis A, Boshell M, et al. Keratin expression in human mammary epithelial cell cultures from normal and malignant tissue: relation to in vivo phenotypes and influence of medium. J Cell Sci. 1989; 94: 403–413. pmid:2483723
- 31. Hass R, Bertram C. Characterization of human breast cancer epithelial cells (HBCEC) derived from long term cultured biopsies. J Exp Clin Canc Res 2009; 28: 127.
- 32. Schmid E, Schiller DL, Grund C, Stadler J, Franke WW. Tissue type specific expression of intermediate filament proteins in a cultured epithelial cell line from bovine mammary gland. J Cell Biol. 1983; 96: 37–50. pmid:6186672
- 33. Hildebrat U, El Fantroussi S, Smidt H, Smoot JC, Tribou EH, Kelly JJ, Noble PA, Stahl DA. Optimization of single-base-pair mismatch discrimination in oligonucleotide microarrays. App Environ Micro. 2002; 69(5): 2848–2856.
- 34. Park CM, Pudrith CB, Skovgaard K, Ren X, Suchyta SP, Stabel JR, Heegard PMH. John’s disease in cattle is associated with enhanced expression of genes encoding IL-5, GATA-3, tissue inhibitors of matrix metalloproteinases 1and 2, and factors promoting apoptosis in peripheral blood mononuclear cells. Vet Immun Immunopath. 2005; 105:221–234.
- 35. Breen GJ, Mader TL, Holt SM, Josey MJ, Rowan KJ. Heat tolerance of Boran and Tuli crossbred steers. J Anim Sci. 1999; 77: 2398–2405. pmid:10492446
- 36. Du J, He-Shuang DI, Liang GUO, Zhong-Hao LI, Ya-Fei CAI, Wang GL. The effect of high temperature on mammary epithelial cells proliferation and apoptosis. Acta zoological Sinica. 2006; 52: 959–965.
- 37. Liu HY, Yang JY, Wu HH, Wu YM, Liu JX. Effects of methionine and its ratio to lysine on expression of casein αs1 in cultured bovine mammary epithelial cells. J Anim Feed Sci. 2007; 16: 330–334.
- 38. Du J, He-Shuang DI, Wang GL. Establishment of a bovine epithelial mammary cell line and its ultrastructural changes when exposed to heat stress. Sheng Wu Gong Cheng XueBao. 2007; 23: 471–476.
- 39. Collier RJ, Collier JL, Rhoads P, Baumgard LH. Invited Review: Genes involved in the bovine heat stress response.J Dairy Sci. 2008; 91: 445–454. pmid:18218730
- 40. Han Hu, Wang JQ, Li FD, Bu DP, Zhou LY. Responses of cultured bovine mammary epithelial cells to heat stress. J Agri Biotech.2011; 19(2): 287–293.
- 41. Agnew LA, Colditz IG. Development of a method of measuring cellular stress in cattle and sheep. Vet Immuno Immunopathol. 2008; 123: 197–204.
- 42. Fehrenbach E, Passek F, Niess AM, Pohla H, Weinstock C, Dickhuth HH, et al. HSP expression in human leukocytes is modulated by endurance exercise. Med Sci Sports Exerc. 2000; 32(3):592–600. pmid:10731000
- 43. Dangi SS, Gupta M, Maurya D, Yadav VP, Panda RP, Singh G, et al. Expression profile of HSP genes during different seasons in goats (Capra hircus). Trop Anim Health Prod.2012; 44(8): 1905–1912. pmid:22535151
- 44. Shimizu T, Ohshima I, Ozawa M, Takahashi S, Tajima A, Shiota M, et al. Heat stress diminishes gonadotropin receptor expression and enhances susceptibility to apoptosis of rat granulosa cells. Reproduct. 2005; 129(4): 463–472.
- 45. Sheil JM, Shepherd SE, Klimo GF, Paterson YJ. Identification of an autologous insulin B chain peptide as a target antigen for H-2Kb-restricted cytotoxic T lymphocytes. J Exp Med. 1992; 175(2): 545–52. pmid:1370687
- 46. Jorge MS, Patrick F, Miguel V, CecõÂlia MA. The stationary-phase morphogenebol A from Escherichia coli is induced by stress during early stages of growth. Mol Microbiol. 1999; 32(4): 789±798. pmid:10361282
- 47. Adnan M, Morton G, Hadi S. Analysis of rpoS and bolA gene expression under various stress-induced environments in planktonic and biofilm phase using 2−ΔΔCT method Mol Cell Biochem. 2011; 357: 275–282. pmid:21630090
- 48. Dimova DK, Stevaux O, Frolov MV, Dyson NJ. Cell cycle-dependent and cell cycle-independent control of transcription by Drosophila E2F/RB pathway. Genes Dev. 2003; 17: 2308–2320. pmid:12975318
- 49. Marta S, Montserrat C, Hugo S, Tapio H, Ernst H, Florenci S. The Cdi/TESK1 kinase is required for Sevenless signaling and epithelial organization in the Drosophila eye. J Cell Sci. 2006; 119:5047–5056. pmid:17118962
- 50. Minchenko AG, Leshchinsky I, Opentanova I, Sang N, Srinivas V, Armstead VE, et al. Hypoxia-inducible factor-1-mediated expression of the 6-xx phosphofructo-2-kinase/fructose-2, 6-bisphosphatase-3 (PFKFB3) gene. Its possible role in the Warburg effect. J Biol Chem. 2002; 277: 6183–6187. pmid:11744734
- 51. Minchenko OH, Ogura T, Opentanova IL, Minchenko DO, Esumi H. Splice isoform of 6-phosphofructo-2-kinase/fructose-2,6-bisphosphatase-4: Expression and hypoxic regulation. Mol Cell Biochem. (2005), 280: 227–234. pmid:16311927
- 52. Bonnet A, LeCao KA, SanCristobal M, Benne F, Robert-Granie C, Law-So G, et al. In vivo gene expression in granulosa cells during pig terminal follicular development. Reprod. 2008; 136: 211–224.
- 53. Tilman G, Catalgol B, Licht A, Ermak G, Pickering AM, Ngo JK, et al. HSP70 mediates dissociation and reassociation of the 26S proteasome during adaptation to oxidative stress. Free Rad Biol Medi. 2011; 51: 1355–1364
- 54. Patsialou A, Wilsker D, Moran E. DNA-binding properties of ARID family proteins. Nucleic Acids Res. 2005; 33(1):66–80. pmid:15640446
- 55. Masuda K, Kimura A, Hanieh H, Nguyen NT, Nakahama T, Chinen I et al. Aryl hydrocarbon receptor negatively regulates LPS-induced IL-6 production through suppression of histamine production in macrophages. Int Immunol. 2011;23(10):637–645 pmid:21930594
- 56. Moya A, Ganot P, Furla P, Sabourault C. The transcriptomic response to thermal stress is immediate, transient and potentiated by ultraviolet radiation in the sea anemone Anemoniavirdis. Mol Ecol. 2012; 21: 1158–1174. pmid:22288383
- 57. Norgaard JV, Theil PK, Sorensen MT, Sejrsen K. Cellular mechanisms in regulating mammary cell turnover during lactation and dry period in dairy cows. J Diary Sci. 2008; 91: 2319–2327.
- 58. Amaral J, Grant JR, Riggs PK, Nedenia BS, Filho Rodrigues EA, Goldammer T et al. A first generation whole genome RH map of the river buffalo with comparison to domestic cattle. BMC Genomics. 2008; 9:631 pmid:19108729
- 59. Haeseleer F, Sokal I, Verlinde CL, Erdjument-Bromage H, Tempst P, Pronin AN, et al. Five members of a novel Ca (2+) binding protein (CABP) subfamily with similarity to calmodulin. J Biol Chem. 2000; 275, 1247–1260. pmid:10625670
- 60. Deak M, Horvath GV, Davletova S, Torok K, Sass L, Vass I, et al. Plants ectopically expressing the iron-binding protein, ferritin, are tolerant to oxidative damage and pathogens. Nature Biotech. 1999; 17: 192–196.
- 61. Calderón-Villalobos LIA, Nill C, Marrocco K, Kretsch T, Schwechheimer C. The evolutionarily conserved Arabidopsis thaliana F-box protein AtFBP7 is required for efficient translation during temperature stress Luz Irina A. Gene. 2007; 392: 106–116. pmid:17240087
- 62. Desarnaud F, Jakovcevski M, Morellini F, Schachner M. Stress downregulates hippocampal expression of the adhesion molecules NCAM and CHL1 in mice by mechanisms independent of DNA methylation of their promoters. Cell Adhes Migrat. 2008; 2(1): 38–44.
- 63. Guerriero V Jr, Raynes DA. Synthesis of heat stress proteins in lymphocytes from livestock. J Anim Sci. 1990; 68: 2779–2783. pmid:2211407
- 64. Malayer JR, Hansen PJ, Buhi WC. Effect of day of the oestrous cycle, side of the reproductive tract and heat shock on in-vitro protein secretion by bovine endometrium. J Reprod Fertil.1988; 84(2): 567–578. pmid:3199376
- 65. Putney DJ, Malayer JR, Gross TS, Thatcher WW, Hansen PI, Drost M. Heat stress-induced alterations in the synthesis and secretion of proteins and prostaglandins by cultured bovine conceptuses and uterine endometrium. Biol Reprod.1988; 39: 717–728. pmid:3196802
- 66. Mishra A, Hooda OK, Singh G, Meur SK. Influence of induced heat stress on HSP70 in buffalo lymphocytes. J Anim Physiol Anim Nutr. 2011; 95(4):540–544.
- 67. Beck SC, Paidas CN, Tan H, Yang J, De Maio A. Depressed expression of the inducible form of HSP 70 (HSP 72) in brain and heart after in vivo heat shock. Am J Physiol. 1995; 269: R608–613. pmid:7573563
- 68. Albers R, van der Pijl A, Bol M, Seinen W, Pieters R. Stress proteins (HSP) and chemical-induced autoimmunity. Toxicol Appl Pharmacol.1996; 140: 70–76. pmid:8806871
- 69. Meza-Herreraab CA, Martíneza L, Aréchigac C, Bañuelosc R, Rincónc RM, Urrutiab J. et al. Circannual identification and quantification of constitutive heat shock proteins (HSP 70) in Goats. J App Anim Res. 2005; 29(1): 9–12.
- 70. Hayashi Y, Iwai T, Toshio K, Tatsuya K, Kenzo O. Translocation of hsp-70 and Protein Synthesis during Continuous Heating at Mild Temperatures in HeLa Cells. Radiat Res. 1991; 125:1. pmid:1986394
- 71. Kim D, Virginia W, Somji S, Garrett SH, Sens MA, Shukla D et al. Expression of hsp 27, hsp 60, hsc 70, and hsp 70 by immortalized human proximal tubule cells (hk-2) following exposure to heat shock, sodium arsenite, or cadmium chloride. J Toxicol Environ Health A. 2001; 63(7): 475–493. pmid:11497330
- 72. Dehbi M, Baturcam E, Eldali A, Ahmed M, Kwaasi A, Chishti MA et al. Hsp-72, a candidate prognostic indicator of heatstroke. Cell Stress Chaperones. 2010; 15:593–603. pmid:20174993
- 73. Ferencz A, Juhasz R, Butnariu M, Deer AK, Varga IS, Nemcsok J. Expression analysis of heat shock genes in the skin, spleen and blood of common carp (Cyprinuscarpio) after cadmium exposure and hypothermia. Acta Biologica Hungarica. 2012; 63(1):15–25. pmid:22453797
- 74. Givisiez PE, Ferro JA, Ferro MI, Kronka SN, Decuypere E,Macari M.Hepatic concentration of heat shock protein 70 kD (Hsp70) in broilers subjected to different thermal treatments. Br Poult Sci.1999; 40: 292–296. pmid:10465398
- 75. Hernandes R, Ferro JA,Gonzales E, Macari M, Bernal FEM, Ferro MIT. Resistance to ascites syndrome, homoeothermic competence and levels of Hsp70 in the heart and lung of broilers. revistabrasileira de zootecnia-brazilian. J Anim Sci. 2002; 31:1442–1450.
- 76. Zulkifli I, Omar AR, Sazili AQ, Rajion MA. Crating and heat stress influence blood parameters and heat shock protein 70 expression in broiler chickens showing short or long tonic immobility reactions. Anim Welfare. 2003; 1–6.
- 77. Taylor P, Zulkifli I, Norma MTC, Israf DA, Omar AR. The effect of early-age food restriction on heat shock protein 70 response in heat-stressed female broiler chickens. Br Poult Sci. 2010; 43.
- 78. Wang S, Edens FW. Stress-induced heat-shock protein synthesis in peripheral leukocytes of turkeys, meleagrisgallopavo. Comp Biochem Physiol. 1993; 106(3): 621–628.
- 79. Lewis S, Handy RD, Cordi UB, Billinghurst Z, Depledge MH. Stress proteins (HSP’s): Methods of Detection and Their Use as an Environmental Biomarker. Ecotoxicology. 2000; 8(5): 351–368.
- 80. Yadav A, Singh KP, Singh MK, Saini N, Palta P, Manik RS et al. Effect of physiologically relevant heat shock on development, apoptosis and expression of some genes in buffalo (Bubalus bubalis) embryos produced in vitro. Reprod Domest Anim.2013; 48(5): 858–65. pmid:23581430
- 81. Thanida S, Nawapen P, Catherine N, Sukanya M, Mongkol T, Threerawat T. In Vitro Culture of Feline Embroys Increases Stress induced Heat Shock Protein 70 and Apoptotic related Genes. J Reprod Dev. 2013; 59(2).
- 82. Wang Q, Sun X, Pritts TA, Wong HR, Hasselgren PO. Induction of the stress response increases interleukin-6 production in the intestinal mucosa of endotoxaemic mice. Clin Sci. 2000; 99:489–496 pmid:11099391
- 83. Takii R, Inouye S, Fujimoto M, Nakamura T, Shinkawa T, Prakasam R et al. Heat Shock Transcription Factor 1 Inhibits Expression of IL-6 through Activating Transcription Factor 3.J Immunol. 2010; 184(2): 1041–1048. pmid:20018623
- 84. Imao M, Masahito N, Hisataka M. Dual effects of heat stress on tumor necrosis factor-a-induced hepatocyte apoptosis in mice. Lab Invest. 2006; 86: 959–967. pmid:16832353
- 85. Tao S, Connor EE, Bubolz JW, Thompson IM, do Amaral BC, Hayen MJ et al. Short communication: Effect of heat stress during the dry period on gene expression in mammary tissue and peripheral blood mononuclear cells. J Dairy Sci. 2013; 96(1): 378–383. pmid:23141830
- 86. Louisa RA, Payton RR, Gondro C, Saxton AM, Nagle KA, Jenkins BW et al. Heat stress effects on the cumulus cells surrounding the bovine oocyte during maturation: altered matrix metallopeptidase 9 and progesterone production. Reproduction. 2013; 146(2):193–207. pmid:23744615
- 87. Jensen K, Günther J, Talbot R, Petzl W, Zerbe H, Schuberth HJ et al. Escherichia coli- and Staphylococcus aureus-induced mastitis differentiallymodulate transcriptional responses in neighbouring uninfected bovine mammary gland quarters BMC Genomics. 2013; 14: 36. pmid:23324411
- 88. Zanello G, Berri M, Dupont J, Sizaret PY, D'Inca R, et al. Saccharomyces cerevisiae Modulates Immune Gene Expressions and Inhibits ETEC-Mediated ERK1/2 and p38 Signaling Pathways in Intestinal Epithelial Cells. PLoS ONE. 2011; 6(4): e18573. pmid:21483702
- 89. Abee T, Wouters JA. Microbial stress response in minimal processing. Int J Food Microbiol. 1999; 50: 65–91. pmid:10488845
- 90. Yu J, Yin P, Liu F, Cheng G, Guo K, Lu A et al. Effect of heat stress on the porcine small intestine: A morphological and gene expression study. Comp Biochem Physiol. 2010; 156: 119–128.
- 91. Groh V, Steinle A, Bauer S, Spies T. Recognition of stress-induced MHC molecules by intestinal epithelial gamma delta T cells. Science. 1998; 279(5357):1737–40. pmid:9497295
- 92. Hulst M, Smits M, Vastenhouw S, de Wit A, Niewold T, van dM J. Transcription networks responsible for early regulation of Salmonella-induced inflammation in the jejunum of pigs. J Inflamm. 2013; 10: 18.
- 93. Wang Z, Cao N, Nantajit D, Fan M, Liu Y, Li J. Mitogen-activated Protein Kinase Phosphatase-1 Represses c-Jun NH2-terminal Kinase-mediated Apoptosis via NF-{kappa}β Regulation. J Biol Chem. 2008; 283(30): 21011–21023. pmid:18508759
- 94. Fernando WGD, Ramarathnam R, Krishnamoorthy AS, Savchuk S. Identification and use of bacterial organic volatiles in biological control of Sclerotiniasclerotiorum. Soil Biol Biochem. 2005; 37: 955–964.
- 95. Vesco DAP, Gasparino E. Production of reactive oxygen species, gene expression, and enzymatic activity in quail subjected to acute heat stress. J Anim Sci. 2013; 91:582–587. pmid:23148249
- 96. Cabello CM, Bair WB, Lamore SD, Ley S, Bause AS, Azimian S, et al. The Cinnamon-derived Michael Acceptor Cinnamic Aldehyde Impairs Melanoma Cell Proliferation, Invasiveness, and Tumor Growth. Free Rad Biol Med. 2009; 46(2): 220–231. pmid:19000754
- 97. Hosack DA, Dennis G, Brad TS, Clifford HL, Lempicki RA. Identifying biological themes within lists of genes with EASE. Genome Biol. 2003; 4:R70. pmid:14519205
- 98. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucl Acid Resol. 2009; 37(1):1–13