宋世文， 盛桂蓮,2)*， 袁俊霞， 肖 博， 胡家銘， 鄧妙璇， 侯新東， 孫國(guó)江， 王林英， 賴旭龍2),

(1)中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，武漢 430078；2)中國(guó)地質(zhì)大學(xué)生物地質(zhì)與環(huán)境地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430078；3)中國(guó)地質(zhì)大學(xué)(武漢)材料與化學(xué)學(xué)院，武漢 430078；4)中國(guó)地質(zhì)大學(xué)(武漢)地球科學(xué)學(xué)院，武漢 430074)

在過(guò)去的十年中，古DNA(ancient DNA, aDNA)領(lǐng)域每年公布的古基因組數(shù)量都在快速增加，古代DNA數(shù)據(jù)集序列的完整性也在逐步提高[1]。古DNA研究已在諸多領(lǐng)域取得了巨大進(jìn)展，例如：家養(yǎng)動(dòng)物的起源與馴化[2-4]、人類起源以及生活方式[5-8]、晚更新世以來(lái)大型動(dòng)物的滅絕原因探究[9-11]等。從古老的生物材料中得到的古DNA序列能提供絕滅物種的遺傳信息[12]，豐富人們對(duì)過(guò)去環(huán)境氣候變化以及人類文明發(fā)展等方面的認(rèn)知。

然而，古代樣品中的DNA保存情況仍然制約著古DNA研究領(lǐng)域的發(fā)展。生物體死亡后，DNA損傷修復(fù)機(jī)制隨之崩壞，水解、氧化及微生物作用將直接造成DNA雙螺旋主鏈的斷裂及分子的高度片段化[12]，古代樣本中通常僅含有極微量并具有高度損傷的內(nèi)源性DNA[13,14]，進(jìn)而阻礙對(duì)物種序列信息的挖掘。目前普遍認(rèn)為古DNA損傷有兩大機(jī)制，一是去嘌呤作用，使DNA鏈容易在高嘌呤位置處斷裂形成小片段，但是其分子機(jī)制目前尚不清楚；另一種機(jī)制是胞嘧啶脫氨基作用，使得胞嘧啶轉(zhuǎn)化為尿嘧啶并最終導(dǎo)致測(cè)序錯(cuò)誤(C→T，G→A)[15]，由于古DNA片段末端容易形成不穩(wěn)定的“單鏈懸垂”結(jié)構(gòu)，脫氨基作用在DNA鏈末端發(fā)生頻率更高[16]。在各種關(guān)于古DNA的研究中，DNA損傷被頻頻提及，古DNA損傷被用于和現(xiàn)代DNA損傷模式進(jìn)行對(duì)比，以驗(yàn)證所得到的序列信息屬于古代個(gè)體而非來(lái)源于現(xiàn)代污染[17]，但不同樣品中古DNA的損傷程度是否存在差異、古DNA損傷程度是否受特定環(huán)境因素影響等研究相對(duì)較少[18,19]。

近年來(lái)，中國(guó)出土的古人類及動(dòng)物遺存中的古DNA研究呈上升之勢(shì)[6,20,21]，特別是中國(guó)東北地區(qū)出土了大量DNA保存較好的古脊椎動(dòng)物化石，其中以猛犸象-披毛犀動(dòng)物群化石(Mammathas-Coelodontafauna)最具代表性，該動(dòng)物群化石對(duì)于第四紀(jì)晚期大型動(dòng)物的地理分布、遷徙路線及絕滅機(jī)制等研究具有重要價(jià)值[22-24]。結(jié)合中國(guó)不同地理位置特定的環(huán)境特征，對(duì)中國(guó)古代生物材料中殘留的古DNA損傷模式及相關(guān)影響因素的探究，是有效發(fā)掘中國(guó)古代生物材料中蘊(yùn)含的遺傳信息、解決多個(gè)生物類群演化疑團(tuán)的必要工作。本文對(duì)采自中國(guó)東北的 33 個(gè)古脊椎動(dòng)物化石及亞化石樣品中進(jìn)行二代測(cè)序(next generation sequencing, NGS)，對(duì)得到的古DNA序列進(jìn)行堿基損傷分析，將各樣品中古DNA的損傷程度與其埋藏時(shí)間、埋藏的地質(zhì)時(shí)期、建庫(kù)方法、材料類型等因素相結(jié)合進(jìn)行系統(tǒng)分析，嘗試確定影響中國(guó)東北地區(qū)古生物樣品DNA損傷的關(guān)鍵因素，為探明該地區(qū)樣品中的古DNA分子特征、有效選取古DNA研究材料提供更多參考。

