基于核心家系全外显子组测序数据探索新生突变与非综合征型唇腭裂的关联

陈 曦，王斯悦，薛恩慈，王雪珩，彭和香，范 梦，王梦莹，武轶群，秦雪英，李 劲，吴 涛△，朱洪平，李 静，周治波，陈大方，胡永华

(1.北京大学公共卫生学院流行病与卫生统计学系，北京 100191；北京大学口腔医学院·口腔医院2.口腔颌面外科，3.儿童口腔科，国家口腔医学中心，国家口腔疾病临床医学研究中心，口腔生物材料和数字诊疗装备国家工程研究中心，口腔数字医学北京市重点实验室，国家卫生健康委员会口腔医学计算机应用工程技术研究中心，国家药品监督管理局口腔生物材料重点实验室，北京 100081)

非综合征型唇裂伴或不伴腭裂(non-syndromic cleft lip with or without cleft palate, NSCL/P)是一类由胚胎发育过程中颌腭组织发育不全或受阻导致的常见出生缺陷[1]。NSCL/P病因复杂，作为一类出生缺陷，遗传病因探索一直是研究热点。既往候选基因研究、全基因组连锁分析和全基因组关联分析(genome-wide association study, GWAS)已经发现了40多个与NSCL/P相关的基因和位点[2-4]，为病因学研究和产前诊断等提供了证据。85%的唇腭裂患者为散发病例，这一流行病学现象提示新生突变(denovomutations, DNM)可能在其发生过程中发挥作用[5]。新生突变是指减数分裂过程中或受精后发生的突变[6]，而非遗传自父母的变异，而核心家系设计能够通过比对双亲与患者同一位点的等位基因直接识别新生突变位点。新生突变是罕见遗传变异的极端情况，平均每人的基因组中约有74个单核苷酸新生突变[7]。由于致病性新生突变在人群中的基因频率很低，常常无法通过全基因组关联研究发现这些位点与疾病的关联，因此需要采用全基因组测序、全外显子组测序等技术手段识别患者的新生突变位点并观察相似表型患者中是否存在一致的新生突变位点。由于新生突变并非遗传自父母，而是在人群中分散存在，因此相关疾病往往呈现散发的特征。此外，由于其所受的自然选择少，因此其致病性更强，这使得新生突变成为导致散发、复杂疾病的主要因素之一[8]。既往研究在多种遗传相关性疾病中发现了具有致病效应的DNM，尤其在神经发育疾病中，如自闭症谱系障碍等[9-10]。此外，先天性心脏病等出生缺陷与DNM是否存在关联也是研究者们关注的重点。近年来利用全基因组测序技术发现了在多种出生缺陷发生中存在致病效应的新生突变位点，如先天性心脏病、歌舞伎面谱综合征等[11-12]，然而非综合征型唇腭裂的新生突变致病效应的相关研究较少见，尚处于起步阶段[5]。因此，本研究基于核心家系设计，利用全外显子组测序数据定位DNM，开展生物信息学分析，以探索非综合征型唇裂伴或不伴腭裂患者中新生突变的发生情况及其潜在致病性。

1 资料与方法

1.1 研究对象

本研究以2016—2018年于北京大学口腔医院募集的22个中国NSCL/P核心家系为研究对象，每个家系由一名患儿和父母双亲组成，共66人。患儿父母均未患唇腭裂。研究开始前获得北京大学生物医学伦理委员会审查批准(IRB00001052-15081)，所有研究对象均签署知情同意书，未成年人由其监护人签署知情同意书。由临床医生负责诊断并进行疾病亚型分类，通过临床检查排除综合征型唇腭裂。

1.2 全外显子组测序检测变异及质量控制

1.2.1遗传变异检测 患儿及其父母每人采集4 mL静脉血，采用盐析法从血细胞中提取DNA样本进行全外显子组测序。测序由武汉华大医学检验所有限公司完成，使用超声波高性能样品处理系统将基因组DNA样本随机打断并进行末端修复。使用Nimble Gen Seq Cap EZ V3(64M)平台进行外显子组捕获，对富集片段进行扩增。扩增产物经Agilent 2100生物分析仪(试剂盒为Agilent DNA 1000 Reagents，Agilent公司，美国)和实时荧光定量PCR(quantitative real-time PCR, qPCR)质量控制合格后上机测序。使用Illumina Hiseq4000测序平台对合格文库进行高通量测序，经Illumina碱基识别软件转化为原始序列数据并以FASTQ文件格式存储，下机数据比对到hg19参考序列，共获得471 695个位点。

