The Expression Regulatory Network in the Lung Tissue of Tibetan Pigs Provides Insight Into Hypoxia-Sensitive Pathways in High-Altitude Hypoxia

To adapt to a low-oxygen environment, Tibetan pigs have developed a series of unique characteristics and can transport oxygen more effectively; however, the regulation of the associated processes in high-altitude animals remains elusive. We performed mRNA-seq and miRNA-seq, and we constructed coexpression regulatory networks of the lung tissues of Tibetan and Landrace pigs. HBB, AGT, COL1A2, and EPHX1 were identified as major regulators of hypoxia-induced genes that regulate blood pressure and circulation, and they were enriched in pathways related to signal transduction and angiogenesis, such as HIF-1, PI3K-Akt, mTOR, and AMPK. HBB may promote the combination of hemoglobin and oxygen as well as angiogenesis for high-altitude adaptation in Tibetan pigs. The expression of MMP2 showed a similar tendency of alveolar septum thickness among the four groups. These results indicated that MMP2 activity may lead to widening of the alveolar wall and septum, alveolar structure damage, and collapse of alveolar space with remarkable fibrosis. These findings provide a perspective on hypoxia-adaptive genes in the lungs in addition to insights into potential candidate genes in Tibetan pigs for further research in the field of high-altitude adaptation.


INTRODUCTION
Tibetans are a unique and geographically isolated pig breed that inhabits the Qinghai-Tibet Plateau, which has an extreme environment with high altitudes Ma et al., 2019). This unique ecological condition is characterized by low air pressure, reduced oxygen content, and high ultraviolet radiation, imposing extreme physiological challenges on domestic animals, and failure to adapt will lead to altitude illness or even death (Cao et al., 2017;Lancuo et al., 2019;Qi et al., 2019). Native high-altitude species have been selected through evolutionary processes to evolve adaptive mechanisms to cope with this harsh environment . Special lung properties of the Tibetan pig, yak, and Tibetan sheep living in the plateau, such as larger lungs, thicker alveolar septa, and more developed capillaries, have been previously reported by Qi and Yang (Yang et al., 2014;Qi et al., 2019). Tibetan pigs exhibit heritable adaptations to high-altitude environments as a result of natural selection. Exposure to hypoxia changes the gene profiles in various cell types and is associated with adaptation to high altitudes (Zhang T. et al., 2019). mRNAs and miRNAs are involved in many biological processes in animals, and not surprisingly, transcriptional analyses have revealed the differential expression of hypoxia regulators that enable adaptation to a hypoxic environment (Ni and Leng, 2016). The hypoxia-inducible factor-1 (HIF-1), vascular endothelial growth factor (VEGF), and mitogen-activated protein kinase (MAPK) signaling pathways are typical hypoxia-associated pathways (Lee et al., 2016;Zhang et al., 2018;Nicolas et al., 2019), and some mRNAs (PHD2, VHL, and FIH-1) and miRNAs (miR-363, miR-421, and miR-204) have been implicated in the regulation of the HIF-1 signaling pathway (Semenza, 2007;Ge et al., 2016;Wang et al., 2016;Xie et al., 2016).
Studies of the molecular mechanisms of livestock adaptation to high altitude have focused on miRNA-mRNA interaction networks. Here, we performed an integrative analysis of the miRNA-mRNA expression profiles in the lungs of high-and lowaltitude pigs (Tibetan pigs and Landrace pigs, respectively) to identify molecular pathways and networks involved in the genetic adaptation of Tibetan pigs to hypoxic conditions.

Ethics Statement
All animal experiments were conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People's Republic of China (Approval number: 2006-398). The procedures for animal care were approved by the Gansu Agricultural University Animal Care and Use Committee of Gansu Agricultural University, and all experiments were conducted in accordance with approved relevant guidelines and regulations.

