How Placenta Promotes the Successful Reproduction in High-Altitude Populations: A Transcriptome Comparison between Adaptation and Acclimatization.

As the best adapted high altitude population, Tibetans feature a relatively high offspring survival rate. Genome-wide studies have identified hundreds of candidate SNPs related to high altitude adaptation of Tibetans, although most of them have unknown functional relevance. To explore the mechanisms behind successful reproduction at high altitudes, we compared the placental transcriptomes of Tibetans, sea level Hans (SLHan), and Han immigrants (ImHan). Among the three populations, placentas from ImHan showed a hyperactive gene expression pattern. Their increased activation demonstrates a hypoxic stress response similar to sea level individuals experiencing hypoxic conditions. Unlike ImHan, Tibetan placentas were characterized by the significant up-regulation of placenta-specific genes, and the activation of autophagy and the tricarboxylic acid (TCA) cycle. Certain conserved hypoxia response functions, including the antioxidant system and angiogenesis, were activated in both ImHan and Tibetans, but mediated by different genes. The coherence of specific transcriptome features linked to possible genetic contribution was observed in Tibetans. Furthermore, we identified a novel Tibetan-specific EPAS1 isoform with a partial deletion at exon six, which may be involved in the adaption to hypoxia through the EPAS1-centred gene network in the placenta. Overall, our results show that the placenta grants successful pregnancies in Tibetans by strengthening the natural functions of the placenta itself. On the other hand, the placenta of ImHan was in an inhabiting time-dependent acclimatization process representing a common hypoxic stress response pattern.

Introduction

The Tibetan populations have inhabited the Qinghai-Tibetan plateau for over ten thousand years (Su et al. 2000; Zhao et al. 2009). Among the various selection pressures of high altitude, hypobaric hypoxia presents the most formidable challenges and has shaped the physiological adaptation of Tibetans. In comparison with other populations residing or migrating to similar high altitudes, anthropological surveys have marked several prominent adaptive phenotypes of Tibetans, including, but not limited to, reduced hemoglobin (Hb) concentrations (Wu et al. 2005), lowered risks of chronic mountain sicknesses (Pei et al. 1989; Leon-Velarde et al. 2005), increased capillary density in muscles (Beall 2007), increased activity of antioxidant capacities (Janocha et al. 2017), and decreased mitochondrial volume densities (Kayser et al. 1996; Gelfi et al. 2004; Marconi et al. 2006). Over the past ten years, genome-wide scans identified numerous candidate genes, including EPAS1, PPARA, and EGLN1, which may contribute to the hypoxia adaptation of Tibetans (Beall et al. 2010; Bigham et al. 2010; Simonson et al. 2010; Yi et al. 2010). However, most adaptation-related single nucleotide polymorphisms (SNPs) are located in noncoding regions, suggesting the functional mechanism of these genetic variations via gene expression-related regulations. In particular, EPAS1 (hypoxia-induced factor 2α, HIF2α) has been repeatedly observed as a key regulator of hypoxia adaptation in Tibetans (Beall et al. 2010; Simonson et al. 2010; Yi et al. 2010). Multiple variants of EPAS1 showed the greatest frequency differentiation between Tibetans and Hans, and a comparison study revealed similarities of EPAS1 region in Tibetans to its sequence in ancient Denisovans (Huerta-Sanchez et al. 2014). Interestingly, EPAS1 showed the highest expression in the lung and the placenta, the two tissues with the most extensive oxygen-related exchanges (Fagerberg et al. 2014).

Successful reproduction is pivotal for the population’s survival and growth in a new ecological system. Especially at high altitudes, birth weight decreases by 100 grams every 1,000 meters of elevation (Jensen and Moore 1997; Giussani et al. 2001); therefore, infant birth weight is an important phenotype related to the postnatal survival rate of newborns (Krampl et al. 2000; Giussani et al. 2001). Strikingly, Tibetan newborns showed comparable birth weights to low altitude populations (Zamudio et al. 1993). Furthermore, compared with non-Tibetans born in Lhasa, Tibetan infants are heavier and with fewer premature births (Moore et al. 2001; Yangzom et al. 2008). Recently, a study based on 1,008 indigenous Tibetan women found that a number of SNPs, including two intronic SNPs of CCDC141, were associated with more pregnancies and a higher survival birth rate (Jeong et al. 2018). However, no allele frequency difference between Tibetans and lowlanders was found in these SNPs. Despite the evolutionary and biomedical significance of pregnancy (Brown et al. 2013), the molecular mechanisms behind the successful pregnancies of Tibetans and their underlying genetic bases remain largely unknown.

The placenta is a highly specialized organ supporting the growth and survival of the fetus during pregnancy. In addition to transporting and exchanging gas, nutrients, and wastes, the placenta releases hormones into both maternal and fetal circulations to regulate metabolism, fetal growth, and parturition (Gude et al. 2004). A genetic study revealed the role of placenta transcriptomes in determining postnatal body sizes (Peng et al. 2018). Since transcriptomic variations result from both genomic and epigenomic modifications, and their tissue-specific patterns are affected by particular environments, studying transcriptomes bridges the underlying genomic variations with the resulting phenotypes (Gustincich et al. 2006). For instance, comparing the transcriptome studies on multiple organs between high-altitude animals and their low-land counterpart species revealed multiple variants involved in hypoxia adaptation (Tang et al. 2015; Qi et al. 2019). Unfortunately, the study of placenta transcriptomes in hypoxia has been limited due to the challenge of sample collection.

In this study, we performed multiple comparative analyses of the placenta transcriptomes in samples collected from Tibetan residents in the Qinghai-Tibetan plateau, Han immigrants living in Qinghai, and sea-level Hans living in Beijing. We identified a set of genes and molecular pathways that may promote high-altitude adaptation and acclimatization. We also observed characteristic expression features in native highlanders and immigrants from low regions. Moreover, by integrating the genomic and transcriptomic profiles, we found several cues connecting the genomic alterations and specific placental expression in Tibetans.

Result

The Hyperactive Transcriptome of Immigrant Hans and the Relatively Stable Transcriptome of Tibetans

In total, 91 placental samples and phenotype data were collected from puerperants from three altitudes: 25 Hans from Beijing (40 m) labeled as sea level Hans (SLHan), 19 Hans from Xining (2,000–2,500 m), the capital of Qinghai Province, labeled as immigrant Hans (ImHan), and 47 Tibetans who had their deliveries at Lhasa (≥ 3,650 m), the capital of Tibet Autonomous Region (supplementary table 1, Supplementary Material online; see Material and Methods). After evaluation, 55 mRNA samples from the chorion, amnion, and decidua of 35 individuals were sequenced and yielded 3.2 billion high-quality sequencing reads ranging from 31 to 99 million reads per sample (supplementary table 1, Supplementary Material online).

