The procedures for animal care were approved by the Animal Welfare Committee of the Tibet Agriculture and Animal Husbandry College, and all experiments were conducted in accordance with approved relevant guidelines and regulations. Tibetan and Yorkshire pigs were raised under identical conditions with same diet at the farm of the Tibet Agriculture and Animal Husbandry College, located in Linzhi city (Tibet, 2,900 m, above mean sea level). All pigs were housed in standard conditions with natural, uncontrolled room temperature and light. Complete formula meal feed was fed three times per day and pigs had ad libitum access to water. At the age of 6 months, three Tibetan pigs and three Yorkshire pigs were randomly selected to slaughter and sample. Liver tissues from each individual were collected, immediately frozen in liquid nitrogen, and stored at −80°C until use.

The total RNA from liver tissues was extracted following the phenol-chloroform method. RNA concentration and integrity were determined using an Agilent 2100 Bioanalyzer and the Agilent RNA 6000 Nano Kit (Agilent Technologies, Palo Alto, CA, USA). RNA-seq libraries were constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (# E7530L, NEB; San Diego, CA, USA) according to the manufacturer’s instructions. Paired-end 150 bp reads were generated using the HiSeq Xten platform.

Data analysis included quality control of raw reads, filtering, alignment, assembly, expression count, and annotation. Contaminated reads were filtered using a Perl script in-house. Clean reads were mapped to the reference genome (Sus scrofa genome v11.1, downloaded from ENSEMBL web server) using HISAT2 v2.1.0 [11]. Reads for each gene in each sample were counted by HTSeq v0.6.0 (https://github.com/simon-anders/htseq). Fragments per kilobase per million mapped reads were then calculated to estimate the expression level of each gene. Deseq (http://bioconductor.org/packages/release/bioc/html/DESeq.html) was used to identify the differentially expressed genes (DEGs), with the criteria of fold change >2 and adjusted p-value <0.05 [12].

Differentially expressed mRNAs were selected for validation by quantitative-polymerase chain reaction (qPCR). Total RNA was extracted from liver tissue using TRIzol reagent, according to the manufacturer’s instructions (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA synthesis was performed with 1 μg of total RNA, following the protocol accompanying the FastQuant RT Kit (Tiangen Biotech (Beijing) Co., LTD, Beijing, China). The qPCR amplification was performed as described previously [1]. Relative expression is presented as mean±standard error of the mean. Differences were tested for statistical significance using the Student’s t-test. A p-value <0.05 was considered as the threshold for statistical significance (* p<0.05; ** p<0.01).

Liver tissue protein was extracted using a mammalian tissue total protein extraction kit (AP0601-50, Beijing Bangfei Biotechnology Co., Ltd., Beijing, China). Quantification was investigated using a protein quantification kit (Dingguo Changsheng, Beijing, China). Protein peptides from each group were labeled using the 8-plex isobaric tags for relative and absolute quantitation (iTRAQ) reagents multiplex kit (ABI, Foster City, CA, USA). For each sample, approximately 200 μg protein was precipitated with 10 volumes of acetone at −20°C overnight. After centrifugation, the protein pellet was dissolved in 60 μL of dissolution buffer. The iTRAQ labeling reagents were added to the peptide samples and reacted at room temperature for 1 h, before investigation of the labeling and extraction efficiency. Samples were pooled, vacuum-dried, and dissolved in solution. The reconstituted peptides were analyzed with a Q-Exactive mass spectrometer (Thermo Fisher Scientific, USA), coupled with a nano high-performance liquid chromatography system (Thermo Fisher Scientific, USA). The detailed steps for enzymolysis and iTRAQ labeling have been described previously [5,9].

Data analysis was performed using Proteome Discover v2.1 (Thermo Fisher Scientific, USA) software. Peptide identification was performed using the SEQUEST search engine, using the Uniprot database (Sus scrofa 49003 entries, 20190102). Proteins were identified according to the following criteria: mass error was set to 15 ppm for precursor ions and 0.02 Da for fragment ions. Trypsin was used as a cleavage enzyme, and two missing cleavages were allowed. The false discovery rate was adjusted to 0.01. Differentially expressed proteins (DEPs) were identified as having a fold change >1.2 and p-value <0.05.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) functional enrichment analyses were performed, using the g:Profiler web server [13]. Visualization of enrichment results was performed using the R packages ‘ggplot2’, ‘forcats’, ‘enrichplot’, ‘wordcloud2’, and ‘GOplot’.

After removing low quality reads, an average of 45.86 million paired-end reads were obtained for each sample. The mapping rate aligning to Sus scrofa 11.1 of all samples was above 96.2%. Reads located on exons accounted for more than 68.71% (Supplementary Table S1). A total of 18,791 genes were identified across the six libraries. A number of genes highly expressed in liver tissue, such as albumin, fibrinogen beta chain precursor, cytochrome c oxidase subunit I, apolipoprotein E precursor, catalase, cytochrome P450 2E1, ATP synthase F0 subunit 6, vitronectin, apolipoprotein A-I (APOA1), and 17-beta-hydroxysteroid dehydrogenase 13, may play important roles in determining the metabolic pattern of liver tissue. Tibetan pigs have adapted to a high-altitude environment, and their livers have adapted to function under hypoxia and low temperatures. The current study aimed to characterize the molecular adaptation pattern by comparing the DEGs between Tibetan and Yorkshire (low-land breed) pigs; after a series of filtering steps, 273 DEGs were identified, which included 187 upregulated and 86 downregulated genes in the Tibetan group relative to the Yorkshire group (Figure 1A; Supplementary Table S2).

To further understand the potential molecular mechanisms involved in high-altitude adaptation, DEGs were subjected to functional enrichment analysis (Supplementary Table S3). Relevant GO classifications were identified (Figure 1B). These included the response to oxygen-containing compounds, the small molecule biosynthetic or metabolic process, the lipid metabolic process, the drug metabolic process, the organic hydroxyl compound metabolic process, the steroid biosynthetic or metabolic process, and the cholesterol biosynthetic process. Only four KEGG pathways were enriched, which consisted of metabolic pathways, the peroxisome proliferator-activated receptor signaling pathway, steroid biosynthesis, and terpenoid backbone biosynthesis (Figure 1C; Supplementary Table S3).

A total of 19,622 peptides and 3,781 proteins were identified, of which 3,390 proteins were quantified. Protein molecular weights ranged from 20 kDa to 50 kDa. Peptide coverage, peptide count, and isoelectric points for these proteins are presented in Supplementary Figure S1. To better understand the function of these expressed protein, proteins were subjected to GO enrichment analysis. The top 10 enriched GO classifications are shown in Supplementary Figure S2, grouped according to the major GO categories of biological process, cellular component, and molecular function. For the biological process category, most proteins participated in the metabolic process (GO:0008152), the organic substance metabolic process (GO:0071704), the cellular metabolic process (GO:0044237), and the primary metabolic process (GO:0044238). For the molecular function category, most proteins were involved in binding (GO:0005488), catalytic activity (GO: 0003824), protein binding (GO:0005515), ion binding (0043167), and organic cyclic compound binding (GO:0097159). Finally, for the cellular component category, most proteins were found in the cell (GO:0005623), the cell parts (GO:0044464), the intracellular parts (GO:0005622), the cytoplasm (GO:0005737), and in intracellular organelles (GO:0043229) (Supplementary Figure S2).

Using 1.20/0.833 and a p-value <0.05 as a threshold to classify proteins as upregulated or downregulated, a total of 257 DEPs were identified in Tibetan pigs relative to Yorkshire pigs. These proteins consisted of 38 upregulated proteins and 219 downregulated proteins (Figure 2A; Supplementary Table S4). The most highly expressed proteins in the Tibetan group were phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 gamma (PIK3C2G), myoglobin (MB), carbonyl reductase (CBR2), solute carrier family 9 member C1 (SLC9C1), PDZ and LIM domain 5 (PDLIM5), MYL3, LOC100736765 (Myosin-7), troponin C (TNNC1), and TNNT2 (Supplementary Table S4). The most highly expressed proteins in the Yorkshire group were ADH4, CAD, LOC100739663 (metallothionein), metallothionein-1D (MT1D), keratin 75 (KRT75), LOC110260333 (3-hydroxy-3-methylglutaryl coenzyme A synthase), keratin 5 (KRT5), and fatty acid synthase (FASN) (Supplementary Table S4).

To further characterize protein expression patterns, DEPs were subjected to functional enrichment analysis. Relevant terms associated with GO biological processes and KEGG pathways are presented in Figures 2B and 2C, respectively, and the detailed information is presented in Supplementary Table S5. Relevant terms were mainly divided into metabolic processes, biosynthetic processes, and catabolic processes. The biosynthetic processes consist of small molecule biosynthetic and cofactor biosynthetic processes. These metabolic processes included carboxylic acid, the aromatic amino acid family, fatty acid metabolic, monocarboxylic acid, organic, oxoacid, carboxylic, coenzyme, drug, cellular amino acid, pyruvate, alpha-amino acid, dicarboxylic acid, and glutamine family amino acid processes. The catabolic processes included the carboxylic acid catabolic process, aromatic amino acid family catabolic process, lipid, organic acid, carboxylic acid, and aromatic amino acid processes (Figure 2B). KEGG pathways enriched by DEPs consisted of metabolic pathways, carbon metabolism, glycolysis/gluconeogenesis, biosynthesis of amino acids, the pentose phosphate pathway, tryptophan metabolism, pyruvate metabolism, arginine biosynthesis, the citrate cycle, and proximal tubule bicarbonate reclamation (Figure 2C).

To validate the RNA-Seq results, 12 DEGs were randomly selected for qPCR analysis (Figure 3). Of these, acyl-CoA synthetase medium chain family member 2B (ACSM2B), dimethylarginine dimethylaminohydrolase 1 (DDAH1), adenine phosphoribosyltransferase (APRT), phosphatidylethanolamine binding protein 1, FASN, proteolipid protein 1, sorbin and SH3 domain containing 2 (SORBS2), and nucleobindin 2 were downregulated, and solute carrier organic anion transporter family member 2A1 (SLCO2A1), prenylcysteine oxidase 1 (PCYOX1), perilipin 2 (PLIN2), and ladinin 1 (LAD1) were upregulated. Figure 3 shows the relative expression levels between the two groups, which are consistent with the RNA-Seq data, indicating that identification and abundance estimates of genes are highly reliable in this study.

Cumulatively, 18,791 genes and 3,390 proteins were investigated and compared; 273 DEGs and 257 DEPs, of which 216 proteins were annotated, were identified. Only 11 genes overlapped between the DEGs and DEPs (Table 1). Of these, nine DEGs showed the same expression pattern. The five upregulated genes at both mRNA and protein levels included retinol dehydrogenase 16 (RDH16), APRT, PCYOX1, SORBS2, and ENSSSCG00000036224 (Table 1). RDH16, retinol dehydrogenase 16, has been shown to be reduced by insulin in HepG2 cells, and to regulate energy balance and adiposity [14]. PCYOX1 hydrolyzes the thioether bond of prenylcysteines in the final step of degradation of prenylated proteins and propagates the oxidative burden of low-density lipoproteins [15]. SORBS2, sorbin and SH3 domain-containing 2, could suppress hepatocellular carcinoma tumorigenesis and metastasis [16]. Four downregulated genes at both mRNA and protein levels included PLIN2, LAD1, kynurenine aminotransferase 1 (KYAT1), and DDAH1 (Table 1). PLIN2 has been shown to induce obesity and progressive fatty liver disease via mechanistically distinct hepatocyte and extra-hepatocyte actions [17]. KYAT1 plays a role in protecting against liver damage from thioacetamide [18]. Dimethylarginine dimethylaminohydrolase 1 functions in protecting against hepatic steatosis induced by a high-fat diet and insulin resistance [19]. In addition, four GO enrichment terms overlapped those from the DEGs (Figure 4A) and those from the DEPs (Figure 4B). These four classifications included the small molecule metabolic process, the small molecule biosynthetic process, the drug metabolic process, and the organic hydroxyl compound metabolic process (Figure 4). This shows that biosynthetic and metabolic patterns of small molecular and organic hydroxyl compound metabolic processes were functional in the high-altitude adapted group. Genes involved in these biological processes are presented in Figure 4. Oxygen content is less in a high-altitude environment. Populations in this region have undergone natural selection, leading to improved oxygen delivery against the challenging environmental stress [20]. Growing evidence also suggests that high-altitude adaptation is derived from multiple molecular mechanisms; not only oxygen delivery but also oxygen utilization by cellular metabolism [21,22]. It was reported that decreased fatty acid oxidation could also be a better strategy of hypoxia adaptation since oxidation of fatty acids generates less ATP than carbohydrates per molecule of oxygen consumed [23]. In this study, proteins involved in the fatty acid metabolic process (including 15-hydroxyprostaglandin dehydrogenase, ATP citrate lyase [ACLY], acyl-CoA thioesterase 4, 3-hydroxybutyrate dehydrogenase 2, phosphoenolpyruvate carboxykinase 2 [PCK2], ACSM2B, aldo-keto reductase family 1 member C1, PCK1, apolipoprotein A4 and ethylmalonyl-CoA decarboxylase 1) were lower in Tibetan pigs compared to the Yorkshire pigs. Another change induced by hypoxia could be the conversion of glucose oxidation to glycolysis to decrease oxygen demand and maintain energy production [24,25]. This condition would promote glucose uptake and glycolysis and suppress mitochondrial glucose oxidation. Oxygen demand and utilization was suppressed in the hypoxia adaptation process through a decreasing mitochondrial density and the expression level of tricarboxylic acid cycle (TCA) enzymes [26]. In this study, protein level of five enzymes of the TCA cycle pathway (isocitrate dehydrogenase 1, ACLY, aconitase 1, PCK2, and PCK1) were lower in the Tibetan pig relative to the Yorkshire pig (a low-land breed)(Supplementary Table S5, Table S4). Furthermore, genes involved in glycolysis (glucokinase and phosphomannomutase 1) were upregulated in the Tibetan group relative to the Yorkshire group. Conversely, genes related to the gluconeogenesis process (PCK1, PCK2, and fructose-bisphosphatase 1) were decreased in the Tibetan group. This may imply that glycolysis is more active in the Tibetan pig to reduce oxygen demand and maintain energy production. In this case, it could be deduced that oxygen utilization efficiency was higher in the Tibetan pig than in the low-land pig after long-term, high-altitude adaptation, which mainly depended on reducing fatty acid metabolism and glucose oxidation, and increasing glycolysis. Undoubtedly, these differences between Tibetan and Yorkshire pigs maybe also derived from their different genetic backgrounds. The experimental animals in this study were raised in a high-land region. If the data from the counterparts of the Tibetan and Yorkshire pigs raised in the low altitude area were integrated in this study, key genes could be recognized precisely in the different environments and breeds.

Results of a similar study [2] on the molecular characterization of liver mRNA between Tibetan and Rongchang pigs were compared as in the present study. They compared gene expression levels using RNA-Seq between Tibetan pigs (high-land model) and Rongchang pigs (low-land model) and identified 490 DEGs. Compared to the present study, 22 DEGs and 10 DEPs were found to overlap (Table 1). Of these 22 DEGs, nine DEGs (ENSSSCG00000006985, nocturnin, SLCO2A1, squalene epoxidase [SQLE], UDP-glucose 6-dehydrogenase, Rho guanine nucleotide exchange factor 16, klotho beta, nei like DNA glycosylase 1, and cortactin binding protein 2) showed a similar expression pattern in both studies. The remaining 13 genes (Fos proto-oncogene, growth arrest and DNA damage inducible beta, insulin like growth factor 1, SQLE, FAD dependent oxidoreductase domain containing 2, TNF alpha induced protein 3, early growth response 1, solute carrier family 1 member 1, ENSSSCG 00000016119, ENSSSCG00000006716, RAR related orphan receptor A, 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3, neural EGFL like 2, and Ras association domain family member 6) showed the opposite pattern between studies. Furthermore, 10 DEPs were found to overlap with this previous study. Four proteins (NADH:ubiquinone oxidoreductase subunit B5, aminoadipate-semialdehyde synthase, collagen type XIV alpha 1 chain, and nerve growth factor receptor) showed a similar expression pattern, while six (phytanoyl-CoA dioxygenase domain containing 1, CXXC motif containing zinc binding protein, nudix hydrolase 6, thiopurine S-methyltransferase, decapping enzyme, scavenger, and UDP-galactose-4-epimerase) showed the opposite expression. The divergence between groups could be explained by different living environments and molecular genetic characterization.