世居平原、高原习服与世居高原三类人群DNA甲基化差异及其功能解析

DNA methylation differences and functional analysis across native lowland, acclimatized, and native highland populations

摘要

背景 人体响应高原环境的分子特征与适应机制，对急慢性高原病的预防、诊疗具有重要意义。已有研究表明，世居高原人群具有明确的遗传基础和独特的表观遗传特征，但在全基因组DNA甲基化层面，目前尚未检索到针对世居平原及习服高原、世居高原人群之间进行系统比较的研究报道。 目的 探讨平高原环境下的群体DNA甲基化差异及其功能解析。 方法 自建高原研究队列，包括世居平原人群(native lowlanders，NLs)、迁居至高原的习服人群(acclimatized newcomers，ANs)和世居高原人群(native highlanders，NHs)，利用全基因组亚-硫酸测序技术(whole-genome bisulfite sequencing，WGBS)获取其甲基化数据，人群之间进行差异甲基化区域分析(differentially methylated regions，DMRs)及表观基因组关联分析(epigenome-wide association study，EWAS)，识别群体之间的差异甲基化信号并进行GO/KEGG富集解析其生物功能。 结果 研究纳入117例NLs(男性，年龄20.15±1.43岁)、162例ANs(男性，年龄22.07±2.65岁)和124例NHs(男性101人、女性23人，年龄22.10±7.34岁)。ANs与NLs组相比，甲基化水平在基因组范围上整体降低(P＜0.01)，差异分析得到显著的440个DMRs及611个差异甲基化位点(differentially methylated CpGs，DMCs)，驱动免疫功能(如T细胞活化、NK细胞毒性通路等)与能量代谢模式(如LDHC基因)发生显著改变，但同时伴随心肌病相关通路的异常激活。NHs与Ans相比，具有1 313个DMRs及2 022个DMCs，相关差异主要富集于钙信号通路及EPAS1/EGLN1等基因区域。NHs整体甲基化水平显著低于NLs，两组人群比较具有247个DMRs及688个DMCs，二者差异主要富集于HIF-1(q＜0.01)与自噬信号通路。 结论 DNA甲基化在人群适应高原环境过程中发挥重要作用。NHs不仅在整体甲基化水平上显著低于NLs，还显示出更多特异性的差异甲基化区域，提示其可能通过精准调控关键基因的甲基化状态，实现对低氧环境与细胞稳态的长期稳定适应。

Abstract

Background Understanding molecular mechanisms of high-altitude adaptation is crucial for combating mountain sickness. While native highlanders exhibit distinct genetic/epigenetic traits, systematic whole-genome methylation comparisons across acclimatization stages are lacking. Objective To investigate the differences in DNA methylation among populations in plain and plateau environments and analyze their functional implications. Methods A self-established high-altitude research cohort was composed of native lowlanders (NLs), acclimatized newcomers (ANs), and native highlanders (NHs). Whole-genome bisulfite sequencing (WGBS) was used to obtain their methylation data. Differential methylated regions (DMRs) analysis and epigenome-wide association studies (EWAS) were conducted among the populations to identify differential methylation signals between the groups. GO and KEGG enrichment analyses were then performed to elucidate the biological functions of these methylation differences. Results This study included 117 NLs (all men, aged 20.15 ± 1.43 years), 162 ANs (all men, aged 22.07 ± 2.65 years), and 124 NHs (101 men and 23 women, aged 22.10 ± 7.34 years). After entering the plateau, the methylation levels of NLs decreased genome-wide (P＜0.01). Compared with NLs, ANs exhibited 440 DMRs and 611 differentially methylated CpGs (DMCs), driving significant changes in immune functions (such as T cell activation and NK cell cytotoxicity pathways) and energy metabolism patterns (e.g., the LDHC gene), but also accompanied by abnormal activation of pathways related to cardiomyopathy. Compared with ANs, NHs had 1 313 DMRs and 2 022 DMCs, with differences mainly enriched in the calcium signaling pathway and gene regions such as EPAS1/EGLN1. The overall methylation level of NHs was significantly lower than that of NLs, and NHs had 247 DMRs and 688 DMCs, with differences mainly enriched in the HIF-1 pathway (q＜0.01) and the autophagy signaling pathway. Conclusion DNA methylation plays an important role in human adaptation to high-altitude environments. In comparison, NHs not only exhibit a significantly lower overall methylation level than NLs, but also display a greater number of specific differentially methylated regions, suggesting that they may achieve long-term stable adaptation to hypoxic environments and cellular homeostasis through precise regulation of the methylation status of key genes.

Graphical abstract

关键词

高原适应 / 习服 / 低氧 / DNA甲基化 / 信号通路 / 青年

Key words

high-altitude adaptation / acclimatization / hypoxia / DNA methylation / signal transduction / young adult

