Introduction

High temperatures (HT) can lead to seasonal infertility in livestock, which is particularly evident in the pig breeding industry. HT manifests as reduced ovulation, conception, and pregnancy maintenance rates1,2. In addition, HT induces an ovarian heat shock protein response and increases apoptotic signal transduction in granulosa cells3,4. Therefore, understanding the effects of HT on porcine reproduction is an increasingly important area of concern.

Mitochondria are one of the main organelles affected by HT, and their integrity and homeostasis are disrupted under HT conditions. Research has shown that HT experienced by oocytes during maturation can persist during embryonic development, alter nuclear structure, damage the cytoskeleton, cause mitochondrial swelling, and induce organelle damage5. Mitochondrial membrane potential, an important marker for detecting the developmental ability of preimplantation embryos, is reduced in HT embryos6. HT also affects embryonic development by damaging mitochondrial DNA transcription, replication, and expression of genes encoding adenosine 5’ triphosphate (ATP) to produce oxidative phosphorylation complexes7. In addition, the mitochondrial electron transport chain (ETC) is sensitive to heat stress signals. However, research on this topic in embryos is limited. In somatic cells, the core subunits of mitochondrial complex I, NDUFS8 and NDUFV2, are downregulated under HT conditions to regulate the oxidative stress response8. HT oxidizes thiol groups in mitochondrial complexes I, II, IV, and V, thereby reducing the mitochondrial ETC activity9. These findings indicate that mitochondria play an important role in early embryonic response to HT.

Life begins with a transition of genetic control from the maternal to the embryonic genome during zygotic genome activation (ZGA). The genome of the fertilized zygote is initially in transcriptional inertia, but the zygotic gene begins to rapidly transcribe and initiates ZGA after several divisions in early mammalian embryos10,11. In many species, ZGA is not a single biological event, but rather a process of gradual transcriptional activation that consists of two main phases: minor and major ZGA12. The onset and duration of ZGA vary considerably among species; however, within the same species, the process is tightly controlled. In mice, ZGA occurs in the 2-cell (2 C) embryos13, whereas in humans and pigs, it occurs during the 4 C period14,15. ZGA is the driving force for early embryonic regulation of its own development, which is essential for the establishment of embryonic developmental potential and pluripotency. Aberrant ZGA activation, influenced by genetic and environmental factors, leads to defects in preimplantation embryonic development16,17,18,19. During the embryonic stage, from oocytes to ZGA, the embryos are susceptible to HT; this leads to a decrease in developmental ability. In addition, late-stage embryos, such as morulae or blastocysts, acquire HT resistance. The mechanism of temperature-dependent developmental changes is believed to involve the accumulation of antioxidants in embryos in response to HT-induced production of reactive oxygen species20. However, the differences in the RNA transcriptome expression profiles between normal and HT conditions during the ZGA stage of early porcine embryonic development remains unclear.

In this study, we selected porcine 4 C embryos as the research object and studied the effects of HT on the transcriptome characteristics of ZGA in early porcine embryos at the transcriptome level. RNA-seq was used to screen for differentially expressed genes (DEGs) using gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG)21,22,23, gene set enrichment analysis (GSEA), and other analysis methods. Our results indicate that HT substantially affects pathways related to genetic information integration, oxidative phosphorylation, mitochondrial function, and the cell cycle in early porcine embryos, and has different mRNA expression profiles for ZGA. Among them, NDUFA6 and CDKN1A served as hub genes for DGEs.

Materials and methods

Unless otherwise specified, all the chemicals were purchased from Millipore (Burlington, Massachusetts, USA). All experimental operations were conducted on a constant temperature heating table set at 38.5℃ to ensure the stability and accuracy of the experimental environment.

Ethics statement

This study was approved and performed in accordance with the guidelines of the Institutional Animal Care and Use Committee (IACUC) of Chungbuk National University Laboratory Animal Center, Cheongju, Korea. All procedures did not involve live animal experiments in our work and all porcine ovaries were purchased from local slaughterhouses. All animal studies were approved and performed according to committee guidelines.

