Differentially expressed transcripts study during pregnancy and postpartum anestrus of yak (Bos grunniens)

Abstract

Background: Yak is the main livestock species in the plateau area, and its reproductive performance is low, usually two years or three years. A very few of yaks recover within a certain period of time after delivery and smoothly enter the next estrous cycle, while most of them enter the postpartum anestrus and show no estrus performance. However, the key biological factors and influencing mechanisms that cause postpartum anestrus in yaks are not clear. Objective: To study the expression of differential transcripts in ovaries of yak during pregnancy and postpartum anestrus. Methods: Each three yaks in pregnancy and anestrus under natural grazing conditions in Haiyan County, Qinghai Province were selected and slaughtered, and their ovaries were collected and sent to Biomarker Technologies. Oxford Nanopore Technologies single-molecule real-time electrical signal sequencing technology was used to perform full-length transcriptome sequencing. Astalavista software was used to identify the types of alternative splicing events in yak estrus and pregnancy, and TAPIS pipeline was used to identify alternative polyadenylation. Results: The results showed that there were 1751 differentially expressed transcripts (DETs) between pregnancy and anestrus in yak, of which 808 were upregulated and 943 were downregulated. GO analysis showed that the biological processes of DETs were mainly reproductive, reproductive and rhythmic processes. KEGG analysis showed that the DET cell junction-related adhesion junction protein (β-catenin) and amino terminal kinase (JNK) were involved in FAs (local adhesion). Phosphatidylinositol-3-kinase (PI3K) is involved in the PI3K/AKT/mTOR signaling pathway. Circadian rhythm output cycle failure (Clock) and brain and muscle tissue aromatic hydrocarbon receptor nuclear transporter-like protein 1 (Bmal1) are involved in circadian rhythm signaling pathway. Conclusion: This study found that β-catenin, JNK, PI3K, Clock and Bmal1 were closely related to postpartum anestrus in yak.

Keywords:

• Yak
• Pregnancy
• Postpartum anestrus
• Differentially Expressed Transcripts

Introduction

The yak (Bos grunniens) is native to the Qinghai-Tibet Plateau and its adjacent pastoral areas, and has important economic value.¹ The reproduction rate of yak is very low.² A very few of yaks recover smoothly into the next estrus cycle within a certain period of time after delivery, and most of them enter the postpartum anestrus due to malnutrition, disease and other reasons, showing no estrus performance,³ more than 90% of postpartum female yaks could not estrus in the same year.⁴ Apathy increases the cost of breeding, slows the turnover of livestock and reduces the production efficiency, which seriously affects the preservation and continuation of the population. However, the key biological factors and influencing mechanisms that cause postpartum anestrus in yaks are not clear.

Transcriptomics refer to a full set of RNA transcripts produced by cells or organisms under specific conditions, including messengers, ribosomes and noncoding RNAs and reflect gene expression profiles related to physiological and pathological conditions.⁵ In recent years, high-throughput gene expression methods have been widely used to discover pathogenesis and therapeutic response markers.⁶ Oxford nanopore sequencing [Oxford Nanopore Technologies (ONT)] is a unique third-generation sequencing technology that uses electrical signals to identify base sequences.⁷ Full-length transcriptome sequencing based on ONT single-molecule real-time sequencing technology requires no interruption of RNA fragments, and full-length cDNA is obtained by reverse transcription. The ultra-long reading of the platform contains a single complete transcript sequence information, and the later analysis does not need to be assembled, and the measurement is obtained. The analysis process of obtaining full-length transcriptome mainly includes three stages: full-length sequence recognition, full-length sequence polish to obtain consistent sequence and consistent sequence de-redundancy. The detailed steps are as follows. (1) The low-quality (length less than 500 bp, Qscore less than 7) sequences and ribosomal RNA sequences in the sequences were filtered from the original sequences, and the full-length sequences were obtained according to the presence of primers at both ends of the sequences. (2) Polishing the full-length sequence obtained in the previous step to obtain a consistent sequence. (3) The obtained consensus sequences were de-redundant according to the results of contig alignment with the reference genome or constructed. The final transcript sequence can be directly used for subsequent analysis of isoform, homologous genes, gene families, SSR, alternative splicing (AS), lncRNA, etc. It guides people to understand this life activity at the center of the central law at a deeper level. In addition, it can also be used to upgrade the annotation of the genome and improve the genome database. Lan et al.⁸ applied RNA-seq technology to compare and analyze the transcriptome data of estrus ovaries of yak and cattle, providing new insights for a comprehensive understanding of the specificity of yak reproduction. Transcriptome study of yak ovary further improves the gene structure of yak and excavates new genes related to reproduction, which provides a basis for depicting the normal transcriptome map of yak ovary and further exploring the reproductive performance of yak.⁹ Huo et al. combined RNA sequencing transcriptomics research with quantitative proteomics analysis to identify genes and proteins related to postpartum ovarian cycle.¹⁰ The Xu et al. study showed the histological characteristics of follicular development in yaks and revealed candidate genes that may play an important role in regulating ovarian activity in seasonal reproduction of yaks.¹¹ However, there are few reports on transcriptome studies of yak ovaries during pregnancy and postpartum anestrus.

