Transcriptome and Metabolome Integration Provides New Insights Into the Regulatory Networks of Tibetan Pig Alveolar Type II Epithelial Cells in Response to Hypoxia

Tibetan pigs show a widespread distribution in plateau environments and exhibit striking physiological and phenotypic differences from others pigs for adaptation to hypoxic conditions. However, the regulation of mRNAs and metabolites as well as their functions in the alveolar type II epithelial (ATII) cells of Tibetan pigs remain undefined. Herein, we carried out integrated metabolomic and transcriptomic profiling of ATII cells between Tibetan pigs and Landrace pigs across environments with different oxygen levels to delineate their signature pathways. We observed that the differentially accumulated metabolites (DAMs) and differentially expressed genes (DEGs) profiles displayed marked synergy of hypoxia-related signature pathways in either Tibetan pigs or Landrace pigs. A total of 1,470 DEGs shared between normoxic (TN, ATII cells of Tibetan pigs were cultured under 21% O2; LN, ATII cells of Landrace pigs were cultured under 21% O2) and hypoxic (TL, ATII cells of Tibetan pigs were cultured under 2% O2; LL, ATII cells of Landrace pigs were cultured under 2% O2) groups and 240 DAMs were identified. Functional enrichment assessment indicated that the hypoxia-related genes and metabolites were primarily involved in glycolysis and aldosterone synthesis and secretion. We subsequently constructed an interaction network of mRNAs and metabolites related to hypoxia, such as guanosine-3′, 5′-cyclic monophosphate, Gly-Tyr, and phenylacetylglycine. These results indicated that mitogen-activated protein kinase (MAPK) signaling, aldosterone synthesis and secretion, and differences in the regulation of MCM and adenosine may play vital roles in the better adaptation of Tibetan pigs to hypoxic environments relative to Landrace pigs. This work provides a new perspective and enhances our understanding of mRNAs and metabolites that are activated in response to hypoxia in the ATII cells of Tibetan pigs.

Tibetan pigs show a widespread distribution in plateau environments and exhibit striking physiological and phenotypic differences from others pigs for adaptation to hypoxic conditions. However, the regulation of mRNAs and metabolites as well as their functions in the alveolar type II epithelial (ATII) cells of Tibetan pigs remain undefined. Herein, we carried out integrated metabolomic and transcriptomic profiling of ATII cells between Tibetan pigs and Landrace pigs across environments with different oxygen levels to delineate their signature pathways. We observed that the differentially accumulated metabolites (DAMs) and differentially expressed genes (DEGs) profiles displayed marked synergy of hypoxia-related signature pathways in either Tibetan pigs or Landrace pigs. A total of 1,470 DEGs shared between normoxic (TN, ATII cells of Tibetan pigs were cultured under 21% O 2 ; LN, ATII cells of Landrace pigs were cultured under 21% O 2 ) and hypoxic (TL, ATII cells of Tibetan pigs were cultured under 2% O 2 ; LL, ATII cells of Landrace pigs were cultured under 2% O 2 ) groups and 240 DAMs were identified. Functional enrichment assessment indicated that the hypoxia-related genes and metabolites were primarily involved in glycolysis and aldosterone synthesis and secretion. We subsequently constructed an interaction network of mRNAs and metabolites related to hypoxia, such as guanosine-3′, 5′-cyclic monophosphate, Gly-Tyr, and phenylacetylglycine. These results indicated that mitogen-activated protein kinase (MAPK) signaling, aldosterone synthesis and secretion, and differences in the regulation of MCM and adenosine may play vital roles in the better adaptation of Tibetan pigs to hypoxic environments relative to Landrace pigs. This work provides a new perspective and enhances our understanding of mRNAs and metabolites that are activated in response to hypoxia in the ATII cells of Tibetan pigs.

