Introduction

Low air pressure and low oxygen partial pressure at high altitude seriously affect the survival and development of human beings and other animals1, 2. Indigenous animals of the Tibetan Plateau exhibit heritable adaptations to this extreme environment due to natural selection; therefore, the Tibetan pig is an ideal animal model for research on the molecular ecology and pathology related to hypoxia3,4,5.

Adaptation to hypoxic conditions is a complex trait that is influenced by many factors6. Several genes exhibiting strongly selection have been identified in Tibetan populations by genome scanning, providing a possible genetic basis for adaptation to high-altitude conditions7,8,9. Likewise, a number of rapidly evolving, positively selected genes have been reported in the Tibetan pig10; however, the mechanisms involved in gene expression at high altitude are not fully understood. The development of RNA-seq and isobaric tags for relative and absolute quantitation (iTRAQ) technologies has enabled the identification of eukaryotic genes associated with complex traits via analysis of transcriptomic and proteomic profiles, with low bias, a large dynamic range, low frequency of false-positive signals, and high reproducibility. The iTRAQ method employs a set of amine-reactive isobaric tags to derivatize peptides at the N-terminus and at lysine side chains, thereby allowing for simultaneous protein identification and quantification via mass spectrometry analysis of peptide fragments (m/z range from 150 to 2,200) and signature ions (m/z values from 113 to 121), respectively. Both RNA-seq and iTRAQ have been widely used to screen for functional genes involved in muscle growth and lipid deposition in pigs and other domestic species11,12,13,14. An integrated analysis of the transcriptome and proteome of pig heart tissues at high and low altitudes would comprehensively profile hypoxic adaptations, and enable us to characterize the genes and proteins that are involved in this process.

Heart tissues were previously used to study hypoxic physiology and pathology15, 16. An increased heart rate, blood pressure, and other high-altitude responses result in changes in cardiac structure and function. In particular, animals populating the Tibetan Plateau exhibit strong cardiac function, and their cardiomyocytes show an ability to promote the expression of genes associated with adaptation to hypoxic conditions to protect cardiac cells and tissue structure, and increase blood circulation. Indeed, Tibetan pigs have adapted to high altitudes with well-developed hearts and lungs, and a blunted erythropoietic response17.

Cellular responses to hypoxia are reliant on controlled transcriptional and post-transcriptional events, in which certain genes are positively regulated and others either remain inactive or are actively repressed. In this study, we performed a comparative analysis of the transcriptomic and proteomic profiles of heart tissues obtained from Tibetan and Yorkshire pigs raised at high (TH and YH, respectively) and low (TL and YL, respectively) altitudes using RNA-seq and iTRAQ technologies. Via these analyses, we identified key genes and molecular mechanisms involved in the high- altitude adaptations of the Tibetan pig.

Results

Overview of RNA transcriptomic profiles of pig heart tissues

A total of 551 M paired-end reads were acquired from the RNA-seq experiment. Sequence alignment was performed using the pig genome (Sscrofa10.2) as a reference. Approximately 83.68% (82.07–85.18%) of the mapped reads were acquired from the RNA-seq experiment of which 76.26% (74.86–77.74%) were mapped to unique genomic locations. Additionally, of the total mapped reads, roughly 89% in each group corresponded to exons (Table S3, Fig. S1).

Gene expression levels were evaluated by counting the number of mapped reads per gene. In total, 18,585 expressed genes were detected in the heart tissues examined; of these, 15,701 were detected in each of the four groups (Fig. S2A). Less than 1% of the genes were expressed at greater than 1,000 FPKM (numbers of fragments per kilobase of exon per million mapped fragments); approximately 3% were expressed between 100 and 1,000 FPKM, and approximately 96% were expressed at less than 100 FPKM (Table S3, Fig. S2B). Expression value distributions were similar among the eight samples tested. Likewise, the FPKM values were evenly distributed amongst the four groups, and there were no outliers (Fig. S2C).

Analysis of differentially expressed genes (DEGs)

Pair-wise comparisons with the strict criteria of |log2 (FC)| > 1 and P < 0.01 were utilized to detect DEGs (Fig. S3A). Specifically, 299 and 242 DEGs were detected in the comparisons of TH to YH (Table S4–1, Fig. S3B) and of TL to YL (Table S4–2, Fig. S3B), respectively, including 26 DEGs that were common to each comparison (Fig. S3C). Meanwhile, we detected 169 and 368 DEGs between TH and TL (Table S4-3, Fig. S3B) and YH and YL (Table S4-4, Fig. S3B), respectively, of which 24 overlapped (Fig. S3C).

Functional analysis of DEGs

