Next Article in Journal
Why and How European Farmers Are Dedicated to Breeding the Dwarf Dahomey Cattle
Previous Article in Journal
Effect of Dietary Rosemary and Ginger Essential Oils on the Growth Performance, Feed Utilization, Meat Nutritive Value, Blood Biochemicals, and Redox Status of Growing NZW Rabbits
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Animals
Volume 12
Issue 3
10.3390/ani12030376
animals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Time Series Ovarian Transcriptome Analyses of the Porcine Estrous Cycle Reveals Gene Expression Changes during Steroid Metabolism and Corpus Luteum Development

by
Yejee Park
,
Yoon-Been Park
,
Seok-Won Lim
,
Byeonghwi Lim
and
Jun-Mo Kim
*
Functional Genomics & Bioinformatics Laboratory, Department of Animal Science and Technology, Chung-Ang University, Anseong 17546, Gyeonggi-do, Korea
*
Author to whom correspondence should be addressed.
Animals 2022, 12(3), 376; https://doi.org/10.3390/ani12030376
Submission received: 19 November 2021 / Revised: 1 February 2022 / Accepted: 2 February 2022 / Published: 4 February 2022
(This article belongs to the Section Animal Genetics and Genomics)
Browse Figures
Review Reports Versions Notes

Abstract

:

Simple Summary

The estrous cycle, which is divided into follicular and luteal phases based on ovulation, is influenced by reproductive hormones which affect reproduction and cause changes in the reproductive system of the pig. As the main reproductive organ, the ovary is involved in ovulation and changes in the corpus luteum. We aimed to identify dynamic changes in gene expression through differentially expressed gene profiling and to provide a comprehensive understanding of the molecular mechanisms that occur in the pig ovary during the estrous cycle. The transcriptome analysis revealed that the dynamic change in gene expression was more activated in the luteal phase than in the follicular phase. Functional analysis revealed that the metestrus and diestrus periods are important in preparation for pregnancy or the next estrous cycle after ovulation.

Abstract

The porcine estrous cycle is influenced by reproductive hormones, which affect porcine reproduction and result in physiological changes in the reproductive organs. The ovary is involved in ovulation, luteinization, corpus luteum development, and luteolysis. Here, we aimed to provide a comprehensive understanding of the gene expression patterns in porcine ovarian transcriptomes during the estrous cycle through differentially expressed genes profiling and description of molecular mechanisms. The transcriptomes of porcine ovary were obtained during the estrous cycle at three-day intervals from day 0 to day 18 using RNA-seq. At seven time points of the estrous cycle, 4414 DEG were identified; these were classified into three clusters according to their expression patterns. During the late metestrus and diestrus periods, the expression in cluster 1 increased rapidly, and steroid biosynthesis was significant in the pathway. Cluster 2 gene expression patterns represented the cytokine–cytokine receptor interaction in significant pathways. In cluster 3, the hedgehog signaling pathway was selected as the significant pathway. Our study exhibited dynamic gene expression changes with these three different patterns of cluster 1, 2, and 3. The results helped identify the functions and related significant genes especially during the late metestrus and diestrus periods in the estrous cycle.

1. Introduction

The estrous cycle is a major regulatory factor of female fertility and influences the development of the reproductive organs, ovulation, and hormone secretion [1]. It affects litter size and the weight of offspring in the porcine industry [2,3]. Pigs are anatomically similar to humans, and genetic and physiological aspects are also similar to those of humans, so many studies and experiments are being conducted on pigs [4].The 21-day porcine estrous cycle is divided into the follicular and luteal phases based on the ovulation process, and involves various hormones such as gonadotropin releasing hormone, follicle stimulating hormone, luteinizing hormone (LH), progesterone, and estrogen [1]. Therefore, the female endocrine secretions in the estrous cycle are strongly related to the functional roles of the ovary.
The mammalian ovary is not merely a female reproductive organ where oocytes are stored and ovulation occurs, but it also regulates a variety of physiologically diverse aspects of reproduction [5]. After ovulation, the ovulated follicle turns into a corpus luteum (CL), which secretes reproduction-related hormones. Furthermore, the ovary produces steroid hormones (estrogen and progesterone) and peptide growth factors, which play key roles in the functioning of the ovary itself. These are also critical to the regulation of the hypothalamic-pituitary-ovarian axis and the development of secondary sex characteristics [1].
An ovarian transcriptome analysis has been conducted in recent years to identify differentially expressed genes (DEG). Some RNA-seq studies involving litter size in sows have been conducted to further porcine prolificacy [2,3]. A previous study of the ovary screened the major functional genes or molecular markers of estrus expression during diestrus and estrus in gilts from two different breeds [6]. In addition, previous research has provided comprehensive transcriptome data on the porcine ovary at the proestrus and estrus stages through RNA-seq technology [7]. Although we previously analyzed the whole estrous cycle using RNA-seq, we focused on the interaction and connectivity between three reproductive tissues (ovary, oviduct, and endometrium), rather than on examining the changes that occurred during the estrous cycle in a particular tissue [8]. Furthermore, few studies have been carried out on the entire gene expression changes during the development of follicles and CL in the porcine estrous cycle [8,9]. In this study, we aimed to provide a comprehensive understanding of the gene expression patterns in porcine ovarian transcriptomes during the entire estrous cycle through DEG profiling, cataloging, and description of molecular mechanisms.

2. Materials and Methods

2.1. Ethics Statement