First, we analyzed the transcriptomic features of different placental tissues among the three groups. Using hierarchical clustering, samples were first grouped by tissue type and then by ethnicity (fig. 1), suggesting differentiated chorion, amnion, and decidua functions. The gender of the newborns appeared not to affect the clustering patterns. Contrarily to our initial expectations, we observed a high similarity between the Tibetan and the SLHan groups but distinct patterns in ImHan for all three tissue types (fig. 1). As all three placental tissues showed consistent population-related features, we chose chorion, the primary nutrition-supplying tissue, to perform the analyses next described. Compared with SLHan and Tibetans, ImHan showed a higher heterogeneity (a lower Pearson correlation coefficients distribution, Pcc, fig. 1) and a significantly higher proportion of highly expressed genes (fig. 1, PImHan-SLHan = 1e-07, PImHan-Tibetan = 6.7e-17). In view of the difference from the relative stable transcriptomes in Tibetans and SLHan, we describe the gene expression of ImHan placenta as hyperactive, implying that these individuals are in the process of acclimatization to high altitude. On the other hand, the similarity in the placental transcriptomic patterns between Tibetans and SLHan may indicate each group’s adaptation or reproduction suitability to their habitual altitudes, that is, high land for Tibetans and sea level for Hans (See Discussion).

Fig. 1.

Diverse expressions of placental genes in Tibetans, ImHan, and SLHan. (A) Hierarchical clustering of all samples. Notably, the dominant clustering was categorized by tissues: amnion (purple), chorion (ivory), and decidua (light blue); secondarily, samples were grouped by populations; the sex of the fetus had no effect (light blue for boy and pink for girl) on the clustering results, although certain sampling biases may have influenced the clustering patterns in amnion and decidua due to more female samples in the ImHan group. Under each type of tissue, ImHan (orange) show a distinct expression pattern from both SLHan (green) and Tibetans (blue). (B) Lower global Pearson correlation coefficients (Pcc) in ImHan, corresponding to a higher heterogeneity than in SLHan and Tibetans. (C) A significantly higher proportion of highly expressed genes was found in ImHan than in SLHan and Tibetans. (D) Hierarchical clustering of the ImHan generations. High similarity of the transcriptomes was seen in the residents of multiple generations (G ≥ 3, upper in pink). In contrast, significant heterogeneity was found among the first and second generations (G1 and G2, lower in gray). (E) Distinct transcriptome patterns of native and immigrant highlanders. The left panel shows the expression heat map generated by STEM across the three altitudes’ populations. Inflammatory and other immune-related responses were significantly enriched in ImHan (up-right, top 5 by P-values), whereas female pregnancy and related transporting genes were highly expressed in Tibetans (bottom-right, top 5 by P-values).

The immigration of the Han people from low-lands to Qinghai province, the east region close to the plateau border, started during the Qin dynasty (over 2,000 years ago) and continued to modern times with several waves. To explore whether hypoxia acclimatization among ImHan is inhabitation time-dependent, we divided ImHan into first-generation, second-generation, and multi-generation (≥ 3 generations) based on the history of puerperant families. As illustrated in fig. 1D, the hierarchical clustering of the global genes well separated the multi-generation individuals (upper group) from the first- and second-generation ones (lower group). Furthermore, despite the similarity among the transcriptomes of the multi-generation samples, significant heterogeneity was found among the individuals of the first- and the second-generation groups (supplementary fig. 1, Supplementary Material online). This finding suggests their drastic responses to hypoxic conditions during pregnancy in these recent immigrants. Such a diverse pattern (i.e., the longer residing history at high altitude, the more stable (less active) placenta transcriptomes) is further evidence that, under hypoxic stress, the reproduction of immigrants from the sea level has been in the process of acclimatization. Also, the higher diversity and hyperactive transcriptomes of the first two generations’ individuals were most likely due to their immediate hypoxia responses.

We also conducted an integrated haplotype score (iHS) test. A larger number of high-iHS loci (iHS > 3, 1,226 vs. 948) and a significantly higher iHS value (P = 3.54e-5, t-test) were seen in the multi-generation ImHan (6 samples) compared with SLHan (10 samples) (supplementary fig. 2, Supplementary Material online), implying an ongoing positive selection process within ImHan. In addition, a total of 391 multi-generation ImHan-specific high-iHS loci were identified (annotated to 342 genes) (supplementary table 2, Supplementary Material online), although no significant functional enrichment (false positive ratio, FDR < 0.1) or influence on the placenta transcriptome pattern (supplementary fig. 2, Supplementary Material online) was observed. These data raised the possibility that the adaptation in the ImHan genomes also occurred accompanying their acclimatization.

To further explore the transcriptome profiles related to the inhabiting altitude, we clustered the expressed genes across the three populations using Short Time-Series Expression Miner (STEM) (Ernst and Bar-Joseph 2006). As shown in the heat map in fig. 1E, along with the increase of the altitude, only a limited number of genes (75) showed their expressions up- or down-regulated. Among these, a reported hypoxia-inducible gene, ANGPTL4, was up-regulated, suggesting a common regulation of the metabolism-related hormone upon hypoxia in all populations (Murata et al. 2009). In contrast to the altitude change, the majority of genes (1,731) were exclusively up-regulated in ImHan or Tibetans, suggesting a dramatic difference in the placental functions upon hypoxia between the two groups. Highly expressed genes in ImHan were significantly overrepresented in inflammatory and other immune-related responses (P = 3.93e-14, fig. 1, up-right). However, the functional overrepresentation analysis (ORA) of the most highly expressed genes in Tibetan placentas highlighted the biological process of female pregnancy (P = 4.49e-14, fig. 1, bottom-right). Overall, such diverse processes upon hypoxic stress during the pregnancies of ImHan and Tibetans may reflect the differences between acclimatization and adaptation.

The Overrepresentation of the Hypoxia Response in ImHan and the Up-Regulation of the Reproduction-Related Genes in Tibetans

We identified 4,293 differentially expressed genes (DEGs) among the three sample groups, accounting for 9.9% of the total gene expression. To reveal genes related to high altitude adaptation in Tibetans and hypoxia acclimatization in ImHan, we performed pairwise comparisons of these two groups with SLHan. This resulted in 1,049 DEGs between SLHan and ImHan (793 up- and 256 down-regulated in ImHan, fig. 2) and 1,487 DEGs between SLHan and Tibetans (1,119 up- and 366 down-regulated in Tibetans, fig. 2). Notably, the gene up-regulation was dominant in two highlander groups compared with SLHan, indicating that both ImHan and Tibetans undergo hypoxic stress in their environments. However, the expression patterns of these two highlander groups appeared to be quite distinct, as 3,891 DEGs were found between them (2,050 up- and 1,841 down-regulated in ImHan, fig. 2).

Fig. 2.

DEGs among SLHan, ImHan, and Tibetans and functional comparison of the most significant DEGs. (A) Numbers of DEGs in pairwise comparison. (B) Significantly overrepresented biological processes (cutoff at top 5) and pathways of the up-regulated DEGs in SLHan (olive), ImHan (orange), and Tibetans (purple). (C) Progressively up-regulated genes upon hypoxia treatment in the HUVECs of SLHan. (D) Venn diagram showing the number of up-regulated genes and overrepresented pathways and processes, generated by kKEGG and GO ORA, in ImHan and HUVECs. (E) DEGs Comparison of Tibetans/SLHan (top), ImHan/SLHan (left), and Tibetans/ImHan (right). (F) Functional enrichment of the Tibetan-specific genes. (G) Comparison between the general placenta’s highly expressed genes and the Tibetan-specific DEGs in chorion (1,053, see text). Most of the 26 placental-specific genes were up-regulated only in Tibetans (bottom). (H) The immunohistochemistry staining of CGB, CSHL1, and PSG6 in syncytiotrophoblast illustrates their up-regulation in Tibetans.

