Transcriptome analysis reveals molecular regulation mechanism of Tibet sheep tolerance to high altitude oxygen environment

Abstract

As one of the most important livestock breeds on the Qinghai-Tibet Plateau, Tibetan sheep are of great importance to the local economy, agriculture and culture. Its adaptive mechanism in low temperature and low oxygen at highland altitudes has not been reported. In this study, transcriptome sequencing was used to analyze the heart, liver, spleen, lung, kidney, and muscle tissue of sheep at low and highland altitudes. LOC101112291, SELENOW, COL3A1, GPX1, TMSB4X and HSF4 were selected as candidate genes for adapting to plateau characteristics in Tibet Sheep. Besides, glutathione metabolism, arachidonic acid metabolism, nucleotide excision repair, regulation of actin cytoskeleton, protein digestion and absorption, thyroid hormone synthesis, relaxation signaling pathways may play important roles in the adaptation to plateau hypoxia, and cold tolerance. Structural analysis also showed that sequencing genes related to the adaptation mechanism of Tibet sheep to highland altitude. This study will lay a certain foundation for Tibet sheep research.

SIMPLE SUMMARY

Tibet sheep are an ancient species in the Qinghai Tibet Plateau. After a long period of domestication. Tibet sheep adapt to the hypoxic environment of the plateau in terms of physiology and morphology. At the same time, Tibet sheep is also one of the major sources of material for herdsmen in tibetan. In this study, six different tissue samples (heart, liver, spleen, lung, kidney, and muscle) of Tibet sheep were analyzed to reveal the underlying mechanisms of different tissues respond to hypothermia condition. The results showed that six key genes and eight important signaling pathways involved in regulating the adaptation of Tibet sheep to the plateau. In addition, there were more alternative splicing (AS) events and single nucleotide polymorphism (SNP) sites in highland altitude Tibet sheep than in lowland altitude sheep, which was also a concern in the highland altitude adaptability of Tibet sheep.

Keywords:

• Tibet sheep
• hypoxia adaptation
• cold resistance
• high UV
• transcriptome

Introduction

The Qinghai Tibet Plateau is known for its harsh environment, such as high altitude, low oxygen, alpine cold and strong ultraviolet rays. These conditions limit the survival of many animals.¹ The low-oxygen environment of the plateau mainly affects the heart and lung function of the human body, and people who can easily cause acute altitude sickness rush to the plateau, even in severe cases. Qinghai-Tibet Plateau is the highest cold region on earth.² Long-term exposure to a cold plateau environment can lead to disorders of local nutrition, function, and cell metabolism, and even cause tissue death and shedding. Therefore, thermogenesis is key to survival during long-time exposure to cold environment.^3–5 Qinghai Tibet Plateau has strong solar ultraviolet radiation and can cause some damage to the body, such as the skin.⁶ In the long-term evolution of organisms, native animals living on the plateau were able to adapt well to the extreme environment of hypoxia at altitude and developed unique adaptability,⁷ such as yaks,^8,9 Tibet pigs^10–13 and Tibet sheep.¹⁴

Tibet sheep is one of three major coarse wool sheep breeds on the Qinghai-Tibet Plateau. Through long-term domestication, Tibet sheep has evolved to adapt to the hypoxic environment of the cold plateau. At present, more than 30 million Tibet sheep live on the 2.5 million km² of the Qinghai Tibet Plateau.¹⁵ Tibet sheep are one of the main sources of production and life for farmers and herdsmen on the Qinghai-Tibet Plateau. Through Tang et al. performed whole genome sequencing of Tibet sheep in the early stages, he did not contrast with sheep at lowland altitudes.¹⁶ Therefore, the molecular mechanism of Tibet sheep adapting to the local alpine and low-oxygen environment is not clear. Studies have shown that the heart,¹⁷ liver,¹⁸ spleen,¹⁹ lungs,²⁰ kidneys,²¹ and muscles²² are the most sensitive tissues under highland altitude hypoxia conditions. Different from other studies, we study the adaptation mechanism of Tibetan sheep to the plateau environment from themselves,²³ sequence the internal organs of Tibetan sheep at different altitudes,²⁴ and our research has studied the adaptation mechanism of Tibetan sheep to high altitude more thoroughly.²⁵ Therefore, in this study, transcriptome sequencing was used to explore how different tissues (heart, liver, spleen, lung, kidney, muscle) respond to highland altitude conditions. Our findings can be used to better understand the adaptive characteristics of Tibetan sheep to highland altitude environments. The key pathways and genes identified in this study can improve sheep’s ability to adapt to extreme environments.

