Comprehensive Analysis of Long Non-coding RNA and mRNA Transcriptomes Related to Hypoxia Adaptation in Tibetan Sheep

Tibetan sheep have lived on the Qinghai-Tibet Plateau for a long time, and after long-term natural selection, they have shown stable genetic adaptability to high-altitude environments. However, little is known about the molecular mechanisms of the long non-coding (lnc)RNAs involved in the adaptation of Tibetan sheep to hypoxia. Here, we collected lung tissues from high-altitude Tibetan sheep and low-altitude Hu sheep for RNA sequencing to study the regulatory mechanisms of the lncRNAs and mRNAs in the adaptation of Tibetan sheep to hypoxia. We identified 254 differentially expressed lncRNAs and 1,502 differentially expressed mRNAs. We found 20 pairs of cis-regulatory relationships between 15 differentially expressed lncRNAs and 14 protein-coding genes and two pairs of trans-regulatory relationships between two differentially expressed lncRNAs and two protein-coding genes. These differentially expressed mRNAs and lncRNA target genes were mainly enriched in pathways related to lipid metabolism and immune function. Interaction network analysis showed that 17 differentially expressed lncRNAs and 15 differentially expressed mRNAs had an interactive relationship. Additionally, we used six differentially expressed lncRNAs and mRNAs to verify the accuracy of the sequencing data via qRT-PCR. Our results provide a comprehensive overview of the expression patterns of the lncRNAs and mRNAs involved in the adaptation of Tibetan sheep to hypoxia, laying a foundation for further analysis of the adaptations of plateau animals.

Introduction

The Qinghai-Tibet Plateau is the highest plateau worldwide. It is known as the “Third Pole” of the earth and is known for its low oxygen, low temperatures, and strong ultraviolet radiation. Along with human migration and settlement, many domestic animals have reproduced for generations in these harsh living conditions and have thus adapted unique and distinctive characteristics at the morphological, physiological and genetic levels (1–3). Scientists have used large-scale omics data to reveal the adaptive genetic mechanisms of domestic animals on the Qinghai-Tibet Plateau and identified a number of plateau-adaptive candidate genes, including EPAS1 (endothelial PAS domain-containing protein 1), EGLN1 (Egl-9 homolog 1) and VEGF (vascular endothelial growth factor), in the HIF (hypoxia inducible factor) hypoxia-inducing pathway (4–7).

Long non-coding RNA (lncRNA) is non-coding RNA that contains >200 nucleotides. LncRNA has important roles in many life activities such as the dosage compensation effect, epigenetic regulation, cell cycle regulation and cell differentiation regulation (8, 9). Studies have shown that lncRNA can bind to HIF-1α and activate its expression, thereby playing an important role in hypoxia-induced tumor cells (10–14). LncRNA may also play an important role in the adaptability of domestic animals to plateaus. Recent studies have shown that lncRNA is involved in the adaptability of Tibetan chickens and yaks to plateaus; however, the regulatory mechanisms of this lncRNA remain largely unknown (15, 16).

To reveal the potential role of lncRNA in the adaptability of domestic animals to plateaus, we selected sheep from different altitudes and collected lung tissue for transcriptomic sequencing. We used an integrated analysis of mRNA and lncRNA data to explore the molecular mechanisms of Tibetan sheep's adaptability to plateaus, find the genes and regulatory pathways related to this adaptability, and provide a theoretical basis for Tibetan sheep production and breeding.

Materials and Methods

Sample Collection

We collected lung tissues from six healthy 18-month-old rams, of which, three were from the Qinghai-Tibet Plateau (Zashijia sheep, Qumalai County, Yushu Tibetan Autonomous Prefecture, Qinghai, China, altitude ~4,800 m), and three were from the plains of China (Hu sheep, Minqin County, Wuwei, Gansu, China, altitude ~1,400 m). All six sheep grazed and received supplemented feed and were euthanized by an intravenous injection of phenobarbital solution (Fatal-Plus, 10 mg/kg body weight, Vortech Pharmaceuticals, MI, USA). Tissues were collected from the middle lobe of the right lung and stored in liquid nitrogen, then used for RNA extraction.

RNA Extraction, Library Construction, and Sequencing

We used TRlzol reagent (Invitrogen, Carlsbad, CA, USA) to extract the total RNA and remove rRNA per the manufacturer's instructions. RNA samples that qualified through quality inspection were used for library construction in accordance with the instructions provided by Illumina (Illumina, CA, USA). The RNA was randomly interrupted into short 200–500-nt fragments, and the first cDNA strand was synthesized using random primers. Buffer, dNTPs (dUTP instead of dTTP), RNase H and DNA polymerase I were added to synthesize the second cDNA strand. The cDNA was purified using the QiaQuick PCR kit, then EB buffer was added to elute the proteins, the ends were repaired, and base A and sequencing adapters were added. Uracil-N-glycosylase was used to degrade the second chain. Agarose gel electrophoresis was used to select the appropriate fragments for PCR amplification. Qualified libraries were sequenced using the HiSeq™ 4000 platform.

Raw Data Filtering, Comparison, and Splicing

