基于基因芯片技术的高原红细胞增多症基因组DNA甲基化差异分析*

罗勇军,陈 郁,蒲 懿,张 莉,董红梅△

(1.陆军军医大学陆军卫勤训练基地军事医学地理学教研室,重庆 400038;2.陆军军医大学基础医学院学员五大队十五队,重庆 400038;3.中国人民解放军联勤保障部队第941医院心肾内科,西宁 810001)

当进入海拔超过2 500 m的高原地区后,人体会在低氧环境作用下出现明显的生理反应,因而医学界将海拔超过2 500 m定义为医学高原[1]。为代偿低氧环境、改善低氧血症,红细胞会明显增生,这种增生在一定程度内能够提高运氧能力和血氧水平、减轻缺氧损伤,是人体习服高原环境的重要机制。然而,当红细胞增生超过一定范围时,就会演变为高原红细胞增多症(high altitude polycythemia,HAPC)[2],其带来的危害已经抵消了原有的有益效应。过度增生的红细胞降低血流速度,而动脉血氧饱和度降低进一步加重低氧血症,可造成全身多器官多系统的缺氧损伤,并伴有相应的临床症状及体征,严重影响高原居民的身体健康,尤其在高原移居人群中更为普遍,造成的危害也更为严重。HAPC是高原环境和人体遗传因素相互作用的结果,其发病机制较为复杂,也尚未完全阐明。因此,探讨HAPC的发病机制是降低HAPC危害的重要研究方向。

DNA甲基化是调控基因表达的重要表观遗传学机制之一,能够显著影响疾病的发生与进展[3]。HAPC的遗传易感因素较为复杂,多个基因参与其中,单个基因的甲基化并不能解释HAPC的发病机制,在全基因组层面的甲基化研究也还处于初步阶段。本研究通过全基因组DNA甲基化芯片技术建立高原移居汉族HAPC特异性甲基化基因表达谱,对比HAPC患者与健康者基因组水平DNA甲基化差异,并借助生物信息学的方法和技术筛选相关差异基因,分析其潜在功能和分子作用网络,为探索HAPC的发病机制,以及可能的靶向治疗和预防方案提供理论依据。

1 资料与方法

1.1 一般资料

按照HAPC的诊断标准[4],在世居中国青海省西宁市的汉族男性人群中筛选了4例HAPC患者(病例组),血红蛋白(hemoglobin,Hb)水平216～253 g/L;另选取同期体检的5例高原汉族健康者作为对照组,Hb水平170～181 g/L,年龄、性别匹配。本研究经陆军军医大学伦理委员会批准(批准编号:2020第001-02号),所有受试者均签署知情同意书。

1.2 方法

1.2.1血样采集

抽取受试者晨起空腹外周静脉血2 mL,经乙二胺四乙酸(EDTA)抗凝后,用于提取基因组DNA。

1.2.2基因组DNA提取与质检

采用美国OMEGA生物技术有限公司的全血基因组提取试剂盒(货号:D3392-02)提取血样中的基因组DNA[5]。提取完成后先用分光光度计定量,并将样品调到标准水平50 ng/μL,取20 μL,然后用0.8%的琼脂糖凝胶电泳质检。样品电泳结果提示基因组DNA主带清晰,通常相对分子质量不小于10×103,没有明显降解,总量5 μg以上,能够进行下游的甲基化芯片检测。

1.2.3基因组DNA甲基化水平检测

采用美国Illumina公司的甲基化芯片Infinium MethylationEPIC BeadChip(850 k芯片)检测受试者基因组DNA甲基化水平[6]。该芯片覆盖了853 307个CpG位点,并全面覆盖基因启动子区、基因编码区、CpG岛、增强子区及99%的RefSeq基因。后续的亚硫酸盐转化,DNA扩增,DNA的片段化、沉淀和重悬,DNA与芯片的杂交,芯片清洗、单碱基延伸、染色及芯片扫描和数据提取在北京博奥晶典生物技术有限公司完成。

1.3 统计学处理

