江苏常熟地区汉族人群47个微单倍型的遗传多态性及遗传结构分析*

潘坤鹏 冯耀森 于文帅 刘宗伟 姚伊人 赵杰 康克莱 张驰 王乐* 吴坚

（1）昆明医科大学法医学院，昆明 650500；2）公安部鉴定中心，法医遗传学公安部重点实验室，现场物证溯源技术国家工程实验室，北京 100038；3）中国人民公安大学侦查学院，北京 100038；4）江苏省常熟市公安局，常熟 215500）

微单倍型是近年来国际法医遗传学界广泛关注的一类新型遗传标记，被定义为200 bp或300 bp的DNA片段内，至少两个单核苷酸多态性（single nucleotide polymorphism，SNP）位点组成的序列多态性遗传标记［1］。由于微单倍型多SNP的特点，往往表现出比单个SNP遗传标记更高的多态性［2］。而且，不同于传统短串联重复序列（short tandem repeat，STR）遗传标记，微单倍型在DNA复制过程中不会因重复序列的复制滑脱现象从而干扰DNA分型［3］，因此，微单倍型被认为有助于混合DNA分析［4］。

随着二代测序（next generation sequencing，NGS）技术的发展，微单倍型遗传标记在法医遗传学领域中的研究逐渐广泛［5-8］。微单倍型在多种场景的应用潜能被研究者们发掘，如祖先推断［9］、混合DNA解析［10］、亲缘关系推断［11］等。为了更好地评估微单倍型的应用效能，研究者们公布了不同微单倍型组合在各自研究人群中的频率和法医学参数等信息［12-13］。但目前基于中国汉族人群微单倍型的群体结构分析仍较少，评估微单倍型遗传多态性或遗传结构的研究多是基于公开访问的千人基因组（1000 Genomes Project，1KG）数据。该数据所包含中国地区样本数据仅301个，由中国南方群体（Southern Han Chinese（CHS），n=105）、中国北京群体（Han Chinese in Beijing，China（CHB），n=103）和中国西双版纳的傣族群体（Chinese Dai in Xishuangbanna（CDX），n=93）三个群体［14-15］组成。同时，低覆盖度的数据很可能导致等位基因的丢失［16］。因此，有必要对中国汉族人群微单倍型深度测序进行遗传多态性和遗传结构分析。

本研究将选用针对案发现场常见的低质量DNA而设计的MHSeqTyper47微单倍型检测试剂盒［16］。该试剂盒已整合的47个微单倍型基因座由Kidd等［17］以及Chen等［18-19］前期报道，其中37个微单倍型在中国汉族人群中被证实混合拆分效能较高［13］。另外10个基因座具有良好的低质量DNA检测能力［16］，并具有理想的有效等位基因数（effective number of alleles，Ae）值。使用试剂盒对223名中国江苏常熟汉族无关个体的47个微单倍型基因分型。统计了取样群体的等位基因频率，描述了江苏常熟汉族人群的遗传多态性特征。本研究还对比了47个微单倍型基于1KG数据的理论效能和实际取样人群数据的差异，并基于Nei氏遗传距离和分化指数探索常熟汉族人群和1KG参考人群的遗传结构。最后基于取样人群的法庭科学参数评估了试剂盒的效能。

1 材料与方法

1.1 DNA提取及定量

223份江苏常熟地区汉族无关个体唾液卡样本由志愿者签署知情同意书后提供，经过公安部物证鉴定中心科研伦理委员会审查后批准，批准号为2019-002。剪取大约25 mm2唾液卡，使用M48磁珠提取纯化试剂盒（Qiagen，德国）提取DNA。使用Qubit®3.0 荧光计对DNA进行定量，将定量后的模板DNA稀释至1 mg/L用于扩增。

1.2 文库制备及测序