We used DAVID and Ingenuity Pathways Analysis (IPA), which have related capabilities but utilize different databases, to obtain a complete functional view of the DEGs. Using DAVID, we observed an overrepresentation of gene ontology (GO) terms related to “response to hormone stimulus” (19 genes, P = 3.15E-05) and “mitochondrion” (36 genes, P = 6.30E-05) in TH vs. YH, and “oxidation reduction” (51 genes, P = 1.77E-16) and “cofactor binding” (32 genes, P = 2.44E-15) in TL vs. YL (Table S5). We also observed an overrepresentation of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways “ribosome” (15 genes, P = 8.52E-10) and “dilated cardiomyopathy” (7 genes, P = 9.08E-03) in TH vs. YH, and “glycolysis/gluconeogenesis” (9 genes, P = 3.02E-04) and “renin-angiotensin system” (5 genes, P = 1.22E-03) in TL vs. YL. Analysis of these two comparison groups also highlighted GOs and pathways, such as “regulation of blood pressure”, “vascular process in circulatory system”, “response to hypoxia”, and “complement and coagulation cascades pathway” (Table S5). Additionally, the mammalian target of rapamycin (mTOR) signaling pathway was detected in the TH vs. YH comparison; hypoxia-induced mTOR activation reduces mTOR targeting to enhance angiogenesis in response to hypoxia18. Moreover, we identified a variety of DEGs in pathways upstream of the mTOR pathway19, in the insulin signaling pathway (seven of the DEGs were designated to encode oxygen-independent regulators).

In the TH vs. TL analysis, the overrepresented KEGG pathways were “contractile fiber” (12 genes, P = 1.68E-08), “muscle contraction” (10 genes, P = 1.22E-03), “dilated cardiomyopathy” (6 genes, P = 5.35E-03), “complement and coagulation cascades” (5 genes, P = 1.01E-02), while in YH vs. YL, they were “oxidation reduction” (32 genes, P = 2.10E-10), “=2.10biosynthetic processl (12 genes, P = 7.72E-09), “glycolysis/gluconeogenesis” (9 genes, P = 1.37E-05), and “histidine metabolism” (6 genes, P = 1.79E-04). Notably, GOs and pathways, such as the “mitochondrion”20, “glucose metabolic process”, “dilated cardiomyopathy”16, and “cardiac muscle contraction”21 pathways, which play important roles in the hypoxia response, were also detected in TH vs. TL.

We also performed metabolic pathway analysis on the DEGs using IPA software. Based on the data obtained via this method, the main canonical pathways that were overrepresented in the TH vs. YH and TL vs. YL comparisons were related to hepatic fibrosis/hepatic stellate cell activation (FN1, CTGF, CCL5, A2M, and FGF1, and COL1A2, CTGF, AGT, and COL3A1; P = 3.85E-02 and P = 2.58E-02, respectively), and glucocorticoid receptor signaling in TH vs. YH (DUSP1, SLPI, CCL5, and A2M, and DUSP1 and AGT; P = 1.82E-01 and P = 3.66E-01, respectively). Thus, depending on the environment, Tibetan pigs modulated the expression of distinct DEGs to adapt to hypoxic conditions. In addition, interesting pathways such as the hypoxia-inducible factor 1 (HIF-1) signaling pathway (EGLN3) were detected in both the TH vs. YH and TH vs. TL comparisons (Table S6).

Of the primary molecular and cellular functions that were significantly overrepresented in the four comparison groups, according to IPA, the most relevant were involved in “cardiovascular system development and function” in TH vs. YH (23 genes, P = 1.33E-09–2.40E-03), “cardiac hypertrophy” in TL vs. YL (5 genes, P = 1.393E-05–1.90E-01), “cardiovascular system development and function” in TH vs. TL (16 genes, P = 4.22E-07–4.69E-03), and “skeletal and muscular system development and function” in YH vs. YL (3 genes, P = 7.05E-06–9.10E-03).

