成年双生子空腹血糖、糖化血红蛋白与全基因组DNA甲基化的相关性研究

王兆年，高文静△，王碧琦，曹卫华，吕 筠，余灿清，逄增昌，丛黎明，汪 华，吴先萍，刘 彧，李立明

(1. 北京大学公共卫生学院流行病与卫生统计学系，北京 100191； 2. 青岛市疾病预防控制中心，山东青岛 266033； 3. 浙江省疾病预防控制中心，杭州 310051； 4. 江苏省疾病预防控制中心，南京 210009； 5. 四川省疾病预防控制中心，成都 610041； 6. 黑龙江省农垦总局疾病预防控制中心，哈尔滨 150090)

据世界糖尿病联盟报告[1]，2019年全球范围内成人糖尿病标化患病率为9.3%，而有研究报道[2]，2013年我国糖尿病患病率为10.9%，高于世界平均水平。糖尿病是一组以高血糖为特征的代谢性疾病，高血糖可对人体肾脏、心脑血管等多个器官系统造成损害[3]。空腹血糖及糖化血红蛋白(glycated haemoglobin, HbA1c)均可反映机体糖代谢状况，并作为糖尿病的诊断指标。美国糖尿病协会建议[4]，抽取静脉血测量的空腹血糖≥7.0 mmol/L或HbA1c≥6.5%的患者应被诊断为糖尿病。空腹血糖是最直接反映血糖代谢状况的指标，也是最早应用于糖尿病诊断的指标。HbA1c与空腹血糖相比，更能反映最近几周内中长期血糖变化，2011年被世界卫生组织推荐为糖尿病的诊断指标之一[5]。另外，血糖升高即便没达到糖尿病的诊断标准，也可对心血管系统造成不利影响，提高心血管系统疾病的发病风险[6]。

DNA甲基化是指在DNA的胞嘧啶-磷酸-鸟嘌呤(cytidine-phosphate-guanosine site，CpG)二核苷酸的胞嘧啶5′端在酶催化作用下共价键结合一个甲基基团。DNA甲基化可通过影响转录酶的结合强度调控基因的表达、细胞的蛋白质合成等代谢过程，因此，DNA甲基化研究可以从组织和机体代谢改变层面探讨疾病的致病机制和过程，为探索复杂疾病的病因及疾病间联系提供线索。空腹血糖是一个短期的血糖代谢指标，可反映当下时点的血糖水平；HbA1c代谢周期则较长，可反映患者机体的中长期血糖代谢的平均水平，因此，空腹血糖与CpG位点甲基化的相关性研究可反映短期内的血糖代谢水平和DNA甲基化间的相关性，HbA1c与CpG位点甲基化的相关性研究更能反映较长时间的血糖代谢水平和DNA甲基化间的相关性。

双生子是一类特殊的人群，同卵双生子可匹配遗传和早期发育等环境因素[7]，在相同样本量的条件下，结论具有更高的准确性[8]，但目前尚缺乏基于中国等东亚双生子人群的研究证据。利用双生子人群进行的2型糖尿病及血糖相关指标也多在10～30对之间[9-11]，纳入的双生子人群较少，难以发现更微小的关联，因此本研究以我国成年双生子人群作为研究对象，探索与空腹血糖、HbA1c相关的DNA甲基化现象，比较我国与欧美人群间的差异，为糖尿病相关的表观遗传学机制提供进一步证据。

1 资料与方法

1.1 研究人群

双生子研究对象来自中国双生子登记系统(Chinese National Twin Registry，CNTR)， 2013年6月至12月和2017年6月至2018年10月于山东青岛、浙江、江苏、四川及黑龙江五地募集。研究对象纳入标准：(1)年龄为30岁及以上；(2)双生子两成员均参与了现场体检及问卷调查；(3)问卷初筛为同卵双生子；(4)双生子两成员均可提供空腹血标本。本研究经过北京大学生物医学伦理委员会审查并批准(批准号：IRB00001052-13022和IRB00001052-14021)，双生子两成员均签署知情同意书。

1.2 现场调查及实验室检测

利用问卷收集一般人口学信息、疾病相关信息，以及研究中可能存在的吸烟、饮酒等混杂因素信息。问卷信息由经过培训的调查员面对面访问获得，其中糖尿病患病情况及用药史均为研究对象自报，患病情况要求为区县级及以上医院诊断。体格检查包括身高、体质量、收缩压、舒张压等指标，由工作人员使用统一的身高测量仪、体成分测量仪和电子血压计测量。空腹血糖和HbA1c利用现场抽取的空腹静脉血，由指定的第三方实验室采用统一的方法检测，其中空腹血糖采用己糖激酶法检测，HbA1c采用高效液相色谱法检测。