Next, we conducted gene ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway ORA using these DEGs. Similar to the bulk expression pattern (fig. 1), the highest number of up-regulated genes for the enriched biological processes and pathways were also seen in ImHan. The most significant entries (adjusted P < 0.05) included inflammatory response, extracellular matrix organization, cell adhesion, ribosome, and angiogenesis (35 in total) (supplementary table 3, Supplementary Material online, fig. 2). As some of these pathways were reported to be directly activated by hypoxia (Kwon et al. 2017), we suspected that the acclimatization of ImHan may share the same hypoxia-responsive pathways as SLHan experiencing hypoxia. To verify this hypothesis, we treated the umbilical vein endothelial cells (HUVECs) derived from one of our SLHan individuals with 1% oxygen; then, we performed mRNA sequencing of the samples in a time series (0, 6, 12, 22, 34, and 48 hours). We first looked at the continuously up-regulated genes along the low oxygen treatment time (fig. 2), as they were likely triggered by the hypoxic environment; 17.5% of these transcripts were also up-regulated in ImHan (fig. 2). In particular, of the eight biological processes and pathways activated during hypoxia in the SLHan HUVECs, six were also identified in ImHan (supplementary table 3, Supplementary Material online, fig. 2). These results suggest that common hypoxia response reactions were shared by the ImHan acclimatization and the SLHan under hypoxic conditions. In addition, the inflammatory and other immune-related pathways were found in the comparisons between ImHan and Tibetans as well as between ImHan and SLHan, suggesting a strong stress response to high altitude conditions occurring during Han immigrants’ reproduction.

Notably, the main biological process overrepresented by Tibetan placental DEGs was female pregnancy, which was the only item identified in the Tibetan/ImHan comparison and the top one in the Tibetan/SLHan comparison (fig. 2). Besides, L-amino acid transport, regulation of Guanosine triphosphatase (GTPase) activity, and cell migration were also up-regulated in Tibetans compared with SLHan (fig. 2). To further explore the genetic adaptation of Tibetans, we first defined the Tibetan-specific DEGs as the genes with no expression difference between SLHan and ImHan, but with significant differences between Tibetans and the two Han groups. After such filtrations, 1,053 Tibetan-specific DEGs were found (fig. 2, supplementary table 4, Supplementary Material online). By using ORA, most of these DEGs were found to be involved in signal transduction, adrenergic signaling, female pregnancy, and transport (fig. 2  supplementary table 4, Supplementary Material online). Interestingly, based on their transcription yields, all the highest expressed Tibetan-specific DEGs were categorized to placenta predominately expressing hormones (supplementary table 4, Supplementary Material online). These hormones include the chorionic somatomammotropin hormone (such as CSH1, CSH2, CSHL1, GH1, and GH2), glycoprotein hormones, alpha polypeptide (CGA), and the pregnancy-specific beta-1-glycoproteins family. Therefore, we further checked 77 genes (supplementary table 5, Supplementary Material online) whose transcripts were at least five times more in the placenta than in any other tissues, based on the multi-tissue RNA-seq data from the HPA project (Uhlen et al. 2015). Among these, 26 were Tibetan-specific DEGs (fig. 3 upper panel, P = 3.06e-15), and 21 were up-regulated in Tibetans (fig. 3, bottom panel). Furthermore, we examined five of the 21 genes at the protein level by immunohistochemistry. Three of these five genes, CGB, CSH1, and PSG7, showed increased staining in Tibetan placental chorions (fig. 2), while the expressions of TAC3 and KISS1 were similar in Tibetans and SLHan (supplementary fig. 3, Supplementary Material online). Based on the above observations, we conclude that, for the reproduction of Tibetans, a major adaptation strategy appears to be the rise of the placenta-related hormones promoting fetal growth and development in a hypoxic environment (See Discussion).

Fig. 3.

Placental expression of EPSA1 in Tibetans. (A) EPAS1 and its interacting genes in the PPI network. A total of 52 genes were extracted from STRING with high confidence (≥ 0.9) as the interacting genes of EPAS1; the edges with |Pcc| ≥ 0.4 are displayed in the network. (B) Pcc distribution between EPAS1 and its interacting genes in the PPI network shows significantly higher Pcc values in Tibetans. (C) Co-factors and regulators of EPAS1 showed strong co-expression with EPAS1 in Tibetans. (D) Most EPAS1 targeted genes were significantly up-regulated in Tibetans. Red indicates the genes specifically expressed in the placenta. (E) Placenta-specific genes showing strong co-expression with EPAS1 in Tibetans. Edges with |Pcc| ≥ 0.4 are displayed in the network. The genes in pink were up-regulated and the genes in green were down-regulated in Tibetans. (F) Splice junction patterns between exons five and six. Light red lines indicate various lengths of splicing. A novel junction 56 bp downstream of the wild-type splicing site of exon six was identified and observed mainly in Tibetans. Stars indicate each of four randomly picked samples with rs150877473 sequenced, gray stars indicate CC, green stars show CG, and red stars illustrate the Tibetan-specific GG genotype (supplementary table 8, Supplementary Material online). (G) Alternative splicing of exon six results in a novel isoform of EPAS1 in Tibetans. (H) Consequence of the alternative splicing on the protein structure of EPAS1. The upper diagram shows the protein structure of wild-type EPAS1, and the bottom diagram predicts the truncated EPAS1 protein structure caused by the alternative splicing. (I) Validation of the alternative splicing event in the HUVECs of SLHan.

Despite the large number of Tibetan-specific DEGs (1,053), only 575 genes could be classified as ImHan-specific DEGs using the same categorizing strategy (fig. 2, supplementary table 4, Supplementary Material online). Specifically, these DEGs were overrepresented in inflammatory response, nuclear-transcribed mRNA catabolic process, and translational initiation-associated functions. The top highly expressed ImHan-specific DEGs were involved in the ribosome-associated pathways (supplementary table 4, Supplementary Material online). However, for the above 77 genes highly expressed in placenta, only two were found among these transcripts. This finding suggests that, in response to the hypoxic stress during pregnancy, ImHan did not significantly increase the expression of placenta-specific genes the way we observed in Tibetans. The distinct up-regulated pathways between the two highlander groups suggest certain intrinsic physiological differences under low oxygen conditions during pregnancy. Specifically, Tibetans enhanced the placental function of female pregnancy while ImHan activated the hypoxia response, which induced various common pathways rather than the tissue-specific functions during their acclimatization to high altitude.

Factors Contributing to the Successful Reproduction of Tibetans at High Altitude

To reveal the cause of the up-regulation pattern observed in Tibetan placentas, we calculated the co-expression indexes between the genes specifically expressed in the placenta and 1,788 transcription factors (TFs) from the DBD database (Wilson et al. 2008). This analysis resulted in 43 TFs showing strong correlations with the placental expression in all three populations (supplementary fig. 4, Supplementary Material online, average Pcc > 0.6). Interestingly, almost all of these TFs were up-regulated in Tibetans (supplementary fig. 4, Supplementary Material online), suggesting that the up-regulation of placental genes was caused by these up-regulated TFs.