The regulation of yak during pregnancy and anestrus is different, but it is related to the time course. At present, it is not clear about the expression of differentially expressed transcripts (DETs) in the ovaries of yak during pregnancy and postpartum anestrus. The DETs in the ovaries of yak during pregnancy and postpartum anestrus were analyzed by transcriptome technology, and the key biological factors of anestrus caused by ovarian stagnation after pregnancy and delivery of yak were sorted out and excavated, which laid a foundation for the study of transcriptome and molecular regulation mechanism of yak anestrus and pregnancy.

Materials and methods

Experimental animals

In September 2019, 16 female yaks calving in April were selected in Haiyan County, Qinghai Province, and estrus was observed twice a day under natural grazing conditions. Male yaks dressed in ‘test cloth’ for natural grazing. The estrous responses of female and male yaks to crawling and hurdles were observed daily. The female yak urinates frequently, and the vagina flows more mucus, even a small amount of blood. The bull climbs up across, which can be determined as estrus. The estrus yak was used as the control group. The estrous yaks were used as the control group until October 2019. After two consecutive rectal examinations, the surface of the ovary was irregular and the corpus luteum of pregnancy appeared on one side, which was determined as pregnancy. The ovarian surface is smooth without large follicles or corpus luteum. If the content of progesterone in plasma is less than 0.5 μg/L, it is determined as postpartum anestrus. After clinical examination, each three yaks during pregnancy and postpartum anestrus were selected for slaughter, and the ovaries were quickly taken out. The ovaries were washed with normal saline and labeled, stored in liquid nitrogen and sent to Baimaike for sequencing.

Screening of differentially expressed genes

During the detection of differentially expressed genes, the genes of samples during pregnancy and postpartum estrus were compared in pairs, and the false discovery rate (FDR) was obtained by correcting the significant P value of the difference. The difference of transcript expression between the two groups was greater than or equal to 2 (Fold Change ≥ 2), and FDR < 0.01 was used as the screening standard.

Variable shear analysis

The variable splicing types of each sample were obtained by Astalavista software, and the software default parameters were used. The above five AS events of transcripts were counted, and the number of predicted AS events in each sample was statistically plotted.

Variable polyadenylation analysis

We use TAPIS pipeline to identify alternative polyadenylation (APA), using software default parameters.

DET screening

The experimental procedure was performed according to the standard protocol provided by ONT, including sample quality detection, library construction, library quality detection and library sequencing. In the detection of DETs, the transcripts of samples during pregnancy and anestrus were compared in pairs, and the Benjamini and Hochberg correction method was used to correct the significant P value (P-value). The corrected Fold Change ≥ 2 and FDR < 0.01 were used as the screening criteria.

DETs GO annotation and enrichment analysis

Biomarker Technologies (BMK) performed GO annotation and enrichment analysis on the screened DETs through the GO database and performed significance analysis on the biological process (BP), molecular function (MF) and cellular component of DETs. With P < 0.01 as the significant threshold, high-frequency annotations with statistical significance were obtained.

KEGG analysis of DETs

BMK performed KEGG annotation and enrichment analysis on the screened DETs through the KEGG database and obtained statistically significant signal transduction pathways with P < 0.01 as the significance threshold.

Results

Analysis of differentially expressed genes

By comparing the differentially expressed genes of yak during pregnancy and anestrus, the number of differentially expressed genes was 841, the number of upregulated genes was 347, and the number of downregulated genes was 494 (Fig. 1).

