基于高通量测序技术探究糖尿病足与健康人组织DNA中的5-羟甲基胞嘧啶变化

马睿瑶，赵景会，储金林，铁 璐，李念东，李琳琳

0 引 言

糖尿病是最常见的疾病之一，糖尿病足(Diabetic foot ulcer，DFU)是糖尿病一种严重的并发症，增加患者的经济负担，严重影响患者的生活，同时也是非创伤性下肢截肢的重要原因之一[1-2]。DNA中的5-羟甲基胞嘧啶(5hmC)主要是由5mC通过10-11位易位(TET)酶氧化形成。最近的研究表明表观遗传学在糖尿病创面中的作用[3]，DNA甲基化对糖尿病足的恢复具有一定作用[4]，同时在健康人和糖尿病或糖尿病血管病变患者中检测到不同水平的5hmC。因此，对糖尿病的另一种并发症——糖尿病足的研究引起极大兴趣。同时有研究表明，糖尿病患者的外周血中5hmC水平高于健康人[3]，糖尿病患者血浆中5hmC水平可用于判断糖尿病患者是否存在并发症[5]，但目前尚未有对糖尿病足溃烂组织中5hmC表达的报道。本研究使用5hmC-Seal技术[6]探究DFU和健康志愿者伤口组织DNA中的全基因组5hmC谱，筛选可能成为DFU的潜在治疗分子靶点。

1 材料与方法

1.1 一般资料选取2021年1月至 12月北京大学第三附属医院伤口治疗中心收治的6例糖尿病足患者组织为DFU组，男女各3例，选择同期年龄匹配的3例健康志愿者溃烂组织为对照组，男1例，女2例。本研究通过北京大学第三附属医院伦理委员会审批(批准号：M2020539)，患者均知情同意且签署知情同意书。

1.2方法

1.2.1 样本制备与DNA提取取20 mg新鲜冻存组织于EP管中, 采用Quick-DNATMMini Prep Plus Kit (ZYMO) 对组织中总DNA进行提取，-20 ℃冰箱稳定保存。按TRIzol©Reagent试剂盒说明操作测定组织中的总DNA。使用DNA进行5hmC文库建立[6]。并在IlluminaNextSeq500平台上进行高通量测序。

1.2.2测序数据的处理及分析通过Fragment Analyzer对每个文库样本的DNA条带大小及定位进行分析，使用Illumina公司的Nextseq 500/550 High Output Kit v2(75 cycles)在NextSeq 500 测序平台上进行双端 38-bp高通量测序，每个样本的测序通量为1.5Gb，条带大小为38 bp，最终得到给定文库中DNA片段的碱基序列。得到数据后，首先使用Trimmomatic软件去除样本中原始FASTQ中低质量的数据，以此提高后续分析效率。在使用Bowtie2软件，使用将原始文件比对到人类hg19基因组上，使用SAMtools进行筛选，以保留唯一匹配值bedtools首先将对端读取进行扩展并转换为BedGraph格式，并将其进行归一化处理，在使用UCSC基因组浏览器中的bedGraphToBigWig转换为bigwig格式，可在Integrated Genomics Viewer中进行可视化处理，其次可以合并峰值，只保留出现在3个样本以上且小于1000 bp的峰值区域。潜在的5hmC富集区域使用MACS识别(参数：MACS14-p1e-3-f BAM-g hs)。X和Y染色体被排除在外，生成文件作为后续分析文件[7]。

1.2.3筛选DhMGs用DEseq2软件包进行差异分析，以P<0.05和|log2FoldChange|≥0.5的条件筛选出DhMGs。