文库制备按照MHSeqTyper47试剂盒（公安部鉴定中心，北京）用户手册进行两步PCR，第一步利用多重PCR扩增子捕获技术，扩增目标区域。第二步以连接有唯一标签和测序接头的标签引物扩增第一步PCR产物，获得测序平台兼容的文库。取1 ng DNA与10 µl反应混合物1、4 µl MH引物、补无核酸酶水（Thermo Fisher Scientific，美国）混合得到20 µl PCR1反应混合物。反应混合物在95℃预变性5 min，然后95℃ 30 s、60℃ 2 min、72℃ 2 min进行22个循环，72℃ 5 min，最后4℃保存。取25 µl反应混合物2、1 µl标签7、1 µl标签5、1 µl PCR1产物与22 µl无核酸酶水混合得到50 µl PCR2反应混合物。反应混合物在95℃预变性5 min，然后95℃ 30 s、60℃ 2 min、72℃ 2 min进行15个循环，72℃ 5 min，最后4℃保存。使用MHSeqTyper47试剂盒中包含的纯化磁珠纯化PCR2产物，纯化后文库在7500 PCR仪器（Thermo Fisher Scientific）上使用KAPA文库定量试剂盒（Roche，瑞士）进行定量。将纯化后文库均一化至6 nmol/L并混合到一起，取5 µl混合文库使用5 µl的0.2 mol/L NaOH（Sigma-Aldrich，美国）变性5 min，并使用HT1（Illumina，美国）终止变性和文库稀释。最后，使用MiSeq®Reagent v2试剂盒（Illumina）对600 µl含有90%的12.5 pmol/L文库和10%的12.5 pmol/L PhiX Control V3（Illumina）的混合物进行测序。

1.3 数据分析

通过微单倍型分析软件MHTyper［20］对测序结果进行分析，使用最低测序深度为10，基因座内等位基因测序深度和基因座内最高等位基因测序深度比值0.1作为阈值过滤可能来源于污染或噪声的测序读取序列。

分析得到的微单倍型序列通过MHTyper软件分配的数字等位基因进行命名，这将兼容下游法医学参数计算及群体遗传学分析软件。命名流程参考Song等［21］的研究，即首先将组成一个微单倍型的SNP按照其在染色体上的位置进行排序，然后根据dbSNP数据库（https://www.ncbi.nlm.nih.gov/SNP/index.html）列举出每个SNP的等位基因，接着将SNP的所有等位基因进行组合得到所有可能的微单倍型等位基因，最后将微单倍型的所有等位基因组合按照字母表排序，并以阿拉伯数字序号作为等位基因的数字命名。

使用数字等位基因通过Alrequin v3.5软件进行Hardy-Weinberg平衡以及连锁不平衡检验［22］，使用PowerStats v1.2电子表格软件计算法庭科学参数［23］。使用python 3编写的脚本计算Ae和累积随机匹配概率（combine matching probability，CMP）。线性回归分析由线上软件SPSSPRO（https://www.spsspro.com/）完成。成对遗传分化指数（fixation index，Fst）与Nei氏遗传距离使用数字等位基因通过GenAlEx v6.5计算并通过MEGA Version 11.0.13绘制分子进化树。

2 结果

2.1 测序参数和样本数据质量

223份样本使用3张PE MiSeq®V2流动池（Illumina）完成测序，每张芯片分别检测73、73和77个样本。3次测序中的平均簇通过率分布在77.90%至90.21%之间（表1），与同类测序相近［24-25］，而3次测序的Q30比例分布在72.49%到84.76%之间。3张芯片均得到了足以用于后续分析的读取序列数。

Table 1 Sequencing metrics of three runs for 223 unrelated individuals from Changshu, Jiangsu

排除人工检查判断为测序失败和污染的样本，最终216份样本的数据用于后续数据分析。单个样本总序列数分布在11 004到402 195之间，其中，91.67%（n=198）的样本总序列数在40 000到200 000之间。单样本平均有效序列比例分布在0.931到0.992。图1a展示不同总序列数区间的样本数目及总序列数与有效序列比例的关系，可以观察到有效序列比例随着总序列数而升高。各基因座216个样本的平均等位基因覆盖比（allele coverage ratio，ACR）分布在0.835到0.912（图1b）。

Fig.1 Evaluation of sequencing data

2.2 等位基因频率