INTRODUCTION
Tibetan pigs are well adapted to high-altitude environments and primarily live in semi-agricultural and semi-pastoral areas with an average elevation of 2,500-4,300 m on the Qinghai-Tibet Plateau in Southwest China, as an ideal model for investigating the biological mechanisms of hypoxic adaptation (Groenen et al., 2012;Shang et al., 2020). In particular, hypoxia poses a major barrier to life, as pulmonary respiratory function is typically injury by a modest decline in oxygen (Nathan et al., 2019;Grimmer et al., 2021). Gas exchange in the lung occurs within alveoli, which are mainly composed of alveolar type I and II epithelial (ATI and ATII) cells. The lung is a unique physiological environment from a metabolic perspective, as it constitutes the barrier between external air and the pulmonary vasculature in the body (Liu and Summer 2019). Metabolomics is expected to provide a reliable assessment of the physiological status of cells and can reveal change in cell metabolic pathways under hypoxia via the systematic study of small-molecule metabolite profiles (Gunda et al., 2018). At present, lung epithelium is derived from pluripotent stem cells, which have been extensively characterized to discover targets for new therapeutics for hypoxic disease (Herriges and Morrisey, 2014). The hypoxic phenotype is characterized by a shift away from primary reliance on oxidative phosphorylation to enhanced glycolysis for the maintenance of homeostasis (Lottes et al., 2014). ATII cells have a highly oxidative metabolic phenotype; they can consume lactate for oxidative ATP production and synthesize the majority of lipid components as a site of lactate consumption to produce pulmonary surfactants (Lottes et al., 2014). Mammalian targets of the cAMP signaling pathway, hypoxia inducible factor (HIF) pathway, AMP-activated kinase pathway, and many other signaling pathways are involved in mediating cellular metabolic homeostasis and the response to hypoxia (Zhang et al., 2018b;Li et al., 2019;Jones et al., 2020). For example, ATII cells are a heterogeneous cell population that contributes to alveolar maintenance and through metabolic and molecular adaptation in the alveolus, to maintain bioenergetic homeostasis after experiencing insults in a hypoxic environment. These cell play roles in processes including fluid transport and homeostasis (Matalon, 1991), pulmonary surfactant production (Guillot et al., 2013), and may serve as progenitors that undergo significant regeneration and transdifferentiation to repopulate ATI cells (Cardoso and Whitsett, 2008;Hoffman and Ingenito, 2012). Recent studies of lung development under hypoxia have generated novel insights into the function and molecular regulatory pathways of different cell lineages in the lung. Although ATII cells may repair lung damage under hypoxia (Lottes et al., 2014), the effects of O 2 limitation on the mRNA and metabolic pathways necessary to maintain cellular energy in ATII cells have not been studied extensively. This project aimed to reveal the network and intercellular signal transduction pathways that contribute to energy homeostasis by identifying DEGs and specific metabolic processes via 2D analyses of the transcriptome and metabolome using primary ATII cells of Tibetan pigs cultured in normoxic and 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.

In vitro Hypoxia Model in ATII Cells
Lung tissues were collected from healthy newborn male piglets (Tibetan pigs and Landrace pigs) at 7 days of age. The lung tissues were perfused with normal saline to flush out any remaining blood and soaked in PBS buffer containing penicillin and streptomycin. ATII primary cells were isolated as described previously (Wang et al., 2016), with minor modifications. ATII cells were cultured in complete medium (DMEM) at 37°C under normoxic conditions, with 21% O 2 , 5% CO 2 , and 79% N 2 . Hypoxic conditions were achieved by culturing cells in a hypoxic atmosphere containing 2% O 2 , 5% CO 2 , and 98% N 2 under hypoxic conditions (Bactal 2 gaz; Air Liquide, France).