Sample Collection
In total, 18 Tibetan male piglets from the highlands (TH group; Gannan Tibetan Autonomous Prefecture, Gansu, representing an altitude of 3,000 m) and 18 Landrace male piglets from the lowlands (LL group; Jingchuan, Gansu, representing an altitude of 1,000 m) with similar weights and non-genetic relationships were selected, and nine piglets from each group migrated to low altitude (TL group; Tibetan pigs at low altitude) or high altitude (LH group; Landrace pigs at high altitude) from their original rearing facility at the age of 1 month. We randomly selected six pigs from each group to collect the left lower lobes of the lung from indigenous and imported adult male pigs at the age of 6 months. These animals (n = 6 in each group) were feed restricted for 12 h and slaughtered in their feeding place. Six samples from each group were immediately stored in stationary Abbreviations: TH, Tibetan male piglets from the highlands; LL, Landrace male piglets from the lowlands; TL, Tibetan male piglets migrated to the lowlands; LH, Landrace male piglets migrated to the highlands. liquid for hematoxylin and eosin (H&E) staining, and three of the six samples were randomly selected and collected within 1 h after the pigs were harvested and stored immediately in liquid nitrogen for subsequent RNA extraction.

Hematoxylin and Eosin Staining
Sections from the left lower lobes of the lung were stained with H&E (Ban et al., 2018;Zhang et al., 2019b), observed under a microscope (Sunny Optical Technology Co. Ltd, Ningbo, China), and then photographed using Image View (Sunny Optical Technology Co. Ltd).

RNA Extraction
Total RNA from the lungs was extracted using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol, and eukaryotic mRNA was enriched by oligo (dT) beads (Epicenter, Madison, WI, USA). RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and verified by 1% gel electrophoresis. All samples presented an RNA integrity number (RIN) > 7.5.

Library Construction and Sequencing for mRNA
After total RNA was extracted, eukaryotic mRNA was enriched by oligo (dT) beads (Epicenter) and reverse-transcribed into cDNA using random primers. mRNA was ligated with proper 5 ′ and 3 ′ adapters. The ligation products were reverse-transcribed by PCR amplification to generate a cDNA library, which was sequenced using an Illumina HiSeq TM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Library Construction and miRNA Sequencing
After total RNA was extracted for miRNA sequencing, 18-30 nt RNA molecules were enriched by polyacrylamide gel electrophoresis (PAGE). A 3 ′ adapter was added to enrich the 36-44 nt RNAs, and the 5 ′ adapter was then connected to the RNA. PCR products of 140-160 bp were amplified by reverse transcription. A cDNA library was generated and sequenced using Illumina HiSeq TM 2500 sequencing (Illumina Inc., San Diego, CA, USA) by Gene Denovo Biotechnology Co., Ltd.