Considering the Tibetan genetic background, we then reviewed the correlations of their genomic variants with the regulation of these TFs. By integrating our in-house collection of Tibetan whole-genome sequencing data and eQTLs from GTEx v7 (Consortium 2013), we found that at least 12 TFs were possibly regulated by the genetic polymorphisms of Tibetans (supplementary table 6, Supplementary Material online). For instance, POU2F3, a co-expressed TF that regulates placenta-secreted lactogen, was significantly up-regulated in Tibetans compared with the Han population (PSLHan-ImHan = 5e-05, FCSLHan-ImHan = 3.11; PImHan-Tibetan = 5e-05, FCImHan-Tibetan = 10.67). Additionally, we found that 13 eQTLs of POU2F3 showed significant allele frequency differences between Tibetans and Hans. Regarding the high expression alleles among these 13 eQTLs, 11 showed increased frequencies in Tibetans. Among these 11 SNPs, which ranged from 7 to 20 kb away from TSS of POU2F3, nine were localized in H3K4me1 modified regions (enriched at active and primed enhancers), specifically in the placenta as shown by the Roadmap Epigenomics (http://www.roadmapepigenomics.org). These results suggest that certain genes specifically expressed in the placenta are up-regulated in Tibetans by their genetic variations.

As POU2F3 is a well-known introgressed gene from Neanderthals (Vernot and Akey 2014), we further compared the above 13 SNPs with the Neandertal genomes. As shown in supplementary table 7, Supplementary Material online, none of the above 11 high-expression-related alleles was found in Neanderthal genomes. Furthermore, the remaining two SNPs were also common polymorphisms in Neanderthals and Africans. Therefore, these high-expression-related alleles upstream of POU2F3 appeared not to be introgressed from Neanderthals.

As a key regulator in the hypoxia adaptation of Tibetans, we next analyzed the co-expression pattern of EPAS1 (also named Hypoxia-Inducible Factor 2 Alpha, HIF2A) and its interacting genes in the placental tissues of the three populations. Consistently with previous findings (Peng et al. 2017), the expressions of EPAS1 and its isoforms in the placental chorion showed no significant difference among the three populations (supplementary fig. 5, Supplementary Material online). However, as shown in fig. 3A and 3B, EPAS1 displayed a strong positive correlation with its partner genes in Tibetans, whereas only a weak correlation was observed in SLHan and ImHan. Furthermore, regarding the known EPAS1 interacting genes, including EPAS1 regulators and coactivators of downstream pathways (EGLN2, VHL, UBC, HIF1AN, EP300, ARNT, and CREBBP) (fig. 3), the co-expression patterns of EPAS1 with these factors were significantly different between the placental tissues of Tibetans and their counterparts from the two Han groups. Notably, Hypoxia-Inducible Factor 1 Alpha (HIF1A) showed only a weak relationship with these interactors, and no significant difference was found among the populations (fig. 3, bottom), which suggests that it is HIF2A (EPAS1) rather than HIF1A that acts as the placental hypoxia-inducible factor during Tibetans’ pregnancies.

Among the target genes of EPAS1, we first identified 1,986 putative targets from the reported EPAS1 ChIP-seq data across four cell lines (human renal proximal cell line HKC8, renal cell carcinoma cell line RCC4, human hepatocellular carcinoma cell line HepG2, and human breast cancer cell line MCF-7) (Schodel et al. 2011; Smythies et al. 2019). Among these genes, 60 highly expressed targets largely varied between Tibetans and SLHan (average fragments per kilobase of transcript per million mapped reads, FPKM > 10 and fold change > 2); 53 of them were up-regulated while seven were down-regulated in Tibetans (fig. 3). Among the 53 up-regulated genes, six belonged to the above 43 TFs with a strong correlation to placental expression (supplementary fig. 4, Supplementary Material online, supplementary table 6, Supplementary Material online); four of them, PSG10P, PSG11, SERPINE2, and GPC3 (labeled in red in fig. 3), were placenta-specific, highlighting the particular regulation role of EPAS1 in Tibetan placentas. This finding was further verified by the Tibetan-specific positive correlation between EPAS1 and the 41 DEGs among the 77 placental highly expressing genes (fig. 3E, supplementary table 5, Supplementary Material online).

Notably, a new 3′ splicing site in the sixth exon of EPAS1 was identified in most Tibetan placental samples (fig. 3), resulting in the skipping of the first 56 bp of this exon. Consequently, a frameshift was generated, leading to a premature termination codon 48 bp downstream of the new splice site in exon six. When translated, this novel isoform was predicted to encode a truncated EPAS1 (207 vs. 819 amino acids) with 16 new amino acids at the C-terminal (fig. 3). As illustrated in fig. 3H, the functional elements of EPAS1, including the basic helix-loop-helix (bHLH) domain and one of the per-Arnt-SIM (PAS-A) domains, were retained in this novel peptide. However, PAS-B, the oxygen-dependent degradation domain (ODD), the N-terminal and the C-terminal transactivation domains (N-TADs and C-TADs, respectively) were deleted. Remarkably, our in-house genome data identified one Tibetan-specific SNP (rs150877473) located just five bp upstream of the sixth exon. As shown in fig. 3G, the RBPmap predicted that this allele might attenuate the splicing acceptor of exon six by affecting its binding by the RNA binding protein SRSF3 (Paz et al. 2014). To validate the effect of this SNP on alternative splicing, we sequenced this region in randomly picked four Tibetan and four Han samples (stars in fig. 3, supplementary table 8, Supplementary Material online). As expected, all the samples with the G allele (the Tibetan dominant one) displayed the truncated splicing form. In the only Tibetan sample with the C allele (red star in fig. 3), no such splicing junction was detected. Similarly, the only SLHan sample with the G allele (green star in fig. 3) showed limited reads supporting this splicing junction. We further validated this alternative splicing in the HUVECs of SLHan and did not detect this splicing junction under the low oxygen treatment (fig. 3).

Taken together, the above results suggest the role of genomic variations in the expression profiles of Tibetan placentas. In particular, EPAS1 appears to be a key factor involved in the regulation of placenta-specific genes, which in turn promote the successful pregnancy of Tibetans at high altitudes.

Pathways Promoting the Successful Reproduction of Highlanders

To gain a comprehensive understanding of the successful reproduction of both Tibetans and ImHan, we next analyzed known hypoxia-targeted pathways, especially those correlated to angiogenesis, metabolism, and oxidative stress (Muz et al. 2015).

The More Efficient Aerobic Metabolism of Tibetans

As listed in supplementary table 8, Supplementary Material online, our gene set enrichment analysis (GSEA) resulted in several activated pathways that shared the same upregulated core gene sets, including Huntington’s disease, oxidative phosphorylation, Alzheimer’s disease, and Parkinson’s disease, in both ImHan and Tibetans. The products of these core genes are located in the electron transport chain (supplementary fig. 6, Supplementary Material online). This up-regulation may reflect the higher efficiency of electron transport and energy generation in the highlanders under hypoxia.

Notably, our GSEA results also indicated the up-regulation of the TCA cycle in Tibetans compared with the two Han groups, as shown by the higher expression of the TCA rate-limiting enzymes including citrate synthase, isocitrate dehydrogenase, and oxoglutarate dehydrogenase (fig. 4). No TCA up-regulation was seen in ImHan, even though, we found a slight up-regulation of the key enzymes involved in glycolysis, such as hexokinase, pyruvate kinase, and phosphofructokinase (supplementary fig. 6 and 6, Supplementary Material online), and an up-regulation of the enzymes inhibiting the access to TCA such as pyruvate dehydrogenase kinase (PDK) and lactate dehydrogenase A (LDHA) (supplementary fig. 6, Supplementary Material online). These findings suggest higher glycolysis levels during pregnancy in Han immigrants, consistently with previous results showing the energy metabolism shift from oxidative phosphorylation to glycolysis in sea level individuals under hypoxia (Kim et al. 2006). The level of placental glycolysis in Tibetans was similar to SLHan, as exemplified by their comparable expressions of PDK and LDHA (supplementary fig. 6, Supplementary Material online). Moreover, Tibetans may accelerate the transfer from glycolysis to TCA by up-regulating pyruvate dehydrogenase E1 (PDHA1), the primary link between the glycolysis and the TCA cycle (supplementary fig. 6, Supplementary Material online). These results show that, despite the slight up-regulation of glycolysis in ImHan under hypoxia, Tibetans may maintain their efficient energy production by up-regulating the TCA cycle and accelerating the transfer from glycolysis to TCA.

Fig. 4.

Differentially activated pathways in highlanders. (A) Activation of the TCA cycle in Tibetans only (left) and activated oxidative phosphorylation in both ImHan and Tibetans (right). Several rate-limiting enzymes from the TCA cycle were up-regulated in Tibetans. (B) Activation of the antioxidant-associated genes in ImHan and Tibetans. The upper row shows the expression of three genes as ROS markers, and the lower row lists the expression of three antioxidant-associated genes. (C) Autophagy was activated in Tibetans only, and proteasome was activated in both ImHan and Tibetans. (D) Highly expressed autophagy-associated genes in Tibetans. The expression of three autophagy marker genes is indicated in the upper row. The lower row shows three autophagy genes mediated by calcium signaling that were highly expressed in Tibetans. (E) A proposed model of the signaling cascades granting successful reproductions in two groups of highlanders. These cascades are illustrated in the bio-processes of energy metabolism (green square), angiogenesis (pink square), and antioxidant-related pathways (blue square, see text). Significantly up-regulated genes and pathways are highlighted with a blue (Tibetans only), pink (ImHan only), or purple (ImHan and Tibetans) frame. The gray frame shows the pathway without significant expression changes. Asterisks mark the genes and pathways that contain Tibetan-specific genetic variants as shown by previous studies as well as our recent results (data not shown).

Enhanced Angiogenesis in Highlanders

An overrepresentation of angiogenesis genes was observed in ImHan (supplementary table 3, Supplementary Material online). For example, vascular endothelial growth factor A (VEGFA) was more expressed in ImHan than in Tibetans and SLHan. However, when we selectively analyzed the expression variability of 274 hypoxia-related genes (including angiogenesis-associated ones) (Simonson et al. 2010), angiogenic signals were also observed in Tibetans. Compared with SLHan, 26 and 35 hypoxia-related genes showed significantly different expressions in ImHan and Tibetans, respectively (supplementary table 3, Supplementary Material online). Of these DEGs, angiogenesis-associated genes were overrepresented in both ImHan (ADM, HMOX1, IL1B, PGF, THBS1, and VEGFA) and Tibetans (ADM, ENG, NOS3, PGF, and THBS1). Interestingly, except for the up-regulation of ADM in both highland populations, different angiogenic genes and regulations were observed in ImHan and Tibetans. For instance, PGF was up-regulated and THBS1 was down-regulated in Tibetans, whereas the regulation of these two genes in ImHan was the opposite (supplementary table 3, Supplementary Material online). Therefore, ImHan and Tibetans appear to increase different angiogenesis-related pathways to support the oxygen supply for the embryo development in hypoxic conditions. In particular, consistently with the high expression of placenta-related genes in Tibetans, the up-regulation of the placenta-specific angiogenic factor PGF further suggests the tissue-specific angiogenesis in native highlanders (see Discussion).

Enhancement of the Antioxidant System in Tibetans

As illustrated in the upper row of fig. 4B, three markers of reactive oxygen species (ROS), HSB1 (hsp29), HSP90B1 (hsp90), and HSP90AA1, were up-regulated in both ImHan and Tibetans. However, these genes were most highly expressed in Tibetans, suggesting that highlanders undergo oxidative stress during pregnancy, especially Tibetans. Two types of up-regulation antagonizing the oxidative state were also observed in these high-altitude residents, including antioxidant system components, such as SOD1, SOD3, and GPX4 (fig. 4, lower row) and the proteasome pathway (fig. 4, right), which acts as the second antioxidant system to degrade damaged, misfolded, unfolded, or pathogenic proteins caused by accumulated ROS (Lefaki et al. 2017).

Notably, Tibetans, unlike ImHan, appeared to have a third strategy to cope with hypoxia: the strong activation of autophagy, as demonstrated by both our GSEA results (supplementary table 9, Supplementary Material online; fig.4C, left) and the significantly higher expression of the autophagy markers MAP1LC3B (LC3II) and SQSTM1 (p62) (fig. 4). Among the various upstream pathways of autophagy, the dramatic up-regulation of CALM1, DAPK1, BECN1, and BCL2 in Tibetans (fig.4D) implies the activation of autophagy possibly initiated by calcium signaling. Furthermore, our genomic data showed significant allele frequency differences between Tibetans and Hans in multiple genes and loci of the calcium-associated pathway (data not shown), suggesting the possibility that certain genetic variations are involved in the Tibetan-specific activation of autophagy.

Unfolded Protein Response in Tibetans

Altered redox homeostasis in the endoplasmic reticulum (ER) may lead to ER stress dictating the fate of cells. In our Tibetan placental samples, we observed the higher expression of key genes involved in ER stress, such as ERN1, XBP1, and HSPA5 (fig. 4; supplementary table 4, Supplementary Material online, supplementary fig. 7, Supplementary Material online). Moreover, we observed the activation of specific processes in response to ER stress. One such process was the autophagy pathway, which may serve as an adaptive response to ER stress (Rashid et al. 2015) and, as mentioned above, it was uniquely activated in Tibetans (fig. 4). ER stress may also be reduced by the unfolded protein response (UPR) network, which was initially activated by ER stress sensors PERK, IRE1α, and ATF6 (Hetz and Papa 2018). The key genes participating in UPR-related pathways, including ATF6, PERK, and ATF4, were primarily up-regulated in Tibetans (supplementary fig. 7, Supplementary Material online). One of the downstream processes of UPR is the N-glycan biosynthesis, which was up-regulated in Tibetans with the most significant P-value in our entire dataset (supplementary table 9, Supplementary Material online, supplementary fig. 7, Supplementary Material online). Moreover, key genes of the JNK signaling, which is downstream of UPR and can modulate autophagy, were also up-regulated in Tibetans (supplementary fig. 7, Supplementary Material online). No significant ER stress-related up-regulation was seen in ImHan, possibly due to the relatively lower altitude of the ImHan habitancy region (2,000–2,500 m).

Integrating the results described above, we propose a possible model showing the signal cascades of highlanders for coping with hypoxia during pregnancy. As shown in fig. 4E, the up-regulation of the TCA cycle in Tibetans and, possibly, the higher glycolysis activity in ImHan lead to the higher energy production of these highlanders. The accumulated ROS activate the antioxidant enzymes and the proteasome pathway in both ImHan and Tibetans. Moreover, the activation of autophagy and N-glycan biosynthesis in Tibetans implies that native highlanders hold another energy-saving strategy: a higher efficiency in recycling the building blocks under hypoxia, promoting cell survival rather than apoptosis (as indicated by the up-regulation of the anti-apoptotic BCL2). Furthermore, the angiogenesis pathway was upregulated in both ImHan and Tibetans, although with rather different regulation mechanisms. The fine-tuning of energy production and utilization, as well as the high efficiency of energy and nutrition transport in the placenta, allow the highlanders, especially Tibetans, to thrive in hypoxic environments with successful reproductions. Compared to ImHan, Tibetans appear to possess certain unique placental regulations in autophagy and angiogenesis, which are largely correlated to their specific genetic polymorphisms.

Discussion

Our comparative study of the placenta transcriptomes from SLHan, ImHan, and Tibetans revealed several features of the native and/or immigrant highlanders, allowing them to adapt to hypoxic environments during reproduction.

Immigrants of Han Are Still in the Process of Acclimatization

Compared with the other two groups, the Han immigrants in Xining have a hyperactive transcriptome, with a global higher gene expression and a greater heterogeneity among individuals (figs 1B and 1C). Their main enriched biological processes such as inflammatory response, ECM, cell adhesion, and angiogenesis, which are shared with sea level residents when they are confronting hypoxic stress (fig. 2B and 2C; supplementary table 3, Supplementary Material online), appear to be mediated by HIF1A (Hartmann et al. 2000; Norman et al. 2000; Semenza 2000, 2010; Solaini et al. 2010; Eltzschig and Carmeliet 2011; Krock et al. 2011; Gilkes et al. 2014). The heterogeneity of ImHan individuals is most likely due to their different family histories, as Han people started to move from the central China to Qinghai 2,000 years ago with several waves of immigration. Except for the sample donors who lived at high altitudes for a few generations only, the detailed immigration times for most of the families are unclear. Therefore, we divided our samples into 1–2 and ≥ 3 generations. The pathways enriched in ImHan and the more drastic response to hypoxic conditions of the first two generations immigrants (fig. 1, supplementary fig. 1, Supplementary Material online) suggest that the Han population from this region of two to three thousand meters (See Materials and Methods) is still in the process of acclimatization. The iHS distribution towards high scores within multi-generation ImHan also suggests a possible ongoing selection process (supplementary fig. 2, Supplementary Material online). No clear functional influence of high-iHS loci was pinpointed, possibly due to several difficulties. One was the relatively small sample size of multi-generation ImHan. Another difficulty in revealing more reliable signals is the complexity of the migration times of multi-generation individuals whose families moved to high altitudes tens to thousands of years ago. Finally, even for the earliest migrants, the selection power would still be subtle, given the scale of the migration for thousands of years. Contrary to our initial expectations for a more significant difference between Tibetans and SLHan, our observation for the similarity in gene expression between these two groups suggests that both populations have lived in their adapted conditions (fig. 1). Certainly, placental samples from the Han population in Lhasa would have been better suited to our analysis. However, it is hard to collect these samples at such high altitudes because most Han puerperants deliver at lower regions due to the higher pregnancy risks and the hypoxia-associated diseases.

More Consumable Nutrients in the Tibetan Fetal Development

Our results suggest that during fetal development, the placentas of Tibetans provide more nutrients to their fetuses based on three features. 1) Tibetan placentas may have a more extensive vasculature. The vascular network delivers oxygen and nutrients, and hypoxia may induce angiogenesis to increase the number of capillaries (Ali et al. 1996; Krock et al. 2011). The exchange of oxygen and nutrients between mother and fetus occurs in the placenta, and a modified capillary branching pattern was observed in residents living at high altitudes (Ali et al. 1996; Krock et al. 2011). Our results showed that general factors for vasodilation and angiogenesis, including NOS3 and its regulator ENG, and placental-specific angiogenesis factors, such as PGF and PSGs, significantly increased in Tibetans (fig. 4; supplementary table 4, Supplementary Material online). Such up-regulation suggests a richer vasculature in Tibetan placentas, allowing a more efficient oxygen and nutrients delivery. This observation is also consistent with the finding that Tibetans have a higher muscle capillary density (Beall 2007). 2) The nutrient transportation system of Tibetans may be more efficient in transferring glucose to the fetus. This appears to be achieved by the up-regulation of glucose transporter type 1 (SLC2A1) and glucose transporter type 11 (SLC2A11) (supplementary table 4, Supplementary Material online) in the placentas of Tibetans, whereas we did not observe any difference in the nutrient transportation factors between ImHan and SLHan. 3) The possibly elevated glucose usage in Tibetan fetuses. Generally speaking, placental growth hormone and chorionic somatomammotropin hormones induce peripheral insulin resistance (Economides et al. 1989; Barbour et al. 2002; Mannik et al. 2012), which consequently results in elevated circulating glucose for the fetus by reducing maternal glucose utilization (Semenza 2010). In our data, we observed multiple signals from the insulin-related pathways in Tibetans. The first was the significantly higher expression of both placental growth hormone and chorionic somatomammotropin hormones (CSH1, CSH2, and CSHL1) in the Tibetan placental samples (supplementary table 4, Supplementary Material online). Furthermore, insulin resistance appears to be more active in Tibetans, as shown by the significant expression changes in genes such as INSR, PRKCZ, SREBF1, GFPT1, PPP1CB, and GYS1 (supplementary table 4, Supplementary Material online). In addition, insulin secretion seems down-regulated in Tibetans, as demonstrated by their lower expression of the placenta-specific insulin growth factor (IGF2). The down-regulation of IGF2 may be due either to its negative regulators, such as the placenta-specific glypican-3 (GPC3), which was specifically up-regulated in Tibetans (supplementary table 4, Supplementary Material online), or to its receptor IGF2R (Harris et al. 2011), which showed polymorphism between Tibetans and Han (data not shown). All these observations point to an increased glucose consumption during the fetal development of Tibetans. Such an increase may be one of the causes of the well-known fact that Tibetans have higher fetal weights than other populations inhabiting similar altitudes (Moore et al. 2001; Tripathy and Gupta 2005). Expression variations in glucose transporting and usage-related genes were not seen in ImHan, suggesting that the activation of angiogenesis may be the main adaptation mechanism for providing oxygen and nutrients for the embryo development in Han immigrants.

