High-Altitude Stress Orchestrates mRNA Expression and Alternative Splicing of Ovarian Follicle Development Genes in Tibetan Sheep

Simple Summary To achieve optimal growth performance and improved fertility in animals living on high plateaus, it is important to understand how high-altitude stress reduces fertility in females. This study analyzed the transcriptome dynamics of Tibetan sheep ovaries under high-altitude stress. High-altitude stress suppressed the expression of follicular development marker genes and impaired the luteinizing hormone/follicle-stimulating hormone signaling pathway. High-altitude stress also increased abnormally spliced isoforms of transcription factors and RNA processing factors. Therefore, high-altitude stress may reduce the fertility of Tibetan sheep by disrupting the normal expression/hormone signaling of follicular development genes. Further work is needed to decipher whether this phenomenon is a unique feature of Tibetan sheep or a general mechanism in animals under high-altitude stress. Abstract High-altitude stress threatens the survival rate of Tibetan sheep and reduces their fertility. However, the molecular basis of this phenomenon remains elusive. Here, we used RNA-seq to elucidate the transcriptome dynamics of high-altitude stress in Tibetan sheep ovaries. In total, 104 genes were characterized as high-altitude stress-related differentially expressed genes (DEGs). In addition, 36 DEGs contributed to ovarian follicle development, and 28 of them were downregulated under high-altitude stress. In particular, high-altitude stress significantly suppressed the expression of two ovarian lymphatic system marker genes: LYVE1 and ADAMTS-1. Network analysis revealed that luteinizing hormone (LH)/follicle-stimulating hormone (FSH) signaling-related genes, such as EGR1, FKBP5, DUSP1, and FOS, were central regulators in the DEG network, and these genes were also suppressed under high-altitude stress. As a post-transcriptional regulation mechanism, alternative splicing (AS) is ubiquitous in Tibetan sheep. High-altitude stress induced 917 differentially alternative splicing (DAS) events. High-altitude stress modulated DAS in an AS-type-specific manner: suppressing skipped exon events but increasing retained intron events. C2H2-type zinc finger transcription factors and RNA processing factors were mainly enriched in DAS. These findings revealed high-altitude stress repressed ovarian development by suppressing the gene expression of LH/FSH hormone signaling genes and inducing intron retention of C2H2-type zinc finger transcription factors.


Introduction
The growth and reproduction of animals are highly sensitive to changes in their surrounding environment [1]. Most animals use reproductive strategies such as seasonal collected from two-year-old nulliparous ewes. For the high-altitude stress sample group, 6 ewes were randomly selected from 42 nulliparous Tibetan sheep ewes at two-year-old age. Similarly, a total of 6 ewes were also randomly selected from 90 two-year-old ewes as the low-altitude sample group. The left ovary of selected ewes was collected, and the blood and stains on the surface of the ovarian tissue were quickly washed twice with 1 × PBS prepared with 4 • C precooled RNase-Free water. The surface liquid was quickly sucked dry with a dust-free paper towel, and the tissues were quickly stored in a liquid nitrogen cylinder for transportation.
The experiment was conducted during the local slaughter season, and the ewes were sampled at government-designated slaughter and processing sites, in strict accordance with the "Regulations on the Administration of Livestock and Poultry Slaughter in Qinghai Province", as well as the relevant regulations on experimental animal welfare and ethical review by the Ethics Committee of the Laboratory Animal Management Committee of Qinghai Academy of Animal Husbandry and Veterinary Sciences (protocol code QHMKY-2020-09 on 24 June 2020).

Library Construction and Sequencing
Total RNAs were isolated from the ovary of Tibetan sheep using the RNAprep pure Tissue Kit (Tiangen, Cat.DP431). RNA-seq libraries were constructed using three micrograms of highly qualified total RNAs (RNA integrity > 8) by TruSeq RNA Library Prep Kit v2 (Illumina, Cat. RS-122-2001) [12]. The RNA-seq libraries were then sequenced on the Illumina Hiseq X Ten platform to generate pair-end 150 bp reads.

Analysis of Alternative Splicing
RNA-seq reads were qualified and mapped using the same strategy as in the section: Analysis of RNA-seq. The output sam files of HISAT2 were converted to bam files and sorted by Samtools (version 1.9) [25]. The sorted bam files were then normalized according to RPKM by deepTools (version 3.5.1) [26]. The gtf files were assembled from the normalized bam files using StringTie (version 2.1.6) [27]. The merged gtf files were generated by TACO (version 0.7.3) [28]. The merged gtf was annotated to the sheep genome using gffcompare (version 0.11.2) [29]. Differential alternative splicing events were identified by rMATS (version 4.1.2) using the criteria: false discovery rate (FDR) < 0.05, |IncLevelDifference| > 0.1 [30].