Figure 1. Differentially expressed genes and transcripts of yak during pregnancy and postpartum anestrus. Note: DEG: Differentially expressed genes; DET: differentially expressed transcripts; upregulated: number of upregulated transcripts or gene; downregulated: downregulate the number of transcripts or gene.

DET AS analysis

The results of AS are shown in Fig. 2. Five AS forms were detected in yak ovarian transcripts: Exon skipping, Intron retention, Alternative 5′ splice site and Alternative 3′ splice site. A total of 2056 AS events were identified during pregnancy, including 854 (41.55%) SE, 528 (25.68%) RI, 342 (16.61%) A3, 290 (14.08%) A5 and 42 (2.04%) MEE. A total of 2487 AS events were identified during anestrus, including 1050 (42.2%) SE, 390 (15.68%) A3, 645 (25.92%) RI, 344 (13.85%) A5 and 58 (2.35%) MEE. SE is the main shear type, and MEE is a rare shear type (Fig. 2).

Figure 2. Pie chart of AS event types in pregnant and postpartum anestrus yaks. (A) Pie charts of AS event types of yak pregnancy genes. (B) Pie charts of AS event types of yak postpartum anestrus genes. A3: Variable transcoding termination site; A5: variable transcription start site; SE: exon jumping; RI: intron retention; MEE: variable exons.

DET analysis

By comparing the DETs of yak during pregnancy and anestrus, a total of 1751 DETs were obtained, with 808 upregulated transcripts and 943 downregulated transcripts (Fig. 1).

Variable polyadenylation analysis of DETs

We used TAPIS pipeline to identify APA and found that a total of 8924 genes with one or more APA sites were identified during yak pregnancy, of which one APA site had the largest number of genes, reaching 3908 (43.73%). The number of genes with two and three APA sites was 2006 (22.43%) and 1077 (12.03%), respectively. Another 829 (9.29%) genes had more than five APA sites. A total of 8899 genes with one or more APA loci were identified in yak during anestrus, of which one APA locus had the largest number of genes, reaching 3869 (43.48%). The number of genes at two and three APA sites was 1954 (21.96%) and 1100 (12.36%), respectively. Another 882 (9.91%) genes had more than five APA sites (Fig. 3).

Figure 3. Distribution of APA loci in genes in pregnant and postpartum anestrus yaks. Note: abscissa: Number of polyadenylation sites; vertical axis: number of genes. (A) Numerical distribution of APA loci of the yak pregnancy gene. (B) Numerical distribution of APA loci of the yak postpartum anestrus gene.

DETs GO annotation and enrichment analysis

The BP enrichment analysis of DETs in yak during pregnancy and anestrus showed that they were mainly enriched in BPs such as reproduction, reproductive processes and rhythmic processes. DETs were distributed in 1,136 GO terms. GO terms related to yak reproduction, reproductive process and rhythm process were selected, and the upregulated DETs were distributed in 51 GO categories. In terms of reproduction, a total of 47 DEGs were involved in reproduction and reproductive processes. Seven DEGs were related to prosodic process terms. The downregulated DETs were distributed in 56 GO categories. In terms of reproduction, a total of 52 DETs were involved in both reproductive and reproductive processes. Six DEGs were related to the rhythm process. The MF enrichment analysis of DETs in yak pregnancy and anestrus showed that a total of 17 functional items were enriched in the GO category. There were 688 DETs associated with binding molecular functional categories, of which 357 were upregulated and 331 were downregulated. There were 359 DETs associated with catalytic activity, of which 142 were upregulated and 217 were downregulated (Fig. 4).

Figure 4. GO annotation classification. Statistical chart of DETs in pregnant and postpartum anestrus yaks. Note: The abscissa is the GO category, the left side of the ordinate is the percentage of transcripts, and the right side is the number of transcripts. This map shows the enrichment of transcripts of each secondary function of GO under the background of DETs and all transcripts, reflecting the status of each secondary function under the two backgrounds. The secondary function with a significant proportion of differences indicates that the enrichment trend of DETs is different from that of all transcripts, and can focus on whether this function is related to differences.

KEGG annotation and enrichment analysis of DETs

KEGG analysis of DETs between yak ovaries during pregnancy and anestrus showed that they were identified in the KEGG database. A total of 905 DETs were enriched in 185 signaling pathways in yak ovarian tissues during pregnancy and anestrus, mainly related to cellular processes, environmental information processing, genetic information processing, metabolism and biological systems. Among them, metabolic and biological systems are 10 subcategories. Genetic information processing and cellular processes are four subcategories. In addition, three sub-categories of environmental information processes are related (Fig. 5)