Materials

Ethics statement

All experiments were conducted in the Key Laboratory of Conservation and Utilization of Animal Genetic Resources on the Qinghai-Tibet Plateau of the Ministry of Education. And the animal experiments were approved by the Animal Experiment Ethics Committee of Mian yang Normal University (2020-1-20).

RNA extraction

In this study, one-year-old male Tibetan sheep at highland altitude (from Ruoergai County, Sichuan Province, 3400 m above sea level), yearling Tibetan sheep (n = 3) and lowland altitude Tibetan (from Hu-Zhu County, Qinghai Province, 2550 m above sea level) sheep (n = 3) were collected as experimental animals. Total RNA was extracted from the tissues (lungs F, liver G, muscle J, spleen P, kidneys S, and heart X) using TRIzol® Reagent (Invitrogen, California, USA) according the instructions and genomic DNA was removed using DNase I (TaKara, Japan). RNA degradation and contamination was monitored on 1% agarose gels. Then RNA quality was determined by 2100 Bio-analyser (Agilent Technologies, Majorbio, China) and quantified using the ND-2000 (Nano Drop Technologies). Only high-quality RNA sample (OD_260/280 = 1.8–2.2, OD_260/230 ≥ 2.0, RIN ≥ 8.0, 28S:18S ≥ 1.0) was used to construct the sequencing library.

Library preparation, and Illumina Hiseq-Xten/nova Seq 6000 sequencing