All animal experiments were approved by the Institutional Animal Care and Use Committee of the National Institute of Animal Science, Republic of Korea (No. 2015-137). All experimental procedures were carried out in accordance with the Guide for Care and Use of Animals in Research and reported according to ARRIVE guidelines (https://arriveguidelines.org, accessed on 20 April 2017).

2.2. Animals and Sampling

We used 22 crossbred (Landrace × Yorkshire) gilts of approximately similar age (6–8 months) and weight (100–120 kg) who had undergone at least two normal duration estrous cycles. The pigs were selected from different pens from a single research farm in the National Institute of Animal Science (NIAS) to follow estrus behavior. Gilts were subject to estrus detection daily in the presence of boars; the first day of detected estrus behavior was designated as day 0. At days 0 (D00; n = 3), 3 (D03; n = 2), 6 (D06; n = 3), 9 (D09; n = 3), 12 (D12; n = 4), 15 (D15; n = 4), and 18 (D18; n = 3) gilts were sacrificed, and whole ovaries were dissected aseptically from the reproductive organs as previously described [10,11] in the euthanized gilts (Figure 1). The collected ovaries were ground whole and used as samples for the study. Ovarian tissues were snap-frozen in liquid nitrogen and stored at −80 °C for RNA extraction. These sampling materials were used in the previous study to integrate different kinds of female reproductive tissues and identify the core regulation module in the integrative gene expression networks [8].

2.3. RNA Extraction, Library Preparation, and Sequencing

Total RNA was extracted from the ovarian tissues using the TRIzol reagent (Invitrogen, Life Technology, Carlsbad, CA, USA) in accordance with the manufacturer’s recommendations. RNA integrity was assessed using electrophoresis in 1% agarose gel, and RNA with a RIN value of 7 or higher was used. The quantity was validated using the NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Individual libraries were generated using the reagents provided in the Illumina TruSeq RNA sample preparation kit as described previously [12]. The sequencing was performed using the constructed cDNA library with the 100-base pair (bp) paired-end method on an Illumina HiSeq 2000.

2.4. Data Processing and DEG Identification

FastQC (v0.11.4) software was used for the quality checking tool of raw reads [13], and the reads were trimmed using Trimmomatic (v0.38) for adapters and low quality [14]. Hisat2 (v2.1.0) was utilized to map the clean reads against the porcine reference genome (Sus scrofa Sscrofa11.1.98, GCA_000003025.6) of the Ensembl genome browser (http://www.ensembl.org/Sus_scrofa/, accessed on 31 July 2020) as the default option of the program [15]. Samtools (v1.9) was used to sort the mapped reads and change the SAM file to a BAM file. The raw counts corresponding to the genes for each library were calculated based on exons in Sus scrofa GTF v98 (Ensembl) as the genomic annotation reference file using featureCounts (Subread package v1.6.3) [16]. The DEG analyses were performed at each time point (D03, D06, D09, D12, D15, and D18) and compared with D00, which resulted in six comparisons (D03 vs. D00, D06 vs. D00, D09 vs. D00, D12 vs. D00, D15 vs. D00, and D18 vs. D00). All DEG analyses for obtained raw counts were performed using the Bioconductor edgeR package v3.28.1 [17]. The genes were removed when all samples had raw counts of ≤10 to reduce statistical bias in the DEG analyses. Normalization for raw counts was performed using the trimmed mean of M-values method [18]. Limma, an R package, was used for performing MDS to identify the similarity among samples [19], and the ggplot2 R package was used for visualization of the MDS and volcano plot [20]. The DEG were identified based on a negative binomial generalized linear model, the p-value was adjusted using the Benjamini–Hochberg correction with a false discovery rate of <0.05, and absolute log2 fold change (FC) ≥ 1 was applied as the cutoff for DEG [21].

2.5. Functional Validation of DEG

Gene clustering analysis and visualization were conducted using Multi Experiment Viewer (MeV v4.9.0) to identify similarly expressed patterns of time-serial DEG [22]. Stringent thresholds (absolute log2 FC ≥ 2) for DEG were selected for the k-medians clustering analysis with 1K iterations. To annotate gene ontology (GO) terms and Kyoto encyclopedia of genes and genomes (KEGG) pathways using DEG corresponding to 7 time points and each cluster, the database for annotation, visualization, and integrated discovery (DAVID) v6.8 was used as a functional annotation tool [23]. Biological process (BP), cellular component (CC), and molecular function (MF) were included in GO terms. Next, REVIGO treemap [24] and ClueGo (v2.5.5) plugin of Cytoscape (v3.7.2) [25] were used to visualize BP, CC, and MF GO terms. A pie graph was composed with the pathways which met a threshold of p-value < 0.05, and the interval of GO tree was from the minimum level of 4 to the maximum level of 8. In addition, pathways with at least five genes were functionally grouped to define term–term interrelations based on their kappa score level (≥0.3). The treemap was modified to show representative terms among several similar terms, and the pie graph shows the proportion of each term in the functional analysis results. Subsequent analyses were performed using the DEG and functional annotation results on D06.

2.6. Gene Modulation in the KEGG Pathway Analysis

Considering the results of the BP terms in each cluster, functionally similar KEGG pathways were selected as a representative term. Next, the modulations of responsible genes (proteins) in the selected representative KEGG pathway through functional analyses were verified using the clusterProfiler R package [26]. Gene modulations were performed to identify how many core genes are expressed in the selected KEGG pathways and which parts of the pathway the gene is involved in. The genes with the greater log2 FC value are shown in dark red, whereas the genes with the lower log2 FC values are shown in dark blue in the selected pathway. Among the genes corresponding to each protein, genes indicating the absolute maximum change were used as a representative value. Additionally, the core enriched genes of the significant pathway were expressed as heatmaps.

3. Results and Discussion

3.1. Integrative RNA-seq Analysis and DEG Profiling during the Estrous Cycle

In the current study, we divided the porcine estrous cycle into seven time points based on the development of follicles (D00), ovulation (D03), luteinization and luteal development (D06–12), luteolysis (D15–18), and hormone secretions and generated RNA-seq data for each time point (Figure 1). For these seven time points, 430 million paired-end sequence reads were produced from 22 samples with an average of 19.5 million reads per sample. The average numbers of the clean reads after trimming were 16,428,136 (D00), 16,024,399 (D03), 17,400,840 (D06), 16,910,350 (D09), 15,876,680 (D12), 16,541,879 (D15), and 17,577,294 (D18). The average unique mapping rate (94.52%) and the average overall mapping rate (98.36%) were determined (Table 1). The RNA-seq data revealed a separation between the late metestrus and diestrus phases (D06, D09, and D12) and other phases (D00, D03, D15, and D18) in the multidimensional scaling (MDS) plot, which may indicate that important changes occurred during the late metestrus and diestrus phases compared to the other phases (Supplementary Figure S1).
Identification of DEG in time series analysis is a key step towards better understanding the functional role of genes during the entire estrous cycle [27]. In our study, a total of 4414 DEG were identified from all time points based on D00 (Figure 2a and Supplementary Table S1). The number of DEG on D06 (n = 3538) and D09 (n = 3163) increased sharply in comparison to that on D03 (n = 965) and then decreased dramatically on D15 (n = 232) and D18 (n = 45) compared to that on D12 (n = 1683). This might be attributed to the fact that the oocyte-releasing follicles in the ovary develop into CL after ovulation (D06), which causes rapid physiological and morphological changes in the ovary [1]. Furthermore, a large proportion of downregulated DEG found on D06 corresponding to the metestrus, D09, and D12 corresponding to the diestrus phase were identified (D06: 55.90%; D09: 53.68%; D12: 52.88%). The overlapped DEG between D06 and D09 was the largest with 1167 genes (Supplementary Figure S2).
A clustering analysis was applied to characterize the DEG according to similar expression patterns during the overall estrous cycle. As a result, the DEG were grouped into three clusters showing a similar pattern with dynamic expression changes (Supplementary Table S2). A total of 832 DEG were included in cluster 1, which was mainly upregulated on D06 and D09 (Figure 2b), whereas 1205 DEG in cluster 2 were involved with the downregulated pattern during the same period (Figure 2c). Cluster 3 included 436 DEG, which represented the only remarkable downregulation on D06 (Figure 2d). These clustered expression patterns indicated the processes of ovulation, luteinization, CL development, and luteolysis during the ovarian estrous cycle. In the process of CL development, apoptosis, angiogenesis, steroidogenesis, signal transduction, translation, cell proliferation, and tissue remodeling occur [28]. These findings show that marked changes in the ovarian transcriptome occur more in the luteal phase rather than in the follicular phase.

3.2. Functional Annotations and Enriched Pathway Analyses in Clusters

The GO and KEGG enrichment analyses were conducted to discover the connection between the estrous cycle and its molecular functions. The treemaps and pie graphs described the enriched terms in BP of GO, and the KEGG pathway terms were organized into bar plots (Figures 3, 5 and 7, and Supplementary Figure S3). For BP terms, a total of three KEGG pathways were selected for each cluster (Supplementary Figures S4–S6). Heatmaps were drawn to visualize specific DEG, based on D06, when the most dynamic change occurred in all three clusters (Figures 4, 6 and 8). This dynamic change in D06 could be attributed to the process of follicles transforming into CL after ovulation [1].

3.2.1. Reproductive Hormones in Luteal Phase

In cluster 1, intestinal absorption and sterol biosynthesis were the most significant representative GO terms (Figure 3a), and similarly, sterol metabolic process and lipid biosynthetic process were the most significant terms in the pie graph (Figure 3b). As indicated in cluster 1, lipid metabolism (steroid biosynthesis), secretory and digestive system (bile secretion, pancreatic secretion, insulin secretion, gastric acid secretion, salivary secretion, and fat digestion and absorption), and amino acid (alanine, aspartate, and glutamate metabolism; valine, leucine, and isoleucine degradation; taurine and hypotaurine metabolism, and tyrosine metabolism)-related terms were the main pathways (Figure 3c). These significant terms were consistent with the changes that occur in hormones involved in the estrous cycle. Steroid biosynthesis related to steroid hormones and the digestive system including bile, pancreas, and insulin are also part of the lipid digestion process. In addition, they are directly and indirectly related to hormone secretion [29]. Several terms appear to display similar functions in the GO and KEGG results (sterol metabolic process, lipid biosynthetic process, and steroid biosynthesis). Based on BP terms, steroid biosynthesis was chosen as the major pathway in cluster 1. There were several KEGG orthologues related to steroid biosynthesis. Most of the upregulated genes shown in Figure 4a were DEG and were related to the enzymes and reactions for synthesizing cholesterol. Steroid biosynthesis in the ovary is very important for proper ovarian function. Previous studies have shown that rat androstenedione synthesis maintains the function of the CL and increases progesterone biosynthesis [30]. It has also been shown that androgens can maintain luteal function in mice lacking the classic progesterone receptor [31]. Therefore, it can be concluded that various steroid hormones are involved in steroid biosynthesis and progesterone formation in the CL. Progesterone especially is related to the morphological change from corpora lutea to corpus luteum [32]. In cluster 1, the fold change in the expression level of 10 steroid biosynthesis DEG (HSD17B7, DHCR7, LSS, SQLE, SC5D, FDFT1, NSDHL, MSMO1, DHCR24, and TM7SF2) were visualized at the time point of D06 through a heatmap (Figure 4b). HSD17B7 participates in post-squalene cholesterol biosynthesis, which is a leading substance in steroids [33]. LSS, SC5D, NSDHL, DHCR24, and TM7SF2 are involved in all stages of cholesterol biosynthesis, along with SREBF1, which is involved in lipid homeostasis and can be presumed to take part in steroid biosynthesis for pregnancy preparation during the late metestrus and diestrus (Figure 4a,b) [34].
Steroid hormones in the ovary are closely related to blood flow in the uterus; estrogen decreases blood flow, whereas progesterone inhibits its action [35]. Decreased blood flow in the uterus and the effect of hormones decreases the quality of oocytes, which causes infertility [36]. Therefore, proper blood supply to the uterus indicates endometrial acceptability and is associated with successful fertilization [37]. Late metestrus to diestrus (D06–D09), corresponding to all clusters, are the periods when CL development is most active, resulting in increased concentrations of progesterone. From this, it can be hypothesized that the expression of genes selected for DEG during this period had consequences associated with the development of the CL at the time and consequently an increase in progesterone concentrations. Progesterone also suppresses the development of new follicles by inhibiting the activity of the hypothalamus and the pituitary gland. Owing to these processes, the recurrence of the estrous cycle is controlled during the luteal phase or gestational phase. Next, in the latter part of diestrus (approximately D12), luteolysis begins if fertilization has not occurred. In addition, this might reflect that the amplitude of LH secretion in the luteal phase is similar to that in the gestational phase [1]. This implies that during the diestrus phase before luteolysis, the ovary forms substantial structures similar to those in the gestational phase for the preparation of pregnancy. As part of the preparation, the CL synthesizes progesterone, which is essential for the establishment and maintenance of early pregnancy. Therefore, the CL secretes progesterone to maintain its structure and continues to secrete progesterone during the maintenance process to repeat and maintain the process. In addition, progesterone prevents the degeneration of CL and maintains its structure during pregnancy [31].

3.2.2. Immune Responses in Luteal Phase

The greatest enriched GO term for cluster 2 was neutrophil chemotaxis (Figure 5a). Immune-related terms, such as regulation of immune response and phagocytosis, constitute a large proportion of enriched GO terms (Figure 5b). KEGG pathways involved in immune disease (rheumatoid arthritis, asthma, autoimmune thyroid disease, allograft rejection, inflammatory bowel disease, graft versus host disease, and hematopoietic cell lineage), infectious disease (Staphylococcus aureus infection, leishmaniasis, and tuberculosis), and signaling molecules and their interactions (cytokine–cytokine receptor interaction, cell adhesion molecules, and neuroactive ligand–receptor interaction) were enriched in cluster 2 (Figure 5c). The KEGG pathway cytokine–cytokine receptor interaction was chosen as the significant pathway in consideration of BP terms in cluster 2. There are several KEGG orthologues related to cytokine–cytokine receptor interaction. One-third of these are expressed, and most downregulated genes were differentially expressed as shown in Figure 6a. Ovulation occurs at approximately D03, and the ovulated oocyte undergoes a chemical reaction with the sperm within two or three days, resulting in the establishment of pregnancy. It appears that immune-related pathways are downregulated to accept cells from other individuals to prepare for pregnancy until luteolysis begins around D12. Sperm influx causes the female immune system to recognize them as a foreign antigen and respond; however, the immune response at that time is downregulated [38,39]. It is postulated that the organs of the same reproductive system are organically connected, so that the immune system of the ovary may also change under the influence of the oviduct that the sperm cells enter. Cytokines stimulate and inhibit the development of steroids by ovarian cells, and particularly, inflammatory cytokines can regulate the number of leukocytes and the function of the endocrine system in the ovary. In addition, cytokines act on follicle growth, activation of leukocytes, ovulation, luteinization, and luteolysis and are also involved in both inhibition and stimulation of follicular responses to gonadotropin [40]. As important mediators of the immune response, cytokines are involved in pregnancy and can stimulate or inhibit cell growth, induce chemotaxis of cells, and regulate cell differentiation and the expression of other cytokines. In the cytokine–cytokine receptor interaction in cluster 2, INHA, IL33, IL12, RB1, GDF5, LTA, CXCR6, IL2RA, CSF2RB, IL17C, TNFRSF9, IL7R, BMP15, CCL22, IL25, CCR4, CCL17, TNFSF15, CX3CR1, TNFSF14, TNFSF8, TNFSF13B, TNFSF13, IL20RA, CSF3R, IL11, CXCL14, and TNFRSF1B were chosen and visualized using a heatmap (Figure 6b). More importantly, CCL17, CCL22, and CCR4 are structurally related molecules that regulate cell trafficking of various types of leukocytes and play a fundamental role in T-cell development, homeostasis, and the immune system [41]. Both TNFSF9 and the IL family also contribute to regulating the function of immune cells and play an anti-tumor role [42]. CXCL14 and CSF3R play major roles in the growth and differentiation of granulocytes [43]. This led us to assume that the downregulation of these genes is closely related to the activation of leukocytes required for ovulation and immune response related to pregnancy preparation during the late metestrus and diestrus periods (D06–D09).

3.2.3. Cell Proliferation and Growth in Luteal Phase

The treemap showed negative regulation of interleukin-6 production as the highly enriched processes in cluster 3 (Figure 7a). Other biological processes enriched in cluster 3 included the morphogenesis-related terms such as negative regulation of the BMP signaling pathway, skeletal system morphogenesis, and embryonic forelimb morphogenesis (Figure 7b). In the KEGG bar plot of cluster 3, embryo development-related pathways (the hedgehog signaling pathway and osteoclast differentiation) and immune system-related terms (complement and coagulation cascades, Staphylococcus aureus infection, and basal cell carcinoma) were observed as the major pathways (Figure 7c). Considering the BP term (negative regulation of interleukin-6 production, negative regulation of BMP signaling pathway, skeletal system morphogenesis, and embryonic forelimb morphogenesis), we focused on the hedgehog signaling pathway among the KEGG pathways (Figure 8a). In this pathway, the displayed DEG were downregulated on D06 (Figure 8b).
The hedgehog signaling pathway affects embryonic development and morphogenesis and is known to be involved in follicular growth and development in mouse and Drosophila ovaries [44]. The current study revealed that PTCH2, GLI1, GLI2, ARRB1, ARRB2, HHIP, LRP2, GRK3, and IHH were significantly downregulated in the hedgehog signaling pathway. Among them, IHH, GLI1, and HHIP were identified as the noticeable DEG in the late metestrus stage because of their high expression levels. IHH is involved in cell proliferation and differentiation, especially bone cell development [45]. The GLI protein, encoded by the GLI gene, acts as a factor for hedgehog signaling, determining the fate of many types of cells and most organs during embryonic development and contributing to cell proliferation and pattern formation [46]. Similarly, HHIP is known as a gene involved in cell proliferation, growth, and morphogenesis [47]. Among the DEG of the hedgehog signaling pathway, studies have shown that transcription factors GLI1 and HHIP are expressed in the CL and act as either transcriptional activators or repressors [48]. This pathway was downregulated just after the follicles ovulated and before the CL began to develop in earnest. This was because the adequate differentiation of ovarian cells and normal follicular development are essential for ovulation and subsequent CL formation [49,50]. Thus, this pathway reflects when both follicular growth and CL development have occurred. Therefore, it can be postulated that the DEG in the hedgehog signaling pathway involved in cell growth and proliferation is downregulated because it affects the development of follicular growth in the ovary.

4. Conclusions

In conclusion, this study aimed to show the dynamic changes in gene expression that take place in the ovary of the nonpregnant gilt during the estrous cycle. According to the results of the transcriptome analysis, the gene expression wave was more active during the luteal phase than in the follicular phase. Our results showed that the late metestrus and diestrus periods exhibited a remarkable expression pattern, and we confirmed that these periods are important to prepare for pregnancy or for the next estrous cycle after ovulation. Therefore, it is necessary to understand the mechanisms of ovarian regulation during this period and study the relationship between pregnancy and estrous cycles that may occur later. Our study suggests a promising basis for the understanding of successful reproduction at the genomic level in the porcine industry and is an informative resource for further studies on the porcine ovary.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ani12030376/s1, Figure S1: Multidimensional scaling plot of total data, Figure S2: Upset plot of DEG in 6 time points, Figure S3: Biological Process (BP) of Gene Ontology (GO) treemaps for each time point (by DEG of Figure 1C), Figure S4: Log2 fold change in genes of the steroid biosynthesis (KEGG) during the entire estrous cycle, Figure S5: Log2 fold change in genes of the cytokine-cytokine receptor interaction (KEGG) during the entire estrous cycle, Figure S6: Log2 fold change in genes of the Hedgehog signaling pathway (KEGG) during the entire estrous cycle, Table S1: Differentially expressed gene (DEG) list for each tested time point of the porcine estrous cycle, Table S2: Differentially expressed gene (DEG) list for each tested time point of the porcine estrous cycle after clustering analysis.

Author Contributions

Y.P., Y.-B.P. and S.-W.L. analyzed the data and prepared the main figures. Y.P. wrote the main manuscript. B.L. edited the manuscript. J.-M.K. designed and supervised the project. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by Korea Institute of Planning and Evaluation for Technology in Food, Agriculture and Forestry (IPET) through High Value-added Food Technology Development Program (or Project), funded by Ministry of Agriculture, Food and Rural Affairs (MAFRA) (321037051WT011).

Institutional Review Board Statement

The study was conducted according to the guidelines of the Institutional Animal Care and Use Committee of the National Institute of Animal Science, Republic of Korea (No. 2015-137).

Informed Consent Statement

Not applicable.

Data Availability Statement

The data generated in this study are deposited in the NCBI SRA data base under accession no. SRP127622.

Acknowledgments

This research was supported by the Chung-Ang University Research Grants in 2020. I would like to offer my special thanks to Rajesh Kumar Pathak and other Bigkim lab members.

Conflicts of Interest

The authors declare that there are no conflicts of interests.

References

  1. Soede, N.; Langendijk, P.; Kemp, B. Reproductive cycles in pigs. Anim. Reprod. Sci. 2011, 124, 251–258. [Google Scholar] [CrossRef]
  2. Zhang, X.; Huang, L.; Wu, T.; Feng, Y.; Ding, Y.; Ye, P.; Yin, Z. Transcriptomic analysis of ovaries from pigs with high and low litter size. PLoS ONE 2015, 10, e0139514. [Google Scholar]
  3. Zhang, F.P.; Ran, X.Q. Changes in Ovary Transcriptome and Alternative splicing at estrus from Xiang pigs with Large and Small Litter Size. bioRxiv 2019, 547810. [Google Scholar] [CrossRef]
  4. Meurens, F.; Summerfield, A.; Nauwynck, H.; Saif, L.; Gerdts, V. The pig: A model for human infectious diseases. Trends Microbiol. 2012, 20, 50–57. [Google Scholar] [CrossRef] [PubMed]
  5. Edson, M.A.; Nagaraja, A.K.; Matzuk, M.M. The mammalian ovary from genesis to revelation. Endocr. Rev. 2009, 30, 624–712. [Google Scholar] [CrossRef]
  6. Chu, Q.; Zhou, B.; Xu, F.; Chen, R.; Shen, C.; Liang, T.; Li, Y.; Schinckel, A.P. Genome-wide differential mRNA expression profiles in follicles of two breeds and at two stages of estrus cycle of gilts. Sci. Rep. 2017, 7, 5052. [Google Scholar]
  7. Yang, S.; Zhou, X.; Pei, Y.; Wang, H.; He, K.; Zhao, A. Identification of differentially expressed genes in porcine ovaries at proestrus and estrus stages using RNA-Seq technique. BioMed Res. Int. 2018, 2018, 8. [Google Scholar] [CrossRef] [PubMed]
  8. Kim, J.-M.; Park, J.-E.; Yoo, I.; Han, J.; Kim, N.; Lim, W.-J.; Cho, E.-S.; Choi, B.; Choi, S.; Kim, T.-H. Integrated transcriptomes throughout swine oestrous cycle reveal dynamic changes in reproductive tissues interacting networks. Sci. Rep. 2018, 8, 5436. [Google Scholar] [CrossRef]
  9. Tanski, D.; Skowronska, A.; Tanska, M.; Lepiarczyk, E.; Skowronski, M.T. The In Vitro Effect of Steroid Hormones, Arachidonic Acid, and Kinases Inhibitors on Aquaporin 1, 2, 5, and 7 Gene Expression in the Porcine Uterine Luminal Epithelial Cells during the Estrous Cycle. Cells 2021, 10, 832. [Google Scholar] [CrossRef] [PubMed]
  10. Murray, F., Jr.; Bazer, F.W.; Wallace, H.; Warnick, A. Quantitative and qualitative variation in the secretion of protein by the porcine uterus during the estrous cycle. Biol. Reprod. 1972, 7, 314–320. [Google Scholar] [CrossRef]
  11. Buhi, W.; Vallet, J.; Bazer, F. De novo synthesis and release of polypeptides from cyclic and early pregnant porcine oviductal tissue in explant culture. J. Exp. Zool. 1989, 252, 79–88. [Google Scholar] [CrossRef] [PubMed]
  12. Srikanth, K.; Kwon, A.; Lee, E.; Chung, H. Characterization of genes and pathways that respond to heat stress in Holstein calves through transcriptome analysis. Cell Stress Chaperones 2017, 22, 29–42. [Google Scholar] [CrossRef] [PubMed]
  13. FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online: https://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed on 19 March 2021).
  14. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef]
  15. Kim, D.; Paggi, J.M.; Park, C.; Bennett, C.; Salzberg, S.L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 2019, 37, 907–915. [Google Scholar] [CrossRef]
  16. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930. [Google Scholar] [CrossRef] [PubMed]
  17. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed]
  18. Robinson, M.D.; Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11, 1–9. [Google Scholar] [CrossRef]
  19. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef]
  20. Wickham, H. ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 2011, 3, 180–185. [Google Scholar] [CrossRef]
  21. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodol.) 1995, 57, 289–300. [Google Scholar] [CrossRef]
  22. Howe, E.; Holton, K.; Nair, S.; Schlauch, D.; Sinha, R.; Quackenbush, J. Mev: Multiexperiment viewer. Biomed. Inform. Cancer Res. 2010, 18, 267–277. [Google Scholar]
  23. Dennis, G.; Sherman, B.T.; Hosack, D.A.; Yang, J.; Gao, W.; Lane, H.C.; Lempicki, R.A. DAVID: Database for annotation, visualization, and integrated discovery. Genome Biol. 2003, 4, R60. [Google Scholar] [CrossRef]
  24. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef]
  25. Bindea, G.; Mlecnik, B.; Hackl, H.; Charoentong, P.; Tosolini, M.; Kirilovsky, A.; Fridman, W.-H.; Pagès, F.; Trajanoski, Z.; Galon, J. ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 2009, 25, 1091–1093. [Google Scholar] [CrossRef]
  26. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R package for comparing biological themes among gene clusters. Omics J. Integr. Biol. 2012, 16, 284–287. [Google Scholar] [CrossRef]
  27. Lim, B.; Kim, S.; Lim, K.-S.; Jeong, C.-G.; Kim, S.-C.; Lee, S.-M.; Park, C.-K.; Te Pas, M.F.; Gho, H.; Kim, T.-H.J.V.R. Integrated time-serial transcriptome networks reveal common innate and tissue-specific adaptive immune responses to PRRSV infection. Vet. Res. 2020, 51, 1–18. [Google Scholar] [CrossRef] [PubMed]
  28. Bharati, J.; Mohan, N.H.; Kumar, S.; Gogoi, J.; Kumar, S.; Jose, B.; Punetha, M.; Borah, S.; Kumar, A.; Sarkar, M. Transcriptome profiling of different developmental stages of corpus luteum during the estrous cycle in pigs. Genomics 2021, 113, 366–379. [Google Scholar] [CrossRef] [PubMed]
  29. Dorfman, R.I.; Ungar, F. Metabolism of Steroid Hormones; Burgess Publishing: New York, NY, USA; London, UK, 1953. [Google Scholar]
  30. Goyeneche, A.A.; Calvo, V.; Gibori, G.; Telleria, C.M. Androstenedione interferes in luteal regression by inhibiting apoptosis and stimulating progesterone production. Biol. Reprod. 2002, 66, 1540–1547. [Google Scholar] [CrossRef]
  31. Christenson, L.K.; Devoto, L. Cholesterol transport and steroidogenesis by the corpus luteum. Reprod. Biol. Endocrinol. 2003, 1, 90. [Google Scholar] [CrossRef]
  32. Cavazos, L.; Anderson, L.; Belt, W.; Henricks, D.; Kraeling, R.; Melampy, R. Fine structure and progesterone levels in the corpus luteum of the pig during the estrous cycle. Biol. Reprod. 1969, 1, 83–106. [Google Scholar] [CrossRef]
  33. Marijanovic, Z.; Laubner, D.; Moller, G.; Gege, C.; Husen, B.; Adamski, J.; Breitling, R. Closing the gap: Identification of human 3-ketosteroid reductase, the last unknown enzyme of mammalian cholesterol biosynthesis. Mol. Endocrinol. 2003, 17, 1715–1725. [Google Scholar] [CrossRef]
  34. Eberlé, D.; Hegarty, B.; Bossard, P.; Ferré, P.; Foufelle, F. SREBP transcription factors: Master regulators of lipid homeostasis. Biochimie 2004, 86, 839–848. [Google Scholar] [CrossRef]
  35. Ford, S.; Reynolds, L.; Farley, D.; Bhatnagar, R.; Van Orden, D. Interaction of ovarian steroids and periarterial al-adrenergic receptors in altering uterine blood flow during the estrous cycle of gilts. Am. J. Obstet. Gynecol. 1984, 150, 480–484. [Google Scholar] [CrossRef]
  36. Nargund, G.; Bourne, T.; Doyle, P.; Parsons, J.; Cheng, W.; Campbell, S.; Collins, W. Associations between ultrasound indices of follicular blood flow, oocyte recovery and preimplantation embryo quality. Hum. Reprod. 1996, 11, 109–113. [Google Scholar] [CrossRef]
  37. Dong, Y.; Cai, Y.; Zhang, Y.; Xing, Y.; Sun, Y. The effect of fertility stress on endometrial and subendometrial blood flow among infertile women. Reprod. Biol. Endocrinol. 2017, 15, 15. [Google Scholar] [CrossRef] [PubMed]
  38. Oertelt-Prigione, S. Immunology and the menstrual cycle. Autoimmun. Rev. 2012, 11, A486–A492. [Google Scholar] [CrossRef] [PubMed]
  39. Alvergne, A.; Tabor, V.H. Is female health cyclical? Evolutionary perspectives on menstruation. Trends Ecol. Evol. 2018, 33, 399–414. [Google Scholar] [CrossRef]
  40. Gougeon, A. Regulation intragonadique de la folliculogenese humaine: Faits et hypotheses. Ann. D’endocrinologie 1994, 55, 63–73. [Google Scholar]
  41. Scheu, S.; Ali, S.; Ruland, C.; Arolt, V.; Alferink, J. The CC chemokines CCL17 and CCL22 and their receptor CCR4 in CNS autoimmunity. Int. J. Mol. Sci. 2017, 18, 2306. [Google Scholar] [CrossRef] [PubMed]
  42. Wu, J.; Wang, Y.; Jiang, Z. TNFSF9 is a Prognostic Biomarker and Correlated with Immune Infiltrates in Pancreatic Cancer. J. Gastrointest. Cancer 2020, 52, 1–10. [Google Scholar] [CrossRef]
  43. O’Meara, T.; Marczyk, M.; Qing, T.; Yaghoobi, V.; Blenman, K.; Cole, K.; Pelekanou, V.; Rimm, D.L.; Pusztai, L. Immunological Differences Between Immune-Rich Estrogen Receptor–Positive and Immune-Rich Triple-Negative Breast Cancers. JCO Precis. Oncol. 2020, 3, 767–779. [Google Scholar] [CrossRef] [PubMed]
  44. Wijgerde, M.; Ooms, M.; Hoogerbrugge, J.W.; Grootegoed, J.A. Hedgehog signaling in mouse ovary: Indian hedgehog and desert hedgehog from granulosa cells induce target gene expression in developing theca cells. Endocrinology 2005, 146, 3558–3566. [Google Scholar] [CrossRef] [PubMed]
  45. Shimoyama, A.; Wada, M.; Ikeda, F.; Hata, K.; Matsubara, T.; Nifuji, A.; Noda, M.; Amano, K.; Yamaguchi, A.; Nishimura, R. Ihh/Gli2 signaling promotes osteoblast differentiation by regulating Runx2 expression and function. Mol. Biol. Cell 2007, 18, 2411–2418. [Google Scholar] [CrossRef]
  46. Ruiz i Altaba, A. Gli proteins encode context-dependent positive and negative functions: Implications for development and disease. Development 1999, 126, 3205–3216. [Google Scholar] [CrossRef]
  47. Bak, M.; Hansen, C.; Henriksen, K.F.; Tommerup, N. The human hedgehog-interacting protein gene: Structure and chromosome mapping to 4q31. 21→ q31. 3. Cytogenet. Genome Res. 2001, 92, 300–303. [Google Scholar] [CrossRef]
  48. Russell, M.C.; Cowan, R.G.; Harman, R.M.; Walker, A.L.; Quirk, S.M. The hedgehog signaling pathway in the mouse ovary. Biol. Reprod. 2007, 77, 226–236. [Google Scholar] [CrossRef]
  49. Maruo, T.; Katayama, K.; Barnea, E.R.; Mochizuki, M. A role for thyroid hormone in the induction of ovulation and corpus luteum function. Horm. Res. Paediatr. 1992, 37, 12–18. [Google Scholar] [CrossRef] [PubMed]
  50. Acosta, T.; Miyamoto, A. Vascular control of ovarian function: Ovulation, corpus luteum formation and regression. Anim. Reprod. Sci. 2004, 82, 127–140. [Google Scholar] [CrossRef]