RNA Extraction and Transcriptome Sequencing
ATII cells were collected at 48 h under 21% O 2 (TN, LN) or 2% O 2 (TL, LL) conditions, and three of nine replicates were selected from each group. RNA was extracted using TRIzol reagent kit (Invitrogen, Carlsbad, CA, United States). Second-strand cDNA synthesis was performed with DNA polymerase I, RNase H, dNTPs (dUTP instead of dTTP) and buffer. PCR amplification, sequencing library generation and subsequent sequencing were carried out using the Illumina HiSeqTM 4,000 platform (or other platforms) by Gene Denovo Biotechnology Co. (Guangzhou, China).
Eight mRNAs were randomly selected and quantified using qRT-PCR to validate the accuracy of the sequencing data. The primers used for qPCR were designed and synthesized by Qingke Biological Company (Xi'an, China) (Supplementary Table S1). To normalize sample differences, the β-actin gene served as an internal control, and the 2 −ΔΔCT method was used to calculate relative expression levels. p＜0.05 was considered significant.

Metabolite Extraction
Six of nine replicates were set for each group using the same protocol. Metabolite extraction, endogenous metabolite identification and data processing were carried out by Genedenovo Tech Co., Ltd. (Guangzhou, China) using liquid chromatography-mass spectroscopy (LC-MS). The samples were freeze dried in the same proportions, after which 1,000 µl of methanol (−20°C) was added to the lyophilized powder, and the samples were centrifuged at 12,000 rpm for 10 min at 4°C. A 450 μl aliquot of the supernatant was then concentrated by drying under vacuum, dissolved in 150 μl of a 2-chlorobenzalanine (4 ppm) 80% methanol solution, and filtered through a 0.22 µm membrane. A 20 µl aliquot of each sample was used for quality control (QC) analysis, and the remaining samples were used for LC-MS analysis.

Metabolome Analysis by Liquid Chromatography-Mass Spectroscopy (LC-MS)
UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (ThermoFisher, Germany) coupled with an Orbitrap Q ExactiveTM HF-X mass spectrometer (Thermo Fisher, Germany) in Gene Denovo Co., Ltd. (Guangzhou, China). Samples were injected onto a Hypesil Gold column (100 × 2.1 mm, 1.9 μm) using a 17 min linear gradient at a flow rate of 0.2 ml/min. The eluents for the positive polarity mode were eluent A (0.1% FA in Water) and eluent B (Methanol). The eluents for the negative polarity mode were eluent A (5 mM ammonium acetate, pH 9.0) and eluent B (Methanol). The solvent gradient was set as follows: 2% B, 1.5 min; 2-100% B, 12.0 min; 100% B, 14.0 min; 100-2% B, 14.1 min; 2% B, 17 min. Q ExactiveTM HF-X mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas flow rate of 40 arb and aux gas flow rate of 10 arb.

Processing and Statistical Analysis of Metabolomics Data
The raw data files generated by UHPLC-MS/MS were processed using the compound discoverer 3.1 (CD3.1, Thermo Fisher) to perform peak alignment, peak picking, and quantitation for each metabolite. The main parameters were set as follows: retention time tolerance, 0.2 min; actual mass tolerance, 5 ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, 100,000. Peak intensities were normalized to the total spectral intensity. The normalized data was used to predict the molecular formula based on additive ions, molecular ion peaks and fragment ions. Then peaks were matched with the mzCloud (https://www.mzcloud.org/), mz Vaultand Mass Listdatabase to obtain the accurate qualitative and relative quantitative results. The identified metabolites were analyzed by multidimensional statistical and further to show the distribution of the original data and the classification of variables, including principal component analysis (PCA) was applied in all samples using R package models (http://www.r-project.org/) (Warnes et al., 2007), partial least squares-discriminant analysis (PLS-DA) was applied in comparison groups using R package ropls (http://www.r-project.org/) (Thevenot 2016), and orthogonal partial least squares discriminant analysis (OPLS-DA) was applied in comparison groups using R package models (http://www.r-project.org/) (Warnes et al., 2007). The characteristics of metabolite expression patterns were used to identify the differential abundance of the metabolites according to the criteria of a p value of the T test＜0.05 and a variable importance in projection (VIP) value ≥ 1, which were indicated differential metabolites between different groups. The differential metabolites were annotated and subjected to enrichment analysis according to the KEGG metabolome database. An FDR≤0.05 defined pathways that were significantly enriched in differential metabolites. Pathway over-representation was evaluated by MSEA using the MetaboAnalyst module.

Integrative Analysis of the Metabolome and Transcriptome
We examined and integrated the DEGs and DAMs for overlapping metabolic pathways performed by two way orthogonal PLS (O2PLS) analysis to identify hypoxia related genes and pathways based on both the metabolome and transcriptome datasets. In addition, metabolome and transcriptome data integration was calculated by using pearson correlation coefficients. The pairs of genes and metabolites in common metabolic pathways (with absolute Pearson correlation coefficient >0.995 and p < 0.05) were visualized using Cytoscape (V3.3.0). The top 50 gene and metabolite pairs were selected for heatmap analysis. Additionally, the top 250 pairs of genes and metabolites (with an absolute pearson correlation >0.5) were employed for metabolite transcript network analysis using R graph packages.

Identification of DEGs in ATII Cells
To obtain an overview of the ATII cells transcriptome in normoxic and hypoxic environments, raw data (Q30 > 93.13%, with clean read ratio of >95.61%) were obtained for RNA-Seq (Supplementary Figure S1). Twelve libraries generated an average of 74,968,018 clean paired reads, and 95.61-98.28% of the clean reads were sequentially mapped to the porcine reference genome by Illumina paired-end sequencing technology (Table 1) Table S2 in Supplementary Material S3). A clustering heatmap analysis of these mRNAs showed excellent repeatability and gene expression profiles among the four groups ( Figure 1A). A total of 1,470 DEGs were identified and screened based on the intersection in ATII cells between normoxic (TN, LN) and hypoxic (TL, LL) groups to explore the potential functions of gene responses in hypoxia, which considered as the most interesting candidates ( Figure 2C, Supplementary Material S4). Furthermore, 54 novel genes were identified in the sequencing data. Eight mRNAs were randomly selected and identified using qRT-PCR to validate the accuracy of the sequencing data, which indicated that the qRT-PCR results were consistent with the mRNA-seq data ( Figure 1D).

Enrichment and Functional Annotation of Differential Gene Expression in Response to Hypoxia
According to GO annotation, DEGs in four groups were classified into three major categories including biological process, cellular component and molecular function (Supplementary Figure S2). Comparing 1,470 DEGs of 2% O 2 (TL, LL) groups with the 21% O 2 (TN, LN) groups revealed that biological process was the most enriched category, with a total of 54 terms, followed by the cellular component and molecular function categories. In addition, cellular process, cell, and binding were the most abundant terms in biological process, cellular component and molecular function, respectively ( Figure 2A). Interestingly, 1,025 DEGs between normal oxygen groups and hypoxia groups were mainly enriched in the cellular process term (GO: 0009987) of biological process. GO: 0005634, GO: 0005488, and GO: 0019219 were associated with the nucleus, binding, and regulation of nucleobase-containing compound metabolic processes were most significantly enriched between the normoxia and hypoxia groups. The top 20 pathways with the most significant enrichment among four groups were identified (Supplementary Figure S2). Multiple DEGs shared between the normoxia (TN, LN) and hypoxia groups (TL, LL) were enriched KEGG pathways categories, among which the most significant enriched category was proteoglycans, followed by MAPK signaling and then vascular smooth muscle contraction ( Figure 2B). In addition, signal transduction of environmental information processing, global and overview maps of metabolism, endocrine system of organismal system, cellular community-eukaryotes of cellular process and folding, sorting and degradation of genetic information processing were also the most enriched pathways, respectively.
Interesting, system development, anatomical structure development, and multicellular organism development terms significant enriched by specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the TN and TL groups, meanwhile, intracellular part, intracellular, and intracellular membrane-bounded organelle significantly enriched by specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the LN and LL groups. Wnt signaling pathway, carbon metabolism and glycolysis/gluconeogenesis were most significantly enriched pathways by specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the TN and TL groups, however, systemic lupus erthematosus, necroptosis, and cell cycle were the most significantly enriched pathways by specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the LN and LL groups.

Differentially Accumulated Metabolites Analysis
Nontarget metabolomics was used to profile the potential impact of hypoxic environment on the metabolome of ATII cells. A PCA score plot showed the distributions of origin data, and the samples were separated accordingly (Supplementary Figure  S3Aa, Bb). The PCA score plots of the TN-TL (PC1 36.8%, PC2 31.3%) group in positive (POS) mode were analyzed, and the score plots of PCA based on the TN-TL (PC1 50%, PC2 29.5%) group in negative (NEG) mode were revealed (Supplementary Figure S3 in Supplementary Material S5). Indeed, the PCA results illustrated that TN, TL, LN, and LL were clearly separated according to the corresponding groups. PLS-DA and OPLS-DA revealed that the TN, TL, LN, and LL groups were clearly distinguished from each other, which illustrated that the differences in metabolite accumulation among the four groups were significant (Figure 3). A total of 74 (37 up-and 37 down-regulated) and 32 (8 up-and 24 downregulated) DAMs were identified between the TN and TL groups

Integrated Analyses of Hypoxia-Related Pathways by Metabolome and Transcriptome
To obtain comprehensive insight and explore the correlation between hypoxia related DEGs and DAMs, the DEGs were analyzed in the form of venn diagrams (Figure 1). KEGG Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 13 | Article 812411 assessment indicated that DEGs between normoxic (TN, LN) and hypoxic (TL, LL) groups and DAMs were involved in the MAPK signaling pathway, vascular smooth muscle contraction and cGMP-PKG signaling pathway; specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) and DAMs identified between TN and TL were mainly significantly enriched in huntington disease, metabolic pathways, and alcoholism pathways, while specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) and DAMs between LN and LL were mainly significantly enriched in carbon metabolism and purine metabolism pathways (Supplementary Material S7). The gene-metabolite correlation pairs with absolute pearson correlation coefficient >0.995 and p < 0.05 were used to build hypoxia related regulatory networks. Through the integration of the hypoxia-related genes and metabolites identified in ATII cells during hypoxia, potentially regulatory for mRNA-metabolite pairs were analyzed based on the filtered genes. Three integrated analysis and coexpression networks were constructed for DEGs between normoxic (TN, LN) and hypoxic (TL, LL) groups and DAMs, specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) and DAMs identified during the ATII cells response to hypoxia of Tibetan pigs (TN and TL) and Landrace pigs (LN and LL) under different oxygen levels, and the top 300 relationship pair network diagrams are presented ( Figure 6). Guanosine-3′, 5′-cyclic monophosphate (com-1603-neg), Gly-Tyr (com-471-pos), and phenylacetylglycine (com-43-neg) were significantly correlated with a large number of different transcripts respectively. On the contrary, ncbi-396769 were significantly correlated with a large number of metabolites. Sulfamethazine (com-3927-pos), Guanosine-3′, 5′-cyclic monophosphate (Com-1603-neg), and kynurenic acid (com-3922-pos) were significantly correlated with a large number of specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the TN and TL, meanwhile, 2,3,4,9-Tetrahydro-1H-β-carboline-3-carboxylic acid (com-2935-pos), YMH (com-1241-pos), and D-Raffinose (com-617-pos) were significantly correlated with a large number of specific DEGs (excluding DEGs shared between normoxia and hypoxia groups) between the LN and LL.