在216份江苏常熟汉族无关个体的数据中，共观察到47个微单倍型的239个不同的等位基因，不同微单倍型的等位基因数量从3个（mh03KK-006等3个微单倍型，表2）到13个（mh01KK-117）不等，各等位基因的频率分布在0.002 3~0.664 4之间。在对1KG中的2 504份无关个体的分型统计中，47个微单倍型共观察到326个等位基因。通过对比，47个微单倍型在本研究数据中有4个等位基因在1KG的2 504份个体数据中未观察到，分别是：mh04KK-013的等位基因24、mh06CP-003的等位基因3、mh13KK-213的等位基因7，以及mh22KK-061的等位基因7，上述等位基因均在杂合子中观察到。

Table 2 Allele frequency distribution of 47 microhaplotypes in 216 unrelated individuals from Changshu, Jiangsu

2.3 对比1KG数据的微单倍型多态性分析

通过计算1KG中东亚（EAS）参考人群中5个群体，即CHS群体、CDX群体、越南胡志明市群体（KHV，n=99）、CHB群体以及日本东京市群体（JPT）的Ae值，216名江苏汉族个体计算的Ae与5个群体数据计算的微单倍型的Ae值分布均表现出较小的差异（F检验结果水平显著，P值均小于0.05，图2a及图S1）。常熟汉族人群的Ae作为自变量，CHB群体的Ae作为因变量的线性回归分析显示，CHB与常熟汉族人群的预测模型拟合程度最高（R2=0.903，图2a及图S1），与CDX群体相关性程度最低（R2=0.614，图S1）。另一方面，共18个微单倍型的Ae值在江苏汉族人群中高于CHB群体。

基于216名江苏汉族无关个体数据与EAS参考人群中的5个群体数据计算的CMP存在较小差异（图2b）。得到47个微单倍型在常熟汉族人群中的CMP达到1.25×10-36。CHS、CDX、KHV、CHB和JPT群体计算的CMP分别为：4.51×10-38、4.10×10-36、6.64×10-37、1.41×10-37和1.48×10-37。

2.4 基于1KG的群体遗传分化分析

进一步计算了江苏常熟汉族人群和1KG中5个参考人群，即欧洲人群（EUR）、EAS、美洲人群（AMR）、南亚人群（SAS）和非洲人群（AFR）之间的遗传分化指数，江苏汉族人群与EAS参考人群之间的Fst值为0.001 6，群体之间存在最小遗传分化（Fst<0.05）。与AFR参考人群之间观察到的分化程度最大Fst=0.065 7，两个人群间存在中等遗传分化。这与EAS参考人群和AFR参考人群的结果相似（Fst=0.063 7）。同样地，与1KG中26个群体之间的遗传分化距离显示（图2c），江苏汉族人群与CHB群体之间的遗传分化指数最低（Fst=0.006 8）。与AFR参考人群中的Mende in Sierra Leone（MSL）群体之间达到最大值（Fst=0.337 2）。

成对Nei氏遗传距离显示（图2d），26个参考人群及江苏汉族216个个体的人群共形成两个主要分支：AFR参考人群的7个群体形成第一个分支；AMR参考人群的4个群体、欧洲EUR参考人群的5个群体以及EAS参考人群和SAS参考人群与江苏汉族人群形成第二个分支。在第二个分支中，江苏汉族人群与CHB群体，CDX群体以及越族EAS参考人群中的KHV群体共同形成一个分支。江苏汉族人群与CHB群体之间的Nei氏遗传距离最小，为0.006 8。

Fig.2 Comparison of 1000 Genomes Project data and 216 unrelated individuals data in Changshu, Jiangsu

2.5 法庭科学参数

