Evaluation of heat stress effects on cellular and transcriptional adaptation of bovine granulosa cells

Background Heat stress is known to affect follicular dynamics, oocyte maturation, and fertilization by impairing steroidogenic ability and viability of bovine granulosa cell (bGCs). The present study explored the physiological and molecular response of bGCs to different heat stress intensities in-vitro. We exposed the primary bGCs to heat stress (HS) at 39 °C, 40 °C and 41 °C along with control samples (38 °C) for 2 h. To evaluate the impact of heat stress on bGCs, several in vitro cellular parameters including cell apoptosis, intracellular reactive oxygen species (ROS) accumulation and HSP70 kinetics were assessed by flow cytometry, florescence microscopy and western blot, respectively. Furthermore, the ELISA was performed to confirm the 17β-estradiol (E2) and progesterone (P4) levels. In addition, the RNA sequencing (RNA-Seq) method was used to get the molecular based response of bGCs to different heat treatments. Results Our findings revealed that the HS significantly decreased the cell viability, E2 and P4 levels in bGCs, whereas, increased the cellular apoptosis and ROS. Moreover, the RNA-Seq experiments showed that all the treatments (39 °C, 40 °C and 41 °C) significantly regulated many differentially expressed genes (DEGs) i.e. BCL2L1, STAR, CYP11A1, CASP3, SOD2, HSPA13, and MAPK8IP1 and pathways associated with heat stress, apoptosis, steroidogenesis, and oxidative stress. Conclusively, our data demonstrated that the impact of 40 °C treatment was comparatively detrimental for cell viability, apoptosis and ROS accumulation. Notably, a similar trend of gene expression was reported by RT-qPCR for RNA-seq data. Conclusions Our study presented a worthy strategy for the first time to characterize the cellular and transcriptomic adaptation of bGCs to heat stress (39, 40 and 41 °C) in-vitro. The results infer that these genes and pathways reported in present study could be useful candidates/indicators for heat stress research in dairy cattle. Moreover, the established model of bGCs to heat stress in the current study provides an appropriate platform to understand the mechanism of how heat-stressed bGCs can affect the quality of oocytes and developing embryo.


Background
Mammalian ovarian follicle, consisting of an oocyte that undergoes series of biological events including ovulation, fertilization, and formation of an embryo is surrounded by granulosa and theca cells producing signals and hormones to enable the oocyte to develop [1]. During follicle development, granulosa cells (GCs) replicate, secrete hormones, and provide a critical microenvironment for follicular growth [2]. Proliferation and differentiation of GCs is essential for normal follicular growth, development of oocyte, ovulation, and luteinization [3,4].
Heat stress is one of the environmental factors that have harmful effects on the function of ovaries [5] and subsequently decreases the developmental ability of oocytes to be fertilized and further develop the competent embryo [6]. It significantly reduced the production of estradiol, and stenedione synthesis by theca cells [7], inhibited proliferation and induced apoptosis in swine granulosa cells [8]. In support of this, heat stress during in vitro fertilization increased polyspermy and decreased fertilization success by disrupting the antipolyspermy system in oocytes [9], suggesting that heat stress during fertilization mainly affects the oocyte and its developmental competence. Mammalian cells are known to respond to a wide range of environmental stressors in a variety of ways including; heat shock response protein expression [10], unfolded protein response (UPR) [11] and oxidative stress response [12] to support cell survival under suboptimal conditions. Cells may use constitutive induced heat shock proteins (HSPs), molecular chaperones in response to heat stress that facilitate the synthesis, folding, assembly, and transportation of stressdenatured proteins [13]. Heat shock 70 kDa protein (HSP70) is a major stress protein induced in mouse GCs by high temperature [9]. Increasing evidence suggests that heat stress induces intracellular ROS concentration [14], resulting in apoptosis of granulosa cells in the mouse [15]. In addition, ROS may subsequently alter the development of bovine embryos during in-vitro oocyte maturation [16].
RNA sequencing (RNA-Seq) has emerged as an innovative method for both mapping and quantifying transcriptome signatures associated with traits [17]. One of the most biologically relevant applications of RNA-Seq is the comparison of mRNA transcriptome across samples from diseased vs. normal individuals, or other specific experimental conditions [18]. The usage of highthroughput RNA sequencing technology has become a powerful tool and a standard method for the measurement and comparison of gene expression levels in a myriad of species and conditions [19]. Therefore, in our study, we employed RNA-Seq to characterize the complete transcriptome of bGC and facilitate the discovery of differentially expressed genes as well as novel genes and pathways under heat stress.
This study was conducted in Beijing, China. Temperature levels were selected for the experiment to treat the granulosa cells, isolated from the ovaries of cattle that were well adapted to the local environment. For instance, we attempted to select experimental temperature levels that were relevant to the physiological body temperatures of cattle under HS in Beijing. During the summer, we collected the data from many dairy farms in Beijing, showing how environmental temperature-humidity index (THI) can affect rectal body temperature (RT). We found that in summer, body temperature may rise to 41°C (Fig. 1). Therefore, we evaluated the effects of the four temperature levels [38 (control), 39, 40, and 41°C] on the physiological traits and transcriptomic gene expression profile in bGCs.
Moreover, while much is now known about the effects of various factors on normal granulosa cells [14,20], to the best of our understanding, no attempt has been taken so far to propose a molecular mechanism or explore gene interactions and molecular pathways related to heat stress response in bGCs at different heat intensities. We hypothesize that, relative to control, bGCs exposed to heat stress will experience alterations both in physiological traits and expression of key genes and pathways required for normal cellular functions. Therefore, the present study was aimed to explore cellular adaptation, generate global gene expression profile of bovine granulosa cells under normal and heat-stressed state and identify molecular pathways significantly regulated in heat-stressed bGCs.