DISCUSSION
A large number of metabolome sequences, providing data on important parameters related to hypoxia, were generated some time ago, and their generation accelerated after research on adaptation to hypoxia in native high-altitude species, including Tibetan sheep, yaks, and Tibetan pigs Ma et al., 2019;Du et al., 2021). The morphology of pulmonary artery cast specimens, the values of blood physiological indexes and the expression network of hypoxia-responsive genes were identified in Tibetan pigs compared to Landrace pigs in our previous studies (Yang et al., 2021a;Yang et al., 2021b). It can be assumed that changes in DEGs and DAMs reflect the hypoxic adaptation mechanism by regulating the expression of genes and metabolism to benefit ATII cell regeneration and transdifferentiation. In the present study, we focused on the hypoxia adaptation of ATII cells in pigs and identified 546 genes and metabolites with differential abundance between normal-altitude and high-altitude of Tibetan pigs and Landrace pigs. GO enrichment analysis revealed that hypoxia- related DEGs were associated with the nucleus, binding, and the regulation of nucleobase-containing compound metabolic processes. Moreover, several hypoxia-related signaling pathways, such as the MAPK signaling pathway, vascular smooth muscle contraction and the cGMP-PKG signaling pathway, which may form a complex cascade of responses to promote the regeneration and transdifferentiation of ATII cells, are active in hypoxic environments to protect pulmonary normal function.