The More Efficient Energy Production of Tibetan Placentas

The aerobic capacity reflects the oxygen transport and utilizing ability, and it has therefore been used as a measure of the functional adaptation to high altitudes (Hochachka et al. 1991). Unlike the glycolysis elevation of ImHan under hypoxia (supplementary fig. 6, Supplementary Material online), Tibetans adopt a more efficient approach to energy production by up-regulating the TCA cycle and the oxidative phosphorylation pathways (fig. 4; supplementary fig. 6 and table 9, Supplementary Material online).

Our results suggest an additional Tibetan-specific strategy for saving energy via recycling building blocks under hypoxia. To balance the ROS accumulated under hypoxia, both ImHan and Tibetans up-regulate the antioxidant enzymes and the proteasome pathways (fig. 4, supplementary table 9, Supplementary Material online). In addition, Tibetans use autophagy (fig. 4, supplementary table 9, Supplementary Material online) to promote cell survival by recycling cellular components and damaged organelles (Bellot et al. 2009; Scherz-Shouval and Elazar 2011; Saito and Nakashima 2013; Filomeni et al. 2015; Gui et al. 2016; Jawhari et al. 2016). We hypothesize that the activation of autophagy in Tibetans is mediated by calcium signaling triggered by the accumulated ROS. It has been shown that ROS may increase the cytoplasmic calcium concentration (Wang and Zheng 2010), which activates calmodulin (CAML1) (Grotemeier et al. 2010). CAML1 further activates death-associated protein kinase (DAPK) and then triggers autophagy by releasing Beclin-1 from BCL2 (Zalckvar et al. 2009a, 2009b). In our data, the genes involved in this pathway, including CALM1, DAPK1, BECN1, and BCL2, are specifically highly expressed in Tibetans (fig. 4), suggesting such an activation from calcium signaling to autophagy is specific to native highlanders only. Such a unique activation, at least partially, derives from the genetics of Tibetans: our other results have shown significant allele frequency differences in multiple genes of the calcium-associated pathways between Tibetans and Hans (data not shown).