Among the related specific functions of heart tissue, the main functions identified by IPA in TH vs. YH were “angiogenesis” (ANXA2, CCL5, CD151, CD44, CTGF, EGLN3, FGF1, FN1, HDAC9, ITGA5, JUNB, NR4A1, PFKM, SLPI, and TGM2; P = 1.33E-09), “binding of endothelial cells” (ANXA2, CCL5, CD151, CD44, FGF1, FN1, and NR4A1; P = 7.07E-09), and “development of blood vessel” (ANXA2, ATF3, CCL5, CD151, CD44, CTGF, EGLN3, FGF1, FN1, HDAC9, ITGA5, JUNB, NR4A1, and TGM2; P = 9.53E-08); in TL vs. YL were “hypertrophy of cardiomyocytes” (AGT, CRYAB, CTGF, and DUSP1; P = 1.39E-05) and “hypertrophy of heart” (AGT, CRYAB, CTGF, and DUSP1; P = 2.44E-04); in TH vs. TL were “formation of endothelial tube” (ATF3, CCL2, CD44, FN1, and PPP1R14B; P = 4.22E-07) and “angiogenesis” (CCL2, CD151, CD44, CRYAB, FGF1, FN1, HSPB8, ITGAV, JUNB, PFKM, and TGM2; P = 5.22E-06); and in YH vs. YL were “quantity of muscle cells” (AGT, CCL2, and FLNC; P = 7.05E-06) and “quantity of smooth muscle cells” (AGT and CCL2; P = 1.03E-04) (Table S6). Notably, we obtained the same pathways via functional analysis of TH vs. TL and YH vs. YL by IPA. Thus, even when inhabiting the same environment, Tibetan and Yorkshire pigs utilize distinct genes to adapt to hypoxia. All of the DEGs that were related to hypoxic adaptation in the four groups are shown in Table 1.

Table 1 Potential key differentially expressed genes (DEGs) and their functions related to hypoxic adaption in the Tibetan pig.
Full size table

Based on functional annotation, 21 DEGs were associated with hypoxia (Table 1), which involved the GO terms of “response to hypoxia”, “blood vessel development”, “response to reactive oxygen species”, “circulatory system process”, “regulation of heart contraction”, and “immune system development”, as well as the HIF-1 signaling pathway, the mitogen-activated protein kinase (MAPK) signaling pathway, pathways in cancer, glycolysis/gluconeogenesis, and the TGF-beta signaling pathway. Three genes, NPPA (natriuretic peptide type A, also known as ANP), NPPB (natriuretic peptide type B, also known as BNP), and FOS (FBJ murine osteosarcoma viral oncogene homolog), exhibited substantial differential expression, with fold changes greater than 14 in TH vs. YH (Table 1). In particular, the expression level of NPPB in TH was 41.96, which was 138.23-fold higher than that in YH and TL, and may therefore exert a positive effect on hypoxic adaptation in these animals. Similarly, NPPA was also highly expressed in this group, with an expression level of 37.60, which was 27.27-fold higher than that in YH and TL. In humans, ANP and BNP constitute the dual natriuretic peptide system of the heart. Numerous studies have clearly shown that hypoxia stimulates atrial natriuretic peptide (ANP) secretion, resulting in cellular adaptation to hypoxia and protection of the ischemic heart22,23,24.

Overview of protein identification and quantification by iTRAQ

A total of 128,390 spectra were obtained from the 8PLEX liquid chromatography tandem-mass spectrometry (LC-MS/MS) analysis (Fig. S4A). Pooling of the samples from both groups resulted in identification of 12,880 peptides at the peptide level, including 2578 proteins that had been originally identified with false discovery rates (FDR) ≤0.01, and with high confidence of the correct peptide sequence assignment. The numbers of proteins identified at distinct molecular weight ranges were as follows: 10–20 kDa (433), 20–30 kDa (443), 30–40 kDa (370), 40–50 kDa (318), 50–60 kDa (246), and 60–70 kDa (138). In total, these proteins accounted for 75.56% of those identified (Fig. S4B). In addition, the majority of the proteins were identified with high peptide coverage; 95.30% of the proteins had less than 50% sequence coverage, while 4.70% had greater than 50% sequence coverage (Fig. S4C). Among the proteins identified, 64.08% were represented by fewer than five peptides (Fig. S4D), indicating good sequence coverage of the proteins identified.

Analysis of differentially expressed proteins (DEPs)

Based on the selection criteria of FC > 1.2 or <0.83, we detected 473, 297, 394, and 297 DEPs in the TH vs. YH, TL vs. YL, TH vs. TL, and YH vs. YL comparisons, respectively. Specifically, 253 up-regulated and 220 down-regulated DEPs were identified in the TH vs. YH comparison, while 115 and 182, 275 and 119, and 179 and 188 up-regulated and down-regulated proteins were detected in TL vs. YL, TH vs. TL, and YH vs. YL, respectively (Table S7, Fig. S5). Notably, there were 87 and 86 overlapping proteins in the TH vs. YH and TL vs. YL, and in TH vs. TL and YH vs. YL comparisons, respectively (Table S8, Fig. 1).

Figure 1
figure 1

Venn diagram of differentially expressed genes (DEGs) among the four comparison groups. TH, Tibetan highland pig; TL, Tibetan lowland pig; YH, Yorkshire highland pig; YL, Yorkshire lowland pig.

Full size image

Validation of DEGs and DEPs