1.3 甲基化检测

利用同一时期收集外周血标本，进行DNA甲基化检测。全基因组DNA甲基化检测均采用Illumina公司生产的芯片检测。研究涉及两种芯片，分别为Illumina Infinium HumanMethylation450 BeadChip(简称450K芯片)和Illumina Infinium MethylationEPIC BeadChip(简称850K芯片)，采用minfi包[12]对数据进行读取，本研究重点分析两种芯片重合的CpG位点。

1.4 甲基化数据读取和质量控制

正式分析前首先对DNA甲基化数据进行标准化及样本和位点的质量控制。利用 R软件WateRmelon程序包[13]，对甲基化数据利用数据驱动的独立标准化(data-driven separate normalization，DASEN)方法进行标准化(控制Ⅰ型、Ⅱ型探针的测量偏倚、背景噪音偏倚和芯片间偏倚)，并检测每个样本中的每个CpG位点是否存在缺失[若CpG位点所在Ⅰ、Ⅱ型探针信号与空白对照对比差异无统计学意义(P>0.01)，则该位点缺失]。位点质量控制包括：(1)剔除缺失率>0.01的位点；(2)剔除本身为单核苷酸多态性(single nucleotide polymorphism, SNP)位点，且在亚洲人群中最小等位基因频率(minor allele frequency, MAF)>0.05的CpG位点；(3)所在基因探针与其他MAF>0.05的SNP位点存在交叉的CpG位点。样本质量控制包括：(1)剔除缺失率>0.01的样本；(2)利用甲基化芯片上的59个SNP位点进行双生子的卵型鉴定，SNP位点相关系数<0.8即可能为异卵的双生子，予以剔除[14]；(3)同一对双生子中一人被剔除，则该对双生子被剔除。最终形成待分析的DNA甲基化数据集。

1.5 血细胞成分估计

由于各个CpG位点的DNA甲基化水平在不同细胞中存在差异，不同个体血液中各个组分占比存在差异，所以需要调整血细胞成分。本研究利用已测得的450K及850K芯片上甲基化数据，采用minfi包中estimateCellCounts函数[15]估计全血及血凝块样本中CD4+T细胞、CD8+T细胞、NK细胞、B细胞、单核细胞、巨噬细胞等有核细胞组分所占的百分比。

1.6 批次效应及其他协变量调整

甲基化检测时，不同的检测平板及实验批次会由于试剂、操作仪器等差异，使得DNA甲基化数据存在一定误差。本研究使用代理变量分析方法(surrogate variable analysis, SVA)调整这些批次效应及其他潜在的混杂因素，在代理变量调整时采用SVA包的sva函数[16]。

1.7 人群分层控制

1.8 统计学分析

本研究所涉及统计分析及图像(包括曼哈顿图及Q-Q图)绘制均使用R 3.5.3软件进行。连续变量采用均数±标准差描述；分类变量采用频数(百分比)描述。若无特殊说明，显著性水平设为P< 0.05。全基因组甲基化与空腹血糖或HbA1c间相关性分析采用混合效应模型(linear mixed effect model，LME)， 分析基于nlme包的lme函数[18]，以空腹血糖或HbA1c分别作为主效应，甲基化水平(β值)作为因变量，将年龄、性别、体重指数(body mass index, BMI)、血压、血细胞组成成分、用SVA包生成的代理变量等连续变量，吸烟、饮酒、是否服用降糖药等分类变量作为协变量纳入固定效应模型，将双生子对编号纳入随机效应模型，截距设置为随机，其他变量为默认设置，进行回归分析，找出分别与空腹血糖或HbA1c相关的CpG位点。采用Bonferroni法对回归结果中的P值进行多重比较校正，显著性水平设为错误发现率(false discovery rate,FDR)< 0.05。

2 结果

2.1 样本质量控制结果及基本信息描述

本研究最终纳入同卵双生子338人(169对)，CpG位点412 459个，其中男性同卵双生子114对、女性55对，平均年龄(48.2±11.9)岁，其他相关变量基本描述如表1，双生子对内血糖及HbA1c水平差异均有统计学意义(P< 0.001)，空腹血糖平均对内差值(1.1±1.8) mmol/L，HbA1c平均对内差值(0.6±0.9)%。

表1 双生子人口学特征及相关变量基本信息