1 材料與方法

1.1 樣品地點(diǎn)及年代信息

本研究用到的二代測(cè)序數(shù)據(jù)來(lái)源于中國(guó)東北第四紀(jì)晚期以來(lái)9種哺乳動(dòng)物共計(jì)33件化石和亞化石樣品，其中15個(gè)樣品序列來(lái)自本研究組已發(fā)表文章[25-28]，18個(gè)樣品序列尚未發(fā)表。33個(gè)樣品的地點(diǎn)信息正如Fig.1所示，包括骨骼化石或亞化石25個(gè)、牙齒化石7個(gè)、角化石1個(gè)。所有樣品年代信息由美國(guó)BETA實(shí)驗(yàn)室(Beta Analytic, USA)通過(guò)放射性碳測(cè)年(radiocarbon dating)方法獲得，具體信息正如Table 1所示。古DNA提取及建庫(kù)方法參照Meyer等已有研究[29]。

Table 1 Sample information table

Fig.1 Schematic diagram of sample distributions The circular size represents numbers of samples.The colors represent different Marine isotope stages.Green: 11 ka→now, Blue: 24 ka→11 ka; Red:60 ka→24 ka; Yellow: >48 000 BP

1.2 原始測(cè)序數(shù)據(jù)的生物信息學(xué)處理

測(cè)序得到的原始文件為FASTQ格式，使用“cutadapt 1.12”[30]分別對(duì)每個(gè)DNA文庫(kù)測(cè)序文件切除接頭序列，并過(guò)濾掉長(zhǎng)度低于30 bp的短片段，再使用“BWA 0.7.15”(Burrows-Wheeler Alignment)軟件[31]中的“aln”和“samse”功能將片段與參考線粒體基因組比對(duì)，舍棄比對(duì)質(zhì)量得分低于30的序列；使用SAMtools v1.3.1軟件[32]中的“view”和“sort”算法，將剩余序列按5′ 端位置在參考基因組上進(jìn)行排序，通過(guò)“rmdup”去除可能存在的PCR重復(fù)序列，用“merge”算法得到每一文庫(kù)的bam格式文件。

1.3 序列損傷檢測(cè)

對(duì)古DNA分子堿基損傷的評(píng)估，使用各樣品所有文庫(kù)合并處理后的bam文件作為輸入文件，使用“mapDamage”軟件包[33]統(tǒng)計(jì)片段末端堿基損傷頻率與每個(gè)樣品的片段分布，將樣品DNA片段的末端堿基替換率(5′ C→T frequency)以及平均片段長(zhǎng)度(average fragment length)作為損傷分析的指標(biāo)[15]。

1.4 統(tǒng)計(jì)學(xué)分析

按照深海氧同位素階段(marine isotope stages，MIS)將樣品年代劃分為不同地質(zhì)時(shí)期[34]：11 ka→現(xiàn)在，MIS 1; 24 ka→11 ka, MIS 2; 60 ka→24 ka, MIS 3。獲得的堿基損傷信息按照樣品埋藏的地質(zhì)時(shí)期、建庫(kù)方法、化石類型分別進(jìn)行分組分析。為減少偶然誤差，每組含有足夠的樣品數(shù)(n≥ 5)，通過(guò)Shapiro wilk test檢驗(yàn)每組數(shù)據(jù)的正態(tài)性，使用ANOVA方差分析檢測(cè)不同樣品類型、不同建庫(kù)方法得到的測(cè)序文庫(kù)的末端堿基替換率和古DNA平均片段大小，以及不同地質(zhì)時(shí)期古DNA平均片段大小的差異。通過(guò)Kruskal Wallis test分析不同地質(zhì)時(shí)期末端堿基替換率的差異性。最后使用線性回歸統(tǒng)計(jì)分析每個(gè)樣品中堿基損傷情況與埋藏時(shí)間的相關(guān)趨勢(shì)。統(tǒng)計(jì)分析采用SPSS 26.0軟件(IBM SPSS Statistics，IBM公司Armonk，紐約)完成，以P值 < 0.05為顯著性標(biāo)準(zhǔn)。