将处理好的芯片放入扫描仪,利用激光激发芯片上的单碱基延伸产物的荧光基团,扫描仪获取由荧光基团发出的荧光,并生成原始数据,记录扫描结果存放的位置。由此所得的数据直接导入GenomeStudio软件进行分析,根据Illumina官方甲基化分析算法获得每个位点的原始信号强度值。然后经过对不同荧光、探针类型引起的偏差校正及归一化、位点过滤,得到高质量CpG位点的归一化后甲基化水平,即β值,用于质控。根据β值计算样品间Pearson相关系数(Pearson correlation coefficient,PCC),以评估受试样品芯片信号质量的重复率(β值为0表示该位点无甲基化,β值为1表示该位点完全甲基化)。质控后的数据,采用R软件包IMA3.1.2进行差异甲基化分析,其采用的方法为Limma中的经验贝叶斯统计[7-8]。同时针对多重假设检验问题,计算错误发现率(false discovery rate,FDR)校正P值,以降低假阳性率。差异甲基化位点的选取标准为:校正P≤0.05(如果所有位点均校正P>0.05,或≤0.05的位点很少,则以校正之前的P≤0.05为标准)。针对差异甲基化位点,使用R脚本对样本进行聚类分析,并使用通路富集工具KOBAS注释这些差异位点映射到的基因,通过与京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)GENES数据库交联链接,挖掘出统计上明显的基因本体论(gene ontology,GO)功能富集分析和KEGG信号通路分析[9],探讨差异甲基化基因的潜在功能及在HAPC发病机制中可能的作用。通过上述生物信息学分析,探讨差异甲基化基因及其作用网络在高原移居汉族HAPC发病机制中可能的作用。所有检验为双侧,以P<0.05为差异有统计学意义。

2 结 果

2.1 甲基化芯片评价

将扫描得到的原始数据通过GenomeStudio软件获得每个位点的原始信号值和测量P值(DetectionP值),二者用于质控以评估数据质量和后续分析。测量P>0.05表明样品的质量较低,需要去除。原始信号值提示所有样品的位点检测率均在99.9%以上,保证了检测范围。质控包括样本独立对照(sample-independent controls,SICs)、样本非独立对照(Sample-dependent controls,SDCs)、位点质控及个体质控。SICs用于评估操作步骤,SDSs用于评估样品质量。位点质控要求位点在95%以上的个体中测量P<0.05,同时去除位于X、Y染色体上的位点。个体质控要求个体在95%以上位点中的测量P<0.05。在经历严格的质控后,获取了病例组和对照组的β值。PCC分析提示病例-对照β值的PCC为0.997 9,表明芯片信号质量很好,所得数据可以用于后续分析,见图1。

图1 病例组与对照组甲基化水平Pearson相关性分析

2.2 差异甲基化分析

按照差异甲基化位点的筛选标准,共筛选得到96 360个差异甲基化位点。其中,根据位点所属区域差异,分为CpG岛(10 699,11.1%)、岛滩区(23 799,24.7%)、岛架区(7 370,7.6%)和其他区域(54 492,56.6%)。根据甲基化水平变化,病例组的高甲基化位点共10 564个(11.0%),低甲基化位点共85 796个(89.0%)。根据对照组与病例组对比结果显示:对照组中大部分CpG位点呈现高甲基化状态,而病例组大部分甲基化水平降低,与对照组相比差异甲基化基因24 054个。根据差异甲基化位点所映射的基因变化,与对照组相比,病例组的高甲基化基因共5 981个,低甲基化基因共18 073个。

2.3 差异甲基化基因谱

在前述筛选出的差异甲基化位点所处的基因中,按照筛选标准获取了10个最明显的高甲基化基因,分别是含黄素单氧化酶3(flavin-containing monooxygenase 3,FMO3)、G蛋白核仁2(G protein nucleolar 2,GNL2)、钙调磷酸酶类似EF-Hand蛋白2(calcineurin like EF-hand protein 2,CHP2)、酰基辅酶A硫酯酶2(acyl-CoA thioesterase 2,ACOT2)、染色体1开放阅读框25(chromosome 1 open reading frame 25,C1orf25)、OCA2黑素体跨膜蛋白(OCA2 melanosomal transmembrane protein,OCA2)、肌动蛋白相关蛋白2/3复合亚基4(actin related protein 2/3 complex subunit 4,ARPC4)、锚蛋白重复域55(ankyrin repeat domain 55,ANKRD55)、GRB10交互作用GYF蛋白2(GRB10 interacting GYF protein 2,GIGYF2)、溶酶体蛋白跨膜4β(lysosomal protein transmembrane 4 beta,LAPTM4B),见表1。筛选的10个最明显低甲基化基因分别是黏附蛋白偶联受体B1(adhesion G protein-coupled receptor B1,ADGRB1)、Rab亲和蛋白3A类似物(不含C2结构域)[rabphilin 3A Like (without C2 domains),RPH3AL]、酰辅酶A硫酯酶1(acyl-CoA thioesterase 2,ACOT1)、ATP结合盒亚家族C成员13(ATP binding cassette subfamily C member 13,ABCC13)、钙粘蛋白22(cadherin 22,CDH22)、染色体1开放阅读框109(chromosome 1 open reading frame 109,C1orf109)、防御素β128(defensin beta 128,DEFB128)、锚蛋白重复域23(ankyrin repeat domain 23,ANKRD23)、酰基辅酶A合成酶家族成员3(Acyl-CoA synthetase family member 3,ACSF3)、钙/钙调蛋白依赖性蛋白激酶1D(calcium/calmodulin dependent protein kinase 1D,CAMK1D),见表2。