基于216名常熟汉族无关个体数据计算的47个微单倍型的基因型分布均符合Hardy-Weinberg平衡（P＞0.05），所有配对基因座LD检验P值范围为0.004 6~0.933 9（P＞0.05），矫正值为0.000 2（公式为0.5/（n×（n-1）），即微单倍型个数乘以微单倍型个数减去1，用0.5除以得到的数值，n为微单倍型个数），均处于连锁平衡中。表3列出了各个微单倍型的法庭科学参数。47个常染色体微单倍型的Ae均处于1.908 6（mh04KK-010）到4.842 6（mh21KK-320）之间，有6个微单倍型的Ae值大于4，平均Ae值为3.194 4。期望杂合度（expected heterozygosity，He）范围为0.476 0（mh04KK-010）到0.793 5（mh21KK-320），平均He为0.674 5；个体识别概率（discrimination power，DP）在0.453 6~0.917 0之间，平均值为0.813 9，非父排除概率（probability of exclusion，PE）在0.160 7~0.609 4之间，平均值为0.390 9，累积非父排除概率（combine probability of exclusion，CPE）为1-3.59×10-11。

Table 3 Forensic parameters of 216 unrelated individuals in Changshu, Jiangsu

Continued to Table 3

3 讨论

目前，MHSeqTyper47试剂盒已经完成了低质量DNA检测中的效能评估［16，26］。然而，微单倍型的筛选和评估更多的是基于人类遗传项目数据［18，27-28］，可能会由于缺少个体基因型或测序数据覆盖度低等因素而限制分析。因此，探索高测序深度下微单倍型在不同人群中的效能将起到重要作用。本研究使用MHSeqTyper47试剂盒对江苏常熟地区汉族人群进行微单倍型遗传多态性检测分析，统计了江苏常熟汉族代表人群频率数据。另外，对比1KG数据，将为评估实际采集样本与基于数据库的体系理论效能差异提供参考。

早期的研究中统计了低质量样本的测序参数，在对微量、降解、模拟案件现场样本和抑制剂影响的文库的3次测序中，簇通过率分布在73.0%到90.3%之间，3次测序的碱基质量大于Q30的比例分布在75.0%到86.2%之间［16］。与本研究检测的人正常唾液卡样对比，并未观察到测序参数受样本质量的影响。然而，测序的簇通过率和Q30将直接影响测序所得序列数［24，29］。本研究中绝大部分样本得到了足够后续分析的序列数（单样本总测序深度在11 004以上）。

江苏常熟地区汉族人群MHSeqTyper47的CMP与EAS参考人群中的5个群体均结果相近（图2b），但明显高于基于1KG 2 504名样本代表的全球水平的8.70×10-44。这可能是因为受人群之间的特有等位基因影响，比较5个参考人群分别的CMP可以证实这一现象，如基于5个参考人群EUR、EAS、AMR、SAS和AFR数据计算所得CMP分别为2.04×10-37、1.64×10-38、1.58×10-39、1.97×10-40和1.45×10-37。常熟地区汉族人群Ae值与1KG中EAS参考人群的5个群体间并没有显著的差异，但结果表明，常熟汉族人群的Ae分布与CHB群体之间相似性最高，这将说明，CHB群体更好的代表常熟汉族人群的多态性。另外，对比现有商业化STR试剂盒，47个微单倍型在常熟地区216名汉族无关个体的CMP低出GoldeneyeTM20 A试剂盒19个常染色体STR在江苏地区10 000名个体中的数值（1.16×10-23）［30］。

为了评估1KG中EAS参考人群与江苏常熟汉族人群的遗传分化程度，成对的Fst显示：常熟汉族人群与EAS参考人群较为相似，而与AFR和EUR参考人群存在较大的遗传分化。群体间Nei氏遗传距离也表现出相似的结果。综上所述，相比于1KG全部群体数据，基于1KG中EAS参考人群数据进行微单倍型筛选或评估更符合中国人群遗传特征。

4 结论

本研究报告了47个微单倍型基因座在江苏常熟地区汉族人群中的等位基因频率及遗传多态性信息，为47个微单倍型在法医学应用中提供了数据基础。另外，比较了1KG参考人群与江苏常熟地区汉族人群的多态性差异，并揭示了47个微单倍型在江苏常熟地区汉族人群中的遗传结构。总的来说，1KG中的东亚人群参考数据更符合江苏常熟地区汉族人群的遗传特征。

附件见本文网络版（http://www.pibb.ac.cn或http://www.cnki.net）：PIBB_20230027_Figure_S1.pdf