2 結(jié)果

2.1 堿基損傷

通過(guò)“mapDamage”軟件包得到各樣品中古DNA片段的堿基損傷結(jié)果，隨機(jī)抽取2個(gè)可視化結(jié)果正如Fig.2所示。樣品5′端C→T的替換率以及DNA的平均片段長(zhǎng)度結(jié)果正如Table 2所示。

Table 2 End replacement rates and average fragment lengths of the specimens

Fig.2 DNA damage characteristics of samples ZDT9 and CADG542 (A)DNA base frequencies shown by different colors: blue-A, green-C, black-G, red-T;(B)Fragment lengths distribution of two different DNA strands, black & white histograms-single-end DNA strands, colored histograms-double-end DNA strands.Results are generated after the ancient DNA sequencing fragments aligning with the reference sequences.The merged bam file is used to compare with the reference genomes.The maximum read length is 70 nt, and the genomic window upstream and downstream of the read point is 10 nt

2.2 末端堿基替換率與片段大小無(wú)明顯相關(guān)性

對(duì)各樣品古DNA片段進(jìn)行末端堿基替換率與平均片段長(zhǎng)度分布情況統(tǒng)計(jì)，結(jié)果正如Fig.3所示。平均片段長(zhǎng)度絕大多數(shù)(97%)處于45～75 bp之間(Fig.3A)。所有樣品堿基錯(cuò)配概率在0.05～0.35區(qū)間內(nèi)，并且91%的樣品(30個(gè))錯(cuò)配率低于0.2，錯(cuò)配率超過(guò)0.2的樣品僅占9%(3個(gè), Fig.3B)，這3個(gè)樣品(CADG449/DLH2/TH3)都出自黑龍江哈爾濱市。此外，片段長(zhǎng)度在45～55、55～65、65～75 bp三個(gè)區(qū)間內(nèi)樣品的末端堿基替換率未發(fā)現(xiàn)明顯區(qū)別(Fig.3C)。

Fig.3 Statistical results of ancient DNA molecules characteristics for 33 samples (A)Number of samples in different average fragment lengths;(B)Number of samples in different 5′ C→T frequencies;(C)Differences of 5′ C→T frequencies in different fragment length ranges.In(C), values are the mean ± SD of multiple(n ≥ 5)determinations from separate experiments.One-way analysis of variance(ANOVA)followed by post hoc LSD test is applied to determine the significant differences in three fragment length ranges.P < 0.05 indicates significant difference in three fragment length ranges

2.3 古DNA末端堿基替換率隨埋藏時(shí)間增加升高

將樣品的年代信息(Table 1)與mapDamage軟件得到的古DNA損傷指標(biāo)(Table 2)相結(jié)合，構(gòu)建古DNA損傷水平隨埋藏時(shí)間變化的趨勢(shì)圖(Fig.4)。運(yùn)用SPSS軟件對(duì)年代與DNA損傷水平進(jìn)行線性回歸分析，結(jié)果表明古DNA的末端堿基替換率與樣品的埋藏時(shí)間顯著相關(guān)，末端堿基替換率隨著埋藏時(shí)間的增長(zhǎng)而升高(Fig.4A)，而平均片段長(zhǎng)度與埋藏時(shí)間無(wú)明顯關(guān)系(Fig.4B)。

2.4 MIS1時(shí)期古DNA末端堿基替換率最低

根據(jù)年代信息將樣品劃分為MIS1(n= 6)、MIS2(n= 5)和MIS3(n= 14)三組，分析組間樣品古DNA損傷水平差異。不同時(shí)期間古DNA末端堿基替換率和平均片段大小的差異水平并不相同(Fig.5)，末端堿基替換率在3個(gè)時(shí)期間具有極顯著的差異，MIS1與MIS3兩組間差異最大(Fig.5A，supplementary information Table 3)，而平均片段大小則未見(jiàn)明顯區(qū)別(Fig.5B)。

2.5 古DNA損傷情況在不同樣品類型間無(wú)明顯差異

選擇MIS2和MIS3時(shí)期的樣品，根據(jù)樣品材料不同將樣品分為牙齒(n= 5)和骨骼(n= 14)兩組，分析不同材料類型對(duì)古DNA損傷程度的影響。結(jié)果顯示不同的材料類型對(duì)古DNA損傷的影響未見(jiàn)顯著差異(Fig.6)。

2.6 不同建庫(kù)方法對(duì)古DNA損傷無(wú)顯著影響

選擇MIS2和MIS3時(shí)期的樣品，按構(gòu)建測(cè)序文庫(kù)方法的不同將所有樣品分為單鏈(n= 9)和雙鏈(n= 9)兩組，比較兩組間堿基損傷水平差異。結(jié)果正如Fig.7所示，不同的文庫(kù)構(gòu)建方法對(duì)古DNA末端堿基替換率和平均片段大小均無(wú)影響(Fig.7A, Fig.7B)。

3 討論

3.1 不同埋葬環(huán)境對(duì)古DNA末端替換率的影響

本研究中的古DNA樣品采集自地理位置相對(duì)較近的區(qū)域，但是相近地區(qū)間的堿基損傷情況仍有差異，主要體現(xiàn)在來(lái)自哈爾濱地區(qū)的樣品(CADG449/DLH2/TH3)表現(xiàn)出較高的末端堿基替換率(Fig.3B)。對(duì)東北各個(gè)樣品出土地的氣候條件進(jìn)行對(duì)比分析，并未發(fā)現(xiàn)在這幾個(gè)地區(qū)間存在較大氣候差異。然而在自然地理方面，相較于東北其他地區(qū)，哈爾濱地區(qū)水系更加發(fā)達(dá)，該地區(qū)的樣品出土點(diǎn)大多靠近松花江段，可能是該地區(qū)三個(gè)樣品末端堿基替換率較高的原因之一。已有研究表明，在化石形成的初期，樣品埋藏環(huán)境中較高的溫度與較低的含水量能夠幫助樣品脫水，從而形成更加有利于DNA保存的生物化石[35,36]。由于古DNA片段末端部分降解后，兩條鏈長(zhǎng)度不一，容易形成單鏈懸垂結(jié)構(gòu)，在含水量較高的環(huán)境中，古DNA懸垂的單鏈可能更容易發(fā)生DNA呼吸作用，該作用可能導(dǎo)致C→T的增高[16]，土壤中較高含水量因此成為造成古DNA損傷的潛在原因之一。

3.2 DNA損傷程度與埋藏時(shí)間的線性關(guān)系

已有研究顯示，埋藏環(huán)境中不同的pH、離子強(qiáng)度、腐殖酸等都會(huì)影響化石中古DNA的保存[14]，古DNA的堿基損傷水平與樣品埋藏時(shí)間的相關(guān)性比較模糊[37,38]。本研究中，來(lái)自哈爾濱地區(qū)的3個(gè)樣品埋藏時(shí)間相對(duì)較長(zhǎng)，都超過(guò)了4萬(wàn)年，其堿基替換率也相對(duì)較高(Fig.3B)，中國(guó)東北地區(qū)古脊椎動(dòng)物樣品古DNA分子，其末端堿基替換率與埋藏時(shí)間有著極顯著的線性關(guān)系(P<0.01)，其片段大小則與埋藏時(shí)間無(wú)明顯關(guān)聯(lián)(Fig.4)。該結(jié)果與Sawyer等人對(duì)歐洲與非洲出土的六萬(wàn)年內(nèi)化石樣品(包括馬、牛、靈長(zhǎng)類)的堿基損傷模式研究結(jié)果一致[15]，即埋藏時(shí)間不會(huì)影響古DNA的平均片段大小，但在古DNA分子的5′端，堿基C→T取代的頻率會(huì)隨著時(shí)間的推移而增加，說(shuō)明埋藏時(shí)間對(duì)古DNA分子的影響主要體現(xiàn)在末端堿基損傷的增加而非片段長(zhǎng)度的變化。

Fig.4 Correlation between DNA damage and different burial times for 33 samples (A)5′ C→T frequencies of ancient DNA fragments;(B)Average fragment lengths of ancient DNA fragments.Each dot represents one sample in this study.The linear relationship between DNA damage level and burial times is calculated by ANOVA and used to determine whether the equation is valid.P < 0.05 indicates a statistically significant correlation between DNA damage level and the time of burial

3.3 不同地質(zhì)年代間損傷情況的差異

在不同的地質(zhì)時(shí)期組間的古DNA分子5′端C→T替換率存在著極顯著差異(Fig.5A)，兩兩比較發(fā)現(xiàn)MIS1/MIS3兩組間差異最大(P<0.01，supplementary information Table 3)，而其他兩個(gè)地質(zhì)時(shí)期間，即MIS1/MIS2、MIS2/MIS3之間的損傷水平未見(jiàn)明顯差異。MIS1時(shí)期的樣品年代大多處于大暖期(megathermal)，氣候條件比較均一；整體而言，MIS2屬于末次冰期寒冷的階段，MIS3屬于溫暖的階段，這兩個(gè)階段的氣候變化都比較頻繁且劇烈[39,40]，屬于這兩個(gè)時(shí)期的樣品可能會(huì)遭遇到相似的氣候環(huán)境，兩個(gè)時(shí)期間溫度條件的差異對(duì)古DNA損傷造成的影響受到一定程度的削弱，這可能是MIS2與MIS3之間古DNA損傷水平無(wú)差異的原因。古DNA損傷水平在MIS1/MIS3兩組間存在著明顯差異，而在MIS1/MIS2卻無(wú)明顯差異，或許說(shuō)明溫度帶來(lái)的影響并沒(méi)有預(yù)期的那么重要，在不同地質(zhì)時(shí)期間樣品埋藏地點(diǎn)的其他因素(含水量、腐殖酸含量等)可能對(duì)DNA的損傷做出了更大的貢獻(xiàn)。

Fig.5 Statistical results of correlation between DNA damage and different geological periods (A)5′ C→T frequencies distribution for samples in different geological periods;(B)Average fragment length distribution for samples in different geological periods.Values are the mean ± SD of multiple(n ≥ 5)determinations from separate experiments.The Kruskal Wallis test and pairwise comparisons are used to determine whether there are significant differences between different geological time periods,**P < 0.05 indicates significant difference in geological periods

3.4 樣品類型及不同建庫(kù)方法對(duì)古DNA損傷的影響

古代樣品材料中DNA的保存情況與材料孔隙大小有很大聯(lián)系[41]。耳骨是顳骨的一部分，是哺乳動(dòng)物體內(nèi)最堅(jiān)硬、密度最大的骨頭，被認(rèn)為是古DNA保存最好的部位，但其對(duì)于形態(tài)學(xué)和鍶同位素鑒定具有重要意義，用于進(jìn)行古DNA測(cè)序的情況較少。除耳骨外，牙齒樣品中古DNA的保存情況相對(duì)較好[42,43]，其他骨骼化石(腿骨、掌骨等)古DNA保存情況相對(duì)較差。牙齒最外層的牙釉質(zhì)是六棱柱狀的釉柱及其間質(zhì)規(guī)則排列形成的致密結(jié)構(gòu)，羥基磷灰石占據(jù)絕大部分，有機(jī)質(zhì)含量極低，化石形成過(guò)程中能夠有效地降低外界因素對(duì)于牙齒內(nèi)DNA的損害；腿骨、掌骨等骨骼有機(jī)質(zhì)相對(duì)較高，所形成的化石更容易出現(xiàn)疏松多孔的結(jié)構(gòu)。為了減小時(shí)間和其他因素(溫度、降水量)帶來(lái)的干擾，本文選擇了MIS2和MIS3時(shí)期的化石，用其分析不同樣品類型及不同建庫(kù)方法對(duì)古DNA損傷的影響，因?yàn)檫@兩個(gè)時(shí)期間的樣品的DNA損傷情況并未顯示出明顯區(qū)別(supplementary information Table 3)。本研究未顯示出骨骼和牙齒樣品間古DNA的損傷情況存在不同(Fig.6)，這在很大程度上受到材料選擇的影響，因?yàn)楣趋阑尤菀仔纬墒杷山Y(jié)構(gòu)，而疏松結(jié)構(gòu)的骨骼化石中通常并不能提取出DNA，本研究進(jìn)行分析的骨骼化石皆是具有致密結(jié)構(gòu)、已經(jīng)得到古DNA序列的樣品，這也說(shuō)明了具有致密結(jié)構(gòu)的骨骼化石對(duì)DNA的保護(hù)效果并不弱于牙齒化石。