Speculation of the Hypoxia Signaling Pathway Affecting the Regulation of ATII Cells
As we all known, HIF-1α protein synthesis can result in the activation of MAPK pathways through interaction with their cognate receptor tyrosine kinases (Shin et al., 2010;Olson and Vliet, 2011). MAPK is a member of the family of enzymes involved in oxygen sensing and plays an important role in the Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 13 | Article 812411 8 mechanism underlying the increased tolerance of the bighead carp heart to acute hypoxia (Zhou et al., 2020). In this study, cAMP (Com_18786_neg) was found to be increased in the hypoxia groups and showed a strong correlation with hypoxiarelated DEGs (above 300); cAMP may activate MAPK signaling to promote cell clone formation and proliferation in alveolar epithelial cells under hypoxia. Similar results have been found in periodontal ligament stem cells (He et al., 2016;Kabiri et al., 2018). Previous studies have shown that the sodium hydrosulfide-mediated inhibition of the cigarette smoke-induced phosphorylation of ERK, JNK and p38 MAPK may be a novel strategy for the treatment of chronic obstructive pulmonary fibrosis through the inhibition of inflammation, epithelial cell injury and apoptosis (Guan et al., 2020;Shologu et al., 2018;Singh et al., 2018). Therefore, we focused our further studies on MAPK signaling. Numerous hypoxia-related DEGs and metabolites of ATII cells identified between the TN and TL groups were enriched in the MAPK pathway; thus, the activation of MAPK pathways may arise under exposure to hypoxia and play critical roles in opposing the inflammatory response and regulating cell proliferation, differentiation, and apoptosis (Dong et al., 2019;Estaras et al., 2021). Through an indepth analysis of the DEGs and the metabolite interaction network, we found that interactors were enriched in the MAPK pathway and could further regulate other pathways to promote the hypoxia-induced transdifferentiation of ATII cells (Xia et al., 2018;Wang et al., 2019). Taken together, these results suggested that MAPK signaling may contribute to suppressing the effect of hypoxia-induced inflammation and alveolar epithelial cell injury.