In addition to autophagy, ROS also trigger the UPR, a pathway to restore homeostasis (Ali et al. 2017). We found the activation of UPR in Tibetans. In particular, one of the UPR downstream pathways, the N-glycan biosynthesis, was specifically up-regulated in Tibetans with our most significant P-value (supplementary table 9, Supplementary Material online). A slight up-regulation of genes from UPR was observed in ImHan (supplementary fig. 7, Supplementary Material online), and no significant genetic variation was found between Tibetans and Hans. Therefore, UPR might be a general response to ROS stimulation, and because of the relatively mild hypoxia environment in Xining (2162–2500 m), ImHan does not have a strong UPR.

Factors Resulting in the Specific Expression Profile of Tibetans

To explore why the placenta-specific genes are up-regulated in Tibetans, we identified dozens of up-regulated TFs that may promote the high expression of these genes (supplementary fig. 4, Supplementary Material online). In particular, using our in-house data, we discovered that some genetic and epigenetic alterations in Tibetans appear to contribute to the up-regulation of these TFs. In addition to the mentioned polymorphism of EPAS1, for example, nine SNPs near the TSS of POU2F3 were identified in the H3K4me1 modified regions, which are enriched in enhancers in the placenta as shown by the Roadmap Epigenomics dataset. Furthermore, our methylation analysis (data not shown) also identified several TFs, including AFF1 and ESRRG, that had significant methylation variations between Tibetans and SLHan, which may also promote the up-regulation of placenta-specific genes in Tibetans.