Animals 12 00376 g001 550
Figure 1. Overview of ovary sampling throughout the porcine estrous cycle. Stages of the estrous cycle in porcine ovary at seven time points (D00, D03, D06, D09, D12, D15, and D18).
Figure 1. Overview of ovary sampling throughout the porcine estrous cycle. Stages of the estrous cycle in porcine ovary at seven time points (D00, D03, D06, D09, D12, D15, and D18).
Animals 12 00376 g001
Animals 12 00376 g002 550
Figure 2. Dynamic changes of differentially expressed genes (DEG) during the porcine estrous cycle. Dynamic views and expression patterns of DEG using a volcano plot and k-medians clustering algorithm in Multi Experiment Viewer (MeV). (a) Comparing the start of estrous cycle (D00) with each time point (D03, D06, D09, D12, D15, and D18). The x and y axes of the volcano plots were scaled using log2 fold changes and −log10 p-value. Significant DEG (false discovery rate < 0.05; log2FC ≤ |1|: total, log2FC ≥ 1: up, log2FC ≤ −1: down; the total number of DEG: n) in the volcano plots are represented by six different colors corresponding to each time point. The line graph reveals the changes in the number of DEG (up, down) for each time point. (b–d) Gene clustering analysis disclosed three clear expression patterns for each time point (N: the number of DEG against each cluster).
Figure 2. Dynamic changes of differentially expressed genes (DEG) during the porcine estrous cycle. Dynamic views and expression patterns of DEG using a volcano plot and k-medians clustering algorithm in Multi Experiment Viewer (MeV). (a) Comparing the start of estrous cycle (D00) with each time point (D03, D06, D09, D12, D15, and D18). The x and y axes of the volcano plots were scaled using log2 fold changes and −log10 p-value. Significant DEG (false discovery rate < 0.05; log2FC ≤ |1|: total, log2FC ≥ 1: up, log2FC ≤ −1: down; the total number of DEG: n) in the volcano plots are represented by six different colors corresponding to each time point. The line graph reveals the changes in the number of DEG (up, down) for each time point. (b–d) Gene clustering analysis disclosed three clear expression patterns for each time point (N: the number of DEG against each cluster).
Animals 12 00376 g002
Animals 12 00376 g003 550
Figure 3. Visualization of functional analyses of cluster 1. (a) A biological process (BP) treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.
Figure 3. Visualization of functional analyses of cluster 1. (a) A biological process (BP) treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.
Animals 12 00376 g003
Animals 12 00376 g004 550
Figure 4. Pathway-based expression information including DEG of cluster 1. (a) Log2FC change in the genes of the steroid biosynthesis (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Figure 4. Pathway-based expression information including DEG of cluster 1. (a) Log2FC change in the genes of the steroid biosynthesis (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Animals 12 00376 g004
Animals 12 00376 g005 550
Figure 5. Visualization of functional analyses of cluster 2. (a) A BP treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of KEGG pathways.
Figure 5. Visualization of functional analyses of cluster 2. (a) A BP treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of KEGG pathways.
Animals 12 00376 g005
Animals 12 00376 g006 550
Figure 6. Pathway based expression information including DEG of cluster 2. (a) Log2FC change in the genes of the cytokine–cytokine receptor interaction (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Figure 6. Pathway based expression information including DEG of cluster 2. (a) Log2FC change in the genes of the cytokine–cytokine receptor interaction (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Animals 12 00376 g006
Animals 12 00376 g007 550
Figure 7. Visualization of functional analyses of cluster 3. (a) A BP treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of KEGG pathways.
Figure 7. Visualization of functional analyses of cluster 3. (a) A BP treemap was constructed using REVIGO. (b) BP pie graph using Cytoscape (** p < 0.05). (c) Bar graph of KEGG pathways.
Animals 12 00376 g007
Animals 12 00376 g008 550
Figure 8. Pathway-based expression information including DEG of cluster 3. (a) Log2FC change in the genes of the hedgehog signaling pathway (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Figure 8. Pathway-based expression information including DEG of cluster 3. (a) Log2FC change in the genes of the hedgehog signaling pathway (KEGG) using clusterProfiler at D06 and DEG in the pathway were marked. (b) Heatmap of selected DEG on comparing D06 with D00.
Animals 12 00376 g008
Table
Table 1. Summary of the RNA-sequencing results and mapped reads alignment of porcine ovarian samples obtained during the estrous cycle.
Table 1. Summary of the RNA-sequencing results and mapped reads alignment of porcine ovarian samples obtained during the estrous cycle.
GroupSampleRaw ReadsClean Reads
Rate (%)		Uniquely Mapped
Reads Rate (%)		Overall Alignment Rate (%)
Day 0D00C-Ovary-119,596,59716,133,880 (82.33)95.1298.36
D00C-Ovary-219,851,79416,852,302 (84.89)94.8098.31
D00C-Ovary-319,119,94016,298,225 (85.24)94.8498.49
Day 3D03C-Ovary-118,523,86615,737,437 (84.96)94.7298.39
D03C-Ovary-319,201,08816,311,360 (84.95)94.8398.46
Day 6D06C-Ovary-119,286,19916,446,174 (85.27)94.9398.56
D06C-Ovary-220,334,81717,479,844 (85.96)94.4698.50
D06C-Ovary-321,255,42218,276,501 (85.99)94.0798.30
Day 9D09C-Ovary-119,901,12317,088,002 (85.86)94.3198.42
D09C-Ovary-221,040,93418,103,889 (86.04)94.1098.45
D09C-Ovary-318,146,29315,539,159 (85.63)94.4798.45
Day 12D12C-Ovary-119,069,93616,168,172 (84.78)94.8498.45
D12C-Ovary-218,537,73415,955,652 (86.07)93.3998.28
D12C-Ovary-318,664,74916,005,643 (85.75)93.5098.00
D12C-Ovary-417,946,02015,377,254 (85.69)94.0598.32
Day 15D15C-Ovary-118,858,26815,953,501 (84.60)94.5798.32
D15C-Ovary-218,255,01915,646,924 (85.71)94.5698.34
D15C-Ovary-319,733,96216,849,829 (85.38)94.5598.28
D15C-Ovary-420,762,17617,717,261 (85.33)94.4798.17
Day 18D18C-Ovary-122,291,46919,003,553 (85.25)94.6798.31
D18C-Ovary-219,663,92016,727,802 (85.07)94.8898.33
D18C-Ovary-319,968,68817,000,527 (85.14)95.2098.36
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Park, Y.; Park, Y.-B.; Lim, S.-W.; Lim, B.; Kim, J.-M. Time Series Ovarian Transcriptome Analyses of the Porcine Estrous Cycle Reveals Gene Expression Changes during Steroid Metabolism and Corpus Luteum Development. Animals 2022, 12, 376. https://doi.org/10.3390/ani12030376

AMA Style

Park Y, Park Y-B, Lim S-W, Lim B, Kim J-M. Time Series Ovarian Transcriptome Analyses of the Porcine Estrous Cycle Reveals Gene Expression Changes during Steroid Metabolism and Corpus Luteum Development. Animals. 2022; 12(3):376. https://doi.org/10.3390/ani12030376

Chicago/Turabian Style

Park, Yejee, Yoon-Been Park, Seok-Won Lim, Byeonghwi Lim, and Jun-Mo Kim. 2022. "Time Series Ovarian Transcriptome Analyses of the Porcine Estrous Cycle Reveals Gene Expression Changes during Steroid Metabolism and Corpus Luteum Development" Animals 12, no. 3: 376. https://doi.org/10.3390/ani12030376

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop