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

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

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

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

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

1 材料與方法

1.1 樣品地點及年代信息

本研究用到的二代測序數(shù)據(jù)來源于中國東北第四紀晚期以來9種哺乳動物共計33件化石和亞化石樣品，其中15個樣品序列來自本研究組已發(fā)表文章[25-28]，18個樣品序列尚未發(fā)表。33個樣品的地點信息正如Fig.1所示，包括骨骼化石或亞化石25個、牙齒化石7個、角化石1個。所有樣品年代信息由美國BETA實驗室(Beta Analytic, USA)通過放射性碳測年(radiocarbon dating)方法獲得，具體信息正如Table 1所示。古DNA提取及建庫方法參照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 原始測序數(shù)據(jù)的生物信息學(xué)處理

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

1.3 序列損傷檢測

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

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

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

2 結(jié)果

2.1 堿基損傷

通過“mapDamage”軟件包得到各樣品中古DNA片段的堿基損傷結(jié)果，隨機抽取2個可視化結(jié)果正如Fig.2所示。樣品5′端C→T的替換率以及DNA的平均片段長度結(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 末端堿基替換率與片段大小無明顯相關(guān)性

對各樣品古DNA片段進行末端堿基替換率與平均片段長度分布情況統(tǒng)計，結(jié)果正如Fig.3所示。平均片段長度絕大多數(shù)(97%)處于45～75 bp之間(Fig.3A)。所有樣品堿基錯配概率在0.05～0.35區(qū)間內(nèi)，并且91%的樣品(30個)錯配率低于0.2，錯配率超過0.2的樣品僅占9%(3個, Fig.3B)，這3個樣品(CADG449/DLH2/TH3)都出自黑龍江哈爾濱市。此外，片段長度在45～55、55～65、65～75 bp三個區(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末端堿基替換率隨埋藏時間增加升高

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

2.4 MIS1時期古DNA末端堿基替換率最低

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

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

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

2.6 不同建庫方法對古DNA損傷無顯著影響

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

3 討論

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

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

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

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

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

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 樣品類型及不同建庫方法對古DNA損傷的影響

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

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

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

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