Among these TFs, EPAS1 has been repeatedly validated as a key regulator of hypoxia adaptation in Tibetans (Beall et al. 2010; Bigham et al. 2010; Simonson et al. 2010). From our results, we also believe that EPAS1 may be a key factor during the fetal development of Tibetans. First, the expression of EPAS1 strongly correlated with the transcription level of its regulators/co-factors, as well as the placenta-specific genes (fig. 3). Furthermore, by analyzing the ChIP-seq datasets from cell lines of various tissues (Schodel et al. 2011; Smythies et al. 2019), we found that many placenta-specific genes are indeed EPAS1 targets, including PSG10P, PSG11, and GPC3 (fig. 3), whose products participate in angiogenesis and glucose regulation. Moreover, we identified a novel alternative splicing of EPAS1 in Tibetans, which most likely resulted from a Tibetan-specific SNP (fig. 3F and 3G). If translated, this new isoform is predicted to hold an enhanced transcription activity because of its similar structure to HIF-1α417, which was shown to be more active than the wild type (Lee et al. 2004). In particular, this short isoform implies a putatively hypoxia-adapted placental function of Tibetans, since the deletion of several oxygen-sensitive domains (fig. 3) may result in its functioning in an oxygen-independent manner. These results suggested that the positive selection of EPAS1 in Tibetans directly contributes to their successful reproduction at high altitudes.

In addition to the genetic variation promoting the up-regulation of placenta-specific genes in Tibetans, our Tibetan genomic data also showed several genes and pathways whose regulations were involved in the Tibetan-specific genetic background. For instance, in the calcium signaling and insulin secretion pathways, multiple SNPs have been found to have significant frequency differences between Hans and Tibetans. The coherence between the transcriptomic profiles and the genomic polymorphism of Tibetans strongly suggests that the Tibetan genetic variations significantly contribute to the hypoxia adaptation of native highlanders.

An interesting question arising from our study is how the transcriptome activities change in indigenous Tibetans that live at lower altitudes. This question requires further investigations due to the insufficient number of samples and, to the best of our knowledge, no such public data is available yet. However, by detecting the hypoxia adaptation-related physiological indicators, Petosi et al. found that indigenous Tibetans exhibited decreased hematocrit and Hb concentration, as well as increased resting ventilation compared to the SLHan when they moved to low land (Petousi et al. 2014). Hence, the Tibetan-specific functions activated in the placenta may be slightly inhibited when they are no longer exposed to a hypoxic environment.

In conclusion, our study provides numerous cues to how highlanders adapt or acclimatize to high altitude environments during fetal development. Some of the hypoxia response functions, such as the proteasome and the antioxidant systems, are conserved in both Han immigrants and native Tibetans. Both ImHan and Tibetans show the activation of certain functions, but with different routes. For instance, angiogenesis is improved by the up-regulation of general angiogenesis-related factors, such as VEGFA, in ImHan, while in Tibetans the up-regulation of placenta-specific factors, such as PGF, is required. Most importantly, some pathways and genes were found in the population-specific response to a hypoxic environment. In Tibetans, these include the activation of the placenta-specific hormone, the TCA cycle, and the autophagy, whereas in ImHan the inflammatory and the ECM pathways are up-regulated. Overall, the Tibetan-specific activation of the TCA cycle and the insulin signaling may allow more energy production and glucose usage in native highlanders. To a certain extent, the activation of the Tibetan-specific characteristics is determined by their unique genetic variations. At the same time, the specific transcriptome features of ImHan appear to be an immediate result of the hypoxic stress response.

Material and Methods

Tissue Collection and Ethnic Statement

Placenta samples were collected from healthy and full-term pregnant (37–42 weeks of gestation) puerperants after either vaginal or cesarean delivery. SLHan samples were from Peking University Third Hospital at Beijing (n = 25), ImHan samples were from QingHai Red Cross Hospital at Xining (n = 19), and Tibetan samples were from the Second People’s Hospital of Tibet Autonomous Region in Lhasa (n = 47). This study was designed and conducted in accordance with the principles of the Declaration of Helsinki. The project was approved by the IRB of Beijing Institute of Genomics. All donors signed their informed consenting forms.

The average altitude of Beijing is 43.5 m. The elevation of Xining city ranges from 2,162 m to 4,877 m with an average of 3,137 m (supplementary fig. 8, Supplementary Material online). Since QingHai Red Cross Hospital locates at downtown region (dark green in the map of supplementary fig. 8, Supplementary Material online) of 2,000–2,500 m, we describe the living altitude of ImHan as this range. The average altitude of Lhasa is 3,650 m and puerperants of the Second People’s Hospital of Tibet Autonomous Region come from Lhasa city or nearby regions. According to the residing altitude of each Tibetan donor (supplementary table 1, Supplementary Material online), we describe the living altitude of this group as >3,650 m in text.