Expression Analysis of mRNAs
High-quality clean raw data were screened by removing lowquality data with fastp (Chen et al., 2018). The short-read alignment tool, Bowtie 2 (Langmead and Salzberg, 2012) was used to map reads to the ribosome RNA (rRNA) database. An index of the reference genome was built, and paired-end clean reads were mapped to Sus scrofa RefSeq (Sus scrofa 11.1) using HISAT 2 (Kim et al., 2015). The mapped reads of each sample were assembled using StringTie v1.3.1 (Pertea et al., 2015(Pertea et al., , 2016 in a reference-based approach. For each transcription region, a fragment per kilobase of transcript per million mapped reads (FPKM) value was calculated to quantify its expression abundance and variations using RSEM software. RNA differential expression analysis was performed with DESeq 2 (Love et al., 2014) software between the two groups. The raw mRNA-seq data (accession number PRJNA687172) were submitted to the Sequence Read Archive (SRA) database of NCBI.

Expression Analysis of miRNAs
Clean reads were obtained by filtering raw reads, and all of them were aligned with small RNAs in the GenBank database (Benson et al., 2013). All the clean reads were aligned with small RNAs in the Rfam database (Griffiths-Jones et al., 2003) to identify and remove rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs. All the clean reads were also aligned with the reference genome and were searched against the miRbase database (Griffiths-Jones et al., 2006) to identify known (Sus scrofa) miRNAs. All the unannotated reads were aligned with the reference genome by HISAT2. 2.4. Novel miRNA candidates were identified according to their genome positions and hairpin structures predicted by mirdeep2 software. The miRNA expression levels were calculated and normalized to transcripts per million (TPM). The raw miRNA-seq data (accession number PRJNA687649) were submitted to the NCBI Sequence Read Archive (SRA) database.

Functional Annotation of DEmRNAs
DEmRNAs were analyzed using Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses using the online tool Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang et al., 2009) to explore their roles, functions, and enrichment in different biological pathways. Gene Ontology (GO) terms and pathways with q < 0.05 were considered significantly enriched by DEmRNAs. The hypoxic DEmRNAs were filtered based on the intersection of our results and published hypoxia-related genes in the HIF-1 signaling pathway. The hypoxia-related genes and target genes of miRNAs were also mapped to GO terms in the GO database and pathways in the KEGG (Kyoto Encyclopedia of Genes and Genomes) database to further elucidate their functions.
Target Prediction and Integrative Analysis of the Hypoxia-Related miRNA-mRNA Regulatory Network We identified mRNAs with a fold change ≥2 and a false discovery rate (FDR) <0.05 as DEmRNAs. To explore more DEmiRNAs, we identified miRNAs with fold change ≥ 2 and p < 0.05 as DEmiRNAs. The potential target genes of DEmiRNAs were predicted using RNAhybrid 89 (version 2.1.2) + svm_light (version 6.01), miRanda (version 3.3a), and TargetScan (version 7.0), and the genes at the intersection of the results from the three software packages were selected as predicted miRNA target genes. Because mRNAs and miRNAs have potential negative regulatory relationships, we assessed the expression correlation between a miRNA and its predicted target by the Pearson correlation coefficient (PCC). Subsequently, the negatively coexpressed miRNA-mRNA pairs with PCC < −0.7 and p < 0.05 were screened to construct miRNA-mRNA networks.
The coexpression network diagram of DEmRNAs and DEmiRNAs was generated using the PCC, and only the relationship pair network diagram of the top 300 is shown. The coexpression network diagram of the 273 hypoxic DEmRNAs is displayed, and the correlation between miRNA and mRNA was required to account for the top 5% of the total correlation. The potential regulatory network was constructed by Cytoscape (Szklarczyk et al., 2015).

Quantitative Real-Time PCR Validation
Total RNA from pulmonary tissues was extracted with a TRIzol reagent kit and reverse-transcribed into cDNA using a FastQuant cDNA first-strand synthesis kit (TianGen, China). SYBR R Premix Ex Taq TM II (TaKaRa, China) was used for real-time fluorescence quantitative analysis. In total, eight DEmRNAs and eight DEmiRNAs were randomly selected to determine sequencing accuracy. The primers used here were designed using Primer 5.0 software and are listed in Supplementary Tables 1, 2 (Supplementary Material 1).
The experimental data were analyzed with the 2 − CT method (Livak and Schmittgen, 2001). Statistical analyses were performed using GraphPad Prism 8.0 (GraphPad Software, San Diego, CA, USA) and SPSS 20.0 (SPSS, Chicago, IL, USA). The comparisons were conducted by one-way analysis of variance (ANOVA), and p < 0.05 was considered statistically significant.

Morphological Structure
H&E staining showed that the lung sections exhibited the following connective tissues with epithelia: pulmonary alveoli, smooth muscle, blood capillaries, bronchial tubes, and alveolar septa ( Figures 1A-H). The sections from the TH group were characterized by smooth muscle hyperplasia and larger alveoli, while those from the LH group were characterized by a thicker alveolar septum. In addition, the analysis showed that MMP2 expression had a similar tendency to the alveolar septum thickness among the four groups ( Figure 1I).

Identification of DEmRNAs in the Lung
In total, 12 cDNA libraries, which included six Tibetan pigs and six Landrace pigs at high and low altitudes, were sequenced from lung tissues (Supplementary Material 2). After quality filtering, 51,193,112,222 clean paired reads were obtained with 99.70-99.79% of clean reads mapped to the porcine reference genome (Table 1). A total of 471 DEmRNAs (247 up-and 224 downregulated) were identified in the TH group compared to the TL group (Figure 2, Supplementary Table 3 in Supplementary Materials 1, 3). Furthermore, 809 novel genes were identified in the sequencing data. Eight mRNAs were randomly selected and detected using qRT-PCR to validate the  accuracy of the sequencing data. Our verification test indicated that the qRT-PCR results were consistent with the mRNA-seq data ( Figure 4A).

Identification of DEmiRNAs in the Lung
A total of 12 cDNA libraries were sequenced from lung tissues.
In the miRNA-seq data, 10,810,538-14,920,316 clean reads were obtained by removing low-quality data and data with sequences shorter than 18 nt and longer than 30 nt, and 94.380-97.30% clean reads were obtained and mapped ( Table 2). A total of 464 DEmiRNAs (324 up-and 140-downregulated) were identified in the TH group compared to the TL group (Figure 3, Supplementary Table 4 in Supplementary Materials 1, 4). Eight miRNAs were randomly selected and detected using qRT-PCR to validate the accuracy of the sequencing data. Our verification test indicated that the qRT-PCR results were consistent with the miRNA-seq data ( Figure 4B).

Functional Analysis of DEmRNAs
GO and KEGG enrichment analyses showed that most DEmRNAs were involved in cellular processes and pathways related to cytokine-cytokine receptor interaction, the PI3K-Akt signaling pathway, and pathways in cancer ( Figure 5). Interestingly, a number of genes were mainly enriched in "response to stimulus (GO: 0050896)" of biological process among the four groups. GO: 0001071 is associated with nucleic acid binding transcription factor activity and was significantly enriched between the TH and LH groups. The top 20 pathways with the most significant enrichment were obtained. KEGG enrichment results revealed that most of these genes were significantly enriched in cancer pathways among Landrace pigs (LH and LL) (breast cancer and transcriptional misregulation in cancers) or high-altitude groups (LH and TH) (proteoglycans in cancer, pathways in cancer, breast cancer). A number of genes were significantly enriched in cytokinecytokine receptor interaction, hematopoietic cell lineage, and African trypanosomiasis among Tibetan pigs and Landrace pigs in the high-or low-altitude groups. Six pathways were significantly enriched among the high-(TH and LH) or lowaltitude (TL and LL) groups, and 15 pathways were significantly enriched between the Tibetan pig (TH and TL) and Landrace pig (LH and LL) groups.

Identification and Prediction Targets of DEmiRNAs
A total of 59,636 target DEmRNAs of 1,630 DEmiRNAs (365 functionally annotated miRNAs, 989 known miRNAs and 276 novel miRNAs) were analyzed (Supplementary Materials 5). In addition, multiple pathways and GO terms were associated with hypoxia traits. The analysis revealed KEGG pathways that were significantly related to genes targeted by DEmiRNAs, and the Wnt signaling pathway, metabolic pathway and hepatocellular carcinoma were the most significantly related (Figure 6). Interestingly, the results showed that the targets were primarily enriched in terms related to hypoxia adaptation. ssc-miR-210, ssc-miR-101, ssc-miR-7136-5p, ssc-miR-10b, ssc-miR-206, ssc-miR-1343, ssc-miR-142-5p, ssc-miR-421-5p, and ssc-miR-4331 were identified as key miRNAs. Functional assessment showed that 100, 56, and 104 putative targets were mainly enriched in the HIF, PI3K-Akt, and MAPK signaling pathways, respectively.

Screening of Differentially Expressed Hypoxia-Related mRNA Target miRNAs and Their Functional Enrichment Analysis
Functional analysis was conducted to understand the pathways and molecular interactions of DEmRNAs and DEmiRNAs. The DEmRNAs were enriched in a number of important pathways related to hypoxia, and we identified 273 significant DEmRNAs involved in hypoxia adaptation among the four groups (Supplementary Material 6). We predicted potential target miRNAs of mRNAs according to the negative regulatory effects of miRNAs on mRNAs, which were further considered veritable miRNA-mRNA pairs. To further reveal the regulatory relationship of node mRNAs and non-coding miRNAs, the resulting potential regulatory networks of miRNA-target genes associated with hypoxia-genes were constructed (Figure 7). The target DEmRNAs of DEmiRNAs were assessed using KEGG and GO analyses. The results indicated that 71.09, 17.00, and 11.90% (total of 273) of the genes were enriched in the biological process (BP), cell component (CC), and molecular function (MF) categories, respectively, in the TH-vs.-TL comparison (p < 0.05). A number of genes were targeted by hub miRNAs, such as novel-m0237-5p, novel-m0173-3p, and novel-m0142-5p, which had 45, 19, and 14 target mRNAs among the four groups, respectively. Furthermore, miR-2465x targeted HIF-1α, while novel-m0087-3p and novel-m0237-5p targeted HIF-3α.

Construction of the Coexpression Network Between DEmRNAs and DEmiRNAs in Response to Hypoxia
To explore the relationship between miRNAs and mRNAs in a hypoxic environment, a coexpression network of DEmRNAs and DEmiRNAs was constructed, and the top 300 relationship pair network diagrams are listed (Figure 8A,  Supplementary Material 7.1). The intersection of differentially expressed hypoxia mRNAs and miRNAs identified from the four group comparisons represented their differential expression in pig lungs with increasing altitude. TAR1-A, GPD1, ST8SIA5, and LENG8 were selected as the most affected mRNAs, and there were strong correlations with a number of miRNAs. Furthermore, a coexpression network of 273 hypoxic DEmRNAs and DEmiRNAs was constructed ( Figure 8B, Supplementary Material 7.2). MEF2C, AKAP6, NTRK2, MAPT, and GPR146 were selected as the most affected mRNAs, and there were strong correlations with a number of miRNAs.

DISCUSSION
A high-altitude environment plays an important role in the adaptation of native species, and it may modify gene transcription and may irreversibly affect specific phenotypes (Zhang et al., 2015;Ni et al., 2019). We used a complete migrant design to evaluate genes interacting with the environment and selected Tibetan pigs and Landrace pigs in both their native altitude environments and as migrants in a non-native environment. Our previous research identified that Tibetan pigs have heavier and wider lungs, thicker alveolar septa, and a denser vascular network than Landrace pigs. The hemoglobin (HGB) and mean corpuscular hemoglobin concentration (MCHC) of high-altitude pigs (Tibetan and Landrace pigs) were significantly higher than those of low-altitude pigs (Yang et al., 2021). We next investigated whether there are gene expression changes specific to Tibetan pigs that are responsible for hypoxic adaptation. Sequencing of multiple pigs from different breeds revealed that certain genomic regions, including genes involved in the hypoxia response, were under selection in Tibetan pigs (Zhang B. et al., 2016;Zhang et al., 2017Zhang et al., , 2019a. We screened for key genes related to hypoxic adaptation through genotype and environment interaction effects via RNA-seq analyses. Several pathways were enriched in DEmRNAs among Tibetan pigs and Landrace pigs at different altitudes, including the VEGF signaling pathway, PI3K-AKT signaling pathway, and mTOR signaling pathway (Ai et al., 2015;Zhang et al., 2017). Moreover, GO enrichment analysis revealed that these DEmRNAs were associated with vascular regulation, regulatory region DNA binding, or extracellular region. The identified hypoxia-related signaling pathways may form a complex cascade of responses that occur in hypoxic conditions in Tibetan pigs to reduce the risk of pulmonary damage. Hypoxia-regulated miRNAs play vital roles in cell survival and have been implicated in the regulation of both upstream and downstream HIF-1 signaling pathways under hypoxic conditions. For example, miR-199a, miR-17-92, and miR-20b induce HIFs (Dai et al., 2015;Chen et al., 2016;Danza et al., 2016). HIF-1 regulates the expression of various genes to protect cells from hypoxic injury through cell apoptosis, glucose metabolism, and mitochondrial function (Bhattarai et al., 2018;Yu et al., 2018Yu et al., , 2020. HIF-1α is a potential therapeutic proangiogenic molecule that regulates the levels of VEGF to elevate interstitial pressure (Zhi et al., 2018;Lin et al., 2019). Several putative target genes (FOXO3, RASGRF1, and CX3CR1) that are regulated by ssc-miR-214, ssc-miR-320, and ssc-miR-101 have been found to be involved in the HIF-1 related signaling pathway. miR-210 is located on human chromosome 11p15.5 and correlates with angiogenesis and VEGF regulation in breast cancer patients (Forkens et al., 2008;Dai et al., 2015;Tang et al., 2018;. In the present study, the expression of miR-210-x and miR-210-z was significantly lower in TH than in TL but not significantly different in the other groups, which may play vital roles in the expression of proteins in homology-dependent repair pathways and nucleotide excision repair pathways to reverse cellular DNA damage in the lungs of Tibetan pigs during hypoxia (Crosby et al., 2009;Hui et al., 2019). HBB is involved in the malaria reference pathway and downregulates IL-6, which is a key gene in the HIF-1 pathway. Comparison of HBB expression between Tibetan pigs and Landrace pigs showed that among the beta globin amino acid substitutions at positions 58, 75, 119 and 137, the replacement of alanine at position 137 with valine and the locus mutation improved the affinity of HGB and O 2 (Zhang B. et al., 2016). The expression of the HBB gene in Tibetan pigs (TH and TL) was significantly higher than that in Landrace pigs (LH and LL), agreeing with a similar trend previously reported by other authors, and there was similar variation in the HGB concentration in Tibetan pigs (Taliercio et al., 2013;Zhang B. et al., 2016;Yang et al., 2021), indicating that hypoxia transcriptionally upregulates HBB to increase HGB in the blood to ensure the transport of blood and nutrients. These findings may (Jang et al., 2014;Cai et al., 2018) explain why Tibetan pigs have better adaption than Landrace pigs in hypoxic environments regardless of altitude.
The PI3K/Akt pathway is an intracellular signaling pathway that is promoted by several biological molecules, including calmodulin, insulin-like growth factor (IGF), and multiple EGFlike domains 6 (MEGF6) (Pompura and Dominguez-Villar, 2018;Ellis and Ma, 2019;Revathidevi and Munirajan, 2019;Zhang et al., 2020). IGF2 is the target gene of miR-506y and ssc-miR-181d-3p. The expression levels of IGF2 and MEGF6 were significantly upregulated in LH compared to LL, but no differences were found in the Tibetan pigs (TH and TL). We hypothesized that these genes may induce the growth, proliferation, and differentiation of tumor cells in the lungs of Landrace pigs living in a hypoxic environment (Mohlin et al., 2013). Activated Akt induces various biological FIGURE 7 | Integrated miRNA-target negative correlation regulatory network. The square nodes represent miRNAs, and the circle nodes represent target genes. The differentially expressed miRNAs/target genes are highlighted as follows: red indicated differential expression between Landrace pigs raised in lowlands (LL) and Landrace pigs raised in highland (LH), green indicates differential expression between LH and Tibetan pigs raised in highland (TH), blue indicates differential expression between TH and Tibetan pigs raised in lowland (TL), and purple indicates differential expression between TL and LL.
processes, including activating mTOR, localizing FOXO to the cytoplasm, and activating cAMP-response element binding protein (CREB) (Zhang et al., 2011;Gaecía-Morales et al., 2017;Marquard and Jücker, 2020). The FOXO signaling pathway was also enriched in a comparison of pigs living at different altitudes. It has been shown that alcohol suppresses P450 oxidoreductase (POR) and glutathione reductase (GSR) gene expression by upregulating miR-214, which induces oxidative stress and plays a crucial role in adaptation to hypoxia (Zhou et al., 2013;Dong et al., 2014;Stefanetti et al., 2018;. FOXO3 is a targeted gene of ssc-miR-214-3p, and the expression of FOXO3 in the TH group was significantly higher than that in the TL group, but not significantly different between the LH and LL groups. Changes in ssc-miR-214-3p expression may inhibit the cell cycle and promote apoptosis, thereby inhibiting cell proliferation through FOXO3 prolyl hydroxylation in hypoxic conditions. The regulation of the expression levels of IGF2, MEGF6, and FOXO3 through miRNAs may lead to the better adaption of Tibetan pigs in hypoxic environments compared to the Landrace breed. Collagens, such as COL1A1, COL1A2, and COL3A1, are widely represented in ECM-receptor interactions and focal adhesion pathways (Gelse et al., 2003), and their expression was significantly higher in the LH group than in the LL group but not significantly different between the TH and TL groups. All of these genes function as mechanoreceptors and may provide a forcetransmitting physical link between the EMC and cytoskeleton, indicating that enhanced expression of COL1A1, COL1A2, and COL3A1 may be another reason for the superior adaption to hypoxic conditions of TH. Our study revealed that high expression of fibroblast growth factors (such as FGF1, FGF2, FGF9) was higher in the native groups (TH and LL) than in the migrated groups (TL and LH), which was alleviated by activating AKT3 (Pompura and Dominguez-Villar, 2018;Revathidevi and Munirajan, 2019). These findings indicated that Tibetan pigs may increase the expression of FGF1 and the cross-sectional area of a blood vessel to increase blood flow in response to hypoxia (Karar and Maity, 2011;Kir et al., 2018;Sajib et al., 2018).
GPR146 may be upregulated by a number of miRNAs (such as miR-8903, miR-11972 and miR-466-x) under hypoxic stimulation and has been suggested to be an important hypoxiainducible gene in recent years. C-peptide inhibits low O 2induced ATP release in erythrocytes as a putative ligand of GPR146, which was consistent with our results (Richards et al., 2014). Ncbi_397391 (MMP2), ncbi_102159047 (FOXC1), ncbi_100738910 (PRRX1), and ncbi_100520318 (TUB) are potentially regulated by novel-m0237-5p. In the present study, MMP2 expression was significantly higher in the LH group than in the LL group, but no significant differences were found between the TH and TL groups. The expression of MMP2 showed a similar tendency to the results of alveolar septum thickness among the four groups, indicating that MMP2 activities may lead to the widening of the alveolar wall and septum as well as alveolar structure damage and collapse of the alveolar space with remarkable fibrosis in Landrace pigs (Tan et al., 2006).

CONCLUSION
The comparisons between Tibetan pigs and Landrace pigs from high or low altitudes revealed genes and regulatory pathways with possible adaptive changes in response to high-altitude hypoxia. We identified several molecular pathways and hypoxia genes showing adaptive changes in the lung, including increased blood circulation and regulation of blood pressure and circulation as well as regulation of HGB concentration and angiogenesis. Integrated analysis of mRNAs and miRNAs demonstrated that a number of hypoxia genes may be regulated by miRNAs and participate in the hypoxic regulation of the lung. For example, novel-m0237-5p may potentially upregulate the expression levels of MMP2, resulting in widened alveolar septum and alveolar structure damage. These results provide a better understanding of the molecular mechanisms regulating the hypoxia response in the lungs of Tibetan pigs and will help to prevent damage to the lungs caused by hypoxia.

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 Ministry of Science and Technology of the People's Republic of China (Approval number: 2006-398). Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
SZ was the overall project leader who provided financial support and experimental conception. YY was involved in data analyses, statistical analyses, language revisions, journal selection, and manuscript submissions and revisions. HY, TJ, and TY contributed to the experimental design and implementation. CG contributed to the supervision and assistance of students in managing animals and collecting and analyzing samples. YL was responsible for the trial implementation, supervision of students collecting and analyzing samples, and manuscript preparation. YC contributed to supervision of sample collection and analysis and manuscript editing. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Natural Science Foundation of China (31760644, 32060730).