Collection of bovine ovaries and granulosa cells isolation
Ovaries from dairy cattle were collected from a local abattoir and transported to the laboratory in thermally insulated bottles containing sterile physiological saline with 100 U/mL Penicillin and 0.1 mg/mL Streptomycin, at 28-30°C within 2 h of harvesting. After washing with warm 0.9% NaCl solution three times and rinsing in 70% warm ethanol for 30 s, ovaries were washed thrice with warm Dulbecco's Phosphate-Buffered Saline (DPBS). For isolation of BGCs, small healthy follicles (with a diameter of 2-6 mm) were selected using an 18-gauge sterile needle (B-Braun, Germany) and transferred into 15 mL conical centrifuge tubes (Corning, NY, USA). The follicular fluid containing cumulus-oocyte complexes (COCs) and granulosa cells was filtered using a filter with a diameter of 70 μm leaving COCs on the filter. The filtrate with granulosa cells was carefully transferred into 15 mL conical centrifuge tubes, centrifuged at 1500×g for 5 min. The supernatant of the follicular fluid was discarded by aspiration, and granulosa cells were washed thrice in phosphate-buffered saline (PBS), pH 7.4. The GCs were then resuspended in DMEM/F-12 (Gibco, Life Technologies Inc., Grand Island, NY, USA) supplemented with 1% penicillin-streptomycin and 10% fetal bovine serum (FBS, Gibco, Life Technologies Inc., Grand Island, NY, USA).
After 48 h of preculture, cells were attached to the bottom of the wells with the confluence of more than 80%; the medium was replaced with the fresh medium of the same composition. GCs were then cultured at temperatures control group (38°C), or heat treatment groups (39, 40 and 41°C) for 2 h, and then the cells were cultured at 38°C for 12 h. The cells and culture media were collected for further analysis immediately after the culture. Following heat treatments, cultured GCs were harvested using 0.25% trypsin-EDTA (Sigma-Aldrich Chemie, Taufkirchen, Germany).

Western blot analysis of HSP70
Western blot analysis was used in all samples to determine the expression of inducible HSP70 under heat stress. Granulosa cells from each group were washed thrice with 0.1% PVA/PBS, lysed in RIPA lysis buffer (Beyotime, Shanghai, China) containing protease inhibitors. Total protein concentration was measured with Protein Assay (Bio-Rad, 500-0002) and a spectrophotometer at 595 nm (Beckman, DU 530). Proteins were denatured at 100°C for 10 min, separated by SDS-PAGE (12% acrylamide gel containing 0.1% SDS) and transferred onto a nitrocellulose membrane (BioTraceNT, Pall Corp., Port Washington, NY, USA). Membranes were then blocked with 5% (w/v) skim milk in Trisbuffered saline (TBS) containing 0.1% Tween 20 (TBST) at 37°C for 1 h. The membranes were incubated overnight at 4°C with primary antibodies against HSP70 and β-actin after three washes in TBST. All primary antibodies were purchased from Cell Signaling Technology (Beverly, MA, USA), and diluted to a concentration of 1: 1000. After washing in TBST thrice, the membranes were incubated at room temperature for 1 h with horseradish peroxidase (HRP)-conjugated secondary antibody (Zhongshan Biotechnology, Beijing, China). Based on the manufacturer's instructions, the protein bands were detected using enhanced chemiluminescence (ECL) detection kit (Tanon, Shanghai, China) and analyzed by densitometry using Image J 1.44p software. The final data exported from Image J was analyzed in Microsoft Excel. Western blot in triplicate was performed for all samples.