以低压、低氧、低温、高紫外线等为典型特征的高原极端环境条件，易引发人体产生多种生理和病理反应，发生急慢性高原病，包括可直接危及生命的肺水肿或脑病^[1]。解析人体适应或习服高原的分子机制，实现对急慢性高原病的有效预防和快速治疗，是当前亟待解决的关键问题。已有研究表明，世居高原人群的适应性具有明确的遗传基础，其在血红蛋白调控、心肺功能^[2]、生殖能力^[3]等多个方面展示出独特的适应优势，是开展高原研究的理想对象。高原适应不仅具有遗传基础，DNA甲基化等表观遗传调控也发挥重要作用，如EPAS1、EGLN1等关键基因附近的遗传变异对于高原适应至关重要^[4]，其附近区域也呈现独特的DNA甲基化模式^[5]。针对高原适应的遗传基础解析，多个大队列人群测序研究已报道众多高原适应相关变异和功能基因^[5-7]，多聚焦于世居高原人群的遗传基础。针对高原适应性的表观遗传调控机制，目前研究并不充分，缺乏利用高分辨率全基因组DNA甲基化测序技术，系统描绘从平原至高原适应过程中DNA甲基化动态变化图谱的研究，现有研究初步提示DNA甲基化参与高原适应，但存在如样本量有限^[4]和检测技术覆盖不全^[8]等明显的局限性。本研究首次基于世居平原人群(native lowlanders，NLs)、习服人群(acclimatized newcomers，ANs)和世居高原人群(native highlanders，NHs)三个群体队列的全基因组DNA甲基化测序(whole-genome bisulfite sequencing，WGBS)数据，通过三组整体DNA甲基化水平比较与两两差异分析，系统解析了不同高原适应阶段人群的DNA甲基化差异及其功能：ANs与NLs的比较，旨在阐明高原暴露初期的表观遗传应激响应；NHs与ANs的比较，目的是区分短期习服与长期适应的表观差异；而NHs与NLs的比较，则聚焦于高原长期适应过程中形成的稳定表观遗传特征。本研究系统阐释了DNA甲基化在人群高原习服与长期适应中的独特调控作用，为深入理解高原适应的分子机制提供了重要理论依据，并为高原相关疾病的生物标志物发掘和潜在干预靶点的探索提供了重要参考。

1 对象与方法

1.1 研究对象

本研究为一项横断面调查研究，已获得解放军总医院医学伦理委员会批准(伦审第S2018-298-02号)。研究对象选取与样本收集工作于2021年7月至2022年1月完成，依据高原暴露背景分为3组。

1.1.1 世居平原组

来自四川省甘孜地区(海拔约1 300 m)，其纳入标准为世居平原(海拔＜500 m)、无超过2周的高原(海拔＞2 500 m)旅居史，且经体检确认健康、无重大疾病史；排除标准包括近期有高原旅行史、患有慢性心肺疾病或代谢性疾病等基础疾病，以及有长期服药史者。

1.1.2 高原移居组

来自西藏阿里地区(海拔约4 300 m)，为首次从平原进驻高原并持续居住20～24个月的健康个体，入高原前无重大疾病史及慢性高原病史；排除标准为累计居住时间＜20个月或居住期间曾离开高原连续超过1个月、患有慢性高原病或入高原前后存在重大疾病或服药史。

1.1.3 高原世居组

同样来自西藏阿里地区(海拔约4 300 m)，其纳入标准为至少三代世居高原、无连续超过6个月的低海拔(海拔＜2 500 m)生活史，经体检确认健康且无慢性高原病；排除标准包括三代内有与低海拔居民通婚史、患有慢性高原病或其他重大慢性疾病，以及有长期(连续＞6个月)低海拔地区生活史者。

1.2 标本采集与数据质控

1.2.1 标本采集

静脉血样品均于清晨受试者空腹状态下采集，经EDTA抗凝后，在4℃条件下2 h内完成离心分装，并保存于-80℃冰箱。遵循华大智造提供的标准流程，在深圳完成标准的文库构建：使用TIANamp Blood DNA Kit提取DNA，使用MGIEasy全基因组甲基化文库制备试剂盒进行末端修复、加接头、重亚硫酸盐转化及PCR扩增等步骤。使用DNBSEQ-T7平台进行150 bp双端全基因组重亚硫酸盐测序，测序的目标覆盖深度为30X​。

1.2.2 数据质控

下机原始数据经FastQC进行质控，使用Bismark等软件将质控后的数据比对至参考基因组(hg38)。为控制性别混杂效应，本研究排除X/Y染色体所有CpG位点并在回归模型中校正性别，甲基化水平采用β值(β=甲基化C读数/ [甲基化C+非甲基化T读数])表示，以样本有效CpG位点的平均β值作为该样本甲基化水平；通过严格质控筛选保留覆盖深度(C+T)＞5的可靠位点，同时剔除跨样本缺失率＞20%的CpG位点。

1.3 统计学分析

1.3.1 整体DNA甲基化水平及差异比较