Experimental Validation
RT-PCR and qRT-PCR experiments were conducted as described in previous studies [12,31], and the results of RT-PCR were visualized in 1% agarose gels. The primers used in this study were listed in Table S1.

High-Altitude Stress Affected Gene Expression in Tibetan Sheep
To elucidate the effect of high-altitude stress on Tibetan sheep ovary development, we performed transcriptome analysis of Tibetan sheep ovaries at different altitudes. Nine RNAseq libraries were constructed and sequenced using total RNA isolated from high-and lowaltitude Tibetan sheep ovaries. Illumina RNA-seq generated a total of 221 million reads, and over 96.56% of the reads were mapped to the sheep genome (Table S2). The ovaries of Tibetan sheep at low altitude (3000 m.a.s.l, Haiyan, Qinghai, China) were used as the control group, and the ovaries of Tibetan sheep at high altitude (4400 m.a.s.l, Yushu, Qinghai, China) were used as the high-altitude stress group.
A total of 18,561 genes (>30 reads per gene) were expressed in the ovaries of Tibetan sheep. In total, 104 genes were identified as differentially expressed genes (DEGs) under high-altitude stress (Table S3). The number of downregulated genes (64) in high-altitude ovaries was 1.60 times higher than that of upregulated genes (40) (Figure 1a, Table S3). To understand whether the repression of gene expression at high altitude is a global event or occurs only in certain gene groups, we analyzed the expression patterns of all expressed sheep genes at different altitudes. The results showed that the mean expression counts of all sheep genes did not show any clear pattern between high-and low-altitude sheep ( Figure S1). Therefore, we believe that high-altitude stress is not a globally regulated expression of Tibetan genes but is more likely to affect the expression of certain gene groups.

High-Altitude Stress Repressed the Expression of Follicle Development Genes
To further understand the function of high-altitude-related DEGs, we performed KEGG and gene ontology (GO) analyses on the DEGs. KEGG analysis revealed that DEGs were mainly enriched in environmental information processing and cancer-related pathways ( Figure S2). The results of GO analysis showed that high-altitude stress-related DEGs were over-represented in several biological processes, including regulation of vascular development, ovulation cycle, ossification, response to cAMP, regulation of cytokine production, response to lipopolysaccharide, positive regulation of cell motility, positive regulation of cell migration, epithelial cell proliferation, cellular responses to growth factor stimulation and cell adhesion ( Figure 1b). Importantly, these processes are closely related to the regulation of ovarian follicle development or ovulation [32][33][34][35][36][37][38]. A total of 36 DEGs were identified from these GO terms, many of which were well-characterized ovarian follicle development or ovulation-regulating genes (Table S4). For example, lymphatic vessel endothelial receptor-1 (LYVE1) is a molecular marker of lymphatic vessels and is used to distinguish lymphatic vessels from blood vessels [39]. LYVE1 has been reported to be involved in the development of the ovarian lymphatic vasculature and is critical for follicular fluid control [33,40]. In the present study, LYVE1 was one of the most significantly repressed genes under high-altitude stress (Table S4), suggesting that high-altitude stress may inhibit the development of the ovarian lymphatic system and follicular fluid through suppressing the expression of LYVE1. Another example is ADAMTS-1 (a disintegrin and metalloproteinase 1 with a thrombospondin motif); ADAMTS-1 is a central regulator of structural remodeling during ovarian follicle growth, and knockout of ADAMTS-1 results in a defective lymphatic network in mouse ovary [32,33,41]. The downregulation of ADAMTS-1 under high-altitude stress unveiled that high-altitude stress may repress the lymphatic development of ovarian follicle via suppressing the expression of ADAMTS-1. Interestingly, we also observed that 28 DEGs from these GO terms were downregulated under high-altitude stress ( Figure 1b, Table S4) and, together with their functional enrich- stress may inhibit the development of the ovarian lymphatic system and follicular fluid through suppressing the expression of LYVE1. Another example is ADAMTS-1 (a disintegrin and metalloproteinase 1 with a thrombospondin motif); ADAMTS-1 is a central regulator of structural remodeling during ovarian follicle growth, and knockout of ADAMTS-1 results in a defective lymphatic network in mouse ovary [32,33,41]. The downregulation of ADAMTS-1 under high-altitude stress unveiled that high-altitude stress may repress the lymphatic development of ovarian follicle via suppressing the expression of ADAMTS-1. Interestingly, we also observed that 28 DEGs from these GO terms were downregulated under high-altitude stress (Figure 1b, Table S4) and, together with their functional enrichment in ovarian follicle development (Figure 1b), we speculate that high-altitude stress may inhibit ovarian follicle development in Tibetan sheep.
Transcription factors are components of hormone signal transduction pathways. To further explore the central regulatory network of high-altitude stress-related DEGs, we used STRING to perform protein-protein interaction network (PPI) analysis on DEGs encoding genes [22]. PPI network analysis revealed that the number of edges for the PPI-related DEG was 31, which is much larger than the expected number of edges (12) with a PPI enrichment p-value of 3.18 × 10 −6 ( Figure 1d). The significant enrichment of the PPI network further highlights that there were more interactions between proteins encoded from DEGs than expected, and that these proteins may be biologically partially connected as a group [22]. To further explore the biologically connected group of these factors, we performed local network clustering using STRING. The results showed that c-fos/v-fos, and protein fosb, and the basic region leucine zipper of hormone ligand binding were significantly enriched biologically groups, with FDR values equal to 4.46 × 10 −5 and 2.40 × 10 −4 , respectively (Figure 1e, Table S6). A total of six genes were identified from these two biologically connected groups, which were FOS, DUSP1, FOSB, NR4A3, FKBP5, and EGR1 ( Figure 1e, Table S6). Interestingly, four of them (FOS, DUSP1, FOSB, and EGR1) were responses to cAMP and they were downregulated under high-altitude stress (Figure 1f, Table S4). cAMP acts as a secondary messages of luteinizing hormone (LH) and follicle-stimulating hormone (FSH) in the ovary [34]. cAMP is not only responsible for activating primordial follicles into follicle growth, but also plays an important role in the growth and development of all other ovarian follicular stages [35]. Therefore, we suggest that altitude stress may inhibit follicle development by impairing follicle-stimulating hormone signaling transduction in Tibetan sheep.