Porcine oocyte collection and in vitro maturation

Fresh ovaries were collected from local slaughterhouses and placed in 37 ℃ constant temperature saline and sent to the laboratory immediately. Aspirate follicles with a diameter of 3–6 mm and select high-quality cumulus-oocyte complexes (COCs). After thrice washing with operating solution, the COCs were transferred to in vitro maturation (IVM) medium and covered with mineral oil. IVM refers to the addition of 10 ng/mL epidermal growth factor, 0.5 µg/mL luteinizing hormone, 0.1 g/L sodium pyruvate, 0.6 mmol/L L-cysteine, 10% (v/v) porcine follicular fluid, and 10 IU/mL follicle-stimulating hormone. Place approximately 30–50 COCs were added to each well and were cultured at 38.5 °C in a humidified 5% CO2 incubator for 44 h. After IVM, COCs were denuded of cumulus cells by continuous pipetting about 30–50 times in a 1 mg/mL hyaluronidase solution.

Parthenogenetic activation and in vitro culture

The denuded oocytes were put in 297 mmol/L mannitol containing 0.05 mM MgSO4, 0.1 mM CaCl2, 0.5 mM HEPES, and 0.01% (w/V) polyvinyl alcohol (PVA), and the pH was adjusted to 7.2. Under these conditions, oocytes were treated with two direct-current pulses of 120 V for 60 µs. To suppress pseudo extrusion of the second polar body, 7.5 µg/mL of cytochalasin B and 4 mg/mL of bovine serum albumin (BSA) were added to bicarbonate-buffered porcine zygote medium 5 (PZM-5), and oocytes were cultured in this medium for 3 h. Add 4 mg/mL BSA to PZM5 to prepare in vitro culture (IVC) medium and wash the activated oocytes thrice. Divide the oocytes into control and HT groups, then culture them in IVC medium at 38.5 °C and 41.5 °C to 4 C stage, respectively.

RNA-seq analysis

Following the manufacturer’s protocol, about 300 samples (4 C) in the lysis buffer were used to construct sequencing libraries with the SMART-Seq® v4 Ultra® Low Input RNA Kit for Illumina. This protocol was composed of first-strand cDNA synthesis, amplify full-length ds cDNA, end repair, dA tailing, ligation of the indexing adapters, and PCR enrichment in sequence to create a final cDNA library. Libraries were quantified by qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using an Agilent Technologies 4200 TapeStation (Agilent Technologies, Waldbronn, Germany). Paired-end sequencing reads were generated using an Illumina NovaSeq sequencing platform (Illumina, San Diego, CA, USA). Filtering out low quality sequences and joint sequences, the Clean Reads were analyzed. Gene- and transcript-level quantification were calculated as the raw read count, fragments per kilobase of transcript per million mapped reads (FPKM), and transcript per million mapped reads (TPM). RNA library sequencing was performed on the Illumina HiseqTM 2500/4000 by Gene Denovo Biotechnology Co., Ltd (Guangzhou, China).

Real-time quantitative PCR

We collected 30 4 C embryos from each group and subsequently extracted mRNA according to the instructions of the Dynabeads mRNA Direct Purification Kit (61012, Thermo Fisher Scientific, Waltham, USA). cDNA was obtained by reverse transcription of mRNA using oligo (dT) 20 primers and SuperScript III reverse transcriptase (Thermo Fisher Scientific). Real-time quantitative PCR (qPCR) was performed using the WizPure qPCR Master Kit and the amplification was performed at 95 °C for 3 min, followed by 40 cycles at 95 °C for 15 s, 60 °C for 25 s, 72 °C for 10 s, and a final extension at 72 °C for 5 min. The mRNA quantification data were analyzed using the 2−ΔΔCt method. Sequences of all primers are provided in Table 1.

Table 1 Summary of PCR primers.
Full size table