To validate the accuracy of the DEGs detected by RNA-seq analysis, we used real-time reverse transcription-quantitative polymerase chain reaction (RT-qPCR) to evaluate the expression levels of eight DEGs: ADHFE1, CCL5, CKM, DECR1, DPP4, EGLN3, FOS, and SLA-5. The expression levels of these genes in each group are shown in Fig. S6. The eight genes selected were differentially expressed among the four comparison groups (TH vs. YH, TH vs.TL, YH vs. YL, and TL vs. YL), and the RNA-seq data were concordant with those obtained by RT-qPCR. Moreover, RT-qPCR analysis confirmed the DEG expression patterns observed in the two pig breeds at each altitude (Fig. S7).

Notably, the results of previous studies suggest that the label-free method is typically more efficient, in regard to the number of proteins identified, than the iTRAQ approach25,26,27. In the current study, a total of 238 of the 351 (67.81%) and 686 (34.69%) proteins identified using the iTRAQ and label-free methods, respectively, were identified by both methods (Fig. S8A). For these proteins, the fold changes observed via each method were plotted, and a correlation with an R2 value of 0.604 was calculated (Fig. S8B). Of the overlapping 238 proteins, only three were found to be differentially expressed in divergent directions (i.e., up-regulation vs. down-regulation) in the comparison of TH vs. YH via iTRAQ and label-free experiments. Indeed, the overall trends of the commonly quantified proteins were in agreement between the two methods.

Functional analysis of DEPs

To better understand the biological functions of DEPs in Tibetan and Yorkshire pigs at different altitudes, the DEPs were further classified to identify pathways in the heart tissues based on GO and KEGG functional annotations (Table S9, Fig. 2). In TH vs. YH and TL vs. YL, the GO terms were mainly associated with “mitochondrion” (involved in 39 proteins and 44 proteins, respectively), “oxidation reduction” (involved in 23 proteins and 23 proteins, respectively), “blood circulation” (contained proteins such as CAV2, ACE, RENBP, AGT, MYL1, COL1A2, EPHX2, ERAP1, STAT1, HBB, TCAP, and MYH6), “NAD or NADH binding” (involved in 5 proteins and 4 proteins, respectively), “response to oxygen levels” (involved in 8 proteins and 7 proteins, respectively), and “iron ion homeostasis” (contained proteins such as TF, TFRC, HPX, FTH1, HP, and CP). In addition, we found interesting GOs in TH vs. YH, such as “gas transport” (HBQ1, CA2, AQP1, and HBB), “regulation of blood pressure” (ACE, RENBP, AGT, COL1A2, EPHX2, ERAP1, and HBB), and “oxygen and reactive oxygen species metabolic process” (CRYAB, AGT, EPHX2, BNIP3, and PARK7).

Figure 2
figure 2

Differentially expressed protein (DEP) pathway analysis. TH, Tibetan highland pig; TL, Tibetan lowland pig; YH, Yorkshire highland pig; YL, Yorkshire lowland pig.

Full size image

It is worth noting that 2–5% of the DEPs detected in the TH vs. YH, TL vs. YL, and TH vs.TL comparisons belonged to the HIF-1 signaling pathway. Likewise, 3–5% of the DEPs identified in the in TH vs. TL and YH vs. YL comparisons were associated with the vascular endothelial growth factor (VEGF) signaling pathway, which regulates angiogenesis and plays important roles in adaptation to hypoxic conditions. Meanwhile, the MAPK signaling pathway was enriched in 3–7% of the DEPs and pathways in cancer were enriched for 3–5% of the DEPs in the four comparison groups. The MAPK signaling pathway participates in HIF activation, and is involved in various cellular functions, including cell proliferation, differentiation, and migration28, while the proteins in pathways involved in cancer influence sustained angiogenesis, cellular proliferation, genomic damage, and inhibition of differentiation, all of which are closely related to cancer growth and development under hypoxic conditions29. In addition, the insulin signaling pathway, hypertrophic cardiomyopathy (HCM), glycolysis/gluconeogenesis, and dilated cardiomyopathy can be involved in hypoxic adaptation. All of the potential DEPs that were related to hypoxic adaptation in the four groups are shown in Table 2.

Table 2 Summary of hypoxia-regulated differentially expressed proteins (DEPs).
Full size table

Integrated analysis of DEGs and DEPs and iTRAQ and RNA-seq data