表1 筛选出的10个最明显高甲基化基因(病例组 vs. 对照组)

表2 筛选出的10个最明显低甲基化基因(病例组 vs.对照组)

2.4 差异甲基化位点所在基因生物功能富集分析

2.4.1GO功能富集分析

按照筛选标准,病例组与对照组的差异甲基化位点所在基因相关的GO条目共19 585个,根据P值筛选了前10个最明显的GO条目。其中,主要富集于生物过程的物质定位、生物过程和细胞过程的正向调控及细胞发育等,以及细胞组分的细胞质及分子功能的结合,见表3。

表3 差异甲基化位点所在基因的GO分析

2.4.2KEGG信号通路分析

病例组与对照组的差异甲基化位点所在基因相关的KEGG信号通路共305个,符合筛选条件的信号通路共42条(P<0.05),并根据P值筛选了前10个最明显的信号通路,见表4。结果提示,差异甲基化位点所在基因涉及的信号通路主要包括代谢通路、癌症通路及丝裂原活化蛋白激酶(MAPK)信号通路。

表4 差异甲基化位点所在基因的KEGG通路分析

3 讨 论

HAPC是高原地区的常见疾病,海拔4 000 m以上地区患病率高达24.0%[10],而在海拔超过5 000 m的喀喇昆仑山脉地区则更高,移居人群如汉族人群可达80%以上[11],严重威胁着高原地区人群的健康。因此,探索HAPC的发病机制是维护高原人群生命安全的重要研究方向。HAPC的发生是高原环境和人体相互作用的结果,涉及多个系统和多个环节,其分子机制不仅与核基因序列、线粒体DNA(mtDNA)序列改变相关,还涵盖表观遗传学的变化,包括DNA甲基化、组蛋白修饰(如乙酰化、苏木化和磷酸化)及非编码RNA调控等多种机制。PENG等[12]发现了缺氧诱导因子(hypoxia inducible factor,HIF)途径的内皮含PAS结构域蛋白1(EPAS1)和脯氨酰羟化酶蛋白(EGLN1)基因可能是高海拔遗传适应的候选基因。FAN等[13]通过全外显子测序发现,在西藏世居人群中磷脂酰肌醇-3激酶催化亚基δ(PIK3CD)和Ⅳ型胶原蛋白α3链(COL4A3)基因上存在差异单核苷酸多态性(SNP)位点。ZHOU等[14]通过全基因组关联研究,发现红细胞生成调节基因小泛素相关修饰蛋白特异性蛋白酶1(SENP1)和癌基因酸性富含亮氨酸的核磷蛋白32家族成员D(ANP32D)在患有慢性高原疾病的患者中表达上调。陈辉等[15]发现,白细胞介素-12受体B1(IL12RB1)基因的SNPsrs393548、rs436857 和rs845380与西藏世居人群HAPC的发生有关。

目前,关于HAPC遗传易感机制的研究多为DNA水平,表观遗传机制研究较少。本研究借助高通量基因芯片技术,首次从全基因组层面分析了中国汉族HAPC患者的基因组DNA甲基化水平,发现HAPC患者基因组DNA大部分呈低甲基化状态,并筛选了以FMO3为代表的高甲基化基因和以ADGRB1为代表的低甲基化基因,借助生物信息学分析富集了相关的信号通路,提示基因组DNA的异常甲基化也是HAPC发生的重要机制。

除了基因本身的变化外,表观遗传学的改变也是人体适应高原的重要机制[16]。DNA甲基化作为表观遗传的主要机制之一,其核心是由S-腺苷甲硫氨酸作为甲基供体,在DNA甲基化转移酶的催化下,提供的甲基被转移到相应碱基的一种DNA共价修饰方式,并不涉及DNA序列的改变[17]。已有研究提示,作为高原适应的主要候选基因EGLN1,其启动子区域的高甲基化与安第斯山印第安人HAPC有关[16];而另一个高原适应候选基因EPAS1,其甲基化水平与安第斯山印第安人的早期发育和终生高海拔暴露及高海拔适应性表型有关[18-19],而血管紧张素Ⅱ 1型受体相关蛋白的配体(Apelin)基因的高甲基化也在高原肺水肿的发生中发挥着重要作用[20]。在本研究中,发现HAPC患者的基因组DNA多呈低甲基化状态,提示其基因表达较为活跃。在长期的低氧环境刺激下,作为低氧反应重要的调控因子HIFs,在与特异的顺式作用元件相结合后,启动了一系列基因的转录,并以此形成信号转导的级联放大效应,建构复杂的分子作用网络,涵盖人体对缺氧的感知和反应,从造血干细胞分化成红系祖细胞再到成熟红细胞的形成,骨髓造血微环境的调节,促红细胞生成素的分泌与反馈,Hb合成所需的铁、维生素B12等原料的吸收合成与利用。上述过程中不可避免的涉及多种表观遗传学变化,而本研究发现了HAPC患者和高原健康对照在多个基因的甲基化有着不同程度的变化,提示异常的甲基化改变在上述过程中发挥着重要的调节作用。