After sequencing, the raw data (raw reads) were obtained, and fastp software was used for quality control to obtain clean reads. Quality control was performed by removing reads containing adapters, reads with an N content > 10%, reads with only A bases, and low-quality reads (the number of bases with a quality value of Q ≤ 20 accounted for the total number of reads above 50). We used Bowtie2 software to compare clean reads to the rRNA database and remove the ribosomal reads. We used HISAT2 software to compare the clean reads to the sheep reference genome (Oar_v4.0), set the parameter to -rna-strandness RF, and set the remaining parameters to the default values. We used Stringtie software to assemble the clean reads compared to the sheep reference genome and used Cuffmerge software to merge each transcript. We used Cuffcompare software to compare the assembled transcripts with known transcript types (NCBI Refseq, Ensembl transcripts, UCSC), identify known mRNAs, and screen new non-coding transcripts (novel ncRNA).

lncRNA Identification and Target Gene Prediction

To screen for new non-coding transcripts (novel ncRNA), we filtered transcripts with ≥200 bp and ≥2 exons. We used CPC2, CNCI and Pfam software to predict the coding ability of the new transcripts and took the intersection of the transcripts without coding potential as candidate lncRNAs. The lncRNA target genes included cis-target and trans-target genes. mRNA 50 kb upstream and downstream of the lncRNA were used as the cis-target, and mRNA whose expression correlation coefficient was >0.9 was used as the trans-target.

Identification of Differentially Expressed mRNAs and lncRNAs

For the FPKM value, we used StringTie software to calculate the mRNA and lncRNA expression levels, with the equation FPKM = total fragments/mapped reads (millions) × exon length (kB). We used DESeq2 software to analyze the differences between the mRNA and lncRNA, standardize the reads counts, calculate the hypothesis test probability (p-value), and perform multiple hypothesis test corrections to obtain the false discovery rate (FDR). From the difference analysis, P < 0.05 and a fold change > 1.5 were used to screen significantly differentially expressed mRNA and lncRNA.

Functional Enrichment Analysis of Differentially Expressed Genes