In total, 18,585 genes and 2,578 proteins were identified in the RNA-seq and iTRAQ analyses, of which 10,001 and 2,170 were annotated, respectively. Using these methods, the expression of 1,402 genes was detected at both the mRNA and protein levels (Fig. S9A). The proteome and transcriptome reflect gene expression at two different levels, and one of the goals of the combined analysis was to achieve data complementation and to obtain more complete information regarding the gene expression profiles of each organism at high and low altitudes. We obtained 473, 297, 394, and 297 DEPs in TH vs. YH, TL vs. YL, TH vs. TL, and YH vs. YL, respectively, of which 378, 220, 319, and 226, respectively, were annotated. A total of 26, 8, 23, and 7 genes were identified as DEGs and DEPs in the TH vs. YH, TL vs. YL, TH vs. TL, and YH vs. YL comparisons, respectively (Fig. S9B). The Pearson correlation coefficients obtained from a log2 function for TH vs. YH, TL vs. YL, TH vs. TL, and YH vs. YL were 0.40, 0.47, 0.60, and 0.42, respectively (Fig. S9C), indicating that only partial correlations were found at the mRNA and protein levels of overall gene expression. In summary, based on these results, we suggest that a substantial degree of post-transcriptional regulatory activity occurs during hypoxic adaptation, which has not been described in previous studies that focused solely on analyses of transcriptional and/or translational regulation.

Discussion

The Tibetan pig, which inhabits high-altitude regions on the Qinghai-Tibet Plateau, has evolved genetic adaptations to these extreme elevations via natural and artificial selection3. Previous studies demonstrated that Tibetan pigs have a distinct suite of phenotypic and physiological characteristics, in regard to pulmonary vascular structures and blood flow30, 31. Sequencing of multiple individuals from various pig breeds revealed that certain genomic regions, including genes involved in the hypoxia response, olfaction, energy metabolism, and drug responses, are under selection in the Tibetan pig10. Moreover, miRNA-seq analysis identified several miRNAs that play regulatory roles in the hypoxic adaptation of these animals4. High-altitude adaptation is not controlled by a single gene, but by multiple evolved genetic adaptations acting in concert32. In the current study, we screened for key genes and proteins related to hypoxic adaptation by comparing Tibetan and Yorkshire pigs raised at different altitudes via RNA-seq and iTRAQ protein sequencing analyses. Subsequent functional enrichment analysis of each DEG and DEP identified via these analyses using GO and KEGG software identified 21 DEGs (A2M, AGT, ANXA2, ATF3, CD44, CCL5, COL3A1, CRYAB, CTGF, DUSP1, EGLN3, FGF1, DPP4, HBB, HBA, NPPA, NPPB, DECR1, TGFB2, PDLIM3, and FOS) and 20 DEPs (MYL7, PTPMT1, PDLIM3, A2M, C6, CRYAB, PDK4, ERK2, ACE, NAPG, CPS1, HSPE1, DECR1, AHSG, ENO3, MRPS26, DCN, COL3A1, COL1A2, and KRT8) that appear to be associated with adaptation to hypoxic conditions in the Tibetan pig. In particular, A2M, COL3A1, CRYAB, DECR1, and PDLIM3 were identified as both DEGs and DEPs.

Under hypoxic conditions, VEGF-mediated induction of MEF2C and its target gene A2M is strongly reduced, as is the expression of HLX and UNC5B. Together, these data suggest that up-regulation of MEF2C/A2M (and HLX/UNC5B) by VEGF is most strongly pronounced under normoxic conditions, and is inversely correlated with hypoxia33. We found that the expression level of A2M in TH and TL hearts was higher than that in YH and YL hearts, respectively, which could explain why Tibetan pigs adapt better to hypoxic conditions than Yorkshire pigs, regardless of altitude. Meanwhile, TH hearts exhibited 4.68-fold higher levels of CRYAB expression than YH and TL hearts. Given that the siRNA mitochondrial pathway is activated in hypoxic cardiomyocytes in the presence of CRYAB, as evidenced by increased production of cytosolic cytochrome c34, our findings indicate that enhanced expression of this protein could comprise another reason for the superior adaption to hypoxic conditions of TH, compared to YH or TL. PDLIM3 (PDZ and LIM domain protein 3, also known as ALP) is necessary to maintain normal myocardial contractility. GO annotations related to this gene include sion. In summary, based nding and structural constituent of muscle”, and cytoskeletal mutations are known causes of genetically based forms of dilated cardiomyopathy. Moreover, the disruption of the gene encoding ALP is associated with right ventricular chamber dilation and dysfunction, directly implicating α-actinin-associated proteins in the onset of cardiomyopathy35. PDLIM3-specific up-regulation in Tibetan pigs might therefore increase cardiovascular permeability, promote vascular endothelial cell proliferation and migration, and increase the regulation of cardiovascular development to prevent the occurrence of cardiomyopathy.