Figure 5. KEGG annotation and enrichment of DETs in pregnant and postpartum anestrus yaks. Note: The ordinate is the name of the KEGG metabolic pathway, and the abscissa is the number of transcripts annotated to the pathway.

Discussion

At present, AS and APA studies based on third-generation sequencing technology are few and mainly focused on humans and a few model species.¹² The information about AS and APA genes in yak during pregnancy and anestrus is almost blank. The results of this study showed that the AS types of yak during pregnancy and anestrus were mainly SE, followed by A3 and RI, and MEE was the least. This is consistent with the results of John Wiley’s study that in AS events, RI is dominant in plants and SE is dominant in animals.¹³ In this study, the number of AS events identified in yak anestrus genes (2748) was more than the number of AS events during pregnancy (2055), indicating that the AS of yak anestrus genes was more active and had higher protein diversity. The number of APA genes in the yak pregnancy genes (9189) was slightly more than the number of APA genes in the anestrus genes (8924), indicating that the APA of genes during yak pregnancy was more active. APA increased the complexity of the transcriptome of yak during pregnancy and affected the localization, stability, translation efficiency and function of RNA. The information of AS events and APA loci identified in this study can further improve the transcriptome annotation information of yak during pregnancy and anestrus, and also provide reference information for related research in other animals.

Yak postpartum anestrus and cell junction damage

The expression of β-catenin related to cell junction in yak postpartum anestrus was lower than that in pregnancy, and the expression of JNK was higher than that in pregnancy. KEGG analysis showed that β-catenin and JNK were involved in focal adhesion (FA) signaling pathway. Follicle is the basic biological unit of ovarian tissue, mainly composed of granulosa cells, oocytes and membrane cells.¹⁴ The integrity of granulosa cells through cell–cell connections helps maintain the microenvironment of oocyte development.^15,16 Cell junctions are usually specific adhesion structures formed by tight junctions, adhesion junctions and desmosomes.^17,18 They have the functions of maintaining intercellular contact and regulating extracellular and intracellular information exchange.¹⁹ Impaired granulosa cell connection can lead to abnormal follicular development environment. This can cause oocyte necrosis and follicular atresia, thus affecting fertility.²⁰ Recent studies have shown that JNK signaling pathway is involved in the regulation of intercellular junctions, including adhesion junctions,²¹ tight junctions²² and gap junctions.²³ It has been found that JNK signaling pathway affects primordial follicle formation by regulating the expression of E3 ubiquitin ligase MDM2 and degrading E-cadherin in syncytial oocytes.²⁴ Studies have shown that JNK phosphorylation in granulosa cells can lead to a decrease in the expression level of β-catenin, destroy the cytoskeleton and adhesion connection of ovarian granulosa cells and cause damage to the adhesion connection of ovarian granulosa cells, resulting in follicular atresia.²⁵ Study found that beta-catenin causes an increase in cellular autophagy.²⁶ In this study, phosphorylation of the differential transcript JNK and decreased expression of β-catenin during anestrus resulted in negative regulation of the FA signaling pathway, impaired granulosa cell junctions and increased autophagy, which were relevant to the generation of anestrus due to follicular atresia and oocyte necrosis in the postpartum anestrus period of yaks. In the future, we can continue to study the relationship between autophagy and yak anestrus and verify the trend change of autophagy-related genes during anestrus and gestation, and indeed the relationship between anestrus and autophagy in yaks.

Yak postpartum anestrus and follicular atresia