Effects of Hypoxia on Transcription Expression and Metabolite Formation in ATII Cells
Multiple lines of evidence have demonstrated that cyclic GMP (cGMP) mediates the action of natriuretic peptides (NPs) and nitric oxide (NO) as intracellular second messengers that regulate a broad array of physiological processes (Tsai and Kass, 2009). The PKG1 isoform-specific activation of established substrates leads to a reduction in the cytosolic calcium concentration and/or a decrease in the sensitivity of myofilaments to Ca2+, resulting in smooth muscle relaxation (Adebola and Usha, 2011;Rainer and Kass, 2016). The cGMP-PKG signaling pathway could contribute to the attenuation of hyperlipidemia by Lei-gong-gen formula granules in rats, and this pathway has long been known to be a crucial therapeutic target in injury (Braun et al., 2018;Lan et al., 2020). As expected, we also discovered that a number of hypoxiarelated DEGs and metabolic pathways were significantly enriched in cGMP-PKG signaling pathway by analyzing the results of KEGG pathway enrichment, for example, adenosine 5′-phosphate (com_27274_pos), adenosine (com_190_pos) and MADS-box transcription enhancer factor 2C (MEF2C), that may provide the cellular protection and serve as a potential therapeutic target for protect cell against injury of hypoxia (Mahneva et al., 2019). Likewise, PPARγ antagonizes the hypoxia-induced activation of hepatic stellate cells by cross-mediating the PI3K/AKT and cGMP/PKG signaling pathways (Zhang et al., 2018a). As a result, the higher expression of cGMP observed in ATII cells in the hypoxia groups (LL, TL) may exert physiological action mediated by two forms of cGMP-dependent protein kinase (PKGs), cGMPregulated phosphodiesterases and cGMP-gated cation channels, among which the PKGs might be the primary mediator (Inserte and Garcia-Dorado, 2015;Olivares-González et al., 2016). Aldosterone is secreted from the adrenal cortex and plays a vital role in water, salt and sodium homeostasis, thereby contributing to blood pressure control, which requires tight regulation of its secretion (Bollag, 2014). Aldosterone production by zona glomerulosa cells is enhanced by hypoxia via the promotion of the hydrolysis and accumulation of cholesterol ester in lipid droplets (Yamashita et al., 2020). Consistent with previous studies, we observed that hypoxia related DEGs and DAMs were greatly enriched in aldosterone synthesis and secretion and that the expression of adenosine was significantly higher in the TL group than in the TN group, which may result in various regulatory mechanisms of cAMP and circulating aldosterone levels, especially plasma potassium levels and the renin-angiotensin system (Bollag, 2014).