Diabetes mellitus and hypertension were diagnosed by hospitals at or above the county level and reported by patients themselves. The situation of taking hypoglycemic drugs was whether the subjects had taken hypoglycemic drugs within 30 days. FPG, fasting plasma glucose; HbA1c, glycated haemoglobin; BMI, body mass index; SBP, systolic blood pressure; DBP, diastolic blood pressure.

2.2 全基因组DNA甲基化与血糖或HbA1c间相关性分析

全基因组DNA甲基化与血糖或HbA1c间相关性分析将338人(169对)同卵双生子以空腹血糖、HbA1c等血糖相关指标作为主效应，甲基化水平(β值)作为因变量，纳入相关协变量用混合效应模型进行分析，找出与血糖相关的CpG位点。图1展示了CpG位点及其显著性(P值)与染色体之间的位置关系，在与空腹血糖相关分析中，在1、2、4、6、7、8号染色体上发现阳性CpG位点(图1A)，在与HbA1c相关分析中，在1、4、6、7、17号染色体上发现阳性CpG位点(图1B)。

经调整年龄、性别、BMI等因素及进行多重比较校正后，发现与空腹血糖相关位点7个，与HbA1c相关CpG位点10个(表2)，其中cg19693031位点在与空腹血糖及HbA1c相关分析中均是P值最小的位点；共有3个CpG位点在与空腹血糖、HbA1c相关分析中被重复发现(表3)。空腹血糖对应的GCF为1.036，HbA1c对应的GCF为1.014，均在1附近(GCF < 1.1)，Q-Q图(图2)中多数CpG位点与参考线基本重合，仅在末端有较大偏离，显示本研究较好地控制了人群分层。

3 讨论

目前，DNA甲基化是一个疾病机制研究中的新领域，目前已有基于欧美人群的研究发现，DNA甲基化与空腹血糖或HbA1c之间存在相关性，发现多个CpG位点(如cg19693031-TXNIP、cg18881723-SLAMF1、cg05201300-ATP6V0E1、cg01676795-POR、cg00936728-FCRL6、cg00574958-CPT1A、cg07805383-未知基因、cg08309687-LINC0069、cg26262157-PFKFB3、cg12655112-EHD3等)的甲基化可能与空腹血糖或HbA1c水平相关[19-21]，但研究结果之间的差异较大，就目前已知的研究结果而言，仅有位于TXNIP基因上的cg19693031 位点在两篇不同研究中被重复发现[19-20]。

本研究利用收集到的338人(169对)同卵双生子进行全基因组DNA甲基化与血糖指标相关分析。经多重校正后，发现与空腹血糖相关的CpG位点7个, 其中4个CpG位点及其所在基因(cg19693031-TXNIP、cg01538969-DHX16、cg04816311-C7orf50、cg06721411-DQX16)与血糖代谢或2型糖尿病间的相关性已被其他文献报道[20,22-24]。本研究新发现3个与空腹血糖存在相关性的位点(cg08501915、ch.8.1820050F、cg26608667)，其中cg08501915位点虽然在以往的研究中没有被发现与血糖之间有相关性，但有研究发现该位点所在PGRMC2基因和血糖代谢之间可能存在相关性[25]，所以cg08501915位点的甲基化，可能和血糖存在实际的相关性；其余两个CpG位点cg26608667、ch.8.1820050F，分别位于ZFAND2A基因和8号染色体的未知基因上，这些CpG位点与血糖代谢之间的关联尚需进一步探讨。

表2 全基因组DNA甲基化与血糖指标相关性分析

Correlation analysis of whole genome DNA methylation and blood glucose related indicators covariates such as age, gender, BMI, smoking, drinking, blood pressure, hypoglycemic drug use, blood cell composition and surrogate variables generated by SVA package were adjusted. SE, standard error; FDR, false discovery rate; FPG, fasting plasma glucose; 1st exon, first exon; Body, gene body; TSS1500, 200-1500 bases upstream of the transcriptional start site; 3′UTR, 3′ untranslated region; 5′UTR, 5′ untranslated region; -, represents the CpG sites does not target on known genes; Island, the CpG site is on CpG island; N_Shore, ≤2 kb far on the 5′ side; N_Shelf, ≤2 kb from CpG island on the 5′ side; OpenSea, ≥4 kb from any CpG island.