The expression of phosphatidylinositol-3-kinase protein (PI3K) in postpartum anestrus of yak was lower than that in pregnancy. KEGG studies found that PI3K was involved in the mammalian target (PI3K/AKT/mTOR) signaling pathway of phosphatidylinositol-3-kinase/serine kinase/rapamycin. PI3K/AKT/mTOR signaling pathway plays an important role in granulosa cell apoptosis and follicular development.²⁷ PI3K signaling pathway plays a regulatory role in controlling follicular activation pathway.²⁸ During the activation and early development of primordial follicles, the rapid growth of oocytes is related to the activation of PI3K signaling pathway.^29,30 PI3K/AKT signaling pathway is the basic signaling pathway in vivo.³¹ It is important for ovarian dormancy³² and initial follicle activation.³³ More and more evidences show that the regulation of mTOR plays an important role in folliculogenesis, oocyte meiotic maturation, ovarian somatic cell proliferation and steroidogenesis, ovarian aging and embryonic development.³⁴ PI3K/AKT/mTOR signaling pathway activates the expression of phosphorylated PI3K protein in ovarian tissue, activates downstream AkT phosphorylation and mTOR protein through cascade reaction, stimulates granulosa cell growth and promotes primordial follicle development.³⁵ PI3K/AKT/mTOR is a signaling pathway associated with apoptosis and autophagy. In this study, decreased PI3K expression during postpartum anestrus caused negative regulation of the PI3K/AKT/mTOR signaling pathway, which promoted apoptosis and activated autophagy, leading to atresia and degeneration of primordial follicles during postpartum anestrus in yaks, resulting in anestrus. Subsequently, inhibitors of the PI3K/AKT/mTOR signaling pathway can be added to detect changes in autophagy-related genes during the deprivation and gestation periods in yaks to determine whether autophagy is caused by this pathway.

Yak reproduction and circadian rhythm regulation

Circadian rhythm movement output cycle failure (Clock) and aromatic hydrocarbon receptor nuclear transporter-like protein 1 (Bmal1) in brain and muscle tissues of yak postpartum anestrus were downregulated compared with pregnancy. Clock and Bmal1 were involved in the circadian rhythm signaling pathway. The circadian rhythm of mammals is produced by the molecular clock encoded by the gene in the body, which regulates the physiology, metabolism and behavior of the body in a cycle of about 24 h.³⁶ BMAL1 and CLOCK are two core transcription factors of the circadian clock system. In mammals, reproductive physiology is regulated by the circadian clock system.³⁷ Studies have shown that the biological clock system regulates mammalian steroid hormone production,³⁸ follicular development³⁹ and ovulation response.⁴⁰ During follicular development, the rhythm of circadian clock gene expression does not exist in primordial follicles and preantral follicles, appears in early antral follicles and becomes more stable in late antral follicles or preovulatory follicles.⁴¹ Downregulation of Bmal1 gene expression in porcine granulosa cells cultured in vitro affected steroid synthesis related gene expression and progesterone secretion.⁴² Knockout of the circadian clock gene Clock can seriously affect the reproductive capacity of mice, resulting in reduced oocyte release and reduced litter size, which seriously damages reproductive activity and leads to reduced reproductive capacity.⁴³ In this study, the expression of Clock and Bmal1 decreased during postpartum anestrus, resulting in negative regulation of circadian rhythm signaling pathway, which had a certain causal relationship with the instability of follicular development before ovulation, decreased ovulation and decreased progesterone hormone synthesized by steroid synthesis in yak postpartum anestrus. In future experiments, Clock and Bmal1 can be validated to further define its role in this pathway.

Conclusion

In this study, we sequenced and analyzed the transcriptomes of ovarian tissues from yaks during pregnancy and postpartum anestrus using third-generation transcriptomics sequencing technology, and found that DET β-catenin decreased and JNK phosphorylated in the FAs of yaks during the postpartum anestrus period, resulting in impaired cell junctions and increased autophagy; downregulation of PI3K in the PI3K/AKT/mTOR signaling pathway led to apoptosis and autophagy; and downregulation of Clock and Bmal1 in the circadian rhythm signaling pathway led to decreased reproductive capacity of the ovaries. The downregulation of Clock and Bmal1 in circadian rhythm signaling pathway leads to the decrease of ovarian fertility, and all of the above results lead to postpartum anestrus in yaks. The present study is helpful to better analyze and understand the physiological changes of yak ovaries during pregnancy and postpartum anestrus, and has certain reference value for the in-depth study of ovarian physiological characteristics during pregnancy and postpartum anestrus, and for the revelation of the key biological factors in the postpartum anestrus period.

Authors’ contributions

Data curation: J. Ma and S. Hu; Formal analysis: J. Wang; Investigation: Yahua Yang and C. Yang; Project administration: Yanmei Yang; Supervision: Y. Zhaxi and Writing—original draft: Y. Li; Writing—review and editing: Y. Li and S. Huo.

Acknowledgements

The authors wish to thank BMK Cloud for their help in the experiments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that supports the findings of this study are available in the supplementary material of this article.

Additional information

Funding

This research was supported by the Fundamental Research Funds for the Central Universities [31920230026] and Lanzhou Science and Technology Plan Project [No. 2022-2-44].