Statistical analysis

Each experiment was repeated at least three times and statistical comparisons were performed using independent sample t-tests. Statistical significance was set at p < 0.05. The gene expression level of the control group was arbitrarily set to 1 and the HT group was expressed relative to the control group and the results are reported as mean ± SEM. Statistical analysis was performed using GraphPad Prism software (version 5.0; GraphPad Inc., San Diego, CA, USA).

Results

Global transcriptome analysis of control and HT oocytes

To investigate the differences in gene expression in response to HT, high-throughput RNA-seq analysis was performed on control and HT 4 C stage embryos. The principal component analysis (PCA) showed good intra-group reproducibility, and the two groups were clearly distinguished from each other, as shown in Fig. 1A. By comparing transcript abundances, 326 differentially expressed genes (DEGs) were identified, of which 141 were up-regulated and 185 were down-regulated, as shown in Fig. 1B and C. To verify the accuracy of the RNA-Seq data, we randomly selected 11 genes for qRT-PCR. All genes showed expression patterns consistent with the RNA-Seq data, verifying the reliability and accuracy of the sequencing analysis, as shown in Fig. 1D.

Fig. 1
figure 1

Global transcriptome analysis of control and HT oocytes. (A) PCA showing the reproducibility between repetitive samples and the difference between groups. (B) Heatmap of mRNA expression profiles in Control and HT embryos, showing changes in a subset of genes in response to HT. The color key (from blue to red) of Z-score value (− 2, − 2) indicates low to high expression levels. (C) DEGs between control and HT embryos were shown in the volcano map. Genes that expressed higher (up-regulated) in HT embryos are shown in red, and genes that were lower expressed (down-regulated) are shown in yellow. (D) Validation of RNA-seq data by qPCR. *p < 0.05; ns means not significantly.

Full size image

HT interferes with genetic information in porcine embryos

To better understand the roles of the DEGs during the HT process, the DEGs were enriched into Gene Ontology (GO) terminology. The enrichment difference bubble chart showed that most genes were enriched in the GO terms contained in the cellular component category, whereas the terms contained in the molecular function category were the most significantly enriched, as shown in Fig. 2A. The top 15 enriched GO terms classified by –log10(Q value) are shown in Fig. 2B. The top enriched terms in molecular function were DNA polymerase activity and nucleotidyltransferase activity, whereas the most enriched terms in biological processes were reverse transcription and DNA replication. In addition, we found that after HT exposure, the expression of GMNC, FBXO4 and RPA3 related to DNA replication, decreased, while the expression of TOP1MT and POLM increased, which together with downregulated POLE2 are related to DNA polymerase activity, as shown in Fig. 2C.

Fig. 2
figure 2

HT interferes with genetic information in porcine embryos. (A) The enrichment difference bubble chart of GO enrichment analysis on DEGs. The abscissa of the image represents the normalized ratio of up-regulated genes to down-regulated genes. Greater than 0 means that more genes in this pathway are up-regulated. (B) The top 15 enriched items analyzed by GO. (C) The heat map for gene expression levels associated with DNA replication and DNA polymerase activity in enriched GO terms.

Full size image

HT interferes with ZGA in porcine embryos

In pigs, ZGA occurs during the 4 C phase, during which the embryo undergoes intracellular maternal mRNA and protein degradation and transcriptional initiation of its own genome. HT led to the enrichment of DEGs in these terms, indicating that porcine embryonic ZGA was disturbed. Therefore, we compared the expression of the porcine ZGA gene24 in the control and HT groups and found 39 ZGA genes among the DEG, as shown in Fig. 3A. Thereafter, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed to explore the pathways affected by these ZGA genes. As shown in Fig. 3B, these ZGA genes are significantly enriched in TGF-beta signaling, HIF-1 signaling, and selenocompound metabolism and cell cycle pathways. Furthermore, GO analysis of these ZGA genes revealed significant enrichment of terms related to cycling binding, as shown in Fig. 3C. The expression levels of these ZGA genes in control and HT embryos are shown in Fig. 3D. A total of 23 ZGA genes, including CDKN1A, BTF3, NDUFA4, and CDC23 were down-regulated in the HT group embryos.