本研究发现与HbA1c相关的CpG位点10个， 其中6个阳性位点(cg19693031-TXNIP、cg04816311-C7orf50、cg01538969-DHX16、cg01676795-POR、cg09029192-TNRC6C、cg16097041-FLAD1)已在其他研究中发现与HbA1c、空腹血糖或糖尿病之间存在相关关系[19-20,22-24]； 4个位点(cg01339781-ZUFSP、cg24667115-BACH2、cg20697417-未知基因、ch.4.1528651F-FRAS1)在之前研究中未发现与糖尿病或血糖相关指标间的相关性，不过BACH2基因上另一个CpG位点cg27644327曾被报道与2型糖尿病及BMI存在相关性[23]，其余3个CpG位点cg01339781、ch.4.1528651F、cg20697417，分别位于ZUFSP、FRAS1基因和1号染色体的未知基因上，其功能尚需进一步探讨。此外，本研究发现的3个与空腹血糖及HbA1c均存在相关性的CpG位点(cg19693031-TXNIP、cg01538969-DHX16、cg04816311-C7orf50)与2型糖尿病及HbA1c等血糖指标间相关性已在先前研究中被多次报道[22-23]。

表3 全基因组DNA甲基化与血糖指标相关分析重叠位点

Correlation analysis of whole genome DNA methylation and blood glucose related indicators covariates such as age, gender, body mass index, smo-king, drinking, blood pressure, hypoglycemic drug use, blood cell composition and surrogate variables generated by SVA package were adjusted.FDR, false discovery rate; FPG, fasting plasma glucose; Body, gene body; 3′UTR, 3′ untranslated region.

本研究发现的TXNIP基因上的cg19693031位点与血糖之间的相关性已在多篇DNA甲基化研究中被报道，该基因表达的蛋白与硫氧还原蛋白相互作用并对其表达和功能起负调控作用[26]，调节氧化还原过程。近年研究发现，TXNIP是细胞内糖转运信号通路的关键节点，该基因的表达产物通过影响GLUT1和GLUT4调节葡萄糖进入细胞[27]。DHX16基因表达产物参与mRNA剪接过程，是细胞代谢和基因表达的重要参与蛋白[28]。C7orf50为位于7号染色体上的开放阅读框，对于其功能的研究较少，尚不清楚该区域的生理功能及其与血糖代谢之间的作用机制。目前已发表的研究[22-23]仅发现DHX16基因及C7orf50区域上的部分CpG位点的甲基化与血糖水平可能存在相关性，但其与血糖转运、代谢之间的联系尚无明确证据。

本研究的优势是相较其他的DNA甲基化相关性研究，研究对象为同卵双生子，同卵双生子共享了全部的遗传物质和早期发育环境，因此具有较高的匹配度。在样本量相当的研究中，双生子研究具有更高的把握度[8]。在利用双生子人群进行的糖尿病及相关指标与DNA甲基化的相关性研究中，本研究具有较大样本量，因此也发现了多个前人研究当中未发现的位点，并且本研究采用更保守的Bonferroni检验，研究结果更具参考价值。

本研究也存在一定局限性：(1)仅对全血样本进行DNA甲基化检测，未对其他组织进行检测。DNA甲基化存在组织特异性，对于血糖相关指标与其他组织中目标位点之间的相关性需要进一步验证。也有研究提示外周血与组织中DNA甲基化具有较高的一致性，由于采血更为便捷，使得外周血成为目前多数甲基化研究的首选材料。(2)仅分析了各个CpG位点的甲基化和空腹血糖、HbA1c等指标间的联系，这些基因转录、翻译的产物的变化及其与人体内血糖代谢状况间的相关性还有待进一步研究验证。(3)为横断面研究，无法确定甲基化与血糖变化在发生时间上先后顺序，因此本研究只能发现CpG位点的甲基化与血糖相关指标间的相关关系，不能研究因果关系。(4)仅分析了单个CpG位点的甲基化与血糖或HbA1c之间的相关性，未考虑多个连续CpG位点组成的区域(CpG岛)的甲基化与HbA1c之间的相关性，可能忽略了多个位点的微小的甲基化与血糖代谢之间的效应联系。(5)仅考虑了单个CpG位点的甲基化和HbA1c及空腹血糖间的相关性，未考虑位点之间及表型之间的交互作用，在后续研究中，可进一步考虑其他影响因素间的交互作用。

综上所述，本研究利用中国双生子人群，全基因组范围内探索与血糖相关指标相关的DNA甲基化位点，发现了与空腹血糖及HbA1c存在相关的DNA甲基化位点。本研究为国内较早地利用双生子人群进行的血糖相关指标与DNA甲基化之间的相关分析，能够为后续研究提供重要的参考，但是本研究所发现的阳性CpG位点的甲基化与近期及远期血糖代谢水平间的相关性，还需要进一步在更大样本量的人群中验证。