在大鼠的HAPC模型中,EGLN1基因的甲基化水平与正常对照并无明显差异[21],而von Hippel-Lindau(VHL)基因启动子的高甲基化会降低VHL基因表达,进而促进HIF-2α和促红细胞生成素(EPO)的表达,诱导HAPC的发生[22]。本研究提示,HAPC患者的FMO3呈现高甲基化状态,而抑制FMO3的表达能够降低全血中胶原依赖的血小板与基质的黏附性,对血栓形成有抑制作用[23]。在HAPC患者中,由于红细胞过度增生易并发血栓,尤其是诱导微循环障碍,加重全身缺氧损伤。因此,推测HAPC患者的高甲基化FMO3状态可能是功能代偿,从血小板角度抑制血栓形成。在低甲基化基因方面,本研究提示HAPC患者的ADGRB1基因差异最为明显。ADGRB1基因主要在脑部表达,在吞噬、炎症、突触形成、抑制血管生成和成肌细胞融合等方面发挥着重要作用[24]。研究发现,ADGRB1基因的高表达能够抑制脑血管生成,与脑部肿瘤发生和肿瘤周围脑水肿呈负相关[25]。此外,研究提示HAPC患者出现了不同程度的脑水肿[26],而本研究中发现的低甲基化ADGRB1基因可能也是通过提高其表达水平以拮抗脑水肿的发生。

在本研究中,借助多种生物信息学分析方式探讨了差异甲基化基因的功能及其潜在的作用网络,探寻其在HAPC发生中的作用。综合分析两组差异甲基化位点,提示病例组和对照组基因组DNA甲基化水平差异明显,而且病例组绝大部分CpG差异甲基化位点数目低于对照组。基因组DNA的高甲基化位点意味着基因表达水平下降,而HAPC患者的高甲基化位点数目较健康者明显下降,表明在长期的低氧环境刺激下,人体内部分基因被激活,其形成的信号放大机制和分子相互作用网络在HAPC的发生、发展中发挥了重要作用。GO分析提示,差异甲基化基因参与了生物过程的正向调控、细胞质、解剖结构发展和解剖结构形态发生等过程。上述富集的功能均可能涉及HAPC的发病机制,如红细胞生成与细胞发育和正向调控有关,细胞质的变化是红细胞在成熟过程中的必然反应,而结合既包含了转录因子和DNA的结合,也包括了蛋白和蛋白之间的结合与相互作用,在红细胞生成和成熟过程中均发挥了重要作用。KEGG通路分析提示,包括代谢通路、癌症通路和MAPK信号通路等在HAPC对照组中呈现高富集。红细胞增生过程中涵盖了代谢通路的调整,也必然涉及能量代谢的改变,而MAPK信号通路在能量代谢过程中发挥着重要作用。此外,红细胞的过度增生,在生物学过程方面类似于癌细胞的过度增殖,故在信号通路富集过程中也涉及癌症通路的变化。上述结果进一步提示,HAPC的发生是机体在应对低氧环境刺激时调控网络失衡的结果,低氧感知、红细胞发育分化与成熟、造血微环境如炎症状态、红细胞生成原料摄取等相关信号通路依然是探索HAPC发病机制的重要方向,如FMO3可影响血栓形成,ADGRB1在炎症免疫、血管生成和肿瘤发生中均发挥着不同程度的效应,这也进一步提示甲基化在HAPC的发生中发挥了重要的调控作用。

综上所述,本研究借助高通量基因芯片技术,从基因组甲基化角度探索了HAPC可能的发生机制,并发现了甲基化异常与HAPC发生有关。虽然纳入的样本量较小,生物学重复也未完成,在结果解释方面有一定的局限性,但是目前Infinium MethylationEPIC BeadChip(850 k芯片)的应用越来越广,研究方法和技术稳定,研究结果仍有一定的借鉴意义,提示甲基化改变也在HAPC发生中发挥着重要作用,异常的甲基化位点可能是HAPC发生的诊断标志物。本研究在基因组层面检测了甲基化水平,并未对某个具体部位或某个基因的具体位置再检测,针对差异甲基化基因,后期还需要借助甲基化特异性PCR进一步扩大样本验证结果。此外,虽然本研究只是甲基化改变在HAPC发生机制研究的开端,但是其中也提供了大量证据,提示HAPC的发病机制研究也要注意炎症免疫、血管生成和血栓形成等信号通路,其相关的表观遗传学调控机制是HAPC早期治疗和预警的客观依据。