Fig. 3
figure 3

HT interferes with ZGA in porcine embryos. (A) Venn diagram shows that 39 DEGs are overlapped with porcine ZGA genes. (B) GO enrichment analysis on the ZGA genes in DEGs showing the top 5 enriched KEGG pathways. (C) GO enrichment analysis on the ZGA genes in DEGs showing the top 5 enriched GO items. (D) The heat map for gene expression levels associated with ZGA.

Full size image

HT interferes with the energy metabolism of porcine embryos

Subsequently, we explored the impact of HT-induced DEGs on signaling pathways in embryos using KEGG pathway enrichment analysis. The top ten enriched pathways are shown in Fig. 4A. We found that oxidative phosphorylation was a significantly enriched pathway, with a lower expression tendency in the HT group. Metabolic reprogramming of the main energy supply pathway, from glycolysis to oxidative phosphorylation, is another key event in the development of porcine preimplantation embryos. As shown in Fig. 4B, the expression levels of the oxidative phosphorylation genes, ATP5MC3, ATP5MC3, NDUFA6, etc. were lower in the HT group.

Fig. 4
figure 4

HT interferes with the energy metabolism of porcine embryos. (A) KEGG enrichment analysis on DEGs showing the top 10 enriched KEGG pathways. (B) The heat map for gene expression levels associated with oxidative phosphorylation in enriched KEGG pathways. (C) GSEA results showing the suppression of the “oxidative phosphorylation”, “inner mitochondrial membrane protein complex”, “ATP synthesis coupled proton transport”, and “mitochondrial protein complex” signaling pathway in HT embryos. (D) The heat map for gene expression levels associated with mitochondrial function.

Full size image

Gene Set Enrichment Analysis (GSEA) captured several gene sets with differential expression during transcriptome changes. GSEA analysis indicated significant down-regulation of the “oxidative phosphorylation”, “inner mitochondrial membrane protein complex”, “ATP synthesis coupled proton transport”, and “mitochondrial protein complex” signaling pathways associated with energy metabolism, as shown in Fig. 4C.

Mitochondria are the main sites of oxidative phosphorylation, and a decrease in oxidative phosphorylation may be due to mitochondrial dysfunction. Based on the GO terms, we identified genes related to mitochondrial function and found that their expression patterns matched this result, as shown in Fig. 4D. Most Genes related to the mitochondria, mitochondrial proton-transporting ATP synthase complex, mitochondrial respiratory chain complex, and ATP binding were down-regulated in the HT group embryos.

NDUFA6 and CDKN1A are the hub genes of the protein interaction network

After determining that HT mainly affected the pathways and genes related to DNA replication and energy metabolism in porcine oocytes, we plotted the regulatory network between the KEGG pathways. The results are shown in Fig. 5A. The cell cycle and oxidative phosphorylation pathways are central to the regulatory network. This further illustrates the important roles of the cell cycle and oxidative phosphorylation pathways in HT-induced porcine embryonic injury. Thereafter, we plotted a protein-protein interaction network diagram of the differentially expressed genes, as shown in Fig. 5B. We found that NDUFA6 and CDKN1A were at the center of the two interaction networks. Interestingly, NDUFA6 and CDKN1A were involved in the oxidative phosphorylation and cell cycle pathways, respectively. Therefore, we believe that NDUFA6 and CDKN1A are hub genes among the differentially expressed genes associated with HT.

Fig. 5
figure 5

NDUFA6andCDKN1Aare the hub genes of the protein interaction network. (A) The string diagram shows the cell cycle and oxidative phosphorylation pathways are in the center of pathway relationship network. (B) The string diagram shows the NDUFA6 and CDKN1A genes are in the center of protein-protein interaction networks.

Full size image

Discussion

HT are considered one of the most challenging issues in livestock management. The reproductive system is more sensitive to temperature than the other parts of the body and the impact of HT on fertility is persistent. Although the animals are no longer exposed to HT conditions, the low fertility rate present from warm months of the year (June to September) can still be detected in autumn (October to November)25. Studies have described the effects of HT on the reproductive function of domestic animals during gamete and early embryonic development. For example, our previous research showed that HT leads to meiotic division failure of oocytes and a decrease in the proportion of embryos developing into blastocysts26,27. The pre-implantation embryos in the early stages are highly sensitive to HT damage, and exposure to temperatures only 1.5–2.0 °C higher than normal body temperature can affect the development of mouse and bovine 1-cell embryos28,29. Furthermore, some studies have exposed embryos to temperature fluctuations for 24 h to mimic the rectal temperature of cows exposed to heat stress. It was found that embryonic development was impaired after 8 days in this environment, but the effect was not significant after only one day of exposure, indicating that embryonic development may be interrupted by short-term severe or long-term mild heat shock30. However, the molecular mechanisms underlying this process need to be elucidated and improved upon. Our data indicated that there are multiple differences in RNA expression during early porcine embryonic development under HT conditions.

First, we compared and analyzed dynamic changes in global gene expression in porcine embryos under normal and HT conditions. Principal component analysis showed good intra group reproducibility between the control and HT group, indicating that the RNA-Seq data were of high quality. Subsequently, we identified 326 DEGs between the two groups and validated their expression, further demonstrating the reliability of the RNA-seq analysis. To elucidate the biological functions of HT during early porcine embryonic development at the molecular level, the DEGs were subjected to functional enrichment analyses using GO. These results indicated that these terms were mainly related to the integration of genetic information. Previous studies using animal models have shown that exposure to HT during fertilization can alter the expression of DNA repair genes during early embryonic development31. HT caused by cryptorchidism can lead to a decrease in the expression of repair genes that play a role in the final stage of DNA repair, such as DNA polymerase beta, which promotes the recruitment of DNA ligase III32. Multiple studies have indicated that GMNC play an important role in initiating DNA replication by interacting with Cdc45 and CDK2/cyclin E. GMNC depletion impairs DNA replication in Xenopus laevis oocytes33. In addition, it has been identified that the F-box protein FBXO4 binds to phospho-cyclin D1Thr286 for targeted degradation in response to DNA damage and MAPK signals34,35,36. The replication protein A (RPA) complex, composed of RPA1, RPA2, and RPA3, is the main single-stranded DNA (ssDNA)-binding heterotrimer complex and is involved in various DNA metabolic pathways, including DNA replication, recombination, repair, and DNA damage checkpoints37,38. The absence of RPA in oocytes triggers canonical DNA damage response mechanisms and pathways such as the activation of ATM and P53. RPA deficiency leads to chromosomal translocation in MI and MII oocytes, ultimately resulting in female infertility39. Furthermore, POLE2 (DNA polymerase epsilon 2) is involved in DNA replication, repair, and cell cycle control40. In our study, the expression of GMNC, POLE2, and RPA3 was significantly downregulated in early porcine embryos treated with HT, which supports this point.

ZGA involves a series of dynamic epigenetic changes, including degradation of maternal factors, activation of the zygotic genome, and chromatin recombination, thereby enabling the embryo to achieve pluripotency41,42. The thermal sensitivity of embryos is phased out, and early stage embryos are more susceptible to HT infection29. ZGA in bovine embryos begins after 8 C, and HT treatment at 42℃ affects the development from 1 C to 8 C stages. However, it had almost no effect at the morula and blastocyst stages. Even if early stage embryos survive and continue to develop, the quality of their blastocysts will significantly decrease29,30,43 which suggesting that the normal progression of ZGA is crucial for the embryo to acquire the ability to resist stress. As for which contribution of sperm or oocyte to embryo heat tolerance is more important, studies have compared the differences between Brahman cattle (embryos are less affected by heat stress) and Holstein cattle. It has been found that after fertilization of Holstein cow oocytes with Brahman or Angus breed bull semen (both types of semen have no difference in heat stress response44), there is no breed × temperature interaction, indicating that oocytes play a more important role than sperm in resisting the harmful effects of heat shock in embryos45. Therefore, we compared the RNA-seq data with the pig ZGA gene library and found that 39 DEGs overlapped with the porcine ZGA gene. This indicates that under HT conditions, ZGA of early porcine embryos was significantly affected.

Mitochondria produce cellular ATP via oxidative phosphorylation (OXPHOS), which is crucial for the maintenance of cellular homeostasis. Mitochondria are one of the main organelles affected by heat stress, and in vitro heat damage can induce electron transfer chain (ETC) dysfunction and inhibit mitochondrial ATP synthesis46. There is a tight regulatory connection between the folding of oxidative proteins within the endoplasmic reticulum (ER) and mitochondrial oxidative phosphorylation (OXPHOS). The accumulation of unfolded proteins induces ER stress, leading to the production of reactive oxygen species (ROS) and the activation of cellular signaling mechanisms to counteract the accumulation of unfolded proteins, known as “unfolded protein response” (UPR). One of the three main sensors of UPR is the activating transcription factor 6 (ATF6). Our previous studies have also shown that AFT6 regulates the homeostasis of the ER, Golgi apparatus, and mitochondria, as well as regulates early porcine embryo damage through AIFM2 mediated apoptosis under HT conditions27. In this study, we demonstrated that the OXPHOS pathway was enriched in early porcine embryos treated with HT, and that genes related to mitochondrial ATP synthesis and ETC function were significantly altered. Proteomic analysis of cows has demonstrated that HT leads to reduced expression of proteins involved in OXPHOS and mitochondrial energy pathways, such as NDUFS8 and NDUFV2, which are the core subunits of complex I8. In addition, complex I inactivates and dissociates under HT conditions, slowing down electron flow and reducing oxygen uptake in the ETC, leading to the depletion of ATP synthesis and an imbalance in mitochondrial division/fusion dynamics, inducing mitochondrial dysfunction47. Under acute stress, changes in Ndufa6 transcripts have been detected in the hippocampi of mice using RNA-Seq analysis48. In addition, analysis of multiple myeloma (MM) using single-cell RNA-seq and bulk RNA-seq data indicated that NDUFA6 is a hub gene for MM and can serve as a novel prognostic biomarker for MM49. For the cell cycle, a wide range of stress stimuli, such as oxidants, irradiation, and genotoxins can induce upregulation of CDKN1A, leading to cell aging mediated by upregulation of CDKN1A, which is also a form of cell cycle arrest50,51,52. Previous studies have shown that CDKN1A expression is detected at all stages of mouse pre-implantation embryonic development, and CDKN1A mRNA gradually increases from the 4-cell stage and reaches its peak at the morula stage before decreasing at the blastocyst stage. The expression level of CDKN1A was significantly upregulated when embryos are exposed to heat stress. Knocking down CDKN1A in heat stressed embryos can alleviate developmental delay and reduced blastocyst formation caused by heat stress53. CDKN1A could also mediate cell cycle arrest during the morula/blastocyst stage and apoptosis after the inner cell mass blastocyst stage. The increase in CDKN1A expression level was related to the arrest before and after the morula stage54,55. These findings indicated that NDUFA6 and CDKN1A play crucial regulatory roles in early porcine embryonic development by actively responding to HT signals.

Conclusion

HT affects mitochondrial function, oxidative phosphorylation pathways, DNA polymerase activity, DNA replication, and cell cycle in early porcine embryos. In addition, the expression patterns of ZGA genes in porcine embryos changed under HT conditions. Hub genes are NDUFA6 and CDKN1A.