Transcriptome Reveals Granulosa Cells Coping through Redox, Inflammatory and Metabolic Mechanisms under Acute Heat Stress

Heat stress affects granulosa cells (GCs) and the ovarian follicular microenvironment, causing poor oocyte developmental competence and fertility. This study aimed to investigate the physical responses and global transcriptomic changes in bovine GCs to acute heat stress (43 °C for 2 h) in vitro. Heat-stressed GCs exhibited transient proliferation senescence and resumed proliferation at 48 h post-stress, while post-stress immediate culture-media change had a relatively positive effect on proliferation resumption. Increased accumulation of reactive oxygen species and apoptosis was observed in the heat-stress group. In spite of the upregulation of inflammatory (CYCS, TLR2, TLR4, IL6, etc.), pro-apoptotic (BAD, BAX, TNFSF9, MAP3K7, TNFRSF6B, FADD, TRADD, RIPK3, etc.) and caspase executioner genes (CASP3, CASP8, CASP9), antioxidants and anti-apoptotic genes (HMOX1, NOS2, CAT, SOD, BCL2L1, GPX4, etc.) were also upregulated in heat-stressed GCs. Progesterone and estrogen hormones, along with steroidogenic gene expression, declined significantly, in spite of the upregulation of genes involved in cholesterol synthesis. Out of 12,385 differentially expressed genes (DEGs), 330 significant DEGs (75 upregulated, 225 downregulated) were subjected to KEGG functional pathway annotation, gene ontology enrichment, STRING network analyses and manual querying of DEGs for meaningful molecular mechanisms. High inflammatory response was found to be responsible for oxidative-stress-mediated apoptosis of GCs and nodes towards the involvement of the NF-κB pathway and repression of the Nrf2 pathway. Downregulation of MDM4, TP53, PIDD1, PARP3, MAPK14 and MYC, and upregulation of STK26, STK33, TGFB2, CDKN1A and CDKN2A, at the interface of the MAPK and p53 signaling pathway, can be attributed to transient cellular senescence and apoptosis in GCs. The background working of the AMPK pathway through upregulation of AKT1, AMPK, SIRT1, PYGM, SLC2A4 and SERBP1 genes, and downregulation of PPARGCIA, IGF2, PPARA, SLC27A3, SLC16A3, TSC1/2, KCNJ2, KCNJ16, etc., evidence the repression of cellular transcriptional activity and energetic homeostasis modifications in response to heat stress. This study presents detailed responses of acute-heat-stressed GCs at physical, transcriptional and pathway levels and presents interesting insights into future studies regarding GC adaptation and their interaction with oocytes and the reproductive system at the ovarian level.