In this study, we identified a set of associated DEGs (CD44, CTGF, DPP4, CRYAB, EGLN3, FGF1, and FOS) with at least one functional enrichment, such as “angiogenesis”, “binding of endothelial cells”, “development of blood vessel”, MAPK signaling, HCM, and pathways in cancer (Table 1). Of these, CTGF (connective tissue growth factor) and FGF1 (fibroblast growth factor 1), which are involved in blood vessel and vascular development, were expressed at significantly higher levels in TH than in YH, suggesting an important role for the cardiovascular system in hypoxia adaptation. In particular, CTCF restricts upstream enhancers from activating VEGF, thereby inhibiting the induction of VEGF and angiogenesis36. Hypoxia is the principal driver of induction of VEGF during both physiological and pathological angiogenesis37. FGF1 is recognized as an angiogenesis factor in vivo, and angiogenesis is known to be induced by hypoxia38. Our findings therefore indicate that Tibetan pigs must increase CTGF and FGF1 expression and modify their cardiovascular response to hypoxia by increasing blood flow and cardiac pumping.

The expression levels of CD44 and PDLIM3 were also significantly higher in the heart tissues of TH than in those of YH, TL, or YL. CD44 (CD44 molecule) transmembrane glycoproteins are cell adhesion molecules associated with aggressiveness and metastasis37,38,39,40, and cycling hypoxia has been shown to increase the population of CD44+/CD24− cells in a metastatic breast cancer cell line41. Therefore, it seems probable that CD44 expression is regulated by hypoxia. The functional effects of the up-regulation of CD44 and its variant isoforms under hypoxic conditions, particularly in the context of hyaluronan levels, should be considered, as cell signaling events that promote anchorage-independent tumor cell growth, survival, migration, and metastasis occur through the binding of hyaluronan with CD4440, 42, 43. In addition, expression of FOS (also known as c-FOS), which was previously reported to be associated with the hypoxia response44, 45, was 17.3-fold higher in TL than in YL. Together, these data suggest that characterization of DEGs that differ between the two pig breeds could provide clues to understanding the distinct hypoxic adaptations of these animals.

According to the GO and pathway analyses, the DEGs CRYAB (crystallin, alpha B), DPP4 (dipeptidyl-peptidase 4), and EGLN3 (egl nine homolog 3) comprise three functional candidate genes that were related to “response to hypoxia”. Indeed, EGLN1 and EGLN3 are strongly induced by hypoxia in most cell types46. his effect is likely important for cellular adaptation to hypoxic conditions, and is responsible for the increased oxygen-mediated HIF-1a degradation observed after long periods of hypoxia47. Indeed, the induction of EGLN3 mRNA expression by hypoxia is particularly significant when compared to that of other hypoxia-responsive genes44, 46, 48, 49. DPP4 is a novel marker of HIF-1 activity in tumors50. In our study, the DPP4 expression level in TH hearts was 0.22-fold higher than that in YH hearts, and 0.31-fold higher than in TL or YL hearts.

In the presence of low atmospheric oxygen, the oxygen levels in highland Yorkshire pigs introduced from the mainland were correspondingly low, and the difference in oxygen levels between capillary blood and mitochondria was small, resulting in hypoxia. A series of genes, such as FLT-1 [fms-related tyrosine kinase 1 (vascular endothelial growth factor/vascular permeability factor receptor)], which was categorized with the GO term “vascular regulation” and HIF-1 signaling pathway in the YH group and is considered the strongest of the VEGF receptors, as well as the HIF-1 signaling pathway, were up-regulated in Yorkshire pigs in responses to hypoxia in the YH group. Hypoxic exercise can result in increased production of VEGF and Flt-1 in muscle tissue; following VEGF protein production, Flt-1 expression at the vascular endothelial cell membrane can be regulated via autocrine or paracrine responses and participate in muscle tissue angiogenesis51. In this study, expression of the Flt-1 gene was the highest in YH and the lowest in TH, possibly because compensatory responses increase vascular permeability in response to a hypoxic environment.

The specific KEGG enrichments in the Yorkshire pig were primarily concentrated in energy metabolism pathways, such as the glycolysis pathway. This was not surprising as, under hypoxic conditions, myocardial contraction requires increased adenosine triphosphate (ATP) consumption. However, while most differentiated cells employ mitochondrial oxidative phosphorylation to generate ATP to maintain cell processes under normoxic conditions, cells exposed to hypoxia typically produce ATP via glycolysis, which produces lower levels of this energy source than aerobic oxidation. As such, under hypoxic conditions, animals must up-regulate glycolysis to meet the consumption requirements of myocardial cells. In particular, the primary DEGs associated with energy metabolism in Yorkshire pigs were ITGA5, PCK1, and PCK2. Glycolysis produces phosphoenolpyruvate (PEP), and the overexpression of PCK1 can lead to increased PEP production and gluconeogenesis52.

