江蘇常熟地區漢族人群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個微單倍型在江蘇常熟地區漢族人群中的遺傳結構?？偟膩碚f，1KG中的東亞人群參考數據更符合江蘇常熟地區漢族人群的遺傳特征。

附件見本文網絡版（http://www.pibb.ac.cn或http://www.cnki.net）：PIBB_20230027_Figure_S1.pdf