为了对不同高原适应状态人群的DNA甲基化模式进行整体比较，对3组人群分别从全基因组、典型高原适应性基因及差异甲基化区域(differentially methylated regions，DMRs)3个层面进行统计分析。个体样本的全基因组DNA甲基化水平以其所有高质量CpG位点的β值的平均值表示，组间整体比较使用Shapiro-Wilk检验对3组内个体DNA甲基化水平分布进行正态性检验，若服从正态分布则采用独立样本t检验，任一组不符合正态分布，则采用Wilcoxon秩和检验。基于基因在文献中被多次独立报道的频率生成的词云图，从中选取了报道频率最高的6个基因作为典型高原适应性基因^[9]，单个基因DNA甲基化水平定义为其注释区域内(hg38)所有CpG位点的平均β值。EPAS1/EGLN1基因DMRs整体比较取3组两两比较显著DMRs的并集，DMRs本身的甲基化水平以其所包含的所有CpG位点的平均β值表示。所有统计分析采用R编程语言环境(v4.4.0)完成，显著性水平设定为P＜0.05，统计图表及可视结果使用ggplot2(v3.4.2)制作。

1.3.2 差异甲基化区域分析DMRs的检测

采用metilene软件(v0.2-4)^[10]。该工具基于二元分割算法(binary segmentation algorithm)结合二维统计检验(two-dimensional Kolmogorov-Smirnov test)，可在全基因组范围内高效识别样本组间显著差异的甲基化区域。分析采用默认参数框架，关键参数设置如下：最小CpG数量设定为5，以增强短区域DMRs的检测灵敏度；最大无信息间隔(maxgap)设定为300 bp(默认值)，避免跨越长距离无CpG位点的区域被错误合并；显著性阈值：DMRs的统计显著性通过二维Kolmogorov-Smirnov检验评估，并采用BH法对P值矫正后得到q，取q≤0.05的区域；甲基化差异阈值：设定平均甲基化水平组间绝对差值(|Δβ|)≥0.1作为DMRs的生物学显著性阈值。基因组区域的功能注释通过R语言环境中的ChIPseeker包(版本1.40.0)完成^[11]，基于参考基因组GRCh38(GENCODE v44基础注释)。输入区域经标准化处理后，采用annotatePeak函数进行多层级注释，包括基因本体距离(转录起始位点上下游2 kb范围)、基因组特征分类(启动子区、外显子/内含子、远端调控区)。

1.3.3 表观基因组关联分析

采用多元线性回归模型进行表观基因组关联分析(epigenome-wide association study，EWAS)鉴别组间显著的差异甲基化CpG位点(differential methylation CpGs，DMC)。因变量(y)为二元分组(如NHs=1，ANs=0)，自变量为CpG位点的甲基化β值，选取年龄(连续型)、性别(二元型)、甲基化数据的前5个主成分(PC1-PC5，控制细胞异质性和技术变异)以及实验批次作为协变量。P值经Benjamini-Hochberg(BH)法校正(q＜0.05)。QQ图对EWAS结果的统计学分布进行验证。使用R包LOLA^[12]基于大规模数据集预测的基因组功能元件注释，对显著 CpG 位点开展功能元件富集分析^[13]。通过从所有测试的 CpG位点进行随机抽样，将得到的显著CpG位点与随机抽样的分布进行比较，得到OR值和P值，并使用Benjamini-Hochberg法进行错误发现率(false discovery rate，FDR)校正，设定显著性阈值q＜0.05。16个基因组元件包括组装间隙和比对伪影(GapArtf)、静息状态(Quies)、异染色质状态(HET)、多梳抑制状态(ReprPC)、乙酰化状态(Acet)、弱活性增强子(EnhWk)、活跃增强子(EnhA)、转录增强子(TxEnh)、弱转录状态(TxWk)、强转录状态(Tx)、外显子关联转录状态(TxEx)、锌指基因区(znf)、DNase I超敏位点(DNase)、二价启动子(BivProm)、启动子侧翼区(PromF)和转录起始位点(transcription start site，TSS)。基因区域富集分析方法同功能元件富集，注释文件来源于UCSC基因组数据库(https：//genome.ucsc.edu/)基于hg38人类参考基因组组装版本。

1.3.4 GO/KEGG通路富集

使用clusterProfiler包(v4.12.0)^[14]进行GO和KEGG通路富集分析。KEGG分析基于2023年1月版本数据库，并限定于人类特异性通路(hsa标识符)，显著性阈值设定为BH校正后q＜0.05。

2 结果

2.1 研究对象基本情况

共纳入403例研究对象：NLs 117例，均为男性，年龄(20.2±1.4)岁；ANs 162例，均为男性，年龄(22.1±2.7)岁；NHs 124例，男101例，女23例，年龄(22.1±7.3)岁。经过WGBS测序，质控后获得14 259 122个可靠CpG位点。

2.2 NHs-ANs-NLs间的整体甲基化水平及差异甲基化比较

2.2.1 整体甲基化水平分析

NLs组的整体甲基化水平显著高于ANs组和NHs组(Wilcoxon检验，P＜0.01)，而ANs组与NHs组间的整体差异无统计学意义(P＞0.05)(图1A)。高原环境显著降低DNA整体甲基化水平，提示高原环境对表观遗传存在全局性影响；短期习服和长期适应相比，整体甲基化水平趋于稳定，表明高原适应可能通过更精细的表观调控实现，而非全基因组甲基化的显著变化。

2.2.2 键高原适应性基因区域的甲基化水平比较

对6个高原适应性关键基因(EGLN1、EPAS1、HBB、HBG2、RHOQ、TRIM67)的甲基化水平进行分析^[15]，结果显示3组人群呈现显著的差异甲基化模式(图1B)。这些基因涉及HIF-1通路响应(EPAS1、EGLN1)、血红蛋白生成(HBB、HGB2)等重要低氧适应的生理过程，表明DNA甲基化可能通过调控这些生理过程影响高原适应。

2.2.3 整体差异甲基化模式比较

整体来看，差异甲基化分析，ANs vs NLs得到显著的DMRs共440个，NHs vs ANs 得到1 313个DMRs，而NHs vs NLs得到247个DMRs，NLs、ANs和NHs3组之间DMRs存在交集，这些DMRs在高原的早期和长期适应中一直存在动态变化，通过对DMRs的注释后3组交集的基因有17个，仅仅占DMRs注释基因总数的1.1％(图1C)，提示短期习服和长期习服分子机制相差较大，在3组比较所得的DMR中，相互之间存在重叠的区域共15个，其中8个位于基因区域(图1E)，7个位于基因间区。研究比较了EPAS1和EGLN1基因区域的DMRs在3组之间的水平，EPAS1基因区域内存在5个DMRs，分布在基因间区、启动子区域和内含子区域。EGLN1基因在第一内含子区域存在一个DMRs，其在世居高原群体中存在明显的高甲基化，而在习服高原和平原人群之间无显著差异(图1D)。

2.3 ANs与NLs的差异甲基化及功能解析

2.3.1 差异甲基化区域及功能解析

相比NLs，ANs的差异区域共440个DMRs，其中98.2%(432个)为低甲基化，仅1.8%(8个)为高甲基化(图2A)。这些DMRs主要位于内含子、基因间区(包括潜在远端调控元件)和启动子区域(图2B)，表明高原习服过程中DNA甲基化主要发挥表观调控作用，远端基因间区的富集提示染色质空间构象(如增强子－启动子互作)可能参与调控过程。对440个DMRs注释基因的功能富集分析显示：GO-BP通路(q＜0.05)显著富集于细胞黏附调控、T细胞活化、白细胞迁移及淋巴细胞分化等生物学过程。KEGG通路主要富集于NK细胞介导的细胞毒性、PD-L1表达与肿瘤免疫逃逸通路、T细胞受体信号及Th17细胞分化等免疫相关通路(图2C)。

2.3.2 EWAS及功能解析

EWAS分析系统比较了ANs和NLs两个群体的DNA甲基化差异位点特征(ANs标识为1，NLs为0)，经多重检验校正后，共鉴定出611个显著DMC(q＜0.05)。这些位点注释至基因LOC107985200、LOC124900240、ANKFN1、TRPM2、JRKL、FAM13C及CDKN3的CpG位点关联信号最强(-log₁₀P＞20)。值得注意的是，LOC107985200、JRKL与TRPM2基因座呈现重复关联信号(图3A)，提示这些基因可能通过表观遗传调控参与高原适应性生理过程。基于所有CpG位点EWAS分析结果绘制的QQ图显示(图3B)，观测P值与期望P值在大于1E-4基本一致，且在分布尾部显著上扬，表明分析发现了潜在的真实甲基化关联信号。DMC主要富集于1号、20号和21号染色体(图3C)，而且在功能区域转录起始位点(TSS，OR=4.8，q=0.003)、双向启动子(BivProm，OR=3.6，q=0.018)呈现显著富集，其OR值均大于3，提示这些调控元件的高甲基化状态与低氧适应存在强正相关，而静止染色质状态(Quies)和多梳抑制区域(ReprPC)呈现最低OR值(OR＜1)，进一步印证活跃调控元件的特异性富集(图3D)。该分布模式提示高原习服相关表观变异主要富集于基因表达调控的核心区域。特别是TSS区域的高显著性与其在染色质开放状态维持中的生物学功能相符，可能通过调控RNA聚合酶Ⅱ的招募效率影响靶基因转录。这些发现为解析低氧适应的表观遗传调控网络提供了重要线索。进一步，GO富集显示其在cAMP依赖性蛋白激酶活性的调控和肌营养不良蛋白相关糖蛋白复合物等通路显著富集(q＜0.05)；KEGG富集显示其在扩张型心肌病、致心律失常性右室心肌病和可卡因成瘾通路显著富集(q＜0.05)(图3E)。

2.4 NHs与ANs的差异甲基化及功能解析

2.4.1 差异甲基化区域及功能解析

与ANs比较，NHs共有1 313个DMRs(上调DMRs：1 230个；下调DMRs：83个)(图4A)，主要富集于启动子区、基因间区及内含子区域(图4B)，提示其可能通过调控基因表达参与高原适应机制。对DMRs注释基因的功能富集分析显示：GO显著富集于负端导向的微管运动活性(q＜0.05)，可能参与细胞骨架重塑可能参与缺氧条件下的细胞定向迁移与分裂；KEGG通路显著富集于钙信号通路(q＜0.05)，反映钙离子稳态在高原适应中的核心作用(图4C)。

2.4.2 EWAS及功能解析

NHs和ANs两个人群(NHs标识为1，ANs为0)共鉴定出2 022个DMC(图5A)，QQ图显示，观测P值与期望P值在le-3以后逐渐上扬(图5B)，显著位点主要富集于2号和16号染色体的功能调控区域，在基因组中呈现出非随机分布(图5C)，显著富集区域包括转录起始位点(TSS)、二价启动子(BivProm)等活跃调控区域(图5D)；而显著缺失区域包括静止染色质状态(Quies)、多梳抑制区域(ReprPC)等非活跃区域。其中EPAS1基因区域呈现最强关联信号(图5A)，提示其可能通过调控HIF-1信号通路在低氧适应中发挥核心作用。值得注意的是，位于EPAS1下游的TMEM247基因亦检出5个显著DMC，相邻的甲基化模式差异表明其可能通过染色质空间互作(如增强子-启动子环)协同调控高原适应性生理过程。KEGG通路显著富集于肌肉细胞中的细胞骨架通路和HIF-1信号通路(P＜0.01)；GO通路显著富集与miRNA的正向调控和RNA聚合酶核心酶结合等通路(q＜0.01)(图5E)。

2.5 NHs与NLs的差异甲基化及功能解析

2.5.1 差异甲基化区域及功能解析

与NLs比较，NHs分析共鉴定出247个DMRs，其中上调DMRs共85个(34.4%)，下调DMRs共162个(65.6%)(图6A)。这些DMRs在基因组功能区域中主要富集于内含子、基因间区及启动子区域(图6B)，表明表观遗传调控可能通过靶向基因的顺式调控元件影响高原适应性生理过程。对DMRs注释基因的功能富集分析显示，KEGG通路显著富集于HIF-1信号通路、自噬通路和扩张性心肌病通路(q＜0.01)；GO显著富集于髓系细胞分化、细胞对缺氧反应和聚集体自噬等(图6C)。

2.5.2 EWAS及功能解析

NLs和NHs的对比(NHs标识为1，NLs为0)系统揭示了NHs在长期高原适应过程中的特异DNA甲基化变异特征(图7A)。结果显示，相较于NLs，NHs在全基因组范围内共鉴定出688个显著DMC(q＜0.05)，QQ图显示P值为1E-3后尾部逐渐上扬(图7B)，1号染色体和16号染色体数量最多(图7C)，最强关联信号集中于EPAS1基因附近，其6个CpG位点形成显著关联峰簇，表明DNA甲基化通过精细调控核心适应基因EPAS1影响低氧适应；同时ANKRD26P1、FAM230C及LOC107985200等基因发现多个CpG位点甲基化变化，提示表观调控在该基因中发挥的重要作用；此外，EGLN1和ANGPT1等缺氧通路基因的强关联信号进一步证实了自然选择压力下表观基因组的定向进化特征(图7A)，可能与世居高原人群的HIF通路活性优化及血管生成调控相关。这些发现共同揭示了NHs在世代适应过程中形成的特异性表观遗传特征。功能元件富集同样显示其在更加能发挥生物学功能的基因区域如转录起始位点(TSS)，启动子(BivProm)区域更加富集，而在发挥生物学功能不活跃的区域如静止染色质区域(Quies)、多梳抑制区域(ReprPC)不富集(图7D)。显著CpG位点KEGG富集结果显著富集于吗啡上瘾、MAPK信号通路、心肌细胞中的肾上腺素能信号通路、自噬和HIF-1信号通路(P＜0.01)，GO富集在蛋白激酶A结合和cAMP依赖蛋白激酶活性(P＜1e-3)等通路(图7E)。

3 讨论

本研究首次系统比较了世居平原、习服高原和世居高原3个人群的全基因组DNA甲基化图谱，以揭示高原环境暴露下平高原相关人群的表观遗传特征及其潜在生物学机制。与徐书华团队^[4]2023年仅32例样本的小规模研究相比，本研究样本量显著扩大(403例)，统计功效更强；与杨剑团队^[8]2025年使用450K芯片(仅覆盖约45万CpG位点)的研究相比，WGBS技术实现了全基因组范围的精细检测，能够捕捉更全面的甲基化变异。在人群设计方面，本研究注重样本同质性：世居平原人群和高原习服人群均主要由青年男性构成，年龄分布相近(NLs：20.2±1.4岁；ANs：22.1±2.7岁)，且具有相似的饮食与生活习惯，从而最小化了年龄、性别及环境因素对甲基化差异的潜在混淆；世居高原人群样本在年龄(22.1±7.3岁)和生活方式方面也保持了较高的同质性。然而，由于高原地区样本采集存在挑战性，为最大化样本量并更全面地代表世居人群的真实特征，本研究纳入了少量女性样本(23例)。这一决策基于以下考量：首先，世居高原人群经过多代自然选择，其高原适应性表观遗传修饰(如EPAS1/EGLN1区域的稳定甲基化模式)具有性别间的稳定性，性别差异可能对整体适应模式影响有限；其次，在统计分析中已采用严格措施控制性别混杂效应，包括将性别作为协变量纳入多元线性回归模型，并排除了X和Y染色体上的所有CpG位点，以确保结果的可靠性。本研究提供的高质量WGBS数据及差异甲基化图谱，不仅揭示了DNA甲基化在高原适应中的重要作用，更为后续寻找可用于评估习服状态的表观遗传标志物、阐释高原病易感性的表观遗传机制等提供了数据基础与候选靶点集合。

研究揭示了DNA甲基化可能介导了关键基因在高原适应中发挥作用。在短期习服阶段(ANs vs NLs)，甲基化改变集中于快速应激保护与代谢免疫调节，如GAS6低甲基化(q=2.33E-23)可能通过抑制凋亡通路提升心肌细胞对低氧缺血再灌注损伤的保护能力^[15]；在离子通道与能量代谢调控方面，TRPM2基因第3内含子在ANs中呈现大幅去甲基化(Δβ=-0.56，q=2.23E-178)，该钙通道激活可能导致胞内Ca²⁺超载及线粒体活性氧簇(ROS)累积，成为高原病防治的新靶点^[16-17]；LDHC基因启动子区甲基化显著下调(q=2.36E-07)，高原习服可能通过其去甲基化促进糖酵解效率^[18]；但MUC4去甲基化(q=9.05E-77)可能通过HIF-α通路激活加重组织纤维化风险，提示高原习服中潜在的病理机制^[19]。经分析比较，长期适应阶段(NHs vs ANs)的特征是对缺氧通路核心基因的精准优化：EPAS1基因(HIF-2α编码基因)在世居高原人群中存在6个显著DMRs，其中启动子区低甲基化(q=3.71E-22)可能通过提高转录因子结合效率提升其转录活性，进而优化下游低氧应答通路，与血红蛋白浓度的精细调控直接关联^[20-23]；EGLN1基因(缺氧感应因子)的第一内含子高甲基化(q=4.90E-21)可能抑制增强子活性，通过钝化HIF通路减少EPO过度表达，优化氧气利用效率，支持世居高原人群低血红蛋白表型^[20-23]；NOS3基因(eNOS编码基因)启动子高甲基化(q=6.81E-24)调控内皮型一氧化氮合酶表达，影响血管舒张功能，帮助维持较高但不过度的NO水平，平衡缺血缺氧保护与潜在负效应^[23]；世居高原人群血液NO水平的调控类似于其血红蛋白的调控，相较于世居平原人群，虽然表现出更高的NO水平，但其NO水平相比于同海拔短期习服人群更低^[24]。而世居高原与习服高原人群的直接比较进一步证实，这种由EPAS1/EGLN1核心轴主导，并得到EGLN3、HIF3A等基因协同的、经过自然选择固定的表观遗传调控，共同优化了低氧感应、血管张力及能量代谢等关键生理过程，标志着从短期的应激反应向长期稳定的适应机制的进化转变。

研究进一步证明，高原习服和适应是两个涉及不同通路的生理过程。进入高原并习服两年后(ANs)，最显著的DNA甲基化变化体现在与免疫系统功能密切相关的通路上。DMRs和EWAS分析均显示，差异甲基化区域和位点显著富集于T细胞活化与分化(如T细胞受体信号通路、Th17细胞分化)、NK细胞活性、白细胞迁移及炎症反应调控等生物学过程和通路。提示表观遗传变化很可能通过抑制过度炎症反应(如TRPM2的去甲基化可能保护细胞免受低氧诱导的钙超载死亡)、调节免疫细胞功能、增强组织屏障防御和促进能量代谢(如LDHC的去甲基化可能增强糖酵解)来帮助机体应对高原低氧压力，达到高原习服状态。然而，这种以免疫调节为主导的短期适应策略，其负面效应也不容忽视：一方面，它提升了机体在高原的生存能力；另一方面，持续的免疫激活状态(如富集到的PD-L1表达通路提示免疫检查点可能参与调节)和心肌相关通路(如扩张型心肌病、致心律失常性右室心肌病)的富集(图3E)，可能潜在地增加了组织损伤、自身免疫风险以及心血管系统负担，这为解释习服高原相较于世居高原人群具有更高的急慢性高原病(如心血管相关疾病)易感性提供了分子线索。相较于习服过继免疫应激通路的激活，长期适应在关键基因以及富集的通路都不同于前者。两个公认与高原适应性相关的关键基因，EPAS1和EGLN1，其邻近区域在世居平原、习服高原两类人群同世居高原的比较中均被鉴定出显著DMRs和强EWAS信号。例如，在世居高原人群中观察到的EPAS1启动子区低甲基化，可能直接促进该基因的基础表达水平，优化其对低氧的感知和应答；而EGLN1内含子区在世居高原人群中的特异性高甲基化状态，考虑是长期自然选择形成的稳定调控标记，可能通过影响该基因的剪接或表达，精细调控HIF通路的活性。这些发现将已知的高原适应遗传基础与其独特的表观遗传特征紧密联系起来，揭示了遗传－表观遗传协同进化在塑造高原适应表型中的核心作用。在生理功能通路富集方面，相较于习服高原人群以免疫为主的富集，世居高原人群的差异甲基化区域/位点富集于各种生理功能通路，如钙信号通路(ANs vs NHs)、HIF-1信号通路(NLs vs NHs)、自噬(NLs vs NHs)、血管功能调节(如NOS3基因的甲基化差异)、能量代谢、血小板功能(GP5/GP6)以及心肌细胞信号传导(如MAPK信号、肾上腺素能信号)。这表明世居高原群体的高原适应不仅仅是应对低氧，同时涉及对呼吸、循环、代谢、细胞稳态维持等核心生理系统的稳态调节，实现整体生理机能的优化。

从整体DNA甲基化水平来看，相较于世居平原群体，习服和世居高原的整体DNA甲基化水平显著下降，这一现象在习服高原人群尤为显著，表明低氧压力可能影响甲基化转移酶DNMT或TET的酶活性，或导致SAM(甲基供体)代谢改变，从而造成广泛的甲基化改变。其差异区域主要富集于基因组的调控功能区，如启动子、基因间区(可能包含远端调控元件)和内含子(可能包含增强子等调控元件)。这与既往部分小规模研究或特定基因观察到的低氧诱导去甲基化趋势基本一致^[7]。全局性的去甲基化可能是机体对高原环境的一种急性应激反应，类似于其他环境压力(如氧化应激^[25]、营养缺乏^[26])下观察到的现象。这种反应可能通过松弛染色质结构、提高基因组结构的可塑性，为后续适应性基因表达的调控变化提供条件。类似地，高原习服人群的全局性低甲基化状态，可能反映其基因组稳定性受到一定挑战，进而增加 DNA 损伤和异常表达的风险；这或许是该类人群急慢性高原病（如高原脑水肿、高原肺水肿）发病率高于世居人群的一个潜在表观遗传机制^[27]。3类人群两两比较的DMRs交集分析清晰地表明，其共有的DMRs及其注释基因比例极低(仅占1.1%)。这再次表明，高原短期习服和长期自然选择适应所依赖的分子机制存在显著不同。短期习服更依赖于快速、可塑性强的免疫和应激反应相关通路的表观变化，而长期适应则通过自然选择固定了在关键适应性基因(如EPAS1、EGLN1)和核心生理通路(如钙信号、HIF信号、自噬)上的稳定表观遗传修饰(如EGLN1内含子区的世居高原群体特异性高甲基化)。

本研究局限性在于目前研究仅提供了表观遗传组的差异信号，缺少多组学整合以及细胞或动物模型的实验验证。未来，可进一步将本文报道的差异甲基化区域/位点与转录组、染色质状态(ATAC-seq，ChIP-seq)、蛋白质组以及已有的基因组变异数据进行多组学整合分析，并针对关键DMRs(如EPAS1启动子、EGLN1内含子、TRPM2内含子、SETDB1相关区域等)，利用细胞或动物模型进行功能验证，旨在明确其甲基化状态对基因表达及下游生理功能的作用。

数据共享声明　本论文相关数据可依据合理理由从作者处获取，Email：[e-mail]。

参考文献

[1] Zhang XK， Yang YX， Shi QH. DNA methylation in adaptation to high-altitude environments and pathogenesis of related diseases［J］. Hum Genomics，2025，19（1）：100.

[2] Moore LG. Measuring high-altitude adaptation［J］. J Appl Physiol （1985），2017，123（5）：1371-1385.

[3] He YX， Guo YB， Zheng WS，et al. Polygenic adaptation leads to a higher reproductive fitness of native Tibetans at high altitude［J］. Curr Biol，2023，33（19）：4037-4051.

[4] Hu H， Petousi N， Glusman G，et al. Evolutionary history of Tibetans inferred from whole-genome sequencing［J］. PLoS Genet，2017，13（4）：e1006675.

[5] Lin ZS， Lu Y， Yu GL，et al. Genome-wide DNA methylation landscape of four Chinese populations and epigenetic variation linked to Tibetan high-altitude adaptation［J］. Sci China Life Sci，2023，66（10）：2354-2369.

[6] Childebayeva A， Goodrich JM， Leon-Velarde F，et al. Genome-wide epigenetic signatures of adaptive developmental plasticity in the Andes［J］. Genome Biol Evol，2021，13（2）：evaa239.

[7] Shi JL， Jia ZL， Sun JX，et al. Structural variants involved in high-altitude adaptation detected using single-molecule long-read sequencing［J］. Nat Commun，2023，14（1）：8282.

[8] Cheng FF， Shen RJ， Zheng ZL，et al. Distinct methylomic signatures of high-altitude acclimatization and adaptation in the Tibetan Plateau［J］. Cell Discov，2025，11（1）：45.

[9] Yu JJ， Non AL， Heinrich EC，et al. Time domains of hypoxia responses and-omics insights［J］. Front Physiol，2022，13：885295.

[10] Jühling F， Kretzmer H， Bernhart SH，et al. Metilene： fast and sensitive calling of differentially methylated regions from bisulfite sequencing data［J］. Genome Res，2016，26（2）：256-262.

[11] Yu GC， Wang LG， He QY. ChIPseeker： an R/Bioconductor package for ChIP peak annotation，comparison and visualization［J］. Bioinformatics，2015，31（14）：2382-2383.

[12] Sheffield NC， Bock C. LOLA： enrichment analysis for genomic region sets and regulatory elements in R and Bioconductor［J］. Bioinformatics，2016，32（4）：587-589.

[13] Vu H， Ernst J. Universal annotation of the human genome through integration of over a thousand epigenomic datasets［J］. Genome Biol，2022，23（1）：9.

[14] Yu GC， Wang LG， Han YY，et al. clusterProfiler： an R package for comparing biological themes among gene clusters［J］. OMICS，2012，16（5）：284-287.

[15] Zhao AZ， Lei WR， Tian JY，et al. Mangiferin attenuates myocardial ischemia reperfusion injury by regulating the GAS6/axl signaling pathway［J］. Phytother Res，2025，39（3）：1388-1402.

[16] Cheng K， Jia QM， Batbatan C，et al. TRPM2-L participates in the interleukin-6 pathway to enhance tumor growth in prostate cancer by hypoxia-inducible factor-1α［J］. J Interferon Cytokine Res，2023，43（11）：495-511.

[17] Ping J， Liu XY， Lu YM，et al. A highland-adaptation variant near MCUR1 reduces its transcription and attenuates erythrogenesis in Tibetans［J］. Cell Genom，2025，5（3）：100782.

[18] An ZF， Wei DB， Wei L，et al. The expression of Ldh-c in the skeletal muscle of plateau pika （Ochotona curzoniae） enhances adaptation to a hypoxic environment［J］. Biol Open，2017，6（9）：1336-1341.

[19] Li ZT， Hu XH， Ni WJ，et al. Hypoxia-inducible factor-1α promotes tubulointerstitial fibrosis via the upregulation of Muc4［J］. Ren Fail，2025，47（1）：2501757.

[20] Tashi T， Scott Reading N， Wuren T，et al. Gain-of-function EGLN1 prolyl hydroxylase （PHD2 D4E： C127S） in combination with EPAS1 （HIF-2α） polymorphism lowers hemoglobin concentration in Tibetan highlanders［J］. J Mol Med （Berl），2017，95（6）：665-670.

[21] Gray OA， Yoo J， Sobreira DR，et al. A pleiotropic hypoxia-sensitive EPAS1 enhancer is disrupted by adaptive alleles in Tibetans［J］. Sci Adv，2022，8（47）：eade1942.

[22] Yasukochi Y， Nishimura T， Ugarte J，et al. Effect of EGLN1 genetic polymorphisms on hemoglobin concentration in Andean highlanders［J/OL］.

[23] Ozden ES， Ozcan MS， Ilhan I，et al. Radiofrequency electromagnetic field ınhibits HIF-1 alpha and activates ENOS signaling to prevent intestinal damage in a model of mesenteric artery ischemia in rats［J］. Int J Med Sci，2025，22（6）：1465-1476.

[24] He YX， Qi XB，Ouzhuluobu，et al. Blunted nitric oxide regulation in Tibetans under high-altitude hypoxia［J］. National Science Review，2018，5（4）：516-529.

[25] Damiano OM， Stevens AJ， Kenwright DN，et al. Chronic inflammation to cancer： the impact of oxidative stress on DNA methylation［J］. Front Biosci （Landmark Ed），2025，30（3）：26142.

[26] Kadayifci FZ， Zheng SS， Pan YX. Molecular mechanisms underlying the link between diet and DNA methylation［J］. Int J Mol Sci，2018，19（12）：4055.

[27] Gatterer H， Villafuerte FC， Ulrich S，et al. Altitude illnesses［J］. Nat Rev Dis Primers，2024，10：43.