The amnion and decidua samples were taken from the membranes outside chorion portion and separated by dissection. Chorion samples were cut as blocks of 1–2 cm. All specimens were washed with PBS twice to remove blood, and were then cut into small pieces either immediately freezing into liquid nitrogen or placing in RNAlater solution (AM7020, Thermo Fisher Scientific, Waltham, MA) and preserved at −80 °C until RNA isolation. All samples and their phenotype data were collected by medical personnel using the same protocol.

Immunohistochemistry

Immunohistochemistry was performed following the instructions provided with the biotin-streptavidin-peroxidase and diaminobenzidine kits (DS-0002, Zhongshan Golden Bridge Corp., Beijing, China), as previously described (Fu et al. 2010). Briefly, after incubation in citrate antigen retrieval solution (pH 6.0) and in normal goat serum, the placental sections (5 μm) were incubated with the primary antibody, β-hCG (ZM-0134, Zhongshan Golden Bridge Corp., Beijing, China), CSHL1 (ab174295, Abcam, Cambridge, UK), and PSG6 (MAB8598, R&D systems, Minneapolis, MN), respectively, overnight at 4 °C. On the following day, after washing the sections were incubated with biotinylated secondary antibody for 30 min. The signals were detected using diaminobenzidine solution.

Hypoxia Treatment and RNA Sequencing of HUVECs

HUVECs were isolated according to the published protocol (Crampton et al. 2007). Briefly, after rinse with pre-warmed PBS, a fresh cord was filled with pre-warmed trypsin in PBS and incubated at 37 °C for 15 min. The cells were suspended followed by centrifuging at 1,000 rpm for 5 minutes. After washing with PBS, cell pellets were spun at 1,000 rpm for 5 minutes and were re-suspended in 10 mL Lonza EGM-2 medium followed by incubation at 37 °C overnight with 5% of CO2. Supernatant was then removed and cells were replaced at 10 cm petri dishes with fresh media for culture at 37 °C. Fibroblasts were repetitively removed by re-suspending and rinse with medium manually under microscope. After several passages, cells were harvested and stained with CD31 antibody (11-0311-82, Abcam, Cambridge, UK). The percentages of endothelial cells were estimated by CD31 staining followed with flow cytometry.

For hypoxic treatment, the HUVECs of the fifth passage were harvested and seeded into six 10 cm petri dishes. These cells were then cultured in a hypoxia station (H35, Don Whitley Scientific, West Yorkshire, UK) supplied with 1% of oxygen. The HUVECs were collected at 0, 6, 12, 22, 36, and 48 hours after hypoxic treatment, followed by RNA extraction, library preparation, and sequencing (details see below).

Placenta RNA Sequencing and Data Processing

Total RNAs were extracted from 100 mg placental tissue of each sample using TRIzol reagent (15596018, Thermo Fisher Scientific, Waltham, MA), and were purified with RNeasy MinElute columns (74204, Qiagen, Hilden, Germany) according to the manufacturers’ protocols. Purity level and concentration of the total RNAs were measured by NanoDrop®ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). RNA integrity numbers (RINs) were estimated using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA).

After enrichment with ployA, the mRNAs of chorion, amniotic, and decidua were subject to strand-specific library preparation with KAPA Stranded mRNA-seq Kit (KK8420, KAPA BIOSYSTEMS, Wilmington, MA) according to manufacturers’ instruction. In brief, double-stranded cDNA fragments were synthesized from mRNA, ligated with strand-specific adapters, and size-selected for library construction. All high-quality strand-specific libraries were sequenced on HiSeq-2000 (Illumina, San Diego, CA).

FastQC and cutadapt were used to filter low quantity reads and adaptors. The rest reads were mapped to the human GRCh38 reference genome using TopHat2 (version 2.1.1) (Trapnell et al. 2009; Kim et al. 2013) with the fr-firststrand library-type. The transcript assembling, abundance estimation, and normalization were performed by Cufflinks (version 2.2.1) (Trapnell et al. 2010). To enlarge the maximum fragments allowed in a bundle before skipping, we set the max-bundle-frags with 10,000,000. Cuffdiff (Trapnell et al. 2013) was used to estimate gene expression levels and to detect significant DEGs with default settings. DAVID version 6.8 (Huang da et al. 2009) was used to obtain the functional overrepresentation of interested gene sets in GO terms as well as the KEGG pathways. Adjusted P values less than 0.05 were considered to be significantly overrepresented in the annotation category.

The expression values of each gene (FPKM > 0 in any of samples) were log transformed for calculating the expression similarity. The expression matrix containing all samples and all expressed genes was converted to the similarity distance based on Pcc. The heat map was generated by pheatmap package with ward.D2 clustering method. STEM (Ernst and Bar-Joseph 2006) was used to identify significant altitudinal expression profiles and associated genes.

Identification of EPAS1 Target Genes

ChIP-seq datasets from GSE120885 (Smythies et al. 2019) and GSE28352 (Schodel et al. 2011) were used to identify potential target genes of EPAS1. The nearest genes of EPAS1 binding peaks (UCSC hg19 reference) identified by ChIPseeker (Yu et al. 2015) were considered as the potential targets of EPAS1.

Identification of High-Frequency eQTLs in TIBETANS

The eQTLs of the 43 placenta-expression-related TFs were obtained from GTEx v7 (https://gtexportal.org/home/). Fst of these eQTLs between Tibetans and Hans was calculated by plink 1.09 (Purcell et al. 2007). Significant differences of allele frequencies were defined as Fst > 0.1.

GSEA

GSEA between sample groups was performed by GSEA v2-2.2.3 (Subramanian et al. 2005) to explore significantly changed KEGG pathways. Gene sets were downloaded from the human C2 KEGG dataset MsigDB v6.1. All positive- or negative-enriched gene sets were ranked by normalized enrichment scores. SLHan and ImHan were combined as the Han population, and ImHan and Tibetans were combined as the highlanders during the process.

Sanger Sequencing

The 256 bp fragments surrounding rs150877473 were amplified and sequenced by the forward primer 5′- CTGGAAGGTAGGGCAACATCA-3′ and the reverse primer 5′- GGTGCTGGATTGGTTCACAC -3′.

Whole-Genome Genotyping

Genomic DNAs were isolated from the amnions by DNA Mini Kit (51304, Qiagen, Hilden, Germany). Whole-genome genotypings were conducted with Illumina Global Screening Array (20030770, Illumina, San Diego, CA). The distribution of B allele frequencies from the genotyping results were used to estimate possible maternal DNA contaminations. And by this test, two samples from the multi-generation ImHan were removed for further analysis.

Scanning of Positive Selection Signals

The iHSs were calculated by selscan 2.0 using unphased mode (Szpiech ZA, unpublished data, https://www.biorxiv.org/content/10.1101/2021.10.22.465497v1, last accessed Oct. 24, 2021). The Neanderthal alleles of POU2F3 were obtained from Altai (Prufer et al. 2014) and Vindija 33.197 (Prufer et al. 2017).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.