Differential Regulation of Glycolysis and MCM of Tibetan Pigs Response to Hypoxia
We subsequently investigated the contribution of the components of HIF-1 signaling pathways, which can elevate interstitial pressure and protect cells from hypoxic injury through glucose metabolism, mitochondrial function and cell apoptosis, as potential therapeutic proangiogenic molecules (Balogh et al., 2019;Nagao et al., 2019;Li et al., 2021). The hypoxic metabolic switch could be regulated by hypoxia inducible factor-1α (HIF-1α), which induces glycolytic enzymes at the transcriptional level (Kierans and Taylor, 2021). Likewise, hexokinase (HK) mediates the first glycolytic enzymatic step and is key rate-limiting enzyme in glycolysis. The gene expression levels of HK1 and HK2 were significantly higher in the hypoxia groups (2% O 2 ) than in the normoxia groups (21% O 2 ) in our study, and the same pattern has been observed in Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 13 | Article 812411 pulmonary artery smooth muscle cells (PASMCs) (Chen et al., 2018a). The increase in glycolytic metabolism under hypoxia (TL and LL) may require the O 2 -independent glycolytic pathway to preferentially produce sufficient ATP to satisfy bioenergetic requirements, which profound effects on cellular physiology and hypoxic immunity in ATII cells under hypoxia (Semenza, 2017;Allen et al., 2020). Similarly, our results indicated that the expression of ENO2, the target gene of HIF, was significantly higher in the hypoxia (TL and LL) groups and promoted glycolysis, glucocorticoid resistance, and cell growth (Figure 7) . Previous research showed that Tibetan pigs living in highaltitude environments have lower hemoglobin levels than those living in low-altitude environments, which indicates a blunted erythropoietic response to hypoxic challenge (Yang et al., 2021a). HMOX2, is involved in heme catabolism, showed significantly higher expression in TN than in TL (although the difference between LN and LL was not significant), which may lead to a more efficient breakdown of heme and help maintain a relatively low hemoglobin level in Tibetan pigs in a hypoxic environment; Tibetan researchers have shown similar results (Landry et al., 2016;Yang et al., 2016). MCM proteins play crucial roles in DNA helicase, replication and cell proliferation functions, which could negatively regulate HIF-1 activity, in addition to mutually regulating MCMs (Aparicio et al., 2012;Hubbi et al., 2013;Joshi et al., 2016). Several DEGs, such as MCM2, MCM3 and MCM4, were enriched in the cell cycle pathway and showed significantly higher levels in the LN group than in the LL group, in accord with the previously reported finding that hypoxia can downregulate MCM2-7 mRNA simultaneously in HCT116 cell lines and endothelial cells to decrease proliferation as a fundamental physiological response to hypoxia (Hubbi et al., 2011;Kondo et al., 2017). However, the different between the TN and TL groups was not significant, which may explain why Tibetan pigs adapt to hypoxic environments better than Landrace pigs.

CONCLUSION
In conclusion, we identified several molecular pathways showing adaptive changes in ATII cells through integrated metabolomic and transcriptomic analyses, revealing that MAPK signaling may affect hypoxia-induced inflammation and alveolar epithelial cell injury and that glycolytic metabolism in hypoxic environments may have profound effects on cellular physiology in ATII cells. These findings provide a novel perspective and direction for the future exploration of the regulatory mechanisms of pulmonary development and function in Tibetan pigs.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA778032.

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).

AUTHOR CONTRIBUTIONS
SZ was the overall project leader who provided financial support and experimental conception. YY was involved in FIGURE 7 | HK, CAMP, ENO2 were high regulated and promote the glycolysis, glucocorticoid resistance, and cell growth processing under hypoxia. The pathway was visualized by the adobe illustrator software.
Frontiers in Genetics | www.frontiersin.org January 2022 | Volume 13 | Article 812411 data analyses, statistical analyses, language revisions, journal selection, and manuscript submissions and revisions. HY, XL, and ZW 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 and YR were responsible for the trial implementation, supervision of students collecting and analyzing samples, and manuscript preparation. YC and TJ contributed to supervision of sample collection and analysis and manuscript editing. All authors contributed to the article and approved the submitted version.

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