We mapped the differential genes to each term in the gene ontology (GO) database (http://www.geneontology.org/), calculated the number of differential genes in each term, and obtained a list of differential genes with specific GO functions. The hypergeometric distribution test was used to calculate the P-value for the significantly enriched GO functions, then the P-value was corrected using the Benjamini-Hochberg multiple test (FDR). GO items with FDR ≤ 0.05 were considered significantly enriched. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to annotate and classify the differentially expressed genes with pathway functions. The KEGG pathway functional enrichment method is similar to the GO functional enrichment analysis.

Construction of the lncRNA-mRNA Network

We used the interaction relationship in the STRING protein interaction database (http://string-db.org) to analyze the differential gene interaction network. The differential gene set was extracted from the database, and the interaction network diagram was constructed and visualized using Cytoscape software.

Validation of lncRNA and mRNA Expression by qRT-PCR

To test the accuracy of the sequencing results, we selected six differentially expressed mRNA and lncRNA for qRT-PCR verification. We used a cDNA synthesis kit to reverse transcribe the extracted total RNA into cDNA. We used Oligo7 software to design primers and perform specific detections in NCBI. We used the TransStart Green qPCR SuperMix and LightCycler 480 II instruments to perform qRT-PCR and ran each sample in triplicate to ensure the accuracy of the quantitative results. The 2^−ΔΔCt method was used to calculate the relative expression of the target genes, and ACTB (beta-actin) was used as the internal reference gene. Supplementary Table 1 lists the primers used in this study.

Statistical Analysis

All data are presented as the mean ± standard deviation and were analyzed using Student's t-test in SPSS software. p < 0.05 was considered statistically significant.

Results

Overview of the Sequencing Data

Using the Illumina HiSeq™ 4000 platform, the six constructed libraries were paired-end sequenced. Table 1 shows the data quality control and genome comparison results. The six libraries produced 499,954,984 raw reads; 498,979,644 clean reads remained after quality control. The average rate of the clean reads was 99.81%. The GC content of the six samples was 47.57–49.25%, which was consistent with the base composition law of Q20 ≥ 97.63% and Q30 ≥ 93.17%. We used HISAT2 software to compare the clean reads to the sheep reference genome. Approximately 92.47% of the reads could be accurately compared, and the match rate was high. Only data aligned to the sheep reference genome were used for further bioinformatics analysis.

TABLE 1

Table 1. Summary of sequenced RNA-seq data.

lncRNA and mRNA Feature Analysis

An average of 17,904 expressed genes and 5,615 lncRNAs were identified from the six libraries. Using CNCI, CPC2, and Pfam to predict the coding ability of the new transcripts yielded 1,740 new lncRNAs (Figure 1A), including 357 sense, 89 antisense, 31 intronic, 37 bidirectional, 941 intergenic, and 285 other lncRNAs (Figure 1B). Most lncRNAs had two or three exons, which was significantly less than the number of exons in the mRNA (Figure 1C). The length distribution range of the lncRNA and mRNA was basically the same, but the lncRNA was longer than the mRNA (Figure 1D). The lncRNA expression level was lower than that of the mRNA (Figure 1E).

FIGURE 1

Figure 1. Characteristics of lncRNAs and mRNAs in the lung tissues of sheep at different altitudes. (A) Screening of candidate lncRNAs by CNCI and CPC2. (B) Data on lncRNA transcript types. (C) Distribution of exon numbers in lncRNAs and mRNAs. (D) Length distribution of lncRNAs and mRNAs. (E) Expression levels of lncRNAs and mRNAs.

Identification and Analysis of Differentially Expressed lncRNA and mRNA

We identified 254 differentially expressed lncRNAs between the Tibetan and Hu sheep, of which, 123 were upregulated, and 131 were downregulated (Figures 2A,C; Supplementary Table 2). The five lncRNAs with the most significant differential expression were MSTRG.19949.1, MSTRG.6422.3, MSTRG.6796.1, MSTRG.16435.3, and MSTRG.19804.1. We identified 1502 differentially expressed mRNAs between the Tibetan and Hu sheep, of which 469 were upregulated and 1,033 were downregulated (Figures 2B,D; Supplementary Table 3). The five most significantly differentially expressed mRNAs were NKIRAS2 (NK-κB inhibitor-interacting Ras-like 2), SEPW1 (selenoprotein W), PDZK1 (PDZ domain-containing 1), PEAK1 (pseudopodium enriched atypical kinase 1), and KIAA1549. Among these, NKIRAS2 is involved in immune response, SEPW1 is involved in oxidative stress, and PDZK1 is involved in fat metabolism. Additionally, the differentially expressed lncRNA and mRNA were combined into one group, and the differentially upregulated and downregulated lncRNA and mRNA were combined into one group (Figures 2E,F). The differentially expressed lncRNA and mRNA had good reproducibility in the groups.

FIGURE 2

Figure 2. Expression analysis of lncRNAs and mRNAs. (A) Numbers of upregulated and downregulated differentially expressed lncRNAs. (B) Numbers of upregulated and downregulated differentially expressed mRNAs. (C) Volcano plots displaying differentially expressed lncRNAs. (D) Volcano plots displaying differentially expressed mRNAs. (E) Heatmap of differentially expressed lncRNAs. (F) Heatmap of differentially expressed mRNAs.

Prediction of lncRNA Target Genes and Analysis of lncRNA-mRNA Interaction

The 254 differentially expressed lncRNAs were predicted for cis- and trans-target genes (Table 2). We found 20 pairs of cis-regulatory relationships between 15 differentially expressed lncRNAs and 14 protein-coding genes, of which, three lncRNAs were upstream of the target gene, and 17 lncRNAs were downstream. There were also two pairs of trans-regulatory relationships between the two differentially expressed lncRNAs and the two protein-coding genes. Among these target genes, SH2D7 (SH2 domain containing) and IL12A (interleukin-12 A) are involved in immune response, LPAR1 (lysophosphatidic acid receptor 1) and FKBP10 (FK506 binding protein 10) are involved in signal transduction, STEAP4 (six-transmembrane epithelial antigen of prostate 4) is involved in oxidative stress, KCNG3 (voltage-gated channel subfamily G member 3) is involved in ion transport, and PHACTR1 (phosphatase and actin regulator 1) is involved in pulse pressure regulation.

TABLE 2

Table 2. Differentially expressed lncRNAs and their targeted mRNAs.

We constructed an interaction network based on these differentially expressed lncRNAs and mRNAs, and 17 differentially expressed lncRNAs and 15 differentially expressed mRNAs had an interaction relationship (Figure 3). Among the differentially expressed lncRNAs, MSTRG.10206.4, MSTRG.10206.3, XR_001435094.1, MSTRG.10206.1 and XR_001022162.2 interacted with two mRNAs, and the rest of the lncRNAs interacted with only one mRNA. Among the differentially expressed mRNAs, MSTRG.10204 interacted with four lncRNAs, MSTRG.10208 interacted with four lncRNAs, and the remaining mRNAs interacted with one lncRNA. Among these, there were more cis-interaction relationships than trans-interaction relations.

FIGURE 3

Figure 3. Co-expression network of differentially expressed lncRNAs and differentially expressed mRNAs.

GO and KEGG Enrichment Analyses

The target genes of the differentially expressed lncRNAs were significantly associated with 38 GO terms, most of which were related to material metabolic processes, including ferric iron transport, trivalent inorganic cation transport, and sulfate transport (Figure 4A; Supplementary Table 4). In KEGG analysis, there were 18 significantly enriched pathways, including the RIG-I-like receptor signaling pathway, taste transduction, and the Toll-like receptor signaling pathway (Figure 4B; Supplementary Table 5). Differentially expressed mRNA was significantly enriched in 215 GO terms, most of which were related to immune function and stress response (Figure 4C; Supplementary Table 6). KEGG analysis was significantly enriched in 37 pathways, including the PPAR signaling pathway involved in fat metabolism and the NF-κβ, NOD-like receptor and Toll-like receptor signaling pathways involved in immune recognition and response (Figure 4D; Supplementary Table 7).

FIGURE 4

Figure 4. GO and KEGG enrichment analyses of differentially expressed lncRNAs and mRNAs. (A) Annotated GO terms of target genes of differentially expressed lncRNAs. (B) Enriched KEGG pathways of target genes of differentially expressed lncRNAs. (C) Annotated GO terms of differentially expressed mRNAs. (D) Enriched KEGG pathways of differentially expressed mRNAs.

qRT-PCR Validation of Differential lncRNA and mRNA Expression

Six differentially expressed lncRNAs and mRNAs were selected for qRT-PCR verification (Figure 5). The qRT-PCR and RNA-seq results showed similar expression trends, confirming the reliability of the transcriptomic sequencing data.

FIGURE 5

Figure 5. Verification of differentially expressed lncRNAs and mRNAs by qRT-PCR.

Discussion

Tibetan sheep have undergone natural selection for tens of thousands of years in the low-oxygen environment of the Qinghai-Tibet Plateau. Specifically, they have developed cardiopulmonary functions and strong low-oxygen adaptability. Under hypoxic stress, different areas of the body show different response mechanisms, forming a unique hypoxia adaptation strategy. In-depth study of the genetic mechanisms of Tibetan sheep in response to high-altitude hypoxic environmental stress can provide a theoretical basis for solving problems with animals entering Tibet, adapting to the high-altitude hypoxic environment and maintaining normal production performance in plains areas. Studying these animals can also help reveal their excellent genetic resources. Thus, we analyzed the complex life processes of Tibetan sheep and the phenotype that allows them to adapt to high-altitude hypoxic conditions at the transcriptomic level. Screening and regulating the genes related to hypoxia adaptation in Tibetan sheep will enable constructing a transcriptional expression regulatory network to analyze the genetic mechanisms of hypoxia adaptation in Tibetan sheep. GO and KEGG analyses of differentially expressed genes revealed the significant enrichment of pathways such as energy metabolism, immune response, oxidative stress, digestion and metabolism, and body temperature regulation. Strong selection on these pathways has also been performed in the study of adaptability of Tibetan pigs and yaks to high-altitude hypoxia (17–19).

Studies have shown that high-altitude hypoxic environments have an important impact on energy metabolism in animals (15, 20, 21). Similarly, we found many differentially expressed genes and signaling pathways related to glucose and lipid metabolism. Among the most significant top 20 differentially expressed genes, five were related to lipid metabolism (PDZK1, APOB, AHSG, DAXX, and C4BPA) (22–26). Additionally, in the KEGG enrichment analysis, the PPAR signaling pathway related to fat metabolism was significantly enriched. Studies on the high-altitude adaptability of Tibetan pigs have also found the significant enrichment of energy metabolism pathways. Fat is the main form of energy storage in animals' bodies. Under hypoxic conditions, the body's fat metabolism changes significantly, and the body intentionally increases its lipid oxidation and phosphorylated ATP synthesis levels with increased heat production to adapt to alpine environments (27–29). Similarly, some genes involved in fat metabolism, such as LPAR1, is also found in lncRNA target genes (30). The PPAR gene is a candidate gene for adaptation to high-altitude hypoxia. It participates in and is regulated by the HIF pathway and is closely related to the production of ATP in animals (31, 32). The activated PPAR gene can enhance fatty acid oxidation and upregulate the mitochondrial β-oxidation process. Therefore, under hypoxic conditions, Tibetan sheep can adapt to the plateau environment by increasing energy metabolism. In addition, there are certain differences in the structure of alveoli among animals inhabiting different altitudes. Type II alveoli are essential for normal lung function and regeneration after hypoxic injury. Respiration depends on surfactants produced and secreted by type II alveoli, mainly composed of phospholipids (33). Some FGF family members (FGF11, FGF14, and FGF4) were included among the differentially expressed genes. Studies have shown that FGF can stimulate the proliferation of type II alveolar cells, which in turn affects lipid homeostasis on the alveolar surface (34). This may also be caused by the enrichment of differential genes in fat metabolism.

The immune system is the body's defense barrier, and exposure to high-altitude and hypoxic environments can change the body's immune functions (35). At an altitude of 3,000 m, the number of circulating dendritic cells in human plasma is significantly reduced (36). At an altitude of 5,000 m, plasmacytoid dendritic cells in human plasma are reduced, and the TNF-α and IL-6 contents are significantly increased (37). We found that several genes, including NKIRAS2, IKBKE, and TRIM7, were involved in immunoregulation and were lncRNA target genes (38–40). Additionally, some immune-related signaling pathways, such as the NF-κβ, Toll-like receptor, and B-cell-receptor signaling pathways, were also enriched. NF-κβ is an important immune response regulator involved in innate immunity and adaptive immune responses. The differentially expressed genes identified in this study, including NKIRAS2, IKBKE (inhibitor of nuclear factor kappa-B kinase subunit epsilon) and TRIM7 (tripartite motif containing 7), are key factors in the NF-κβ signaling pathway. After recognizing pathogenic microorganisms, Toll-like receptors activate NF-κβ through the MyD88 or TRIF-dependent signaling pathways, then initiate innate and adaptive immunity (41, 42). After activation, NF-κB regulates the expression of cytokines such as IL-1, IL-6, and IL-8 and releases them outside the cell to exert an early immune response effect (43, 44). This study found that the expression levels of genes such as NF-κB2, TLR2 (toll-like receptor 2), IL1A, IL1R2, and IL18 were significantly lower in Tibetan sheep than in Hu sheep. The downregulation of TLR2 expression in the lungs of Tibetan sheep may reduce the binding of IL1A, IL1R2, and IL18, thereby reducing downstream TNF-α activity (45). Decreased TNF-α activity may weaken NF-κB-induced kinase activity (NIK). Correspondingly, the decrease in NIK activity may lead to a decrease in the kinase activity inhibited by NF-κB, and the inhibition of NF-κB kinase activity leads to a decrease in the activity of the I-κB kinase (IKK) complex (46). This cascade reaction attenuates the phosphorylation of NF-κB complex inhibitor I-κB, thereby slowing its dissociation from NF-κB and reducing the level of NF-κB, which may also cause the downregulation of TNF-α expression in the lungs of Tibetan sheep (47). Our results indicate that the NF-κβ signaling pathway and its key factors play important roles in the immunoregulation of Tibetan sheep.

Animals living in high-altitude areas for a long time, under low pressure, low oxygen, low temperature and strong ultraviolet stimulation, will experience varying degrees of oxidative stress and produce large amounts of reactive oxygen and nitrogen species (48, 49). The selenoprotein, SEPW1, has glutathione-dependent antioxidant activity (50). GPX1 (glutathione peroxidase-1) is a selenoprotein that can oxidize glutathione, remove excess hydroxyl free radicals and peroxides, and protect cells from oxidative stress (51). The lncRNA target gene STEAP4 is an iron-copper oxidoreductase that uses the reducing power provided by NADPH to reduce Fe³⁺ and Cu²⁺ to Fe²⁺ and Cu⁺, which is essential for maintaining many cellular processes (52). Tibetan sheep have been living in a hypoxic environment for many years, necessitating physiological and genetic adaptation to oxidative stress. Therefore, compared with Tibetan sheep, Hu sheep are more sensitive to oxidative stress and have higher antioxidant activity. Compared with the level in Hu sheep, GPX1 is expressed at a significantly lower level in Tibetan sheep, which is consistent with the lower GPX activity of Tibetans living on plateaus (53). In contrast, we found that the expression of CAT, a gene related to antioxidant enzyme activity, was significantly higher in Tibetan sheep than in Hu sheep.

Compared with Hu sheep, Tibetan sheep have always been under traditional grazing management and have a strong ability to adapt to the low-oxygen environment of the plateau. The Qinghai-Tibet Plateau is extremely deficient in forage in winter and spring, and Tibetan sheep can maintain normal reproduction even if their nutritional intake is severely insufficient, indicating that Tibetan sheep have developed unique digestive and metabolic functions. We found two related genes involved in bile acid metabolism, SLC26A10 and IL12A, among the lncRNA target genes. SLC26A10, one of the main carriers involved in the enterohepatic circulation of bile acids, ingests plasma-bound bile salts into liver cells in a sodium-dependent manner (54, 55). IL12A has an immunomodulatory effect, is significantly involved in the susceptibility to primary biliary cholangitis, and can be used as a molecular target for clinical diagnoses (56, 57). Bile acid is an important component of bile. A compound with an amphiphilic molecular structure produced by cholesterol metabolism in the liver plays an important role in lipid digestion, absorption and metabolism. This suggests that Tibetan sheep can digest and absorb lipids rationally through bile acids as a defense against the harm caused by the high-altitude hypoxic environment.

Tibetan sheep have strong cold tolerance and can grow and develop normally on the Qinghai-Tibet Plateau more than 3,000 m above sea level. A series of genes involved in the regulation of body temperature were discovered in this study, such as UCP3 (uncoupling protein 3) and HTR4 (5-hydroxytryptamine receptor 4). Studies have shown that UCP3 can significantly increase the oxidative respiration rate of Tibetan pig fat cells, indicating that it is an important gene in enabling Tibetan pigs to resist cold and heat (58). Another study also showed that the cold-resistance mechanism of pigs does not simply involve muscle tremor and fever, but activation of the UCP3 protein (59). By upregulating the expression of the UCP3 gene, Tibetan pigs promote the browning of subcutaneous white fat, increase fat production, and maintain body temperature balance. Moreover, the expression level of UCP3 in Tibetan sheep was shown to be significantly higher than that in Hu sheep, indicating that Tibetan sheep can also mediate subcutaneous fat browning and increase the body's heat production through UCP3. Studies have shown that 5-HT participates in the regulation of animal body temperature (60). When an animal is in a high-temperature environment, the 5-HT neuroendocrine system is activated. 5-HT binds to 5-HTR and activates the downstream cAMP and cGMP signaling pathways, resulting in increased activity of heat-sensitive nerves and then the inhibition of heat production (61). In addition to participating in body temperature regulation, HTR4 also plays an important role in animal gastrointestinal sensitivity and food intake (62). In contrast to the case of Hu sheep, the Qinghai-Tibet Plateau forages that Tibetan sheep live on are characterized by high fiber content, meaning that these animals have a long history of nutritional stress. These factors affect the food intake and ruminant activities of Tibetan sheep, and then affect the calories produced by metabolism to maintain body temperature stability. During this period, HTR4 may have played an important regulatory role.

Conclusion

Here, we systematically identified the expression profiles of lncRNA and mRNA involved in the adaptation process of Tibetan sheep to hypoxia. Functional enrichment and interaction network analysis results showed that lncRNA and mRNA may participate in the adaptation of Tibetan sheep to hypoxia via lipid metabolism. These results provide valuable resources for studying lncRNA and mRNA involved in the adaptation of animals to plateau hypoxia and can help clarify the molecular mechanisms of the adaptation of animals to plateaus.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee of Lanzhou Institute of Husbandry and Pharmaceutical Science of Chinese Academy of Agricultural Sciences (Approval No. NKMYD201805; Approval Date: 18 October 2018).

Author Contributions

ZL, JLiu, and BY conceived and designed the study. CY, JLi, TG, CN, and YY collected the samples. ZL performed the experiments, analyzed the data, and wrote the paper. ZL, JLi, and BY contributed to revisions of the manuscript. All authors read and approved the manuscript.

Funding

This research work was supported by the Chinese Academy of Agricultural Sciences of Technology Innovation Project (25-LZIHPS-07 and CAAS-ZDRW202106), the National Natural Science Foundation of China (32002170), and the Special Fund of the Chinese Central Government for Basic Scientific Research Operations in Commonweal Research Institutes (1610322020004).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the Guangzhou Genedenovo Biotechnology Co., Ltd. for assisting in RNA sequencing.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2021.801278/full#supplementary-material

References

1. York JM, Scadeng M, McCracken KG, Milsom WK. Respiratory mechanics and morphology of Tibetan and Andean high-altitude geese with divergent life histories. J Exp Biol. (2018) 221:jeb170738. doi: 10.1242/jeb.170738

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Storz JF, Runck AM, Sabatino SJ, Kelly JK, Ferrand N, Moriyama H, et al. Evolutionary and functional insights into the mechanism underlying high-altitude adaptation of deer mouse hemoglobin. Proc Natl Acad Sci USA. (2009) 106:14450–5. doi: 10.1073/pnas.0905224106

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yang J, Jin Z-B, Chen J, Huang X-F, Li X-M, Liang Y-B, et al. Genetic signatures of high-altitude adaptation in Tibetans. Proc Natl Acad Sci USA. (2017) 114:4189–94. doi: 10.1073/pnas.1617042114

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Simonson TS, Yang Y, Huff CD, Yun H, Qin G, Witherspoon DJ, et al. Genetic evidence for high-altitude adaptation in Tibet. Science. (2010) 329:72–5. doi: 10.1126/science.1189406

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Liu X, Zhang Y, Li Y, Pan J, Wang D, Chen W, et al. EPAS1 gain-of-function mutation contributes to high-altitude adaptation in Tibetan horses. Mol Biol Evol. (2019) 36:2591–603. doi: 10.1093/molbev/msz158

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lorenzo FR, Huff C, Myllymäki M, Olenchock B, Swierczek S, Tashi T, et al. A genetic mechanism for Tibetan high-altitude adaptation. Nat Genet. (2014) 46:951–6. doi: 10.1038/ng.3067

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cai W, Liu S, Liu Z, Hou S, Lv Q, Cui H, et al. Downregulation of lung miR-203a-3p expression by high-altitude hypoxia enhances VEGF/Notch signaling. Aging. (2020) 12:4247–67. doi: 10.18632/aging.102878

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Zhu J, Fu H, Wu Y, Zheng X. Function of lncRNAs and approaches to lncRNA-protein interactions. Sci China Life Sci. (2013) 56:876–85. doi: 10.1007/s11427-013-4553-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Mercer TR, Mattick JS. Structure and function of long noncoding RNAs in epigenetic regulation. Nat Struct Mol Biol. (2013) 20:300–7. doi: 10.1038/nsmb.2480

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liu M, Zhong J, Zeng Z, Huang K, Ye Z, Deng S, et al. Hypoxia-induced feedback of HIF-1α and lncRNA-CF129 contributes to pancreatic cancer progression through stabilization of p53 protein. Theranostics. (2019) 9:4795–810. doi: 10.7150/thno.30988

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ma C-N, Wo L-L, Wang D-F, Zhou C-X, Li J-C, Zhang X, et al. Hypoxia activated long non-coding RNA HABON regulates the growth and proliferation of hepatocarcinoma cells by binding to and antagonizing HIF-1 alpha. RNA Biol. (2021) 18:1791–806. doi: 10.1080/15476286.2020.1871215

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hua Q, Mi B, Xu F, Wen J, Zhao L, Liu J, et al. Hypoxia-induced lncRNA-AC020978 promotes proliferation and glycolytic metabolism of non-small cell lung cancer by regulating PKM2/HIF-1α axis. Theranostics. (2020) 10:4762–78. doi: 10.7150/thno.43839

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Zhang Z, Fang E, Rong Y, Han H, Gong Q, Xiao Y, et al. Hypoxia-induced lncRNA CASC9 enhances glycolysis and the epithelial-mesenchymal transition of pancreatic cancer by a positive feedback loop with AKT/HIF-1α signaling. Am J Cancer Res. (2021) 11:123–37.

Google Scholar

14. Wang X, Wang Y, Li L, Xue X, Xie H, Shi H, et al. A lncRNA coordinates with Ezh2 to inhibit HIF-1α transcription and suppress cancer cell adaption to hypoxia. Oncogene. (2020) 39:1860–74. doi: 10.1038/s41388-019-1123-9

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang Y, Su W, Zhang B, Ling Y, Kim WK, Zhang H. Comprehensive analysis of coding and non-coding RNA transcriptomes related to hypoxic adaptation in Tibetan chickens. J Anim Sci Biotechnol. (2021) 12:60. doi: 10.1186/s40104-021-00582-2

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Wang J, Chai Z, Deng L, Wang J, Wang H, Tang Y, et al. Detection and integrated analysis of lncRNA and mRNA relevant to plateau adaptation of Yak. Reprod Domest Anim. (2020) 55:1461–9. doi: 10.1111/rda.13767

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhang B, Chamba Y, Shang P, Wang Z, Ma J, Wang L, et al. Comparative transcriptomic and proteomic analyses provide insights into the key genes involved in high-altitude adaptation in the Tibetan pig. Sci Rep. (2017) 7:3654. doi: 10.1038/s41598-017-03976-3

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Wang T, Guo Y, Liu S, Zhang C, Cui T, Ding K, et al. KLF4, a key regulator of a transitive triplet, acts on the TGF-β signaling pathway and contributes to high-altitude adaptation of Tibetan pigs. Front Genet. (2021) 12:628192. doi: 10.3389/fgene.2021.628192

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ge Q, Guo Y, Zheng W, Zhao S, Cai Y, Qi X. Molecular mechanisms detected in yak lung tissue via transcriptome-wide analysis provide insights into adaptation to high altitudes. Sci Rep. (2021) 11:7786. doi: 10.1038/s41598-021-87420-7

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Feng S, Ma J, Long K, Zhang J, Qiu W, Li Y, et al. Comparative microRNA transcriptomes in domestic goats reveal acclimatization to high altitude. Front Genet. (2020) 11:809. doi: 10.3389/fgene.2020.00809

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Xin JW, Chai ZX, Zhang CF, Yang YM, Zhang Q, Zhu Y, et al. Transcriptome analysis identified long non-coding RNAs involved in the adaption of yak to high-altitude environments. R Soc Open Sci. (2020) 7:200625. doi: 10.1098/rsos.200625

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Trigatti BL. SR-B1 and PDZK1: partners in HDL regulation. Curr Opin Lipidol. (2017) 28:201–8. doi: 10.1097/MOL.0000000000000396

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Sirwi A, Hussain MM. Lipid transfer proteins in the assembly of apoB-containing lipoproteins. J Lipid Res. (2018) 59:1094–102. doi: 10.1194/jlr.R083451

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Stefan N, Hennige AM, Staiger H, Machann J, Schick F, Kröber SM, et al. Alpha2-heremans-schmid glycoprotein/fetuin-A is associated with insulin resistance and fat accumulation in the liver in humans. Diabetes Care. (2006) 29:853–7. doi: 10.2337/diacare.29.04.06.dc05-1938

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Li TP, Sun SW, Xiong GZ, Qiu F, Yang DM, Sun SY, et al. Direct interaction of daxx and androgen receptor is required for their regulatory activity in cholesterol biosynthesis. Pharmacology. (2021) 106:29–36. doi: 10.1159/000506488

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jiang P, Wang XN, Wang Y, Li AN, Xiao H, Guo PC, et al. Difference analysis of C4BPA gene expression in mammary tissue of dairy cows. Chin J Vet Sci Jum. (2016) 36:1032–5+43. doi: 10.16303/j.cnki.1005-4545.2016.06.28

CrossRef Full Text | Google Scholar

27. Cheviron ZA, Bachman GC, Connaty AD, McClelland GB, Storz JF. Regulatory changes contribute to the adaptive enhancement of thermogenic capacity in high-altitude deer mice. Proc Natl Acad Sci USA. (2012) 109:8635–40. doi: 10.1073/pnas.1120523109

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gangwar A, Paul S, Ahmad Y, Bhargava K. Intermittent hypoxia modulates redox homeostasis, lipid metabolism associated inflammatory processes and redox post-translational modifications: benefits at high altitude. Sci Rep. (2020) 10:7899. doi: 10.1038/s41598-020-64848-x

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Meir JU, York JM, Chua BA, Jardine W, Hawkes LA, Milsom WK. Reduced metabolism supports hypoxic flight in the high-flying bar-headed goose. Elife. (2019) 8:e44986. doi: 10.7554/eLife.44986

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang X, Li M, Yin N, Zhang J. The expression regulation and biological function of autotaxin. Cells. (2021) 10:939. doi: 10.3390/cells10040939

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Finck BN, Bernal-Mizrachi C, Han DH, Coleman T, Sambandam N, LaRiviere LL, et al. A potential link between muscle peroxisome proliferator- activated receptor-alpha signaling and obesity-related diabetes. Cell Metab. (2005) 1:133–44. doi: 10.1016/j.cmet.2005.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Lefebvre P, Chinetti G, Fruchart JC, Staels B. Sorting out the roles of PPAR alpha in energy metabolism and vascular homeostasis. J Clin Investig. (2006) 116:571–80. doi: 10.1172/JCI27989

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Sever N, Miličić G, Bodnar NO, Wu X, Rapoport TA. Mechanism of lamellar body formation by lung surfactant protein B. Molecular Cell. (2021) 81:49–66. doi: 10.1016/j.molcel.2020.10.042

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Mason RJ. Biology of alveolar type II cells. Respirology. (2006) 11:12–5. doi: 10.1111/j.1440-1843.2006.00800.x

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Feuerecker M, Crucian BE, Quintens R, Buchheim J-I, Salam AP, Rybka A, et al. Immune sensitization during 1 year in the Antarctic high-altitude Concordia environment. Allergy. (2019) 74:64–77. doi: 10.1111/all.13545

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Rohm I, Aderhold N, Ratka J, Goebel B, Franz M, Pistulli R, et al. Hypobaric hypoxia in 3000 m altitude leads to a significant decrease in circulating plasmacytoid dendritic cells in humans. Clin Hemorheol Microcirc. (2016) 63:257–65. doi: 10.3233/CH-152035

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Yilmaz A, Ratka J, Rohm I, Pistulli R, Goebel B, Asadi Y, et al. Decrease in circulating plasmacytoid dendritic cells during short-term systemic normobaric hypoxia. Eur J Clin Invest. (2016) 46:115–22. doi: 10.1111/eci.12416

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Sarais F, Rebl H, Verleih M, Ostermann S, Krasnov A, Köllner B, et al. Characterisation of the teleostean κB-Ras family: the two members NKIRAS1 and NKIRAS2 from rainbow trout influence the activity of NF-κB in opposite ways. Fish Shellfish Immunol. (2020) 106:1004–13. doi: 10.1016/j.fsi.2020.08.052

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li Q, Verma IM. NF-kappaB regulation in the immune system. Nat Rev Immunol. (2002) 2:725–34. doi: 10.1038/nri910

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Lu M, Zhu X, Yang Z, Zhang W, Sun Z, Ji Q, et al. E3 ubiquitin ligase tripartite motif 7 positively regulates the TLR4-mediated immune response via its E3 ligase domain in macrophages. Mol Immunol. (2019) 109:126–33. doi: 10.1016/j.molimm.2019.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Kawai T, Akira S. Signaling to NF-kappaB by Toll-like receptors. Trends Mol Med. (2007) 13:460–9. doi: 10.1016/j.molmed.2007.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Anthoney N, Foldi I, Hidalgo A. Toll and Toll-like receptor signalling in development. Development. (2018) 145:156018. doi: 10.1242/dev.156018

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Fields JK, Günther S, Sundberg EJ. Structural basis of IL-1 family cytokine signaling. Front Immunol. (2019) 10:1412. doi: 10.3389/fimmu.2019.01412

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Qin Y, Li H, Qiao J. TLR2/MyD88/NF-κB signalling pathway regulates IL-8 production in porcine alveolar macrophages infected with porcine circovirus 2. J Gen Virol. (2016) 97:445–52. doi: 10.1099/jgv.0.000345

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Lu YC, Yeh WC, Ohashi PS. LPS/TLR4 signal transduction pathway. Cytokine. (2008) 42:145–51. doi: 10.1016/j.cyto.2008.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Kawai T, Akira S. TLR signaling. Semin Immunol. (2007) 19:24–32. doi: 10.1016/j.smim.2006.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Takeda K, Akira S. TLR signaling pathways. Semin Immunol. (2004) 16:3–9. doi: 10.1016/j.smim.2003.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gaur P, Prasad S, Kumar B, Sharma SK, Vats P. High-altitude hypoxia induced reactive oxygen species generation, signaling, and mitigation approaches. Int J Biometeorol. (2021) 65:601–15. doi: 10.1007/s00484-020-02037-1

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Dosek A, Ohno H, Acs Z, Taylor AW, Radak Z. High altitude and oxidative stress. Respir Physiol Neurobiol. (2007) 158:128–31. doi: 10.1016/j.resp.2007.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Jeong Dw, Kim TS, Chung YW, Lee BJ, Kim IY. Selenoprotein W is a glutathione-dependent antioxidant in vivo. FEBS Lett. (2002) 517:225–8. doi: 10.1016/S0014-5793(02)02628-5

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Lei XG, Cheng W-H, McClung JP. Metabolic regulation and function of glutathione peroxidase-1. Annu Rev Nutr. (2007) 27:41–61. doi: 10.1146/annurev.nutr.27.061406.093716

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Scarl RT, Lawrence CM, Gordon HM, Nunemaker CS. STEAP4: its emerging role in metabolism and homeostasis of cellular iron and copper. J Endocrinol. (2017) 234:123–34. doi: 10.1530/JOE-16-0594

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Imai H, Kashiwazaki H, Suzuki T, Kabuto M, Himeno S, Watanabe C, et al. Selenium levels and glutathione peroxidase activities in blood in an andean high-altitude population. J Nutr Sci Vitaminol. (1995) 41:349–61. doi: 10.3177/jnsv.41.349

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Jani M, Beéry E, Heslop T, Tóth B, Jagota B, Kis E, et al. Kinetic characterization of bile salt transport by human NTCP (SLC10A1). Toxicol In Vitro. (2018) 46:189–93. doi: 10.1016/j.tiv.2017.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Claro da Silva T, Polli JE, Swaan PW. The solute carrier family 10 (SLC10): beyond bile acid transport. Mol Aspects Med. (2013) 34:252–69. doi: 10.1016/j.mam.2012.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Hirschfield GM, Liu X, Xu C, Lu Y, Xie G, Lu Y, et al. Primary biliary cirrhosis associated with HLA, IL12A, and IL12RB2 variants. N Engl J Med. (2009) 360:2544–55. doi: 10.1056/NEJMoa0810440

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Carbone M, Mells GF, Alexander GJ, Westbrook RH, Heneghan MA, Sandford RN, et al. Calcineurin inhibitors and the IL12A locus influence risk of recurrent primary biliary cirrhosis after liver transplantation. Am J Transplant. (2013) 13:1110–1. doi: 10.1111/ajt.12132

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Lin J, Cao C, Tao C, Ye R, Dong M, Zheng Q, et al. Cold adaptation in pigs depends on UCP3 in beige adipocytes. J Mol Cell Biol. (2017) 9:364–75. doi: 10.1093/jmcb/mjx018

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Li M, Tian S, Jin L, Zhou G, Li Y, Zhang Y, et al. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars. Nat Genet. (2013) 45:1431–8. doi: 10.1038/ng.2811

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Cox B, Lee TF. Further evidence for a physiological role for hypothalamic dopamine in thermoregulation in the rat. J Physiol. (1980) 300:7–17. doi: 10.1113/jphysiol.1980.sp013147

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Boulant JA. Role of the preoptic-anterior hypothalamus in thermoregulation and fever. Clin Infect Dis. (2000) 31 (Suppl. 5):157–61. doi: 10.1086/317521

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Lucas G. Serotonin receptors, type 4: a new hope? Curr Drug Targets. (2009) 10:1085–95. doi: 10.2174/138945009789735200

PubMed Abstract | CrossRef Full Text | Google Scholar