RNA purification, reverse transcription, library construction and sequencing were performed at Shanghai Major Bio-pharm Biotechnology Co., Ltd. (Shanghai, China) according to the manufacturer’s instructions (Illumina, San Diego, CA). The transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA) using 1 μg of total RNA. Shortly, messenger RNA was isolated according to the poly-A selection method by oligo(dT) beads and then fragmented by fragmentation buffer firstly. Secondly double-stranded cDNA was synthesized using a Super Script double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end-repair, phosphorylation and ‘A’ base addition according to Illumina’s library construction protocol. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina Nova-Seq 6000 sequencer (2 × 150 bp read length). Raw data from this study have been deposited into CNGB Sequence Archive²⁶ of China national GeneBank Data Base²⁷ with accession number CNP0003239 (https://db.cngb.org/mycngbdb/submissions/project).

Read mapping

The raw paired end reads were trimmed and quality controlled by Seq-Prep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. Clean reads were aligned to reference genome (Ovis_aries https://www.ncbi.nlm.nih.gov/genome/83?genome_assembly_id=351950) with orientation model through HISAT2 (daehwankimlab.github.io)) software.²⁸ Mapped reads of each sample were assembled by String Tie (https://ccb.jhu.edu/software/stringtie/index.shtml) t = 6 follow the method inference.²⁹ Finally, we verified the reliability of the sequencing data by RT-qPCR and homogenized the data by the 2^−ΔΔCt method.

Differential expression analysis and functional enrichment

To identify DEGs (differential expression genes) between two different samples/altitude, the expression level of each gene was calculated according to the transcripts per million reads (TPM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/)³⁰ was used to quantify gene abundances. Essentially, differential expression analysis was performed using the DESeq2/DEG-seq/edgeR/Limma/NOIseq,^31–33 DEGs with |log₂(foldchange)|≥1 and P-adjust ≤ 0.05 (DESeq2/edgeR/Limma)/P-adjust ≤ 0.001 (DEG-seq)/Prob > 0.8 (NOI-seq) were considered to be significantly differential expressed genes. In addition, functional-enrichment analysis including GO (Gene Ontology, http://www.geneontology.org) and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at P-adjust ≤ 0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do).³⁴

Alternative splice events Identification

All the alternative splice events that occurred in our samples were identified by using the recently releases program RMATS (http://rnaseq-mats.sourceforge.net/index.html).³⁵ Only the isoforms that were similar to the reference or comprised novel splice junctions were considered, and the splicing differences were detected as exon inclusion, exclusion, alternative 5′, 3′, and intron retention events.

Results

Transcriptome assembly and analysis

Transcriptome sequencing was performed to analysis the expression model of the lung (F1, F2), liver (G1, G2), muscle (J1, J2), spleen (P1, P2), kidneys (S1, S2) and heart (X1, X2) in sheep from highland and lowland altitudes. After quality filtering, for example, each group in Table 1 reads Q20 > 97%, Q30 > 93.75%, and the error rate is less than 0.1%. (1) Raw data: the total data amount of raw sequencing (the number of raw reads multiplied by the read length). (2) Clean reading: total number of entries for the sequencing data after quality control. (3) Clean base: the total amount of data after quality control sequencing. (4) Error rate (%): The basic average error rate corresponding to the quality data is generally below 0.1%. (5) Based on quality control (QC) sequencing data. Q20 and Q30 refer to the percentage of total bases with a sequencing mass of 99% and 99.9% or higher, respectively. (6) GC content (%): The sum of G and C bases, corresponding to quality control data, as a percentage of the total base. The total number of transcripts was 87,761, and the maximum length of transcripts over 1800 was 54,761(Fig. 1(a)). The exact match obtained with the constant transcript ‘=’ complete match of intron chains has the highest number of transcripts 57,169, but ‘C’ assembled transcripts contains a minimum of 3 (Fig. 1(b)).

Figure 1. Transcriptome assembly and analysis. (a) Histogram of transcript length distribution. Distribution of different transcript lengths, abscissa represents the transcript length range, the ordinate is the number of transcripts within the transcript length range. (b) Pie chart of the new transcripts. c: complete match of intron chain. i: transfrag falling entirely within a reference intron. j: potentially novel isoform(fragment), at least one splice junction is share d with a reference transcript. o: generic exon overlap with a reference transcript. u: unknown, intergenic transcript. x: exon overlap with reference on the opposite strand. =: complete match of intron chain.

Table 1. Raw data and quality control data of each sample group.

Sample	Raw reads	Raw bases	Clean reads	Clean bases	Error rate (%)	Q20(%)	Q30(%)	GC content (%)
F1_2	70184828	10597909028	69724696	10362681678	0.0242	98.25	95.01	51.08
F1_3	60158728	9083967928	59662896	8867426228	0.0249	97.94	94.42	51.22
F1_1	60191384	9088898984	59773368	8886074238	0.0245	98.11	94.72	51.2
F2_1	57395128	8666664328	56917040	8339123671	0.0256	97.66	93.71	52.53
F2_2	58236352	8793689152	57699692	8457662081	0.0261	97.46	93.25	52.24
F2_3	58646456	8855614856	58073102	8435054579	0.0261	97.43	93.25	52.38
G1_2	48770454	7364338554	48426560	7118203119	0.0241	98.28	95.13	48.95
G1_3	45658614	6894450714	45379100	6691399512	0.024	98.34	95.24	48.78
G1_1	51939330	7842838830	51620544	7630930166	0.024	98.36	95.28	49.58
G2_1	54630018	8249132718	54181782	7998990305	0.0259	97.55	93.37	50.84
G2_2	54665348	8254467548	54187948	7967178683	0.0254	97.71	93.9	50.55
G2_3	48132530	7268012030	47701138	7014135620	0.0258	97.59	93.53	50.58
J1_2	56976690	8603480190	56500314	8357112519	0.0242	98.21	95.1	54.75
J1_3	68167982	10293365282	67676550	10032051439	0.024	98.31	95.23	52.96
J1_1	53105410	8018916910	52618554	7755810803	0.0246	97.99	94.74	55.04
J2_2	51703320	7807201320	51265796	7541916611	0.0254	97.72	93.99	53.78
J2_1	51334980	7751581980	50909520	7534422144	0.0258	97.61	93.53	52.84
J2_3	51449836	7768925236	51011914	7548933721	0.0256	97.69	93.73	52.18
P1_2	63387084	9571449684	62923692	9294107690	0.0244	98.16	94.87	52.69
P1_3	59557648	8993204848	59128530	8742652218	0.0243	98.19	94.93	52.45
P1_1	70396280	10629838280	69880888	10309661616	0.0243	98.2	94.96	52.05
P2_1	56969340	8602370340	56462914	8328895521	0.0257	97.59	93.61	52.99
P2_2	57241976	8643538376	56739414	8346526089	0.0255	97.67	93.83	52.71
P2_3	55596290	8395039790	55130710	8096530489	0.0254	97.7	93.95	52.99
S1_2	60189588	9088627788	59757166	8749824226	0.0243	98.19	94.95	50.78
S1_3	74429028	11238783228	73892794	10848322490	0.0243	98.21	94.97	50.59
S1_1	71248472	10758519272	70693568	10382801789	0.0243	98.23	95	50.7
S2_1	55108082	8321320382	54639638	8010905482	0.0261	97.44	93.2	50.82
S2_2	58902252	8894240052	58243376	8426398064	0.0263	97.3	93.15	51.35
S2_3	55933384	8445940984	55456436	8159529556	0.0257	97.63	93.62	52.35
X1_1	49744378	7511401078	49374424	7250928453	0.0247	98.08	94.55	47.43
X1_2	50084590	7562773090	49815014	7355655965	0.0242	98.29	94.97	47.04
X1_3	50703408	7656214608	50353904	7370170407	0.0243	98.23	95	47.85
X2_1	56123218	8474605918	55613822	8220169376	0.0256	97.68	93.7	51.8
X2_2	51321192	7749499992	50809132	7422988153	0.0259	97.53	93.47	51.5
X2_3	47409426	7158823326	46847682	6840897550	0.0265	97.23	92.91	48.45

Gene annotation and classification

As shown in Table 2, six databases (NR, Swiss-Prot, PFAM, COG, GO, and KEGG) were used to explore the biological functions. All genes and transcripts were annotated and classified according to biological function, of which 23,581 (92.77%), 22,231 (87.4%), 18,206 (71.6%), 22,147 (87.13%), 19,409 (76.36%) and 16,914 (66.54%). The similarities are highest in the NR database.

Table 2. Gene annotation and classification.

GROUP	Expression Gene number(percent)	Expression_Transcript number (percent)	All gene number(percent)	All_Transcript number(percent)
GO	19,409(0.7636)	35,066(0.7361)	23,343(0.6943)	39,774(0.6957)
KEGG	16,914(0.6654)	33,059(0.694)	20,401(0.6068)	37,363(0.6535)
COG	22,147(0.8713)	43,253(0.908)	26,176(0.7785)	48,402(0.8466)
NR	23,581(0.9277)	45,355(0.9521)	28,505(0.8478)	51,482(0.9005)
Swiss-Prot	22,231(0.8746)	43,576(0.9147)	26,321(0.7828)	48,805(0.8537)
Pfam	18,206(0.7162)	37,425(0.7856)	20,031(0.5958)	40,278(0.7045)
Total_annotation	23,608(0.9288)	45,400(0.953)	28,620(0.8512)	51,641(0.9033)
Total	25,419(1.0)	47,638(1.0)	33,623(1)	57,172(1)

Transcriptome data are annotated and classified through 6 databases (NR, Swiss-Prot, PFAM, COG, GO, and KEGG), all genes and transcripts are annotated and classified through six databases, and the number and ratio of transcripts and genes in each gene pool.

Subsequently, the genes were mapped to the GO and COG databases, respectively. In the GO database, genes were divided into three broad categories: biological processes, cell composition, and molecular function. Biological processes have the most genes. Interestingly, muscle (J) tissue expresses the largest number of genes in each function (Fig. 2(a)).

Figure 2. Gene annotation and classification. (a) COG classification chart A-RNA processing and modification. B – Chromatin structure and dynamics. C – Energy production and conversion. D – Cell cycle control, cell division, chromosome partitioning. E – Amino acid transport and metabolism. F – Nucleotide transport and metabolism. G – Carbohydrate transport and metabolism. H – Coenzyme transport and metabolism. I – Lipid transport and metabolism. J – Translation, ribosomal structure and biogenesis. K – Transcription. L – Replication, recombination and repair. M – Cell wall/membrane/envelope biogenesis. N – Cell motility. O – Posttranslational modification, protein turnover, chaperones. P – Inorganic ion transport and metabolism. Q – Secondary metabolites biosynthesis, transport and catabolism. S – Function unknown. T – Signal transduction mechanisms. U – Intracellular trafficking, secretion, and vesicular transport. V – Defense mechanisms. Y – Nuclear structure. Z – Cytoskeleton. (b). GO classification map.

Genes were clustered into 21 functional classes in the COG databases (Fig. 2(b)). Genes are also classified through the KEGG database according to the pathways or functions involved, namely metabolism, genetic information processing, environmental information processing, cellular processes, organism systems and human disease (Supplementary Fig. S1–S3). In all control groups (2 represent lowland altitude, 1 represent highland altitude), most of the genes belonged to lipid metabolism, immune system, and signal transduction.

Figure 3. Gene expression analysis. (a) Correlation analysis of samples of different altitudes. (b) The abscissa represents the sample name/group name of the project, and the vertical coordinate represents the expression value of the target gene/transcript in each sample/group. (c) Circles of different colors indicate the number of genes/transcripts expressed in the sample/group, values indicate the number of common and unique genes/transcripts between different samples/groups, and intersecting regions of the circles represent the total number of genes/transcripts in each sample/group. (d) The lower right corner of the figure is the sample name, the upper left corner is the sample cluster, and the different colors indicate the correlation between the two samples. After the sample passes the dimensional analysis, there are relative coordinate points on the main assembly. The closer each sample point is, the greater the similarity between samples. The horizontal axis represents the contribution of principal component 1 (PC1) to the differential sample in the 3D plot, and the vertical axis represents the contribution of principal component 2 (PC2) to the differential sample in the 2D plot. (e) Number of gene expressions in the same tissue at different altitudes.

Gene expression analysis

Quantify gene/transcript expression levels of genes/transcripts, then convert them to Transcripts Per Million (TPM) to obtain normalized gene/transcript expression levels, and finally perform statistical scoring analysis on each set of transcripts. The results showed that these genes are most in the muscle compared to other tissues (Fig. 3(a)). Subsequent intersections of genes co-expressed in each tissue indicated that a total of 9473 genes were expressed in six tissues (Fig. 3(b)). Gene expression are similar in different altitudes different altitudes and in the same tissue (Fig. 3(c)). Different samples show the same tissue with equal correlated values in different groups by Principal Component Analysis (PCA) (Fig. 3(d)). More importantly, there were more gene expressions in the highland altitude treatment group than lowland altitude (Fig. 3(e)). The above shows that the expression amount is similar in the same tissue of sheep at different altitudes, and the number of genes expressed is different.

Differential expressed genes in Tibet Sheep at different elevation

When p < .05 and |log₂FC| > = 1, the gene considered to be differentially expressed. Differentially expressed genes are selected in various sheep tissues. 533 up-regulated genes and 452 down-regulated genes were found in the lungs (F1 vs F2), 239 up-regulated genes and 205 down-regulated genes in liver (G1 vs G2), 2338 up regulated genes and 1489 down regulated genes in muscle (J1 vs J2), 72 up-regulated genes and 87 down regulated genes in the spleen (P1 vs P2), 180 up-regulated genes and 200 down-regulated genes in the kidney (S1 vs S2), 517 up-regulated genes and 366 down-regulated genes in the heart (X1 vs X2) were screened, respectively (Fig. 4). Most DEGs are expressed at high land altitude, and interestingly, differential genes in the liver are very abundant at both high and low land.

Figure 4. Histogram of differential expressed genes in different tissues. The abscissa represents a different groups of difference comparison. The ordinate represents the corresponding number of up and down-regulated genes/transcripts. Red represents up and green represents down.

Key genes and pathways for resistance to plateau conditions

Six genes (LOC101112291, SELENOW, COL3A1, TMSB4X, HSF4) were identified by gene set analysis of the same tissues (heart, liver, spleen, lung, kidney, and muscle) of sheep at different altitudes (Fig. 5(a)). Our goal is to identify signaling pathways and candidate genes associated with hypoxic adaptation, cold resistance, and anti-ultraviolet. Functional annotation of differentially expressed genes by KEGG. The top 20 enrichment KEGG pathways were shown in Fig. 5(b,c). KEGG subpopulation analysis showed that differentially expressed genes were enriched in 8 signaling pathways (relaxation signaling pathway, protein digestion and absorption, platelet activation, thyroid hormone synthesis, Regulation of actin cytoskeleton, arachidonic acid metabolism, glutathione metabolism, nucleotide excision repair), and human disease-related pathways were excluded from the analysis. The relaxing signaling pathway is a pathway of significant enrichment, while COL3A1 is a gene involved in most cellular pathways. All 8 pathways are enriched, 2 genes are involved in 3 pathways, and 4 genes are involved in 8 pathways, which are shown in Table 3. In this study single gene involving two or more pathways were considered key genes. The main reasons for the adaptation of Tibet sheep to plateau hypoxia, alpine, and UV strength may be the genes GPX1, TMSB4X and COL3A1 (Type III collagen A chain) expression, which plays an important role in the adaptation mechanism of the Qinghai-Tibet Sheep Plateau.

Figure 5. Key genes and pathways for resistance to plateau conditions. (a) Screening of key genes in response to sheep at different altitudes. The different colored circles represent a set of genes from a gene set. These values represent the number of genes is common or unique between two gene sets or three gene sets, the sum of all numbers in the circle indicates the total number of genes in that gene set, and crossed region of the circle represents the common number of genes between each gene set. (b,c) Annotation of KEGG function of key genes.

Table 3. Key genes and their associated pathways and functions.

Gene Name	pathway ID	Pathway
GPX1	map00480	Glutathione metabolism
	map00590	Arachidonic acid metabolism
	map04918	Thyroid hormone synthesis
	map05016	Huntington disease
	map05014	Amyotrophic lateral sclerosis (ALS)
COL3A	map05146	Amoebiasis
	map04974	Protein digestion and absorption
	map04933	AGE-RAGE signaling pathway in diabetic complications
	map04926	Relaxin signaling pathway
	map04611	Platelet activation
TMSB4	map04810	Regulation of actin cytoskeleton
HSF4	No	No
SELENOW	No	No
LOC101112291	map03420	Nucleotide excision repair

RT-qPCR analysis

Three genes were selected from 6 DEGs, which had the most obvious effect on the adaptability of Plateau. Although the results showed differences in RT-qPCR results for individual genes compared to transcriptome sequencing results, the trend in results was consistent. Thus, the RNA-Seq results are reliable, and the comparisons of RT-qPCR with RNA-Seq is shown in Fig. 6.

Figure 6. RT-qPCR Validation. (a–c) On the left showed differences in the expression of the same gene in the same tissue in sheep at different high land. The results of the RT-qPCR analysis on the right showed that there were differences in the expression of the same gene in the same sheep tissue at different altitudes (three replicates per group). ‘*’ represents p < .05, the difference is significant.

Identification of alternative splicing genes

In most eukaryotic cells, alternative splicing (AS) is ubiquitous in gene expression. They perform different functions in different cell types and tissues, and have their own specific expressions and variations. We suspect that this may also be related to the mechanism of adaptation at highland altitudes in Tibet sheep. In this study, we analyzed the results of high-quality sequencing and found many AS events. RMATs (http://Rnaseq-mats.sourceforge.net/index.html) are alternative analysis for RNA-seq data applications that not only classify alternative splicing events but also analyze alternating splicing events to discover differences between different samples. In this regard, the alternative splicing events analyzed by RMAT_S can be divided into the following five categories (Fig. 7(a)). These include exon jumps (Skipped exon, SE), first exon alternative splicing (Alternative 5’splice site A5SS), last exon alternative splicing (Alternative 3’splice site A3SS), exon selective jumps (Mutually exclusion exon MXE), intron were stranded (Retained intron RI). The specific pattern as shown in (Fig. 7(b)). SE (70.42%) was the most frequent event, MXE (10.16%) was the second most frequent event, and the third was the A3SS event with the lowest frequency RI (3.36%). These results shows that the frequency of each group alternating stitching events in each group is significantly different, and the frequency of SE events is the highest (Fig. 7(c)).

Figure 7. Identification of alternative splicing genes. (a) the Mode and mechanism of alternative splicing. (b) Proportion of each alternate splicing event. (c) Frequency of alternating splicing events in sheep tissue at different altitudes.

Transcriptome analysis of SNP features

Transcriptome sequencing obtained 128274 SNPs (single nucleotide polymorphism) sites, where the number of transversion types (including C/T, T/C A/G, G/A) in Fig. 8(a) is higher than the conversion type (Fig. 8(a)), and the conversion types include A/T, A/C, T/G, C/G, T/A, G/T, G/C (Fig. 8(b)). In addition to this, G/C type events are the most of conversion types in each group, and G/A events are the most of the conversion types. In addition, highland altitude tissues (group1) received more SNP events than lowland altitude tissues (group 2).

Figure 8. Single nucleotide polymorphism occurrence events and their distribution. (a) the distribution of conversion type SNP events in each group. (b) the distribution of transversion type SNP events in each group.

Discussion

The aim of our study was to identify key genes and important pathways for Tibet sheep to adapt to cold plateau by comparing transcriptome information of sheep living at different altitude. The analysis of the adaptation mechanism of the Tibet sheep lays the foundation for the development and utilization of Tibet sheep genetic resource. Animal tissues changes during the development stage to adapt to the environment, such as the plateau pika changing the thickness of the middle membrane to adapt low oxygen. Thus, adaptation to highland altitude tends to form a unique genetic mechanism.³⁶

Transcriptome sequencing is an excellent tool for exploring the potential mechanisms by which Tibet sheep adapt to hypoxia, cold resistance, and UV resistance in Tibet sheep. Differentially expressed genes were identified by GO functional enrichment analysis, and it was found that the changes in gene expression were closely related to biological functions. In addition, RNA-seq has been used in many fields, such as gene expression profiling, alternative splicing variants.³⁷ Transcriptome sequencing analysis can be used to study the altered expression of known genes and to speculate on the function of genes that have not yet been studied. It can therefore reveal the molecular mechanisms and mechanisms of action of certain regulatory genes, or can be used to analyze differences in gene expression levels in different tissues or different physiological states in the same tissue.

In this study, transcriptome sequencing was performed on sheep liver, kidney, lung, spleen, muscle, and heart. Screening of key genes of tibet sheep adapted to hypoxia, UV resistance and low-temperature resistance. Xist (Loc101112291) may influence the molecular mechanism of sheep litter size yields through methylation progresses.³⁸ Glutathione peroxidase 1 (GPX1) plays an important role in antioxidant stress, inflammatory signaling, sugar and lipid metabolism, follicle, and metabolic diseases. It is an important enzyme in the body’s antioxidant system. It can not only remove reactive oxygen species from the body and repair damaged cells, but also block further damage from reactive oxygen species to antibodies to oxidative stress in the organism. In addition, GPX1 is widely present in antioxidant enzymes and metalloproteases in organisms. Superoxide dismutase (SOD), GPX1 and catalase (CAT) form the first-line resistance in radiation tissues. SOD converts super oxygen cations into hydrogen peroxide, which is then detoxified into water by GPX and CAT.^39–43 As can be seen from Fig. 7, the expression level of GPX1 in Tibet sheep lower than that at lowland altitude areas, so we infer that GPX1 may play an important role in the UV resistance of Tibet sheep on the plateau. SELENOW (Selenium protein W) plays a role in antioxidants and supports brain function.^44,45 COL3A1, as a type III collagen, is mainly found in the extracellular matrix, which is the most abundant protein in animals, providing support, connection, nutrition, and protection for parenchymal cells, and is involved in maintaining the integrity of smooth muscle. COL3A1 is widely involved in life activities, and studies have shown that COL3A1 gene expression is affected by body weight and genotype and has an important effect on intramuscular collagen content.^46–49 TMSB4X (Thymosin) is a small-molecule protein that is spherical and modulator of actin. In addition to being a component of the cytoskeleton, TMSB4X is expressed in most cells and is involved in many biological activities such as cell migration and adhesion, wound healing, intracellular signal transduction, anti-inflammatory, cardiac development, retinal neonatal tubules, cytoskeletal network, protein catabolism, apoptosis, cell survival and stem cell formation. Actin polymerization and depolymerization are essential for a wide variety of cells, including cell migration, chemotaxis, maintenance of cellular traits and cell division.^50–53 The results showed that TMSB4X was expressed most in Tibet sheep lung tissue, suggesting that TMSB4X may play an important role in the adaptation of Tibet sheep to plateau hypoxia. HSF4 (heat shock transcription factor4) regulates the expression of non-heat shock protein and participates in cell development.⁵⁴ Our experiment was conducted through transcriptome sequencing of different tissues (heart, liver, spleen, lung, kidney, and muscle). Muscle (J) tissue has the largest number of genes, and a reference transcription database has been established to facilitate the application of RNA SEQ in Tibet or other sheep. As an important regulator in response to hypoxia, hypoxia inducible factors(HIFs) are key regulators of gene expression during hypoxia. HIF1α is also a transcription factor for numerous target genes which mediate multiple biological functions, such as angiogenesis, hematopoiesis and the maintenance of vascular tone to provide or replenish tissues with blood and oxygen. HIF1α further mediates cellular response to hypoxia by regulating glucose uptake and anaerobic respiration in oxygen-depleted environments.⁵⁵ HIF3α was considered as negative mediator of HIF-regulated genes, which negatively affects gene expression by competing with HIF-1α. In this study the expression of HIF1α and HIF3α is differ in different land sheep, but was no significant.⁵⁶

The GO database⁵⁷ is suitable for various of species, it can describe genes and proteins. In the DEGs (differentially expressed genes), the results of GO enrichment analysis in this study show that cellular components are enriched in both molecular function and biological process. Many categories are related to biological regulation, cell modalities and binding, which indicates the differences in these aspects of sheep at different altitude sheep.

KEGG pathway analysis can be used to obtain some of the signaling pathways involved in DEGs, as the biological function of these genes can be further understood through the analysis of gene metabolic pathways.⁵⁸ This study shows that the relaxation signaling pathway, protein digestion and absorption, platelet activation, thyroid hormone synthesis, actin cytoskeletal regulation, arachidonic acid metabolism, glutathione metabolism, nucleotide excision repair may become the key pathways for Tibet sheep to adapt to the plateaus, hypoxia, and UV resistant. The relaxation signaling pathway relaxes the smooth muscles of the womb and softens the cervical and pelvic ligaments. Relaxation signaling pathway may also lower acute heart failure and blood pressure to reduce pimples.⁵⁹ The digestion and absorption of proteins can produce many functional peptides, which play an important role in health.⁶⁰ Platelet activation mediates intracellular signaling pathways that regulate a wide range of biological functions in the body. Platelet activation promotes arachidonic acid metabolism, and its receptor PAFR (platelet activation factor receptor) regulates the body’s broad biological function by coupling to G-protein to participate in multiple signaling pathways within cells.^61,62 Thyroid hormone synthesis is responsible for the metabolism of thermogenic substance and the promotion of growth, in addition to being an important mechanism for regulating the levels of active hormone in cells.^63,64 Regulation of actin cytoskeleton plays an important role in many biological procession, cell division, cell enlargement and cell signal transduction.^65,66 Arachidonic acid metabolism plays an important role in blood pressure regulation and inhibition of multiplexes and promotes the occurrence of inflammatory factors.⁶⁷ Glutathione metabolism depletes free radicals produced by oxidative stress in the body, protects cell membranes and ensures health.⁶⁸ Nucleotide excision repair (NER): repairs pyrimidine dimer due to ultraviolet bonding (pyrimidine dimer) DNA crosslinking and bulky adducts.⁶⁹ These pathways maybe pathways for Tibet sheep to adapt to hypoxia, alpine, and UV resistance.

Sites of SNP/Indel are widely distributed in transcriptomes, SNPs and indel sites and their characteristics were analyzed by sheep transcriptome at different altitudes. the number of conversions is much higher than the inverted type, this phenomenon was called conversion bias.⁷⁰ It may be related to evolution selection, rather than random selection.⁷¹ In addition, the distribution characteristics of AS events are also noteworthy, studies have shown that different RNA splice variants of genes increase in response to hypoxia, such as VEGF clippers are expressed in response to hypoxia.⁷² In addition, when dealing with strong ultraviolet light, the RNA splice variants expression of CYPI7B1 increases, and it leads to reduced enzyme activity and (OH)₂D₃ synthesis.⁷³ For environmental adaptation, the RNA of some genetic viruses will adapt to the environment through mutations in order to adapt to new hosts.⁷⁴ These results of this study will provide data and a theoretical basis for subsequent studies of Tibet sheep.

Conlusion

In this study, transcriptome sequencing was used to analyze the heart, liver, spleen, lung, kidney, and muscle tissue in sheep at different altitudes. LOC101112291, SELENOW, COL3A1, TMSB4X and HSF4 were selected as candidate genes for plateau adaptation in Tibet sheep, as well as pathways such as glutathione metabolism, arachidonic acid metabolism, nucleotide excision repair, regulation of actin cytoskeleton, protein digestion and absorption, thyroid hormone synthesis, relaxation signaling may play an important role in adapting to plateau environment. Structural analysis of genes showed that the highland altitude tissue received more SNP events than the lowland altitude tissue, which may also be related to the adaptation mechanism of Tibetan sheep to highland altitude. This study will lay a certain foundation for Tibet sheep research.

Author contributions

Jingling Wang and An Li designed the study and wrote the manuscript; Yanyan Li, Yaq Lin, Yong Wang, Qilin Dai, Shizhang Du and Ru Yi preformed research; Zhoucuo Qi and An Li analyzed data. All authors read and approved this manuscript.

Acknowledgements

The authors acknowledge assistance from the goat lipid metabolism research lab.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data support this fundings of this study are openly available in CNGB Sequence Archive(CNSA) of China National GeneBank DataBase with accession number CNP0003239.

Additional information

Funding

This work was supported by Sichuan Science and Technology Program (2019YJ0533), National Natural Sciences Foundation of China (No. 32072723), and Fundamental Research Funds for the Central Universities, Southwest Minzu University (2021057), and Scientific Research Foundation of Mian Yang Normal University (QD2020A15).