Granulosa cells are the integral component of the ovaries, where they tend to nourish the developing oocyte inside the follicle and regulate the endocrine signals of the reproduction at the ovarian axis [20,21]. Heat stress is shown to markedly decrease the survival rate of bovine granulosa cells (GCs), increase oxidative stress, limit cell proliferation and transition processes, inhibit steroidogenesis and promote GC apoptosis [22][23][24]. Impairment of GCs through various external or internal stress stimuli may compromise the ovarian functions and development competence of oocytes [25,26]. In a previous study, higher ambient temperatures and associated heat stress were shown to be associated with a decline in ovarian reserves and symptoms of reproductive aging in women [27]. In particular, the pre-ovulation stage of the ovarian follicle remains highly susceptible to the adverse effects of heat stress causing oxidative stress and damage to the developing oocytes [28][29][30]. Heat stress causes energetic metabolism alterations in the body [31] and resulting high nonesterified fatty acids (NEFAs), ketones and inflammatory cytokines alter the biochemical profile of the ovarian follicles [18,30,[32][33][34]. These biochemical and inflammatory changes cause dysfunction of GCs and lead to developmental competence of the oocytes [30,35,36].
Transcriptomic-level studies involving global RNA-sequencing may give essential biological insights into the cellular mechanisms, metabolic level changes and apoptotic and antioxidant pathways [23]. Earlier studies have reported a senescence in proliferation and reported variable ability of camelid GCs to recover their proliferation potential following acute and chronic heat stress exposure [35,36]. In our previous study about GCs [23], there was a disparity among different temperature ranges where 41 • C presented less lethality than 40 • C. Therefore, it was imperative to evaluate the physiological and transcriptomelevel changes in acute heat-stressed GCs. This study will help further our understanding of basic molecular phenomena underlying heat-stress responses in the mammalian body and help us to explore intervention avenues for welfare and general and reproductive health management. Henceforth, in this study we selected 43 • C at the rate of 2 h [37] to inflict acute heat stress to GCs. To the best of our knowledge, there is a lack of systematic transcriptomic research regarding granulosa cell culture in response to acute heat stress (43 • C). We hypothesized that GCs under acute heat stress would exhibit physical insults and undergo extensive changes in biological pathways funneled towards the cellular homeostasis mechanisms.

Granulosa Cell Culture, Heat Treatment and Cell Proliferation Assay
Cell collection, culture and treatment methods followed our previous study [38]. In brief, follicular contents were collected from small follicles (3-8 mm) of healthy cyclic bovine ovaries through an 18-gauge syringe and sieved through a 40 µm filter (Corning Inc., Corning, NY, USA). Follicular fluid was centrifuged, supernatant discarded, and GC pellets were washed twice. GCs were seeded on sterile culture plates (Corning Inc., Corning, NY, USA). Culture media consisting of DMEM/F12 medium containing 10% FBS (both from Thermo Fisher Scientific, Waltham, MA, USA) and 1% antibiotic were cultured under 38 • C and 5% CO 2 in an incubator. GCs were cultured for 48 h and afterwards subjected to 2 h of acute heat stress (43 • C), while the control group remained at 38 • C.
Cells (2 × 10 4 in each well) were cultured in 96-well plates and optical density (OD) was measured using a Cell Counting Kit-8 (CCK-8; Dojindo Laboratories, Kumamoto, Japan) according to the manufacturer's instructions. In brief, 10 µL of CCK-8 was added to each well and incubated for 2 h. Post-treatment absorbance of cells was measured at a 450 nm wavelength by a plate reader, followed by measurements at the intervals of 6,12,24,48,72,96 and 144 h. For the immediate CCK-8 assay (0 h after heat stress), CCK-8 solution was added just before the start of heat stress treatment.

Reactive Oxygen Species, Apoptosis and Hormone Measurements
Reactive oxygen species, apoptosis and hormone measurement methods followed our previous study [38]. Cells were cultured in 96-well black plates and fluorescence OD was measured using a 6-carboxy-2 ,7 -dichloro-dihydro-fluorescein diacetate kit (DCFDA kit, Abcam, Cambridge, MA, USA) according to the manufacturer's instructions. Briefly, cells (1 × 10 4 in each well) were initially cultured, the medium was removed and wells were washed once with 100 µL/well of 1 × buffer (supplied in kit) followed by the addition of 100 µL/well of the 20 µM DCFDA solution. Cells were incubated for 40 min at 38 • C in the dark. After that, DCFDA solution was removed and cells were supplemented with 100 µL/well complete culture medium and incubated for 8 h (heat stress group, 2 h at 43 • C followed by 6 h at 38 • C). After incubation, fluorescence OD was immediately measured on an Infinite M200 PRO (Tecan Deutschland GmbH, Crailsheim, BW, Germany) fluorescence plate reader at excitation/emission = 485/535 nm wavelengths.
The qualitative apoptosis rate of GCs was measured by a fluorescence microscope with an Annexin V-FITC kit (Nanjing Jiancheng Bio Inst., Nanjing, China). After postheat-stress recovery at 38 • C, cells were harvested by 0.25% trypsin and washed twice with warm PBS and then added to 500 µL 1 × binding buffer. Subsequently, 5 µL FITC and PI were added to each replicate. After 20 min incubation at room temperature, slides were prepared and observed under a fluorescence microscope (Axio Imager A2; ZEISS Microscopy, Oberkochen, BW, Germany). Fluorescence microscope pictures at 40× of early apoptotic (green) and late apoptotic (red) cells were recorded. Four replicate slides of each were counted for distinct red and green features.
After 6 h of post-heat-stress recovery at 38 • C, the culture media from both groups were collected and immediately stored at −80 • C for the determination of P4 and E2 concentration through a commercial ELISA kit (Cusabio Technology LLC, Wuhan, Hubei, China) according to manufacturer's protocols. The absorbance was measured by an Infinite M200 PRO (Tecan Deutschland GmbH, Crailsheim, BW, Germany) plate reader at a 450 nm wavelength.

Statistical Analysis of Physiological Parameters
The data from at least six replicates each for cell proliferation, ROS and apoptosis measurements were analyzed using Graphpad Prism version 9.0.0 for windows (GraphPad Software, Inc., San Diego, CA, USA). Analysis of variance was carried out and means were compared through Tukey's honestly significant difference (HSD) test at a 5% level of significance (α = 0.05). All the data presented in the figures are expressed as mean ± standard error (S.E).

RNA Sequencing
The RNA was isolated from GCs according to the manufacturer's instructions of TRIzol Reagent method [39]. RNA quality and concentration were assessed using an Equalbit RNA BR Assay Kit (Invitrogen, Carlsbad, CA, USA) and Nanodrop 2000 (Thermo, Waltham, MA, USA). RNA integrity was detected by 1% agarose gel electrophoresis, and it could be used for library construction with 28S/18S > 1. For the RNA-seq library, a total of 2 µg RNA was used for purification and fragmentation using a NEBNext Poly(A) mRNA Magnetic Isolation Module (Cat No. E7490S, New England Biolabs (UK) Ltd., Hitchin, Herts, UK), followed by the cDNA library with the NEBNext Ultra RNA Library Prep Kit for Illumina (Cat No. E7530S, New England Biolabs (UK) Ltd., Hitchin, Herts, UK). Libraries were quantified by Equalbit RNA BR Assay Kit (Invitrogen, CA, USA) and pooled equimolarly, and finally submitted for sequencing by the NovaSeq 6000 System (Illumina, Inc., San Diego, CA, USA), which generated 150 base paired-end reads.

Differential Expression, Validation and Functional Analyses of Differentially Expressed Genes
The quality of the sequencing reads was evaluated by FastQC software (v0.11.9), and global trimming was carried using Fastp (v0.20.0) [40]. All clean reads were mapped to the bovine genome of version ARS-UCD1.2 using the software package STAR (v2.7.5c), and the Picard query (http://broadinstitute.github.io/picard/, accessed on 20 October 2021) was used to mark duplicates. Related statistical summary analysis of RNA-seq data sets of treated and control cells is given in Table 1. We investigated the counts of gene expression through RNA-SeQC software (v2.3.6) [40]. Principal component analysis and clustering structure were performed using the psych and hcluster R packages. For differential expression gene (DEG) screening, the quantile-adjusted conditional maximum likelihood (qCML) was performed using edgeR in the R package [41] with criteria fold change ≥ 1.5 and 0.05 for the alpha of false discovery rate (FDR; 3 samples in treatment group versus 3 samples in control group). For RNA-seq data validation, cDNA libraries were constructed using a PrimeScript TM RT reagent kit with gDNA eraser (Perfect Real Time, Takara). Quantitative real-time PCR (RT-PCR) was carried out according to the manufacturer's protocol of the Sybr Premix EX Taq Kit (Takara, Kyoto, China) using applied Biosystem ® StepOnePlus™ (Applied biosystems, CA, USA). Bovine β-Actin and B2M (beta-2-microglobulin) were used as internal control genes. Primers for nine randomly selected significant and non-significant genes were designed using Primer-BLAST (www.ncbi.nlm.nih.gov/tools/primer-blast/ index.cgi?LINK_LOC=BlastHome, accessed on 9 February 2022) and double-checked by Oligo v7 (Molecular Biology Insights, Inc., Cascade, CO, USA) [42]. All the primers were synthesized by TsingKe Biological Co., Ltd. (Beijing, China) and are presented in Supplementary Table S1. RNA expression levels relative to the control genes were calculated as 2-Ct according to previous research [43]. To compare the results of the RNA-seq analysis and RT-qPCR, fold change values were log2 transformed.
The DEGs were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and enrichment analysis of Gene Ontology (GO). Online DAVID software (https: //david.ncifcrf.gov/, accessed on 9 February 2022) was used for functional annotation DEGs of the GO terms including biological process (BP), cellular component (CC) and molecular function (MF) enrichment analysis. STRING version 11.5 (https://cn.stringdb.org/, accessed on 8 February 2022) was used for the protein-protein interaction (PPI) network analysis of DEGs.

Influence of Heat Stress on Physical Parameters of Bovine Granulosa Cells
Cultured GCs were exposed to heat stress treatment (43 • C) in vitro, while the control group remained at 38 • C. Cells in the control group maintained steady proliferation activity ( Figure 1A). No change in cell viability was observed until 24 h after heat stress exposure ( Figure 1B), while the group with medium change at 0 h after heat stress ( Figure 1C) did show a certain degree of proliferation activity. Both control and heat stress groups achieved confluences at 72 h after treatment. GCs were subjected to post-treatment exposure recovery for 6 h at 38 • C, after which GC ROS and apoptosis were estimated. Compared to the control group, a significant (p < 0.05) increase in ROS level was observed in the heat stress treatment group ( Figure 2D). Heat-stressed GCs had a significantly higher (p < 0.05) apoptotic rate compared to the control group ( Figure 2C), as shown in representative fluorescent microphotographs of early apoptotic and late apoptotic events in the control group ( Figure 2A) and heat stress ( Figure 2B), respectively. Similarly, P4 and E2 concentration was significantly (p < 0.05) decreased in the culture media of the treatment group ( Figure 2E,F, respectively). Comparison of physical parameters of bovine granulosa cells (GCs) exposed to heat stress (43 • C for 2 h) versus control (38 • C). Representative fluorescence microscope pictures (200×) of late apoptotic (red) and early apoptotic (green) cells after staining with FITC/PI dye for control (A) and heat stress treatment (B) groups. Means comparison of apoptotic rate (sum of red and green events) of GCs under control and heat stress (C). Reactive oxygen species (ROS) through fluorescence OD value (measured at 485/535 nm wavelength) of GCs stained with 2 ,7 -dichlorofluorescine diacetate (DCFDA) is shown on the y-axis, and the temperature treatments are indicated on the x-axis (D). Progesterone (P4) and estrogen (E2) concentration change comparison among control and heat stress groups (E,F), respectively. All data are represented as means ± S.E. of at least three independent cultures with at least three additional replicates for each culture. Means without common letters (a, b) are significantly different (p < 0.05).

Role of Acute Heat Stress on GCs through Deep Sequencing
RNA-seq read counts of heat stress and control groups (3 replicates each) were distinctly different from each other as determined through principal component analysis and dendrogram clustering structure (Supplementary Figure S1). Whole-genome expression profiling with RNA-sequencing analysis revealed that heat stress had a significant effect on the gene expression in GCs. Out of 12,385 differentially expressed genes (DEGs) listed in Table S2, a total of 330 genes were identified as significant DEGs with a fold change of ≥1.5 and FDR < 0.05 (Table S3). Among these DEGs, 75 were upregulated and 255 were downregulated genes ( Figure 3A,B). The log-transformed relative expression fold change of nine randomly selected DEGs in the heat stress and control groups generated from RT-qPCR were in line with the results of RNA-seq data ( Figure 3C). The top 20 upregulated and downregulated significant DEGs are listed below in Table 2. The majority of DEGs are involved in cellular stress response, structural remodeling and cellular reprogramming.  Among the several hundred genes up-or downregulated as a result of in vitro acute heat stress (Supplementary Table S2), an effort was made to filter out genes related to the heat shock protein family, apoptosis, antioxidant support, energetic support (Figure 4), and partial clues of major signaling pathways on the basis of extensive manual querying of these pathways in related literature ( Figure 5).

Pathway Analysis of Differentially Expressed Genes
Detailed results of the KEGG-based pathway analysis of the DEGs are given in Supplementary Table S4. The illustration of all canonical pathways (p < 0.05) is given in Figure 6. Additionally, Supplementary Table S4 was checked for the top enriched genes in each pathway, and these are listed along with enrichment ratios (in all pathways) in Table 3.    Based on Figures 4 and 5, and a manual search of given pathways in the stress-associated literature and molecular signaling mechanisms given in those studies, Supplementary Table S2 was checked for gene regulation and flowcharts were drawn. Flow charts consisting of various integral cellular pathways and associated genes involved in initiation, intermediate signaling, direct or indirect triggering, target and outcome, and inhibitory mechanisms are given (Figure 7).

Functional Annotation of Differentially Expressed Genes
Further, DEGs were checked for gene ontology (GO) term enrichment. All significant (p > 0.05) GO term details given as biological process (BP), cellular component (CC) and molecular function (MF) are presented in Supplementary Table S5. The top enriched genes in each of the functions of BP, CC and MF terms along with their enrichment ratios are listed in Table 4. Genes TEK, VEGFA, ACTA2 and INSL3 have high enrichment ratios in more than one GO term classification, while ACTA2 is the only upregulated DEG in the given DEGs. Table 4. Detailed annotation of gene ontology terms (GO) of significant differentially expressed genes (DEGs) in control versus acute heat stress treatment groups was carried out. Top 10 DEGs enriched in various terms of biological process (BP) and top 5 enriched DEGs in cellular component (CC) and molecular function (MF) are detailed along with their term enrichment ratio and statistical detailing.

Protein-Protein Interaction (PPI) Networks of DEGs under Acute Heat Stress
Function protein association network analysis of proteins coding DEGs comprising various k-means clustering groups (STRING) with medium confidence score (0.4) are presented in Figure 8. The total number of nodes was 278, total number of edges (interactions) was 221, average local clustering coefficient was 0.324 and PPI enrichment p-value stood at 4.93 × 10 −7 .
Additional full STRING network analysis with the highest confidence score (0.9) was also carried out to narrow down the important interaction nodes (Figure 9), while details are given in Supplementary Table S6. TP53, RPTOR, FKBP1B, GNB5 and AKT1S1 were the important functional nodes with the highest score. SERPIND1, MASP2, AMY2B, PYGM, MAOB, ISG15, SEMA6A and PLXNA2 can be regarded as important nodes based on their centrality and edges.

Discussion
Heat stress causes various adverse modifications at physiological and metabolic levels in cows [15,16]. However, the occurrence of atypical negative energy balance in postpartum cows during heat stress is very important, which is characterized by blunt non-esterified fatty acid (NEFA) response [15,31] (keep it like this and remove 47, no need of 47 here, I cannot change these (please do it by yourself), otherwise whole file references and citations will be disturbed) [15,31]. This phenomenon, when accompanied with low body condition scores, clinical ketosis, decreased feed intake, fatty liver and other morbidities, leads to impairment of the reproduction process [44][45][46]. Summer-time heat stress in the presence of negative energy balance and abnormal postpartum conditions [44][45][46] can cause changes in the blood metabolic profile extending effects towards ovarian follicles, developing oocytes and embryos [47,48]. Primary and pre-antral ovarian follicles are particularly susceptible to the adverse effects of heat stress [18,28]; the compromise of GCs in the follicles has a negative impact on oocytes [25,49], leading to a decline in their developmental competence [50]. Hence, we investigated the physical, biochemical and transcriptome-level changes in acute-heat-stressed GCs.
There was a significant decline in cell proliferation potential (cell senescence) in the heat stress group compared to the control group. Previous studies support this finding, where acute-heat-stressed (45 • C) GCs re-attained proliferation potential 48 h posttreatment [35][36][37]. A 24 h delay in cell confluence achievement (at 72 h) observed in our study can be attributed to different experimental methods and the species difference in these previous studies, as camels are thought to be relatively thermo-tolerant [35]. However, the trend in proliferation activity recovery of our GCs follows the results reported by Alemu et al. and Fu et al. [22,37], where they reported post-heat-stress recovery of GC proliferation activity at 48 h. GCs without heat stress treatment proliferated steadily, while the heat stress group with immediate culture media change showed relatively earlier positive proliferation change compared to the heat stress group with culture medium change at 48 h. This phenomenon is interesting and need further investigation; a preliminary hypothesis may be drawn that heat stress is an energy-intensive process, as supported by our previous review of the literature studies [16,51]. Energy balance has a direct relationship with follicular growth and development, ovulation and fertility in postpartum dairy cows [52]. Therefore, careful monitoring of feeding with enough potential to maintain a positive energy balance appears to be the most important intervention avenue in dairy cows exposed to heat stress [15]. NEBAL is the most basic cause for poor reproductive efficiency [13][14][15][16], emphasizing the importance of maintaining positive energy balance manifolds in the case of additional heat stressors [18]. Going back to cell-specific responses to acute heat stress, downregulation of TP53, MDM4, MYC, PARP3, MAPK8 and RAS2 and upregulation of MAP3K7, CDKN2A, SMAD6 and TGFB2 is important. These genes present strong evidence of the involvement of the p53 signaling pathway working in conjunction with other major pathways affecting acute-heat-stressed GC cell proliferation and causing high apoptosis through complex signaling mechanisms [53,54]. Together with other pathways, the p53 signaling pathway can be regarded as one of the major contributing pathways towards the initial high rate of apoptosis and later resumption of cell proliferation at 48 h in heat-stressed GCs [55,56].
Acute heat stress in vitro triggered accumulation of intracellular ROS and caused an increase in apoptosis of GCs in our study; these results are in accordance with the previous studies [23,57]. Previous studies show that post-heat-stress intracellular ROS accumulation increases in a time-dependent manner [22], while our study shows a substantial intracellular ROS accumulation within 10 h of heat stress initiation. ROS production is a normal phenomenon occurring in the cells due to fundamental biological processes [58]. Higher ROS accumulation in the presence of stress stimuli and mitochondrial damage [22] impairs cell antioxidant response elements (AREs) and causes oxidative stress [23,59]. Besides the intracellular high ROS in acute-heat-stressed GCs in our study, other evidence like upregulated sets of BAX, BAK and BID, and HMOX1, NOXA1, NQO1, BCL2L11, SOD, CAT and GPX4, strongly evidence the presence of oxidative stress and reciprocal cellular mechanisms of a higher antioxidant response against it. Previous studies involving cellular stress associates intercellular ROS with oxidative stress [60] and apoptosis in cells [61,62]. In our study, this high rate of apoptosis was typically characterized by upregulation of FADD, TRADD, MAP3K7, RIPK3 and caspases. ROS-mediated oxidative stress induces bovine GC apoptosis and compromises fertility in cows [63,64]. However, pathways and mechanisms underlying heat stress causing extreme oxidative stress and blunt ARE response causing apoptotic fates of the cells are not completely understood.
Heat-stressed GCs decreased the production of P4 and E2 hormones, as demonstrated previously [22,23,55]. Lower levels of these hormones in heat stress cause poor estrus behavior through E2 and low luteal support and fertility through P4 levels [18]. This decline in steroidogenic activity of GCs can essentially be attributed to the downregulation of STAR, CYP11A1 [23] and mitochondrial damage due to ROS [22]. Gene CYP11A1 was also downregulated in our study; however, STAR was mildly (Log 2 (FC) = 0.051) upregulated.
This phenomenon is interesting as many of its orthologues like STARD13, STARD8, STARD3 and STARD4 were downregulated in the heat stress group. We further searched other genes related to P4 synthesis and found that SRD5A1, HSD17B11 and SREBF1 were also downregulated; therefore, we propose a complex transcription mechanism underlying the decreased P4 secretion under heat stress. This claim may be further augmented by our finding about the majority of genes involved in cholesterol synthesis and homeostasis (such as HMGCS1, HMGCR, ACTA1, ACTA2, ABCA1, SLC25A1, SAOT1, AMPK, SRD5A3, DHDDS NS FDPS) being upregulated in heat stress. A very similar dilemma was observed with the specific genes involved in E2 synthesis: where CYP19A1 was upregulated, CYP17A1 was not detected at all, and HSD3B1 and HSD17B1 were downregulated. We could not identify deeper and mechanistic causes underlying these phenomena; however, our results present a further line of investigation. The established phenomenon of blunt NEFA response in heat-stressed cows is characterized by low fat mobilization, as explained by previous studies [15,31,47]; upregulated genes involved in fatty acids and cholesterol metabolism indicate that fatty acid metabolism is an adaptation of heat-stressed GCs. It can be inferred that fat mobilization is disturbed at organism level in heat stress, but cellular metabolism of fats remains enhanced. Therefore, supplemental fats can help to avert the heat-stressed typical NEBAL, as supported by previous studies [65].
As shown earlier in this study, the majority of genes primarily associated with cell damage (CYCS), inflammation (TLR2, TLR4, IL6, TNFRSF6B, TNFSF9), apoptosis initiation (BAD, BAX, BAD, BAK1, BID) and execution (FADD, TRADD and caspases) are upregulated in heat stress, while conversely a large proportion of antioxidant and anti-apoptotic genes (HMOX1, NOXA1, NQO1, SIRIT1, BCL2L11, SOD, CAT, GPX4) are also upregulated. These variable and complicated phenomena indicate evidence of intricate "biological fight" in heat-stressed GCs, and essentially explain the observation of post-heat-stress transient senescence of GCs found in this and a previous study [35].
In the event of cellular stress, the Nrf2 pathway is involved in complex upregulatory functions related to cell metabolism in conjunction with NF-κB (nuclear factor kappa light chain enhancer of activated B cells) [66,67], PI3K/AKT/mTOR [68][69][70] and AMPK signaling pathways [71]. Evidence of high ROS, upregulation of inflammatory cytokines, ERK1/2, AKT1, KEAP1 and NRF1 indicates induction of the Nrf2 pathway, which is augmented by antioxidant and anti-apoptotic genes. Genes GSK3A, GSK3B and serine threonine kinases like STK26 and STK33 also evidence the activity of the Nrf2 pathway. The upregulation of CDKN2A is also interesting, and augments the repression of transcription activity by the Nrf2 pathway in heat-stressed GCs [72]. Gene SLC2A4 encoding GLUT4 was upregulated in heat stress; however, we propose that the Nrf2 pathway cannot be solely implicated in the GC bioenergetics-based protective responses against thermally driven oxidative stress. The complex role of the transforming growth factor-β (TGF-β) pathway is also important in this context, and is involved in signaling mechanisms related to upregulation of NF-κB and inducing apoptosis through cell cycle arrest [73,74]. The most important evidence of cell cycle arrest in heat-stressed GCs can be explained by our finding of the downregulation of MYC, TP53 and MDM4 and upregulated CDKN2A and CDKN1A, which may lead to apoptosis through p53 signaling pathway involvement [53,54]. Similarly, there was evidence of non-canonical induction of the TGF-β signaling pathway as indicated by downregulation of SMAD6 and upregulation of RHOJ, MAPK1 and AKT1 in this study [75]. Regarding the MAPK signaling, upregulation of MAP3K7 (TAK1) and MAPK14 (p38) [54] and downregulation of MAPK8 are of particular importance [76], as evidenced in this study, and may contribute to the apoptosis of heat-stressed GCs. In the event of cellular stress, the Nrf2 and NF-κB pathways are involved in complex regulatory functions [66,67] related to cell metabolism in conjunction with AMPK signaling pathways [71]. Regarding the AMPK signaling pathway in this study, it is characterized by downregulation of ATP gene families and PPARGC1A, and upregulation of SIRIT1, GLUT4, HMGCS1 and SERBP1 [71,77].
GCs play an integral role in ovaries through their steroidogenic roles and oocyte development role [20,21]. Heat stress and its adverse effects on GCs have been implicated as one of the causes of fertility decline in mammals [18,23,27]. In the previous section of discussion, we elaborated the effects of acute heat stress on the physical parameters of GCs. At the same time, we tried to elaborate brief insights about the pattern of major biological pathways and their regulation. Out of 12,385 true RNA reads from the global transcriptomic comparison of control and acute-heat-stressed GCs, 330 significant DEGs were determined through post-analyses. The majority of DEGs were downregulated (225), indicating adverse effects of acute heat stress on the transcriptome of GCs. Going through the list of the top 20 up-and downregulated DEGs, RBM3 appeared as the most significant upregulated DEG; this gene, being a cancer biomarker, is also shown to promote DNA integrity and promote cellular homeostasis through complex signaling [78]. , calcium transport and membrane-channel-related genes (EFCAB2, SMOC1, CACNB2), zinc finger proteins and F-box-only protein families of genes explains a major part of the upregulated DEGs. Similarly, the majority of the top 20 downregulated DEGs relates to cell-structural disturbance, lack of cell proliferation potential and genes related to the energetic homeostasis of cells. The heat-shock protein family response in the presence of upregulated HSF1 was also relatively blunted in response to acute heat stress, with only a few members such as HSPA12A, HSPA9, HSPB6 and HSPBP1 being upregulated genes. The details of pro-apoptotic, anti-apoptotic and antioxidant genes have been discussed earlier, where there is evidence of strong ongoing antioxidant response in GCs after acute heat stress. KEGG functional pathway analysis shows an enrichment of pathways, such as VEGF and TGF-β signaling, showing the disturbance of cellular transcription and proliferation. Various amino-acid-metabolism-associated pathway enrichments show the evidence of protein breakdown and the decline in cellular energetics homeostasis. The remaining KEGG-enriched pathway nodes suggest intercellular structure breakdown and disturbance. These observations are augmented by the downregulated status of the entire list of top DEGs enriched in all KEGG pathways. The majority of GO annotation enrichments were observed in the BP enrichment terms, where, except ACTA2, all DEGs with the highest enrichment ratios were downregulated in response to acute heat stress. Genes ACAT2, TLR2 and SLC2A4 are the only three upregulated DEGs in the top 20 different GO term enrichment processes. Earlier studies involving acute heat stress in bovine skeletal muscles [68] and GCs [38] reported altered metabolism and insulin signaling, lack of anabolic activities and extensive reshuffling in the bioenergetics mechanisms involving multiple signaling pathways. Heat stress and associated oxidative stress are shown to lead to upregulated inflammatory response [79,80], while inflammation is energy intensive and leads to metabolic imbalance [81,82]. Similarly, in a previous study involving acute and chronic heat-stressed GCs, there were evident changes in the cell morphology, proliferation potential, and wound healing potential [35], where cells lost their tone and became rounded, and evidence of granules and cytoplasmic fragmentations were visible. Therefore, among these DEGs, TLR2 can be taken as a representative of inflammatory response and cytokine signaling upregulation, SLC2A4 as a representative of cellular energetic support mechanisms and ACTA2 as the representative of the resort of maintaining the cellular structural integrity in response to heat stress.

Conclusions
Upon the infliction of in vitro acute heat stress, GCs suffer from decreased cell viability, high oxidative stress and increase in the apoptosis rate. The steroidogenic potential of GCs declines, in spite of the evidence of normal cellular cholesterol metabolism. A large number of antioxidant and anti-apoptotic genes remain upregulated, in spite of strong evidence of upregulated pro-apoptotic and apoptotic genes. There is extensive interplay among all major cell signaling pathways, where the evidence suggests repression of cell transcriptional and proliferation activity, and the cellular "resort" of averting the effects of heat stress through remodeling of cellular structural proteins and energetic homeostasis. The evident role of the inflammation-metabolism nexus appears to be central in the acuteheat-stressed mediated GCs' ultimate fate, and hence this study indirectly emphasizes this intervention nexus for reproductive management under heat stress. Our RNA-seq approach gives novel insights into the biological mechanisms and transcriptional activity of GCs exposed to acute heat stress and opens gates of further research in these type of cell-stress stimuli interactions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/cells11091443/s1, Figure S1: RNA-seq read counts of heat stress and control groups (3 replicates each) are distinctly different from each other as shown in principal component analysis and dendrogram clustering structure. Table S1: List of primers used for RT-qPCR validation of genes found in RNA sequencing. Table S2: List of all differentially expressed gene reads along with their statistics parameters, found in control versus heat-stressed granulosa cells. Table S3: List of all significant (p < 0.05) differentially expressed genes (DEGs) in response to heat stress in granulosa cells. Table S4: Kyoto encyclopedia of genes and genomes (KEGG)-based pathway analysis of differentially expressed genes (DEGs) in response to heat-stressed granulosa cells. All the significantly (p < 0.05) regulated pathways and DEG enrichment ratios are illustrated. Table S5: Detailed annotation of gene ontology (GO) terms (p > 0.05) of significant differentially expressed genes in control versus acute heat-stress treatment groups. GO terms including biological process (BP), cellular component (CC), and molecular function (MF) enrichment ratios are given. Table S6: List containing differentially expressed genes and their statistical parameters involved in PPI interaction network analysis at the confidence score of 0.9.