In this study, many genes were found to exhibit consistency in regard to changes in transcript and protein levels; however, there were several cases in which inconsistent results were obtained by RNA-seq and iTRAQ analysis, suggesting that post-transcriptional regulation plays an important role in the adaptive response to hypoxia in Tibetan pigs. Both transcriptomic and proteomic data are important for deciphering the molecular processes involved in this process. Previous studies have shown that the abundance of mRNA transcripts is does not completely correlate with expression level of certain genes53. Indeed, we also found only partial correlations between the mRNA and protein expression levels of certain genes54. Thus, transcriptomic analyses do not fully represent protein expression55. As such, analysis of transcriptomic or proteomic data alone is insufficient for accurate characterization of functional mechanisms; only by combining these approaches can we obtain a more comprehensive understanding of the biological functions of organisms.

In conclusion, we identified important hypoxia-adapted genes and proteins in Tibetan pigs, and by combining RNA-seq and iTRAQ data obtained from pig heart tissues, we elucidated the regulatory relationship between DEGs and DEPs. A combination of transcriptomic and proteomic data revealed several key candidate regulators (A2M, PDLIM3, CRYAB, and ACE) and pathways (HIF-1 signaling pathway, cardiovascular development, VEGF signaling pathway, and pathways in cancer) that might play high-priority roles in the hypoxic adaptation of Tibetan pigs. These results provide new insights into the molecular mechanisms involved in hypoxia-adaptation regulatory networks, and a greater understanding of human hypoxic diseases. However, further studies are required to confirm the DEGs and DEPs involved in regulating hypoxia identified in this work.

Materials and Methods

Sample preparation and extraction

The experimental design comprised four treatment (comparison) groups (n = 8 per group): Tibetan and Yorkshire pigs raised in the highlands [Linzhi, Tibet, 3000 m above mean sea level (AMSL)] and lowlands (Beijing, China, 100 m AMSL). The animals were slaughtered at the age of six months, and heart tissue samples were collected, immediately frozen in liquid nitrogen, and stored at −80 °C (Table S1). The procedures for animal care were approved by the Animal Welfare Committee of the State Key Laboratory for Agro-biotechnology of the China Agricultural University (Approval number XK257), and all experiments were conducted in accordance with approved relevant guidelines and regulations.

RNA isolation, library preparation, and sequencing

Total RNA was isolated using Trizol® reagent (Invitrogen, Waltham, MA, USA)55. The integrity and concentration and purity of each sample was evaluated via 1% agarose gel electrophoresis and using a NanoDrop™ 2000 Biophotometer (Thermo Fisher Scientific, Waltham, MA, USA), respectively. Samples were then reverse transcribed into cDNA using Superscript II reverse transcriptase (Invitrogen) and random hexamer primers. For RNA-seq analyses, RNA pools were created using equal quantities of RNA from four individuals; two biological replicates were included in each group. RNA-seq libraries were constructed according to the manuals provided by Illumina, Inc. (San Diego, CA, USA), and were sequenced using the HiSeq 2000 platform to generate 100-bp paired-end reads. All RNA sequencing data are deposited in the Gene Expression Omnibus under accession number GSE92981.

Mapping and annotation of sequencing reads

Raw RNA-seq reads were arranged using CLC Genomics Workbench 4.8 software (CLC Bio, Aarhus, Denmark). All reads for which the quality of more than half of the bases was less than 10, as well as all reads that contained more than two Ns or were smaller than 20 bp in length, were eliminated from subsequent analyses. After removal of the adapters, the remaining clean reads were aligned to the whole pig genome (Sscrofa10.2.72) (ftp://ftp.ensembl.org/pub/release-72/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa10.2.72.dna.toplevel.fa.gz) using TopHat (version 2.0.9) software56; two mismatches were allowed per 100-bp read for each alignment. Finally, BAM files generated using SAMtools57 were used for subsequent analysis. Quality control and read statistics were determined using FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Quantification and differential gene analysis for RNA-seq

FPKM values obtained using Cufflink software (version 2.1.1) were used as values for normalized gene expression58. Differential expression analyses of four comparison groups (TH vs. YH, TH vs.TL, YH vs. YL, and TL vs. YL) were performed using the DESeq R package (1.10.1)59. The resulting P-values were adjusted using the Benjamini-Hochberg method. Results were expressed as the fold change (FC) of the average expression of case groups to that of the respective control group. Among the four comparison groups, DEGs were identified as those genes for which log2 (FC) > 1 and P < 0.01.

Protein isolation, enzymolysis, and iTRAQ labeling

Heart tissues were ground into powder in liquid nitrogen using lysis buffer (Roche). The resulting samples were then ultrasonically disrupted for extraction of total protein. After centrifugation at 10,000 × g for 30 min at 4ly, supernatants were collected, and protein concentrations were determined using an enhanced BCA (bicinchoninic acid) Protein Assay Kit (P0010; Beyotime Biotechnologies, Ltd., Beijing, China), according to the manufacturer’s instructions. The protein samples (200 μg) were mixed with dl-dithiothreitol and alkylated with iodoacetamide, and then treated with trypsin overnight at a trypsin-to-protein ratio of 1:100.

Protein peptides (15 μg) from each group were labeled using an 8plex iTRAQ reagent multiplex kit (SCIEX, Framingham, MA, USA). The samples were labeled as 113 (TH1), 114 (TH2), 115 (YH1), 116 (YH2), 117 (TL1), 118 (TL2), 119 (YL1), and 121 (YL2). The labeled samples were pooled and further fractionated offline using an ÄKTA Purifier 100 (GE Healthcare) with a strong cation exchange column (PolySULFOETHYL A™; PolyLC Inc., Columbia, MD, USA). The retained peptides were eluted with Buffer A [10 mM KH2PO4 in 25% ACN (acetonitrile), pH 3.0] and Buffer B (10 mM KH2PO4 and 500 mM KCl in 25% ACN, pH 3.0) with a flow rate of 1.0 mL/min.

LC-MS/MS analysis

Eluted fractions were lyophilized using a centrifugal speed vacuum concentrator (CentriVap® Complete Vacuum Concentrator; Labconco, Kansas City, MO, USA) and dissolved in 0.1% formic acid. Equivalent amounts of peptides from each fraction were mixed and then subjected to reversed-phase nanoflow LC-MS/MS analysis using a high-performance liquid chromatography system (EASY-nLC™, Thermo Fisher Scientific) connected to a hybrid quadrupole/time of-flight mass spectrometer equipped with a nano-electrospray ion source. The peptides were separated on a C18 analytical reverse-phase column with mixtures of Solution A (0.1% formic acid in water) and Solution B (0.1% formic acid in ACN). A full MS scan was conducted using a Q Exactive™ mass spectrometer (Thermo Fisher Scientific).

Database search and protein identification and quantification

For peptide identification and quantification, MS/MS data were searched against the “Sus_scrofa_35257_20151120_uniprot.fasta” file using Mascot 2.2 and Proteome Discoverer™ 1.4 software (Thermo Fisher Scientific). A unique protein with at least two unique peptides that had an FDR <0.0160 was used for data analysis. Protein quantification was based on the total intensity of the assigned peptides. The average of eight labeled sample mixes was used as a reference, and was based on the weighted average of the intensity of report ions in each peptide identified. The final protein ratios were normalized to the median average protein content of the 8plex samples. FC > 1.2 or FC < 0.83 was set as the threshold for identifying differentially expressed proteins.

GO and KEGG enrichment analysis of DEGs and DEPs

DEGs and DEPs were classified by GO and KEGG using DAVID online software (https://david.ncifcrf.gov/)61. For these analyses, an official gene symbol for each DEG or DEP was uploaded, and the species with the maximum number of annotations was used. The GO terms used were BP, CC, and MF, and KEGG pathways with corrected P-values < 0.05 were considered significantly enriched. In addition, IPA software (http://www.ingenuity.com/; Ingenuity Systems, Redwood City, CA, USA) was used to compare the DEGs among the four treatment groups. Accession numbers for these genes were imported into IPA, and the “Core Analysis” function was used to analyze the genes in the context of networks, biological functions, and canonical pathways. Detailed information concerning IPA analyses was published previously60, 62, 63. GO annotation and KEGG pathway analysis of the DEPs was also performed using KOBAS 2.0 (http://kobas.cbi.pku.edu.cn/)64 online software.

Verification of RNA-seq and iTRAQ data

RT-qPCR primers were designed to span the exon-exon boundaries of eight genes selected by RNA-seq; the characteristics of the primers are listed in Table S2. RT-PCR analysis was performed using a SYBR® Green I PCR Master Mix Kit (FP204; Tiangen Biotech Co. Ltd., Beijing, China) on a CFX96™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA), according to the manufacturer’s instructions.

The results for TH vs. YH obtained by iTRAQ analysis were validated via the label-free technique, as described previously. Label-free quantification is a mass spectrometry method used to determine the relative amounts of proteins of interest in two or more biological samples. Unlike other methods of protein quantification, label-free quantification does not use a stable isotope-containing compound to chemically bind to, and thus, label the protein63, 65, 66.