Determination of estradiol and progesterone by ELISA
All of the culture media was collected from controlled and heat treated groups and then estimated the levels of P 4 and E 2. The concentrations of P 4 and E 2 were determined using an estrogen and progesterone enzymelinked immunosorbent (ELISA) kits (ENZO life sciences, Germany) according to the manufacturer's instruction.

Determination of intracellular ROS production
About 2 × 10 4 granulosa cells were cultured in 96-well plates. After growing to a confluence of more than 80%, GCs were incubated at 38, 39, 40 and 41°C for 2 h. Following incubation, cells were stained with 10 μmol/L H2DCFDA fluorescent probe (6-carboxy-2′,7′-dichlorodihydrofluorescein diacetate) (Invitrogen, Carlsbad, CA, USA) for 30 min at 38°C in the dark. GCs samples were then washed once in 0.1% PVA/DPBS, and the images were captured immediately under a fluorescence microscope (Olympus, Tokyo, Japan) equipped with a Cool-SNAP HQ CCD camera (Photometrics/Roper Scientific, Inc., Tucson, AZ, USA). Image J 1.44p software was used to analyze the fluorescence intensity.

Estimation of granulosa cells apoptosis
Bovine GCs were harvested by enzymatic digestion using trypsin and washed three times with preheated PBS. Using the FITC-Annexin V/dead cell apoptosis kit (Life Technologies Inc., Grand Island, NY, USA), APC annexin V/PI double staining was performed to evaluate the granulosa cell apoptosis according to the manufacturer's instructions before being analyzed by flow cytometry. The data were analyzed by Flowjo software (version Win64-10.4.0).

Estimation of cell viability
Cultured and heat treated GCs were trypsinized, collected, and washed with warm PBS. The GCs were then passing through APC annexin V/PI double staining by utilizing the FITC-Annexin V/dead cell apoptosis kit (Life Technologies Inc., Grand Island, NY, USA) to assess the cell viability and apoptosis. Samples were washed with 1× annexin-binding buffer for 5 min in accordance with the manufacturer's instructions and incubated in 490 μL 1× annexin-binding buffer supplemented with 10 μL annexin V conjugate at room temperature in the dark for 15 min. A laser-scanning confocal microscope (TCS SP8, Leica, Germany) was used to determine the number of early apoptotic and dead cells.

RNA extraction for RNA-Seq
The RNA was isolated from bovine granulosa cells using RNA kit (Tiangen, Beijing, China) according to the manufacturer's instructions. RNA samples were treated with RNase free DNase I to avoid DNA contamination. RNA degradation and contamination were detected by 1% agarose gels. The RNA concentration was assessed using NanoPhotometer spectrophotometer (Implen, CA, USA). The extracted RNA was stored at − 80°C and all 12 samples (three from each group) were sent to the company (Gene Denovo Biotechnology Co. Guangzhou, China) for RNA-Seq analysis.