1.2.4DhMGs调控网络的构建与分析利用数据分析软件STRING(https://string-db.org/)对248个DhMGs进行相互作用分析(置信度≥0.9，互作最大值=0)，并建立蛋白互作网络 (PPI)。通过Cytoscape3.9.0(https://cytoscape.org/)软件构建基因调控网络，采用力导向算法让节点的分布更加合理，得到更佳的视觉效果。再使用CytoHubba插件，选用MCC方法，从构建的PPI网络中选取关键枢纽基因(Hub基因)。

1.2.5基因功能分析通过在线DAVID网站 (https://david.ncifcrf.gov/)对Hub基因进行GO功能注释和KEGG通路富集，分析这些基因的生物学过程、细胞组成和分子功能以及相关通路。

1.2.6荧光实时定量PCR (qRT-PCR)根据RNA快速提取试剂盒(EZB)说明书提取伤口组织总RNA，Nanophotometer P330测定RNA浓度。取1μg RNA，参照Vazyme公司逆转录试剂进行逆转录，将反转得到的cDNA 与反应液加至 96 孔板，根据EZB 2×SYBR Green qPCR Master Mix在Bio-Rad Q-PCR 仪中进行实时荧光定量PCR扩增。反应条件为：95 ℃预变性5分钟，95 ℃变性10秒，60 ℃退火/延伸30秒，共40个循环，以β-actin为内参，依据目的基因相对表达量RQ=2-△△Ct法计算基因mRNA的相对表达量，引物序列如表1所示。

表 1 qRT-PCR引物列表

1.2.7免疫组化石蜡切片脱蜡至水：依次将切片放入二甲苯-乙醇-蒸馏水中洗涤。抗原修复：组织切片置于盛满EDTA抗原修复缓冲液(pH 9.0)的修复盒中于微波炉内进行抗原修复，中火5 min加热，停火5min保温再转中低火10min，此过程中应防止缓冲液过度蒸发，切勿干片。自然冷却后将玻片置于PBS(pH 7.4)中在脱色摇床上晃动洗涤3次，每次5 min。阻断内源性过氧化物酶：切片放入3%双氧水溶液，室温避光孵育25 min，将玻片置于PBS(pH 7.4)中在脱色摇床上晃动洗涤3次，每次5 min。血清封闭：在组化圈内滴加3%BSA 均匀覆盖组织，室温封闭30 min。敷一抗：轻轻甩掉封闭液，在切片上滴加PBS按一定比例配好的一抗，切片平放于湿盒内4°C孵育过夜。敷二抗：玻片置于PBS(pH 7.4)中在脱色摇床上晃动洗涤3次，每次5 min。切片稍甩干后在圈内滴加与一抗相应种属的二抗(HRP标记)覆盖组织，室温孵育50min。DAB显色：玻片置于PBS(pH 7.4)中在脱色摇床上晃动洗涤3次，每次5 min。切片稍甩干后在圈内滴加新鲜配制的DAB显色液，显微镜下控制显色时间，阳性为棕黄色，自来水冲洗切片终止显色。复染细胞核：苏木精复染3 min左右，自来水洗，苏木精分化液分化数秒，自来水冲洗，苏木精返蓝液返蓝，流水冲洗。脱水封片：将切片依次放入酒精-二甲苯中脱水透明，将切片从二甲苯拿出来稍晾干，中性树胶封片。显微镜镜检，图像采集分析。分别加入一抗MAPK11(1∶200)和CDC42(1∶200)，二抗为HRP山羊抗兔(1∶200)。

2 结 果

2.1 5hmC-Seal测序结果测序数据分析显示，DFU组的整体羟甲基化水平高于对照组。主成分分析图结果显示，利用变化的5hmC可以将DFU组和对照组区分开。见图1。

a:5hmC在不同基因组特征中的分布图； b:主成分分析图

2.2DhMGs的筛选以P<0.05和|log2FoldChange|>0.5的条件下筛选出DhMGs，共有166个上调基因和262个下调基因，其中包括NDUFB6、GALR1、GDAP1等已知的与糖尿病并发症相关的基因。DhMGs主要分布在内含子、基因间和启动子区域。见图2。

2.3DhMGs调控网络的构建与Hub基因筛选涉及35个节点，见图3。节点的颜色与Degree值密切相关，随着节点增大，由蓝变红所代表的Degree值也随之升高。边的粗细程度则代表Combine score，随着边由细变粗，Combine score值也随之变大。边的颜色与coexpression值相关，随着颜色由黄变为紫色，coexpression值也随之变大，见图3。通过CytoHubba分析可得到10个Hub DhMGs，分别是MAPK14、MAPK11、CDC42、MAP3K1、NCOA1、RAP1B、AGO4、RUNX1、HIF1A、PHC3，见表2。

a：DhMRs分布火山图以红色(DhMRs在糖足中表达上调)或绿色(DhMRs在糖足中表达下调)表示，黑点表示没有差异表达；b: DhMRs基因组特征分布图(以橘色(DhMRs在糖足中表达上调)或绿色(DhMRs在糖足中表达下调)表示，橘色表示上调与下调DhMRs总和)

图 3 35个DhMRs网络图

表 2 糖尿病足溃烂组对比健康对照组Hub基因表

2.4GO富集分析对10个Hub基因进行富集分析，在生物过程中(蓝色)涉及到VEGFR信号通路、Ras蛋白信号转导、肌细胞分化的正向调节以及与免疫相关的白细胞介素12分泌的正向调节等；在细胞成分中(红色)，主要与细胞质、胞液等相关；在分子功能中(绿色)，主要有细分化调控、血管内皮生长因子受体信号通路、髓样细胞分化调控等功能，见图4。

图 4 GO富集分析

2.5KEGG通路分析KEGG分析显示Hub基因主要富集于白细胞跨内皮细胞迁移、有丝分裂原活化蛋白激酶信号通路、血小板激活、T细胞受体信号通路和血管内皮生长因子信号通路等，见图5。其中血小板激活包括RAP1B、MAPK11、MAPK14；血管内皮生长因子信号通路包括CDC42、MAPK11、MAPK14；有丝分裂原活化蛋白激酶信号通路包括RAP1B、CDC42、MAPK11、MAP3K1、MAPK14。

图 5 KEGG富集分析

2.6基因表达验证DFU组MAPK11与 CDC42 mRNA(1.75±0.34、1.21±0.12)明显高于对照组(1.02±0.14、0.95±0.10)，免疫组化结果表明MAPK11与CDC42表达水平在DFU中更高(P<0.05)，见图6。

图 6 组织基因的蛋白表达水平(免疫组化)

3 讨 论

DFU是糖尿病一种严重的并发症，近一半的患者寿命不足5年，死亡率高于某些肿瘤[8]，近年来，随着对糖尿病创面的关注不断增加，也有文献对糖尿病创面愈合的表观遗传机制进行了研究，但对于DFU的5hmC表观遗传学研究尚不多见[9]。因此，利用5hmC对DFU具有潜在的研究价值。

首先本研究发现，5hmC在DFU与健康志愿者组织DNA中的5hmC富集有明显差异，另外，5hmC标志物可以将两组区分开。表观遗传修饰已被确定为发生1型糖尿病和2型糖尿病的原因之一，关于糖尿病伤口愈合的表观遗传机制的研究已经出现在文献中，为我们研究提供理论基础[9-11]。

本研究筛选出428个DhMGs，其中166个基因在DFU患者中发生高度羟甲基化，并且大部分DhMGs分布在启动子、内含子和外显子等基因功能区。有研究表明，基因体位置发生高度羟甲基化会促进基因转录[6]，所以，对这些DhMGs的深入分析有助于我们探究DFU的发病机制。进一步地，通过STRING对428个DhMGs进行PPI分析，使用Cytoscape建立基因调控网络图，我们筛选出MAPK14、MAPK11、CDC42、MAP3K1、NCOA1、RAP1B、AGO4、RUNX1等10个Hub基因。我们对MAPK11与CDC42 mRNA水平和蛋白水平进行验证，发现与5hmC表达水平一致。

MAPK14、MAPK11分别称为称为p38-α酶和p38-β酶，均属于p38 MAPK家族，在细胞受到促炎细胞因子或生理应激等外界刺激下引起的细胞级联反应中发挥重要作用，激活转录因子。既往研究中，使用MAPK抑制剂，糖尿病小鼠伤口愈合能力增强[12]，同时还发现在胃溃疡患者中，MAPK通路与也被激活[13]。另外，文献证明，使用全基因组关联研究(GWAS)发现 MAPK14基因与糖尿病足溃烂相关[14]。与MAPK14同家族的MAPK11，目前虽无MAPK11与糖尿病足溃烂愈合的明确文献，根据本研究结果提示MAPK11或也可能影响糖尿病足溃烂的愈合。CDC42作为Rho GTPases家族的成员，参与多种细胞功能的调节，包括细胞周期进程、存活、转录、肌动蛋白细胞骨架组织和膜运输。最近研究表明，CDC42参与胰岛素的合成过程，与胰岛素抵抗和糖尿病肾病的发病机制相关，被认为是糖尿病进展过程中的关键成员[15]。这与本研究结果一致，CDC42在启动子区域富集5hmC表达升高，而在启动子区域富集5hmC可以促进基因的转录[16]，通过改善CDC42的表达有望成为糖尿病足溃烂治疗的方法。

最后，本研究通过 DAVID数据库对10个Hub基因进行 GO 功能和 KEGG通路分析发现，富集到VEGFR信号通路、MAPK信号通路、血小板激活等这些与伤口愈合相关的通路。文献表明糖尿病患者伤口愈合延迟的一个主要原因是缺氧诱导因子-1α/血管内皮生长因子(HIF-1α/VEGF)轴的功能受损，从而导致对缺氧的反应受损的新血管形成，我们富集出VEGFR信号通路也证明这点[17-18]。成纤维细胞是皮肤的基本细胞类型，高度参与伤口再生过程[19]，有研究对GEO数据关于糖尿病伤口(GSE49566 和GSE78891)测序结果分析，其富集到最高的通路为MAPK通路，同时在 MAPKAPK3、HSPA2 和 TGFBR1 中确定了糖尿病伤口愈合的调节作用[20]。

有趣的是我们还富集到中性粒细胞通路、白细胞介素12分泌通路、白细胞跨内皮细胞迁移和T细胞受体信号通路等与免疫相关的通路，提示这些基因与炎症及免疫系统相关，最近文献表明免疫在感染性疾病中占有重要地位[21-22]，我们的研究也佐证了这点。研究表明，多形核白细胞(PMNS)，也称为中性粒细胞，是一组相对较大的白细胞，最为第一道免疫防线的要作用[23]，参与稳态、伤口愈合和组织修复[24]，其与淋巴细胞的比值可以反映炎性的程度[25]，中性粒细胞的招募失调是糖尿病中性粒细胞起到重足溃烂伤口难愈合的原因之一[23, 26]。白细胞迁移可能与向血管损伤部位输送高反应性促炎和促凝血介质有关[27]。白细胞介素12作为免疫微环境重要的细胞因子在感染性疾病中被研究[28]，它由M1巨噬细胞分泌，期其分泌的异常与巨噬细胞的状态相关，而巨噬细胞的招募与极化对糖尿病伤口溃烂极为重要[29-30], 是其免疫异常的主要原因之一[31]。

综上所述，本研究通过对DFU患者及健康对照组的组织DNA中5hmC进行测序，发现大量的DhMGs，通过生物学网络中节点的连接和关系来分析网络特性，得到DFU相关的基因，富集到与免疫与伤口愈合的相关通路，最后得到MAPK11和CDC42可能在DFU发生和发展中起到重要的作用，有望成为DFU的治疗靶点。