基于核心家系全外顯子組測序數據探索新生突變與非綜合征型唇腭裂的關聯

陳 曦，王斯悅，薛恩慈，王雪珩，彭和香，范 夢，王夢瑩，武軼群，秦雪英，李 勁，吳 濤△，朱洪平，李 靜，周治波，陳大方，胡永華

(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]?；贓xAC數據庫的全外顯子數據庫建立的模型“基因雜合功能喪失變異耐受度”(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等基因缺乏研究證據支持，尚需后續研究的驗證。