Library construction for RNA-Seq
Three samples from each group were selected for library preparation. For the RNA sample preparations, A total amount of 2 μg RNA per sample was used as input material. Using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (#E7530L, NEB, USA),sequencing libraries were generated following the manufacturer's recommendations and index codes were added to assign sequences to each sample. Briefly, using Oligo (dT) magnetic beads mRNA was purified from total RNA. Fragmentation was done in NEBNext First-Strand Synthesis Reaction Buffer (5×) using divalent cations at high temperature. First strand cDNA was synthesized using random hexamer primer, and RNase H. DNA polymerase I, RNase H, dNTP, and buffer were used to synthesize second-strand cDNA. Then the fragments cDNA were purified with QiaQuick PCR extraction kit, end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified and sequenced by Gene Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq™ 2500 and generated 150 bp paired-end reads.

Bioinformatics and statistical analysis of RNA-Seq
Raw reads generated by Illumina Hiseq™2500 were initially processed to get clean reads through the following three steps. i) Removing reads with adaptors contamination; ii) Discarding reads containing more than 10% of unknown nucleotides (N); iii) Removing low quality reads containing more than 50% of low quality (Q-value≤20) bases using Next Generation Sequencing (NGS) Quality Control Toolkit version 2.3.3. The filtered reads of each sample was individually mapped to 48370 reference mRNAs from the Bos taurus reference genome (UMD3.1) obtained from Ensembl (ftp: //ftp. Ensembl. org/pub/release-73/fasta/bos_taurus/dna/) HISAT2 software version 2.0.1 (http://ccb.jhu.edu/software/hisat2). The transcripts were then assembled and quantified using StringTie software version 1.2.2 (http://ccb.jhu. edu/software/stringtie). Using StringTie, transcript files generated were added to a single-merged transcriptome annotation to merge transcripts from different replicas of a group into a comprehensive set of transcripts, and then merge the transcripts from multiple groups into a finally a comprehensive set of transcripts for further downstream differential expression analysis. Differentially expressed genes (DEGs) and transcripts were identified among different sample groups using Ballgown. Ballgown was used as a pipeline package in the R programming language version 3.2.2 (https://www.r-project. org) and the Bioconductor software was used for plotting raw data, normalization and downstream statistical modeling. Gene expression values were calculated by counting the number of fragments per kilobase of transcript per million mapped fragments (FPKM), and Cuffdiff was applied to measure significant differences among the four groups. The result was sorted in Microsoft Excel. DEGs were subjected to Gene Ontology (GO) enrichment and Kyoto Encyclopedia of genes and genomes (KEGG) pathway analyses using the Molecule Annotation System http://david.abcc.ncifcrf.gov/) [21]. Using STRING software (version 10), a network was constructed with the genes involved in the significant pathways in order to generate a protein-protein interaction (PPI) network and to predict physical/functional PPIs. The heat map was constructed using ggplot two packages in R (version 3.2).

Quantitative reverse transcription PCR (RT-qPCR) validation for RNA-Seq analysis
RT-qPCR was conducted to confirm the results of RNA-Seq. Total RNA was extracted from three biological replicates of control and heat-treated granulosa cells as described above and was reverse transcribed using first strand cDNA synthesis Kit (Thermo Fisher Scientific, Germany) with oligo (dT) 18 primers according to the manufacturer protocols. The expression levels were checked for 15 genes. Primer3 web version 4.0.0 (http:// bioinfo.ut.ee/primer3/) and Primer blast (http://www. ncbi.nlm.nih.gov/tools/primer-blast/) were used for designing gene-specific primers and are shown in (Additional file 5). The RT-qPCR was conducted using iTaq™ Universal SYBR® Green Supermix (Bio-Rad Laboratories GmbH, Germany) in Applied Biosystem® StepOnePlus™ (Applied biosystems, CA, USA). A reaction volume of 20 μL with 7.4 μL of ddH 2 O, 0.3 μL of forwarding primer, 0.3 μL of reverse primer, 10 μL of 1× SYBR Green master mix (Bio-Rad Laboratories GmbH, Germany), and 2 μL of cDNA template was used. The Light Cycler 480 instrument (Roche, Germany) was used for performing qPCR. The second derivative maximum method was employed for data acquiring and subjected for further analysis. Using GAPDH as a reference gene, the 2 −ΔΔCT method was used to calculate gene expression levels, [21].

Statistical analysis
Data are expressed as mean values ± SEM. Statistical analysis was carried out using SPSS 16.0. The difference between control and heat-treated groups for cell apoptosis, cell viability, steroidogenesis, ROS accumulation as well as the RT-qPCR results were analyzed using Oneway ANOVA followed by multiple comparisons post hoc test. Differences were considered to be statistically significant at P < 0.05.

Heat stress induces HSP70 expression in bovine granulosa cells
Bovine granulosa cells were heat treated at different temperature levels (Control, 39, 40 and 41°C) for 2 h duration to investigate the effect of heat stress on the expression level of HSP70 in bGCs. We performed western blot and RT-qPCR to check relative abundance of HSP70 both on mRNA and protein level. Our results showed that the expression of HSP70 between control and heat stressed group (39°C) was not significantly different. However, the expression of HSP70 was significantly up-regulated in bGCs under heat stress at 40°C and 41°C post-treatment (Fig. 2a, b).

Heat stress exposure elevates bovine granulosa cell apoptosis
The bGCs apoptotic rate was estimated by flow cytometry (FCM) and fluorescence microscopy. It was found that the apoptotic rate (early + late apoptosis) of the GCs was significantly higher (P < 0.05) in the heattreated groups (Fig. 3a, b). During cell culture, bGCs were exposed to heat stress for 2 h with a range of temperatures (Control, 39, 40 and 41°C). Following heat stress exposure, the cell apoptotic rate was increased in temperature-dependent fashion. As shown in Fig. 3a, b, the apoptotic rate (46%) of bGCs was significantly (P < 0.05) higher at 40°C than other treatments. However, with respect to 40°C, the apoptotic rate of GCs was lower at 41°C (35.4%). Heat treatment of 39°C did not change apoptotic rate (9%) significantly than control group(3.96%). Similar effect of heat stress was noticed for cell viability. Significantly (P < 0.05) lower cell viability was found at 40°C (45.3%) as compared to control (96%) and 39°C (82.2%). No significant difference were noticed between 40°C (45.3%) and 41°C (59.4.3%) as shown in Fig. 3a, b. Fluorescence microscopy was also conducted to estimate the apoptotic rate and viability of GCs and found that the relative fluorescence emissions were higher when GCs were exposed to 40°C than the control group. However, the 39°C treatment group did not show any significant (P < 0.05) difference with the control group. Likewise, fluorescence microscopy showed that after 40°C, the apoptotic rate decreased significantly (P < 0.05) in the heat stress group at 41°C (Fig.  3c, d, e, f, g).

Effects of heat stress on E 2 and P 4 secretion by bovine granulosa cells
The concentration of E 2 (Fig. 4a) in the heat-treated groups (40 and 41°C) was significantly lower (P < 0.05) than the control and 39°C group in the culture media. However, heat treated group at 39°C did not show a significant difference with the control group. Furthermore, a significant (P < 0.05) difference was noticed between treatment groups of 40 and 41°C. Similar secretion pattern was also observed for P 4 with the significant difference (P < 0.05) between control and heat-treated groups (40 and 41°C) (Fig. 4b). However, no significant difference was noticed between 40 and 41°C treated groups.

Heat stress enhanced intracellular ROS accumulation in bovine granulosa cells
Following exposure to heat stress at 40°C, an increasing level of intracellular ROS accumulation was observed in granulosa cells compared to other cell cultured groups (Fig. 5c). However, there were no significant differences in ROS accumulation at 39°C (Fig. 5a, b, e). Furthermore, the relative fluorescence emissions were significantly higher (P < 0.05) when bGCs were exposed to 41°C than the control group but lower than 40°C (Fig. 5d).
RNA sequencing data analysis for identifying differentially expressed genes among three groups (control vs. 39°C, control vs. 40°C and control vs. 41°C) In this study, an attempt was made to obtain a global picture of the in-vitro heat stress response by investigating the transcriptome profile of bGCs. Differentially expressed genes (DEGs) of bGCs were identified via RNA-Seq to analyze genome-wide transcriptional expression differences among the three groups. Under the criteria of |fold change| (|FC|) > 1. 5 Fig. 6a). Results revealed the highest number of DEGs in Control vs. 40°C, whereas the least DEG number was detected in control vs. 39°C group. These results indicate a strong inducement of genes in 40°C bGCs cultured group.
Heat stress resulted in the activation of differentially expressed heat shock factors, apoptotic, steroidogenic and oxidative stress-related genes Among the several hundred genes induced or repressed as a result of in-vitro heat stress, an effort was made to filter out genes related to; heat shock protein family, apoptosis; steroidogenesis and oxidative stress ( Table 2). A heatmap and hierarchical clustering of the top 45   Superscripts (a, b, c) show significant difference, P < 0.05 significant (P < 0.05) differentially expressed genes with |FC| > 1.5 and P < 0.05 demonstrates the relatedness of samples, as shown in Fig. 6c.

Pathway analysis of differentially expressed genes in response to heat stress
For an enhanced understanding of signaling pathways regulated by the heat treatments, the identified DEGs in the three comparisons were subjected to pathway analysis using KEGG.

Control vs. 39°C cultured group
A total of 25 canonical pathways enriched by genes differentially expressed in granulosa cells in this comparison (Additional file 2), of which 18 were significantly (P < 0.05) regulated (Fig. 7a, Table 3) while the rest of the eight did not meet the criteria for significance (P < 0.05). The number of up and down-regulated DEGs involved in these 18 KEGG pathways are presented in Fig. 7a. Furthermore, among significantly regulated pathways, 15 were directly related to the qualitative traits of bGCs under heat stress and the genes distributed in each pathway are enlisted in (Additional file 2).

Control vs. 40°C cultured group
Total of 18 canonical pathways was enriched in response to heat stress; out of them, 12 were significantly (P < 0.05) regulated (Additional file 2, Table 3). The 13 pathways which have key roles in the apoptosis, oxidative stress, antioxidant and steroidogenesis regulation of the bGCs were selected and shown in Fig.7b based on up and down-regulated genes. Between the comparisons of Control vs. 39°C and Control vs. 40°C, seven commonly shared pathways were reported. Furthermore, our findings concealed that most of the DEGs among these pathways were up-regulated (Additional file 2). With the increase in heat stress in Control vs. 40°C comparison. Glutathione metabolism pathways were up-regulated to combat stress by regulating antioxidant genes (SOD1, SOD2, etc.) (Fig. 7b).

Control vs. 41°C cultured group
Out of 28 canonical KEGG enriched pathways in control vs. 41°C comparison, 23 reached a significant level (P < 0.05) and are shown based on up and down-regulated genes (Additional file 2, Fig. 7c, Table 3). Furthermore, 14 pathways were involved in the regulation of apoptosis, oxidative stress, antioxidant, and steroidogenesis regulation of the bGCs under heat stress (Additional file 2). Taking all comparisons, it was found that five pathways (Protein processing in the endoplasmic reticulum, FoxO signaling pathway, Apoptosis, p53 signaling pathway, and Pathways in cancer) were shared in all the three comparisons.

Commonly shared genes among all the pathways of the three comparisons
The total of 142, 321 and 294 significantly (P < 0.05) DEGs were documented in the three comparisons of Control vs. 39°C, Control vs. 40°C, and Control vs. 41°C, respectively. Out of these DEGs, 55 genes were commonly shared among the three comparisons. Furthermore, 58, 201, and 179 DEGs were found to be unique genes for Control vs. 39°C, Control vs. 40°C, and Control vs. 41°C, respectively (Additional file 3, Fig.  6b).

Regulation of signaling pathways under heat stress affecting bGCs functions
Heat stress significantly regulated pathways are affecting bGCs physiological attributes i.e. promote the inhibition of cell growth, steroidogenesis, and induction of apoptosis by the accumulation of ROS, etc. These pathways include (MAPK signaling pathway, FoxO signaling pathway, Apoptosis, ovarian steroidogenesis, Protein processing in the endoplasmic reticulum, and Glutathione metabolism. Genes belonging to these canonical pathways were differentially expressed (Fig. 8) in response to HS.

Functional annotation cluster and gene ontology analysis
Detailed annotation of molecular gene function, biological process and cellular distribution of differentially expressed genes (DEGs; > = 1.5 fold change) identified by gene ontology (GO) descriptions in response to heat stress in-vitro cultured bGCs was accomplished in order to explore biological significance.

Control vs. 39°C comparison
A total of 58, 24, and 16 biological processes (BP), cellular components (CC) and molecular functions (MF) respectively, were found to be affected by heat stress.  (Fig. 9a, b, c)     representative genes were performed on the same samples (Additional file 5). Gene expression profiling of granulosa cells showed that some HSP family-related genes were active during cellular heat response ( Table  2). The expressions of HSP family genes such as HSPA13, HMOX1, apoptotic-related genes (CASP3, BAX and BCL2L1), steroidogenic genes (CYP11A1, STAR) antioxidant activity-related genes (SOD2, CAT, GSTA3) and genes related to oxidative stress (FOXO3 and MAP-K8IP1) were significantly (P < 0.05) regulated across all heat-treated granulosa cells compared to the control group. The results showed that all the genes had similar expression trends as detected in the RNA-Seq. This consistency between RT-qPCR and RNA-Seq revealed the reliability of our RNA-Seq data (Additional file 6).

Discussion
Environmental factors, particularly temperature, have a significant impact on animal breeding and reproduction [22]. Heat stress can be defined as a condition that occurs when an animal can not dissipate body heat adequately to maintain thermal equilibrium [23,24]. Heat is proteotoxic stress and causes denatured proteins that, by forming aggregates, can become cytotoxic [25]. The ovarian follicle's granulosa cells play a crucial role in oocyte nourishing, secreting hormones that create functional bidirectional crosstalk with the oocyte [26]. A brief overview of the current study and mechanisms of regulating heat stress response related to follicular function within the bovine ovary is shown in Fig. 10.
In the present study, bGCs were exposed to different levels of in-vitro heat stress and found that heat stress involves compromising the physiological functions of bGCs by increasing intracellular accumulation of ROS, inducing apoptosis and reducing the synthesis of E2 and P4 [7,9,15]. For more understanding, we conducted the transcriptomic study of in-vitro cultured bGCs exposed to heat stress at 39, 40, and 41°C. Amongst the several hundred genes induced or repressed due to heat stress in-vitro, an attempt was taken to screen out genes associated with, heat shock protein family, apoptosis, steroidogenesis and oxidative stress ( Table 2). As expected, the whole set of genes of heat shock family viz., HSPA8, HSPA14, HSP90AA1, HMOX 1, etc. were up-regulated in bGCs at most of the heat stress points. The expression of these genes was more at 41°C of heat stress as compared to other treated groups (39, 40°C). Our findings are supported by previous studies that led to the induction of HSP genes by heat stress [27,28]. Similar to our study, HSPs induction was reported in various cell/ tissue types such as leukocytes/lymphocytes [29][30][31], bovine endometrial tissue, bovine conceptuses [32,33], bovine granulosa cells [34] bovine MECs [22], buffalo lymphocytes [35] due to heat stress. It has been reported that heat stress causes an increase in HSPs in virtually all the vertebrates, including mice [36,37] domestic goats [38], humans [39,40], juvenile Hamadryads baboons [41], common carp [42], domestic chickens [43][44][45][46] and domestic turkey [47]. Our results showed an increased accumulation of inducible HSP70 in heat stressed groups at both the protein and mRNA levels, thus supporting the idea that HSP70 can act as a reliable thermal stress biomarker [42,48]. Likewise, several apoptosis-related genes like BCL2 associated X, apoptosis regulator (BAX), caspase 3, apoptosis-related cysteine peptidase (CASP3) and (CASP6), etc. were also found to be significantly (P < 0.05) up-regulated under heat stress that signal through the apoptosis signaling pathway. Up-regulation of apoptotic genes could lead to disruption of the potential of mitochondrion transmembrane, resulting in the release of cytochrome c leading to apoptosis induction [49]. Data on the induced expression of apoptotic genes at 40°C suggest that the cellular mechanism may not provide protection for bGCs against heat-induced apoptosis while the rate of apoptosis decreased at 41°C of heat stress due to over-expression of HSP70, HSP90 and HSP60 protein levels likely helped bGCs activate self-protection mechanisms and cope with hyperthermia through clearance of damaged proteins. Our results are in line with some previous reports where the MAPK mediated induction of HSP70 at high temperature could play a crucial role in inhibiting caspase-3 and BAX activation [50,51]. Therefore, we suggest that the induction of HSP70 occurs to reduce apoptosis of granulosa cells induced by heat stress. This is a first study that unveiled the effect of heat stress with different intensities on apoptosis-related gene expression and on the cellular defensive mechanism in bGCs. Heat stress results in intracellular ROS accumulation, causing oxidative stress [52] and apoptosis [53], which subsequently lead to a decline in fertility [54,55]. Moreover, the current study also shows for the first time induction of ROS at different intensities of heat stress in bGCs. Compared to the control, significant (P < 0.05) ROS accumulation was evident at 40°C and 41°C of heat stress, but at 39°C the ROS induction was not significant. We found a decline in ROS levels in bGCs by increasing treatment temperature from 40°C to 41°C. This may be due to the cells being able to activate their antioxidant systems at a higher temperature of 41°C, by regulating genes, i.e. superoxide dismutase 1, 2 (SOD1, SOD2,), glutathione-disulfide reductase (GSR) and glutathione S-transferase, alpha 3 (GSTA3) to protect cells against oxidative stress. In addition, the high expression of HMOX1 gene was observed in the culture of human melanoma cells, confirming the induction of cellular oxidative stress during harmful insults [56]. Similar to our results, the activation of forkhead box O3 (FoxO3) and kelch-like ECH associated protein 1 (KEAP1) under heat stress protects cells from oxidative stress by upregulating antioxidant enzymes superoxide dismutase 2 (SOD2) and catalase (CAT) [57][58][59]. In Saccharomyces cerevisiae and quail, genes from the glutathione peroxidase family were also shown to be induced under heat stress [60,61]. Based on these facts, it is reasonable to suggest that the up-regulation expression of SOD2 and CAT may inhibit ROS biosynthesis through the regulation of KEAP1and FOXO3 in ovarian granulosa cells.
Moreover, the regulation of genes related to steroidogenesis, i.e. steroidogenic acute regulatory protein (STAR) and Cytochrome P450, family 11, subfamily A, polypeptide 1 (CYP11A1) was also affected by heat Fig. 8 Regulation of signaling pathways under heat stress affecting bGCs functions: A network map of pathways significantly (P < 0.05) enriched after heat stress. The nodes are the pathways, and edges connect the genes involved in the pathway shock. Previously it was reported that heat stress could inhibit estradiol biosynthesis in bGCs and impair hormone balance [62]. Positive regulation of P450 aromatase family genes (CYP11A1) in the ovarian follicle promotes estrogen biosynthesis [63]. In our study, the mRNA expression of CYP11A1 decreased in GCs by down-regulation of ovarian steroidogenesis signaling pathway after heat treatment that resulted in a decreased level of E 2 in the culture medium. Based on this confirmation, we can postulate that the down-regulation of CYP11A1 may inhibit estrogen biosynthesis in ovarian granulosa cells. In addition, progesterone is also one of the fundamental steroid hormones for bovine estrous cycle regulation, and its biosynthesis is attributed to the increased expression of STAR and CYP11A1 [64][65][66]. Previously it was reported that under heat stress, the mRNA expression of CYP11A1 and STAR decreased, but the P4 level has no significant (P < 0.05) difference between the control and heat treatment group [9]. An over-secretion of ovarian hormones in porcine ovarian granulosa cells were reported under high temperature [67]. Our findings are in line with the previous studies where heat stress attenuates estrogenic activity in rat granulosa cells by diminishing the expression of gonadotropin receptor [68]. We also found a lower mRNA expression of the genes CYP11A1 and STAR in heattreated groups. This is the first study to establish the impact of different intensities of thermal stress on the synthesis of steroid hormones and gene expression profile in bGCs. These findings have provided evidence to suggest the varied expression profile in bGCs during heat stress of apoptotic, steroidogenesis and oxidative stress-related genes. In the current study, the RT-qPCR analysis thus validated the transcriptional expression profile of HSPs, apoptotic genes, steroidogenesis, and oxidative stress-related genes as observed by RNA-Seq analysis. Our research can further be extended to understand bovine oocyte modulation and embryo development in response to the environmental heat load.

Conclusion
In the present study, we demonstrated for the first time a worthy strategy to characterize the cellular and transcriptomic adaptation of bovine granulosa cells to different heat stress intensities (39°C, 40°C and 41°C) in-vitro. Furthermore, our data suggested that 40°C heat treatment is comparatively detrimental for bovine granulosa cell functions. The study identified several heat-responsive genes from different functional classes and their associated pathways related to heat stress chaperons, cell death, and apoptosis, hormonal synthesis, oxidative stress, etc. known to be affected by heat stress. The results infer that these genes and pathways reported in the present study could be useful candidates/indicators for heat stress research in dairy cattle. Moreover, the established model of bGCs to heat stress in the current study provides an appropriate platform to understand the mechanism of how heat-stressed bGCs can affect the quality of oocytes and developing embryo. Fig. 10 Research review: Mechanisms of regulating heat stress response related to follicular function within bovine ovary. Upregulated genes caspase-3, SOD, BCL-2, BAX, and HSPs (HSP70, HSPA13, HMOX1) were involved in the regulating mechanism of bGCs via induced or inhibited cell apoptosis. Under heat stress, Down-regulated genes CAT, FOXO3 were involved in production of reactive oxygen species (ROS). Likewise, downregulation of STAR, and CYP11A1 were involved in the secretion of E 2 and P 4 . Moreover, the decline of E 2 , and enhancing of ROS in turn, might enhance the possibility of GC apoptosis and follicle function