Fig.6 Statistical results of correlation between DNA damage and different fossil types (A)5′ C→T frequencies of ancient DNA from tooth and bone materials;(B)Average fragment lengths of ancient DNA from tooth and bone materials.Values are the mean ± SD of multiple(n ≥ 5)determinations from separate experiments, One-way analysis of variance(ANOVA)is employed to determine the significant differences between fossil teeth and bones.P < 0.05 indicates significant difference between fossil teeth and bones

相比于構(gòu)建雙鏈建庫(kù)，單鏈建庫(kù)方法中增加了尿嘧啶-DNA糖基化酶與核酸內(nèi)切酶 VIII去除樣品中脫氨胞嘧啶的實(shí)驗(yàn)步驟，被認(rèn)為能夠有效降低測(cè)序讀取出的DNA損傷水平[15,18,44]，但本文結(jié)果顯示，不同的建庫(kù)方法對(duì)顯示出的古DNA末端堿基替換率幾乎無(wú)影響，平均片段大小同樣與建庫(kù)方法無(wú)關(guān)(Fig.7)。古DNA樣品的獨(dú)特性或許對(duì)結(jié)果造成了一定影響，此外尿嘧啶的去除會(huì)形成更多的DNA短片段(< 30 bp)，而短片段在古DNA材料中占據(jù)相當(dāng)大的比例，本文在數(shù)據(jù)處理中過(guò)濾了小于30 bp的DNA片段，這可能是造成本研究結(jié)果中DNA平均片段長(zhǎng)度無(wú)明顯變化的主要原因。

Fig.7 Statistical results of correlation between DNA damage and different sequencing library construction methods (A)5′ C→T frequencies of ancient DNA obtained via single-stranded and double-stranded library construction methods;(B)Average fragment lengths of ancient DNA obtained via single-stranded and double-stranded library construction methods.Values are the mean ± SD of multiple(n ≥ 9)determinations from separate experiments.One-way analysis of variance(ANOVA)is employed to determine the significant differences between single-and double-stranded libraries.P < 0.05 indicates significant difference between Single-and double-stranded libraries

本研究對(duì)中國(guó)東北地區(qū)的古脊椎動(dòng)物樣品進(jìn)行了DNA損傷分析，結(jié)果顯示，樣品中古DNA分子末端堿基替換率與埋藏環(huán)境的含水量、樣品埋藏時(shí)間呈正相關(guān)關(guān)系；不同的MIS時(shí)期之間，古DNA末端堿基替換率有顯著差異；樣品古DNA的平均片段大小與所研究的因素均無(wú)明顯關(guān)系。本文的研究為古DNA研究的采樣地點(diǎn)和樣品部位的選擇提供了一定參照。本研究中的樣品均來(lái)自中國(guó)東北地區(qū)，中國(guó)地域廣闊且化石資源豐富，區(qū)域之間氣候的差異較大，不同區(qū)域間古代樣品DNA損傷情況尚不明了。對(duì)中國(guó)各地區(qū)間、各地質(zhì)時(shí)期間的古代樣品DNA損傷情況的研究可以為樣品的選擇和優(yōu)先處理提供有價(jià)值的參考。