Alternative Splicing Is Widespread in Tibetan Sheep
Alternative splicing (AS) is an important post-transcriptional regulation mechanism that regulates mRNA fate and diversifies the proteome [12]. Alternative splicing in Tibetan sheep has not been characterized, and whether high-altitude stress affects alternative splicing in livestock remains elusive. To identify alternative splicing in Tibetan sheep and decipher the relationship between AS and high-altitude stress, we performed alternative splicing analysis using RNA-seq by bioinformatics software. Among the 18,561 expressed genes (>30 reads per gene), about 52.01% (9654/18,561) of Tibetan sheep genes underwent alternative splicing events (40,773 AS events). Next, we classified expressed genes into four types based on the number of AS events per gene: none (AS events = 0), low (AS events = 1), moderate (AS events = 2-10), and high (AS events = 2-10) events > 10 ( Figure 3a). AS was prevalent in Tibetan sheep, as more than 70% (6938/9654) of AS-containing genes had two or more AS events (Figure 3a). Therefore, these results suggest that alternative splicing is widespread in Tibetan sheep.  Wilcoxon signed-rank test was applied to evaluate the statistical significance of difference between two samples, and p values are indicated above (* denotes p ≤ 0.05). ns: represents no significance. (f) KEGG pathway analysis of RI type DAS genes. The left vertical axis is the specific name of the KEGG pathway, and the horizontal axis is the significance. The bar of different color represents the classification name of the pathway.

High-Altitude Stress Induced Differential Alternative Splicing in an AS-Type-Specific Pattern
Alternative splicing events can be classified into five types: skipped exon (SE), retained intron (RI), alternative 5 splice site (A5SS), alternative 3 splice site (A3SS), and mutually exclusive exon (MXE) [30]. In Tibetan sheep, SE was the most dominant AS type (63.8%), and RI was the second largest AS type (10.7%). The percentages of the other three AS types (A5SS, A3SS, and MXE) were much closer to each other, and they each accounted for about 8% of the total number of AS events (Figure 3b, Table S7). Importantly, 917 AS events corresponding to 671 genes were identified as high-altitude stress-induced differential alternative splicing (DAS) in Tibetan sheep, accounting for 6.94% (670/9654) of AS-containing genes (Figure 3c, Tables S8 and S9). In particular, the distribution of AS types from DAS was different from that of total AS (Figure 3b,d). Unlike SE in total AS, the most dominant AS type was RI in DAS, and the percentage of RI in DAS (35.6%) was 3.33 times higher than that in total AS (10.7%) (Figure 3b,d). The percentages of A5SS and A3SS in DAS (A5SS: 16.5%, A3SS: 16.4%) almost doubled compared to total AS (A5SS: 7.9%, A3SS: 8.7%) (Figure 3b,d). In contrast, the percentages of both SE and MXE were significantly decreased in DAS compared with total AS (SE: AS vs. DAS, 63.8% vs. 27.3%; MXE: AS vs. DAS: 8.9% vs. 4.4%) (Figure 3b,d). Overall, high-altitude stress induced RI, A5SS, and A3SS-type AS events but suppressed SE and MXE (Figure 3b,d). Thus, these results suggest that high-altitude stress induces DAS in an AS-type-specific pattern.
To explore whether high-altitude stress-related DAS have global expression patterns, we performed expression analysis of different types of DAS with all expressed genes. The expression of genes containing ES and MXE types of DAS was significantly reduced compared with all expressed genes (Figure 3e). In contrast, there was no global expression pattern change for high-altitude stress associated with RI, A5SS, and A3SS ( Figure 3e). Therefore, these results suggest that the expression of ES-and MXE-type DAS genes may be regulated at the transcriptional level, while the expression of RI, A5SS, and A3SS genes may be regulated at the post-transcriptional level. Overall, DAS associated with high-altitude stress may affect gene expression in an AS-type-dependent manner.