1.2.2质量控制 对测序数据进行质量控制：(1)剔除测序深度小于6或大于500的位点；(2)剔除基因型质量(genotyping quality, GQ)小于20的样本位点；(3)剔除基因型缺失率大于5%的位点；(4)剔除性染色体上的位点；(5)剔除IBD(identity by descent)检验非亲生亲子关系的家系。位点和样本的质控过程采用VCF tools 0.1.16(https://vcftools.github.io/index.html)和plink 1.9(http://www.cog-genomics.org/plink2/)进行。

1.3 新生突变的识别及生物信息学统计分析

1.3.1新生突变的识别及注释 研究采用GATK 4.2.0.0(https://gatk.broadinstitute.org/hc/en-us)的PossibleDeNovo模块识别新生突变，该模块通过比对患者与父母在同一位点上的碱基类型判断该位点在减数分裂过程中是否发生了突变，即是否为新生突变位点。采用SnpEff(http://pcingola.github.io/SnpEff/index.html)软件对测序结果进行注释，注释内容包括变异类型、所在区域、在ExAC等公开数据库中的弱势等位基因频率(minor allele frequency，MAF)以及polyphen2、sift(sorting intolerant from tolerate)、CADD(combined annotation dependent depletion)分值等功能预测指标。

1.3.2富集分析 采用R软件包denovolyzeR对DNM进行富集分析[13]，估计存在多个DNM的基因数量和包含超预期DNM数量的基因,前者比较实际数据中存在多个DNM的基因数量与通过5 000次置换获得的DNM经验分布，后者比较实际数据中每个基因存在的DNM数量与预期数量。预期DNM的数量计算方法是，依据功能注释将DNM分为“syn(synonymous)” “mis(missense)” “non(nonsense)” “splice(canonical splice site)” 和“frameshift”5类，计算功能丧失(loss-of-function, lof, lof = non + splice + frameshift)突变和引起蛋白质改变的突变(protein-altering, prot, prot=mis+lof)的数量，从而估计各类别或全部类别的预期DNM数量，并比较是否与实际观察到的数量存在差异。DNM预期数量的计算方法如下：

expected(DNMs)=∑2Pin,

其中Pi为每个功能类别DNM的概率，n为先证者数量。采用泊松(Poisson)检验对预期与实际DNM数量进行比较。由于本次分析包含的基因总数为19 618，因此采用Bonferroni 校正后的P值为1.3×10-6[0.05/(2×19 618)]。

1.3.3蛋白质交互作用分析 用STRING 11(https://string-db.org/)进行蛋白交互作用分析及功能蛋白关联分析。STRING数据库收集蛋白质相互作用的证据，主要包含遗传学证据、实验室证据和文本挖掘证据等，并将这些证据分为基因临接、基因融合、基因共现、共表达、遗传数据、通路信息、文本挖掘7种类型，在每种类型上根据该类型证据的强度对蛋白质间交互作用进行评分，然后进行标化获得用于量化蛋白质间总体交互作用可信度的得分，得分范围为0～1[14]。纳入部分既往研究证据较为充足的NSCL/P基因(EYA1、PAX1、MACF1、MSX1、IRF6、PVRL1、TP63、SUMO1、BMP4、DLX4、MTHFR、CRISPLD2、CLPTM1)和样本中存在能够引起蛋白质氨基酸改变的DNM的基因，以总评分大于0.400(中等可信)作为判断交互作用的标准。计算纳入的全部基因(本研究中发现存在DNM的基因和既往研究与NSCL/P显著关联的基因)之间存在的交互作用数，并检验其是否高于随机一组相同数量基因之间存在的交互作用数，进行GO(gene ontology)及KEGG(Kyoto encyclopedia of genes and genomes)富集分析。

2 结果

2.1 研究对象及测序情况

研究共纳入22名NSCL/P患者及其生物学父母双亲，其中8例患者(男性4例)为非综合征型单纯唇裂(NSCLO)，14例患者(男性11例)为非综合征型唇裂合并腭裂(NSCLP)，对患者及其父母进行全外显子组测序，质控后共计获得339 908个位点，均位于常染色体。

2.2 新生突变

22个NSCL/P核心家系中共识别出345个高置信度DNM。经过注释，内含子区及基因间区位点278个，外显子区位点67个。外显子区位点包含同义突变位点20个及影响蛋白质编码的突变47个(包括错义突变44个，无义突变1个，剪接位点突变2个)，后者分布于46个基因且有37个位点收录于dbSNP数据库。

2.2.1富集分析 由于样本DNM所在的基因中有7个基因在参考数据库中无对应预期值，因此共有60个基因纳入分析。与预期相比，样本中发现的同义突变、错义突变、引起蛋白质改变的DNM及总DNM数均有统计学意义(P<1.3×10-6)，无义突变、剪接位点突变及功能缺失性突变DNM的数量则与预期值无显著差异(表1)。在已知样本基因集中存在的总DNM数及各类别DNM数(n)的情况下，由置换检验获得预期存在多个DNM的基因的数量(expMean)和预期DNM数的最大值(expMax，表2)。样本数据中仅有一个基因ANKRD36C内存在2个DNM，且均为错义突变(rs768682466、rs202176708)。HMCN2包含的功能缺失性DNM数量显著高于预期，4个基因(ADGRL2、ANKRD36C、DIPK2A、HMCN2)所含的引起蛋白质改变DNM数量显著高于预期(表3)。

表1 各功能类别新生突变富集度

表2 包含多个新生突变的基因数量

表3 存在显著富集新生突变的基因

2.2.2蛋白质交互作用 样本DNM中无义突变、错义突变及剪接位点突变共有47个，存在于46个基因，其中有37个位点在dbSNP中有报道。图1展示了已知与NSCL/P关联的13个基因和携带无义突变、错义突变及剪接位点突变DNM的46个基因构成的交互作用网络[15]，线条颜色代表不同的交互作用证据类型，其中“databases”和“experiments”代表蛋白质间交互作用已知，证据等级较高。该交互作用网络中共存在34对显著的交互作用(expected=12, PPI enrichmentP=1.58×10-7)，即观察到的交互作用数量显著高于在全基因组随机选取相同数量的基因所能观察到的交互作用数量。6个包含DNM的基因PRSS3、GSC、SLIT1、RGPD4、PPM1J、MUC5B与已知NSCL/P基因存在蛋白质交互作用，共6组交互作用(表4)。其中，RGPD4与SUMO1交互作用得分为0.868，二者间的交互作用具有高置信度(>0.700)，接近极高置信度(0.900)。此外，MUC5B与TP63存在交互作用，且在MUC家族的MUC6、MUC16及MUC4中也发现了DNM位点。NSCL/P基因SUMO1和BMP4均与2个DNM基因存在蛋白质交互作用。在STRING数据库中对NSCL/P基因和DNM基因进行GO富集分析与KEGG富集分析，未发现具有显著性的富集通路。与NSCL/P基因存在蛋白质交互作用的DNM基因中，位点chr1:113257689和rs751787967(chr9:33795603)来自同一病例。gnomAD东亚人群的基因频率数据如表5所示，rs832357、rs751787967两位点的MAF>0.01，且样本与gnomAD东亚人群的基因频率近似。rs1295794649、rs61734162及rs1470361138在gnomAD东亚人群数据中未检测到弱势等位基因出现，为罕见变异。而位点chr1:113257689未在gnomAD、ExAC及千人基因组计划等人类遗传数据库中收录。

图1 STRING数据库蛋白质交互作用网络分析

表4 STRING数据库蛋白质交互作用得分

表5 与NSCL/P存在交互作用的基因中的新生突变

3 讨论

本研究基于病例-双亲核心家系设计，通过比对NSCL/P患者与双亲在同一位点的等位基因直接识别新生突变，探索其致病效应。经过新生突变富集分析，共发现ADGRL2、ANKRD36C、DIPK2A和HMCN24个基因所含的引起蛋白质改变的DNM数量高于预期值，蛋白质交互作用中发现了6个基因(PPM1J、RGPD4、PRSS3、SLIT1、MUC5B、GSC)存在DNM，并且其编码的蛋白质与已知NSCL/P基因编码的蛋白质存在交互作用，这些基因可能是NSCL/P遗传病因学研究重要的候选基因。在GWAS catalog中对上述基因的DNM位点与NSCL/P的关系进行检索，各位点既往均未见研究报道。在OMIM及GeneCards中对上述10个基因进行检索，发现既往有动物实验研究报道了ADGRL2和GSC编码的蛋白质在胚胎发育或颌面部发育中发挥作用[16-17]，提示二者的新发突变可能与NSCL/P具有较强的关系，其他基因未检索到与唇腭裂发病相关的生物学或基因组学研究证据。

本研究在ADGRL2(adhesion G protein-coupled receptor L2，又名LPHN2或latrophilin2)基因外显子区域共发现2个DNM位点，其中chr1:82417643位于该基因的第8号外显子区域，为错义突变，rs1182776340(chr1:82451002)为同义突变，尚未见对上述位点功能的报道。该基因编码的蛋白质属于黏附类G蛋白偶联受体(G protein-coupled receptors, GPCRs)，这类受体参与细胞结构的构成和多项细胞生命活动，如细胞黏附和迁移、物质运输、细胞分化和凋亡等[18-19]。基于ExAC数据库的全外显子数据库建立的模型“基因杂合功能丧失变异耐受度”(the probability of intolerance to heterozygous pLoF variation, pLI)能够评估基因存在功能丧失突变(即无义突变、移码突变和剪接位点突变)对基因表达的影响程度，pLI≥0.9为lof突变不耐受基因，即lof突变会明显影响基因表达，pLI ≤0.1为突变耐受基因，即lof突变对基因表达影响较小。单倍剂量不足得分则评估一对等位基因中的一个发生突变对基因表达结果的影响。ADGRL2的pLI 得分为1.00[20]，单倍剂量不足得分为0.57%[21]，表明其对功能丧失变异和剂量变化敏感度高。既往研究表明ADGRL家族是一种钙非依赖性受体，对α-latrotoxin的亲和力较低，因此被认为可以调节胞吐作用[22]。动物实验发现ADGRL2纯合缺失小鼠在胎儿期死亡，而杂合子小鼠肌张力降低[23]。既往与ADGRL2直接相关的疾病表型报道较少见，一项全外显子组测序研究首先报道了其与菱脑融合(rhombencephalosynapsis, RES)的关联，并且在胚胎细胞中检测到ADGRL2在多个组织中的活跃表达[16]。RES是一种罕见的发育畸形，多数为常染色体隐性遗传病，患者症状包括发作性气喘、肌张力减低、共济失调、唇裂、腭裂等[24]。既往研究显示，一些综合征型唇腭裂的致病基因同时与NSCL/P存在关联，如van der Woude综合征的致病基因IRF6与NSCL/P的关联在多个独立研究中得到了证实[2,25-27]。ADGRL2基因与综合征型唇腭裂关联提示其可能为非综合征型唇腭裂的重要候选基因。

GSC(Goosecoid)属于同源盒基因家族，其编码的Goosecoid同源蛋白质是一种高度保守的同源结构域转录因子，调节早期发育和胚胎发生，从原肠胚到器官发育过程均有表达[28-29]，在颅面发育中发挥重要作用[17]。一项在中国人群中开展的GWAS研究纳入了7 404 名NSCLP患者和16 059名对照，新发现14个具有全基因组显著性的SNPs，其中rs1243572和rs1243573均位于GSC基因附近[30]。有研究表明，在哺乳动物胚胎的早期发育中，BMP4诱导的EVX1直接抑制GSC表达，GSC和EVX1构成了基因调节网络的核心，能够控制原条样胚胎干细胞的分化方向，通过TGFβ信号调节原条的形成[31]。GWAS研究发现了BMP4基因内多个与NSCL/P关联的位点，BMP4调控GSC表达的证据及二者在胚胎发育中的作用提示GSC内可能存在NSCL/P的致病位点。动物实验发现，GSC基因敲除小鼠胚胎中，多种神经嵴衍生物(如第一鳃弓和鳃裂)及随后分化的组织中(如下颌骨和耳道)出现了转录表达的滞后[28]，小鼠上颌的凸起减少，出现多种颌面部畸形特征[32]。另一项研究发现GSC纯合失活小鼠颅面缺损且出生后24 h内死亡[33]。Parry等[28]在3个SAMS(OMIM：602471)散发病例中通过全外显子组测序发现了GSC基因的纯合截断突变。动物实验及病例研究表明，GSC基因在胚胎和面部发育过程中起到重要作用，提示其可能包含颌面部出生缺陷类疾病的重要致病位点。

除此之外，本研究RGPD4基因与SUMO1间存在强交互作用(0.868)，提示二者在功能上可能存在一定的联系。动物实验表明HMCN2蛋白在脊椎动物发育过程中发挥重要作用[34]。虽然既往研究尚未发现明显支持ADGRL2和GSC以外的8种基因与唇腭裂的发病或颌面部发育存在联系的证据，但本研究结果提示，这些基因可作为候选基因进行进一步的研究和验证。本研究的不足之处在于样本量较小，未对新生突变位点进行Sanger测序确认，后续仍需在大样本中对阳性基因进一步验证，并开展动物实验等对其生物机制进行探索。

综上所述，本研究通过全外显子组测序在22个NSCL/P核心家系中确定新生突变，并筛选可能的致病位点及其所在的基因，经过富集分析和交互作用分析确定候选基因，通过既往生物学机制研究及遗传关联研究证据进一步挖掘候选基因影响唇腭裂发病的可能性。本研究发现的ADGRL2和GSC基因在既往研究中也得到了较好的支持，二者与综合征型唇腭裂或胚胎期颌面部发育的关联提示了发生在这两个基因上的新生突变很有可能具有致病性。而RGPD4、PPM1J等基因缺乏研究证据支持，尚需后续研究的验证。