C 2 H 2 -Type Zinc Finger Transcription Factors Enriched in High-Altitude Stress-Related DAS
To decipher the function of high-altitude stress-related DAS, KEGG analysis of DAS genes was performed. The results showed that high-altitude stress-related DAS genes mainly acted on RNA splicing, transcription factors, chromosomes and related proteins, DNA repair and recombination, and ribosome biogenesis (Figure 3f, Table S10).
To understand the role of transcription factors on altitude stress, we analyzed the relationship between TF and DAS. About 11.33% of TFs (76/671) had DAS events associated with high-altitude stress (Figure 4a, Table S11). Surprisingly, 73.78% of the DAS-related TFs (56/76) were C 2 H 2 -type zinc finger (ZNF) TFs (Figure 4b, Table S11). Therefore, we speculate that C 2 H 2 -type zinc finger (ZNF) TFs may be regulated by high-altitude-induced DAS.

High-Altitude Stress-Associated Retained Introns Were Enriched in RNA Processing Factors
As retained intron is the major AS type induced by high-altitude stress, we next analyzed whether different types of DAS have distinct enriched biological groups by STRING's protein-protein interaction network (PPI) analysis. PPI analysis revealed no significant associated biological groups in high-altitude stress-related A3SS and MXE (PPI enrichment p-values: 0.101 and 0.58). High-altitude stress-related ES and A5SS had significant PPI-enriched p-values (PPI-enriched p-values: 0.0138 and 0.00529), but they had no significant biological association groups. However, high-altitude stress-associated RI did have significant PPI enrichment and enrichment of biological groups (Figure 4c). PPI network analysis revealed that the high-altitude stress-related RI had a number of edges of 186, which was much larger than the expected number of edges (112) for a PPI enrichment p-value of 9.69 × 10 11 . Local network cluster enrichment of high-altitude stress-related RIs revealed that RNA splicing-related terms, including the thrap3/bclaf1 family and mRNA 3-terminal processing, were biologically significantly enriched, with FDR values of 3.50 × 10 −3 and 7.40 × 10 −3 , respectively (Figure 4d, Table S12). A total of nine genes related to RNA processing were identified, namely SRRM2, CLK1, CLK4, TRA2A, RBM39, ACIN1, PCF11, CLP1, and DDX39B, including splicing and polyadenylation factors. CLK4 is a serine-and arginine-rich (SR) protein-encoding gene that regulates pre-mRNA splicing [45]. According to RT-PCR results, the ratio of RI transcripts to normal transcripts increased from 1.62 at low altitude to 6.0 at high altitude (Figure 4e,f). Thus, high-altitude stress induced the RI of CLK4. Overall, these results suggest that high-altitude stress may affect RNA processing by inducing DAS.

Discussion
High-altitude stress is the main environmental stress for the growth and production of Tibetan sheep. The molecular basis of high-altitude stress response mechanism of Tibetan sheep remains elusive. In the present study, comparative transcriptome analysis was performed using RNA-seq on ovarian tissues of Tibetan sheep exposed to high altitude and low altitude. A total of 104 genes were differentially expressed in high-and lowaltitude Tibetan sheep ovaries (Table S3). Of these, 61.54% (64/104) of DEGs (64/104) were downregulated under high-altitude stress, and the suppression of gene expression under high-altitude stress did not show a global gene expression pattern but may only occur in specific gene groups ( Figure S1). Gene ontology analysis revealed that high-altitude stress-associated DEGs were over-represented in biological processes related to ovarian follicle development or ovulation regulation (Figure 1b).
More importantly, network analysis revealed that hormone-responsive c-fos/v-fos family transcription factors were significantly enriched in the high-altitude stress-related gene expression regulatory network (Table S6, Figure 1e,f). These regulators were FOS, DUSP1, FOSB, and EGR1 (Table S6). In mice, EGR1 significantly induced ADAMTS-1 expression, and an EGR1 binding site (EBS) was found in the ADAMTS-1 promoter [46]. Furthermore, EGR1 regulates ADAMTS-1 expression in a dose-dependent manner [46]. In the present study, the expressions of EGR1 and ADAMTS-1 were both significantly downregulated in Tibetan sheep ovaries under high-altitude stress (Table S4, Figure 2a,f). Taken together, EGR1 and ADAMTS-1 positively regulate follicular development, and EGR1 is an important component of FSH signaling [43,47]. Therefore, high-altitude stress may inhibit follicular development by impairing the FSH-EGR1-ADAMTS-1 signaling pathway. The expression level of TNFα-induced protein 6 (TNFAIP6) is highly increased in follicular cells, and mutation of TNFAIP6 causes infertility in mice [44,48]. The expression of TNFAIP6 in granulosa cells is regulated by the cAMP response element [44]. In this study, high-altitude stress suppressed the expression of TNFAIP6 (Table S4). As the expression of TNFAIP6 is induced by LH, and TNFAIP6 is responsive to cAMP [44], we speculate that high-altitude stress may also inhibit follicular development by impairing the LH-cAMP-TNFAIP6 signaling pathway. Overall, the results of this study suggest that hormonal signaling plays a key role in Tibetan sheep ovarian tissue response to high-altitude stress.
Alternative splicing is one of the important post-transcriptional regulatory mechanisms in response to environmental stress [12]. Abnormal expression of splicing factors and aberrant splicing are closely related to the development of ovarian cancer [49]. In this study, high-altitude stress induced 917 AS events in Tibetan sheep (Figure 3c). In particular, RI but not ES was the predominant AS type in high-altitude stress-related DAS (Figure 3d). Taken together with RI results from aberrant splicing [50], these results indicated that high-altitude stress induced aberrant splicing in Tibetan sheep ovaries. Cdc2-like kinase (CLK4) plays a role in regulating pre-mRNA splicing [45]. In this study, high-altitude stress induced RI and, thus, increased the aberrant splicing isoform of CLK4 (Figure 4e,f). Taken together with mutation of CLK4 reducing cell proliferation in ovaries [45], we speculated that high-altitude stress may repress cell proliferation by increasing abnormally spliced isoform of CLK4. Overall, our results unveiled that high-altitude stress indeed affects the follicular development at both the transcriptional and post-transcriptional levels.

Conclusions
Tibetan sheep live in the harsh environment of Qinghai-Tibet Plateau and high-altitude stress significantly reduces their fertility. Understanding the molecular mechanism of Tibetan sheep responding to high-altitude stress will provide important genetic resources for the research and application of high-altitude response in sheep. In this study, we characterized 104 high-altitude stress-related genes, and follicular-development-related genes were significantly downregulated under high-altitude stress. More importantly, LH and FSH signaling-related transcription factors were significantly enriched in the differentially expressed genes induced by high-altitude stress. In addition, 917 AS events were characterized in response to high-altitude stress and an increase in abnormally spliced isoforms. High-altitude stress induces RI but inhibits ES-type AS events. Taken together, our findings provide evidence that high-altitude stress represses follicle development genes at both the transcriptional and post-transcriptional levels.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani12202812/s1, Figure S1: The expression counts of all expressed genes in ovary of Tibetan sheep from high altitude and low altitude; Figure S2: KEGG enrichment of high-altitude stress-associated differential expressed genes; Table S1: List of primers used in this study; Table S2: Summary of statistical analysis of RNA-seq data; Table S3: List of differentially expressed genes between high-and low-altitude Tibetan sheep; Table S4: List of key differentially expressed genes under high-altitude stress in Tibetan sheep; Table S5: List of differentially expressed transcription factor genes under high-altitude stress; Table S6: List of PPI enriched clusters for high-altitude stress-associated DEGs in Tibetan sheep; Table S7: Summary of statistical analysis alternative splicing events; Table S8: List of RI type differentially alternative splicing events under high-altitude stress; Table S9: List of SE-type differentially alternative splicing events under high-altitude stress; Table S10: KEGG pathway analysis of high-altitude stress-associated DAS in Tibetan sheep; Table S11: List of DAS-associated transcription factors for under high-altitude stress; Table S12: List of PPI-enriched clusters for high-altitude stress-associated RI in Tibetan sheep.