XIAO Bo SHENG Gui-Lian,2* YUAN Jun-Xia WANG Si-Ren HU Jia-Ming CHEN Shun-Gang JI Hai-Long HOU Xin-Dong LAI Xu-Long

(1 School of Environmental Studies, China University of Geosciences (Wuhan) Wuhan 430078 mr.shaw@outlook.com* Corresponding author: glsheng@cug.edu.cn)

(2 State Key Laboratory of Biogeology and Environmental Geology, China University of Geosciences (Wuhan) Wuhan 430078)

(3 Faculty of Materials Science and Chemistry, China University of Geosciences (Wuhan) Wuhan 430078)

(4 Daqing Museum Daqing, Heilongjiang 163319)

(5 Paleontological Fossil Conservation Center, Qinggang County Qinggang, Heilongjiang 151600)

(6 School of Earth Sciences, China University of Geosciences (Wuhan) Wuhan 430074)

Abstract The deer resources in China are abundant, with seven species in the sub-family Cervinae distributing in various areas. The intraspecific phylogeny of Cervinae has been widely explored, while the evolutionary relationship among different species requires further efforts, in which only few molecular studies on ancient materials have been performed. In this study, we carried out ancient DNA research on two Cervinae subfossils from northeastern China, dating of 3800 and 5100 aBP. Through ancient DNA extraction, double-stranded sequencing libraries construction, next-generation sequencing and bioinformatics data analysis, we reconstructed two mitochondria sequences with lengths of 16475 bp (GenBank accession number: MT784751,sequence integrity: 99.83%) and 16167 bp (GenBank accession numberh: MT784752, sequence integrity: 97.96%), respectively. Based on the mitochondrial homologous sequences of the extant Cervinae species in GenBank, we constructed a phylogenetic tree. The results show that: 1) both the average length and the C-to-T substitution frequencies at 5’- end of the NGS short reads indicate the data are from ancient specimens; 2) the two ancient individuals clustered with Cervus elaphus in the phylogenetic tree, and were molecularly identified as C. elaphus; 3) the two ancient samples from Heilongjiang are phylogenetically close to the extant C. elaphus alxaicus, but far from the extant C. elaphus xanthopygus. Combining the dates of the samples, we suggest that these two samples represent a population of ancient C. elaphus in Heilongjiang, which was not the direct maternal ancestor of the extant C. elaphus xanthopygus.

Key words Heilongjiang, China; Cervinae; ancient DNA; molecular identification; phylogenetic analysis

鹿亞科(Cervinae)又稱真鹿亞科，共分為4個(gè)屬：鹿屬(Cervus)、黇鹿屬(Dama)、花鹿屬(Axis)和麋鹿屬(Elaphurus) (Sheng, 1992)。我國(guó)鹿類資源較為豐富，鹿亞科動(dòng)物就包括鹿屬的梅花鹿(Cervus nippon)、馬鹿(C. elaphus)、坡鹿(C. eldi)、水鹿(C. unicolor)和白唇鹿(C. albirostris), 花鹿屬的豚鹿(Axis porcinus)和麋鹿屬的麋鹿(Elaphurus davidianus)共7個(gè)種(Tu et al., 2012)。中國(guó)現(xiàn)存的鹿種能表現(xiàn)出該亞科動(dòng)物幾乎全部系統(tǒng)發(fā)生的階段，因此中國(guó)也被稱為“鹿類動(dòng)物進(jìn)化的展覽場(chǎng)”(Sheng, 1992)。

過(guò)去，人們對(duì)中國(guó)鹿類動(dòng)物起源和系統(tǒng)進(jìn)化及分類研究主要依據(jù)形態(tài)學(xué)特征，例如角的有無(wú)和形態(tài)、體型和體重、頭骨和炮骨、淚窩和腺體的有無(wú)及特征、齒式和毛色特征、行為特征等，其中顱骨衍生物(牙齒、角等)是形態(tài)學(xué)系統(tǒng)分類的重要特征之一(Groves and Grubb, 1987; Sheng, 1992; Cai and Yin, 1992; Dong and Li, 2009; Dong et al.,2018)?；诨男螒B(tài)學(xué)研究是對(duì)古代鹿亞科動(dòng)物研究的方法之一，盛和林(1992)認(rèn)為鹿類動(dòng)物的化石記錄雖能較真實(shí)地反映其進(jìn)化歷史，但由于化石的分類基準(zhǔn)大都使用角或容易變異的骨的形態(tài)，因此容易造成異議。Geist (1998)指出，偶蹄類的比較形態(tài)學(xué)在分類和進(jìn)化研究中并不總是可靠的，因?yàn)檫@些動(dòng)物的形態(tài)學(xué)特征隨著環(huán)境條件變化而產(chǎn)生很大差異，所以結(jié)果不能被普遍接受。Emerson and Tate (1993)指出鹿類動(dòng)物形態(tài)學(xué)研究所依據(jù)的角形和皮毛顏色等特征容易受到趨同或平行進(jìn)化的影響。因此鹿科動(dòng)物的形態(tài)學(xué)分類仍然存在較大爭(zhēng)議。近40年來(lái)，很多學(xué)者在傳統(tǒng)分類基礎(chǔ)上引入了細(xì)胞學(xué)(Wang and Du, 1983; Neitzel, 1987)和分子生物學(xué)以及計(jì)算機(jī)技術(shù)對(duì)鹿類動(dòng)物進(jìn)行分類研究，從細(xì)節(jié)上對(duì)其進(jìn)行了補(bǔ)充和修改。然而，涉及豚鹿、麋鹿和馬鹿的系統(tǒng)分類地位，不同學(xué)者的觀點(diǎn)尚存爭(zhēng)議(Emerson and Tate, 1993; Grubb, 1993; Geist, 1998; Polziehn and Strobeck, 2002; Liu et al., 2003)。

古DNA是殘存在古生物化石或亞化石、考古材料中的生物大分子，能從分子水平反映生物個(gè)體之間的遺傳組成差異。發(fā)展利用古代鹿科動(dòng)物的古DNA線粒體甚至核基因組信息，并結(jié)合現(xiàn)生鹿科動(dòng)物的遺傳信息和動(dòng)物譜系地理學(xué)，綜合分析研究中國(guó)鹿科動(dòng)物系統(tǒng)發(fā)育是一個(gè)良好的發(fā)展方向(Liu et al., 2017)。本研究對(duì)采自黑龍江省肇東市和賓縣的兩個(gè)鹿科動(dòng)物遺骨進(jìn)行古DNA提取、DNA雙鏈文庫(kù)構(gòu)建和高通量測(cè)序，通過(guò)生物信息學(xué)方法對(duì)序列數(shù)據(jù)進(jìn)行處理，從分子水平對(duì)亞化石材料進(jìn)行物種鑒定，結(jié)合現(xiàn)生部分鹿科動(dòng)物基因序列，比較其同源序列之間的差異，通過(guò)系統(tǒng)發(fā)育分析探討馬鹿亞種間的進(jìn)化關(guān)系。

1 材料與方法

1.1 實(shí)驗(yàn)材料

兩個(gè)實(shí)驗(yàn)材料皆為黑龍江省大慶博物館館藏標(biāo)本。據(jù)博物館資料記載，CADG522于2008年發(fā)現(xiàn)于黑龍江省肇東市太平鄉(xiāng)段松花江古河道；CADG527于2003年發(fā)現(xiàn)于黑龍江省賓縣鳥河鄉(xiāng)的泥溝。樣品年代由美國(guó)BETA實(shí)驗(yàn)室通過(guò)放射碳同位素測(cè)定方法獲得，相關(guān)信息如表1、圖1所示。

表1 研究樣品信息Table 1 Samples information

圖1 產(chǎn)自黑龍江肇東太平鄉(xiāng)(A)和賓縣鳥河鄉(xiāng)(B)的鹿類研究樣品Fig. 1 Samples of Cervus sp. from Taiping, Zhaodong (A) and Niaohe, Binxian (B), Heilongjiang

1.2 古DNA提取

首先用4.5%的次氯酸溶液擦拭骨骼樣品表面，再用無(wú)水乙醇擦拭兩次并風(fēng)干。使用切割機(jī)切取約300 mg骨骼樣本，用已滅菌研缽研磨成粉末。按Rohland and Hofreiter(2007a, b)的方法，使用含EDTA (0.465 mol/L)和蛋白酶K (0.4 g/L)的裂解液充分裂解樣品，再使用DNA純化試劑盒(QIAquick PCR Purification Kit)結(jié)合超濾管(Millipore)完成樣品古DNA的提取。

1.3 古DNA雙鏈文庫(kù)構(gòu)建與測(cè)序

古DNA雙鏈文庫(kù)構(gòu)建包括4個(gè)步驟(Matthias and Martin, 2010): 1) 末端修復(fù)(bluntend repair), 使用T4多聚核苷酸激酶(T4 polynucleotide kinase)和T4 DNA聚合酶(T4 polymerase)將DNA分子片段末端修復(fù)平齊；2) 接頭連接(adapter ligation), 用快速連接酶(quick ligase)將接頭(adapter)片段連接至修復(fù)平齊的雙鏈DNA分子末端；3) 接頭補(bǔ)齊(adapter fill-in), 使用Bst DNA聚合酶將連接接頭的黏性末端補(bǔ)齊為平末端；4) 文庫(kù)分子標(biāo)記(indexing), 通過(guò)添加不同的小片段序列對(duì)DNA分子進(jìn)行特異性標(biāo)記，以在混合文庫(kù)測(cè)序時(shí)對(duì)單個(gè)文庫(kù)進(jìn)行區(qū)分，該過(guò)程主要使用的酶為Q5高保真DNA聚合酶(high-fidelity DNA polymerase)。最后對(duì)構(gòu)建的文庫(kù)進(jìn)行一次磁珠純化即可得到用于測(cè)序的古DNA雙鏈文庫(kù)。構(gòu)建完成的文庫(kù)送北京泛生子基因測(cè)序公司進(jìn)行雙向測(cè)序。本研究對(duì)兩個(gè)樣品分兩個(gè)批次各構(gòu)建兩個(gè)文庫(kù)(522-1, 522-2及527-1, 527-2), 對(duì)第一批次兩個(gè)文庫(kù)進(jìn)行預(yù)測(cè)序得到測(cè)序數(shù)據(jù)522-1Y和527-1Y, 深度測(cè)序得到測(cè)序數(shù)據(jù)522-1S和527-1S。第二批次文庫(kù)直接進(jìn)行深度測(cè)序得到測(cè)序數(shù)據(jù)522-2S和527-2S。

1.4 測(cè)序數(shù)據(jù)處理與分析

測(cè)序公司提供的測(cè)序序列文件為FASTQ格式，對(duì)每一樣品的每一文庫(kù)，分別使用“cutadapt-1.12” (Martin, 2011)切除接頭序列并過(guò)濾掉低于30 bp的短片段，然后用軟件“BWA-0.7.15” (Burrows-Wheeler Alignment) (Li and Durbin, 2010)中的“aln”和“samse”將片段與參考基因組比對(duì)，舍棄比對(duì)質(zhì)量得分低于30的序列；運(yùn)用SAMtools v1.3.1軟件(Li et al., 2009)中的“view”和“sort”算法將剩余序列按5’端位置在參考基因組上進(jìn)行排序，用“rmdup”去除可能存在的PCR重復(fù)序列，用“merge”算法得到每一文庫(kù)的單一bam格式文件。將同一樣品的不同文庫(kù)bam文件合并，再運(yùn)行“rmdup”去除可能存在的重復(fù)序列，得到同一樣品的合并bam文件，最終用angsd-0.916軟件(Korneliussen et al., 2014)“-doFasta”導(dǎo)出匹配得到的各樣品基因序列FASTA格式文件。

對(duì)古DNA分子堿基損傷的評(píng)估，使用各樣品不同文庫(kù)處理合并、并去掉重復(fù)后的bam文件作為輸入文件，使用“mapDamage”軟件包(Aurelien et al., 2001)統(tǒng)計(jì)片段末端5’端和3’端堿基損傷。

序列匹配時(shí)，先以麋鹿線粒體基因組序列(GenBank收錄號(hào)：NC018358)作為參考序列對(duì)測(cè)序數(shù)據(jù)進(jìn)行初步比對(duì)，將所得線粒體序列數(shù)據(jù)導(dǎo)出，使用BLAST (https://blast.ncbi.nlm.nih.gov)在線檢索，得到樣品序列的物種信息，初步判定實(shí)驗(yàn)樣品的種屬類型。本研究中對(duì)兩個(gè)樣品的第一批次得到的第一個(gè)文庫(kù)，先后使用麋鹿(GenBank登錄號(hào)：NC018358)和馬鹿阿拉善亞種(GenBank登錄號(hào)：KU942399)的線粒體基因組作為參考基因組對(duì)預(yù)測(cè)序數(shù)據(jù)(522-1Y和527-1Y)進(jìn)行匹配分析。得到第二批次的各兩個(gè)文庫(kù)的深度測(cè)序結(jié)果(522-1S, 522-2S和527-1S, 527-2S)后，只選用在第一個(gè)文庫(kù)匹配分析中匹配度較高的參考基因組用作合并后的參考基因組，對(duì)兩個(gè)樣品各文庫(kù)分別進(jìn)行數(shù)據(jù)過(guò)濾和合并處理，最終得到本研究的線粒體古DNA序列。

1.5 系統(tǒng)發(fā)育分析

利用Bioedit 7.0.9.0 (Hall, 1999)對(duì)測(cè)得的序列進(jìn)行排列對(duì)比并輔以人工校對(duì)，以所得兩個(gè)樣品的古DNA序列與36個(gè)鹿亞科動(dòng)物和一個(gè)外類群駝鹿(Alces alces)的同源序列(15109 bp, 不包含D-loop區(qū))在CIPRES網(wǎng)站(http://www.phylo.org)使用PartitionFinder2(Cognato and Vogler, 2001)和RAxML V8.0 (Alexandros, 2014)通過(guò)500次自舉檢驗(yàn)估計(jì)節(jié)點(diǎn)構(gòu)建系統(tǒng)發(fā)育樹，堿基替換模型(model for substitution)選擇GTRCAT。系統(tǒng)發(fā)育樹的展示輸出通過(guò)Figtree v1.4.2 (Rambaut, 2009)完成。分析中所用序列信息見(jiàn)表2, 序列選取兼顧鹿亞科在中國(guó)分布的4屬9種。同時(shí)，與樣品聚類較近的中國(guó)馬鹿5亞種和梅花鹿選取了當(dāng)前GenBank中所收錄的所有高質(zhì)量序列。

表2 系統(tǒng)發(fā)育分析所使用序列信息Table 2 Sequence information used in phylogenetic analysis

2 結(jié)果

2.1 古代線粒體DNA片段末端堿基特性

運(yùn)用mapDamage軟件對(duì)兩個(gè)樣品所得線粒體DNA二代測(cè)序序列進(jìn)行分析，末端堿基替換統(tǒng)計(jì)結(jié)果如圖2所示。CADG522和CADG527在5’端的“C→T”堿基替換頻率分別為6.883%和6.897%, 3’端的“G→A”堿基替換頻率分別為6.105%和9.123%。

圖2 樣品CADG522 (A)和CADG527 (B)線粒體DNA片段末端核苷酸錯(cuò)配統(tǒng)計(jì)Fig. 2 Statistical analysis of nucleotide misincorporation observed on reads of samples CADG522 and CADG527

2.2 二代測(cè)序數(shù)據(jù)及序列比對(duì)結(jié)果

先后使用麋鹿(GenBank收錄號(hào)：NC018358, 全長(zhǎng)16355 bp)和馬鹿阿拉善亞種(GenBank收錄號(hào)：KU942399, 全長(zhǎng)16503 bp)的線粒體基因組作為參考基因組，對(duì)兩個(gè)樣品第一批次的文庫(kù)(522-1和527-1)的預(yù)測(cè)序數(shù)據(jù)(522-1Y和527-1Y)進(jìn)行匹配分析，結(jié)果顯示兩個(gè)樣品均與阿拉善馬鹿匹配度較高。最終以阿拉善馬鹿線粒體基因組作為參考基因組對(duì)文庫(kù)深度測(cè)序數(shù)據(jù)(522-1S, 522-2S和527-1S, 527-2S)進(jìn)行過(guò)濾、匹配和合并處理，分別得到兩個(gè)樣品16475 bp (GenBank收錄號(hào)：MT784751, 序列完整度：99.83%)和16167 bp (GenBank收錄號(hào)：MT784752, 序列完整度：97.96%)古DNA線粒體基因序列。文庫(kù)測(cè)序信息如表3所示。

表3 二代測(cè)序數(shù)據(jù)信息Table 3 Next-generation sequencing data

2.3 系統(tǒng)發(fā)育樹

基于本研究所得兩條古DNA線粒體序列和從GenBank檢索所得36條鹿亞科動(dòng)物的線粒體基因組序列，以駝鹿(A. alces)為外類群(outgroup)構(gòu)建的分子系統(tǒng)發(fā)育樹(圖3)顯示，兩個(gè)古代樣品(CADG522和CADG527)都與馬鹿阿拉善亞種聚為一支，馬鹿阿拉善亞種、天山亞種、甘肅亞種和東北亞種一起與梅花鹿互為姐妹群，且與后者一起共同構(gòu)成馬鹿塔里木亞種的姐妹群。白唇鹿兩個(gè)個(gè)體聚為一支處于水鹿大分枝內(nèi)部。坡鹿與3個(gè)麋鹿聚為同一分枝，往樹根依次為黇鹿屬的黇鹿和花鹿屬的花鹿和豚鹿。

圖3 鹿亞科線粒體基因組15109 bp同源序列ML法系統(tǒng)發(fā)育樹Fig. 3 ML phylogenetic tree of Cervinae species based on 15109 bp mitochondrial genome homologous sequences

3 討論

3.1 古代樣品線粒體DNA真實(shí)性鑒定

本研究所得古DNA序列真實(shí)性可從如下方面論證：首先，采集樣品時(shí)嚴(yán)格按照古DNA研究規(guī)范，樣品的處理、提取及建庫(kù)都在古DNA專用實(shí)驗(yàn)室的獨(dú)立環(huán)境進(jìn)行，實(shí)驗(yàn)者嚴(yán)格遵守古DNA實(shí)驗(yàn)防污染規(guī)則操作。其次，實(shí)驗(yàn)過(guò)程中，提取、建庫(kù)全程設(shè)置有空白對(duì)照及正對(duì)照，對(duì)照實(shí)驗(yàn)組顯示無(wú)異常，排除試劑、操作者及樣品之間的交叉污染。再次，本研究所得序列經(jīng)BLAST查詢，所得推薦高可信度序列均為鹿科動(dòng)物線粒體基因序列；本實(shí)驗(yàn)室從未進(jìn)行現(xiàn)代鹿科動(dòng)物分子生物學(xué)實(shí)驗(yàn)，所得序列確定來(lái)源于古代樣品。最后，本研究各文庫(kù)所得線粒體DNA片段平均長(zhǎng)度較短，約45～75 bp, 符合古DNA片段特征(Hofreiter et al., 2014)。同時(shí)，經(jīng)二代測(cè)序所得兩個(gè)樣品的DNA片段的5’端堿基“C→T”替換頻率分別為6.883%和6.897%, 符合Sawyer et al. (2012)關(guān)于古代線粒體DNA片段5’端堿基“C→T”替換頻率與樣品年代的相關(guān)性研究結(jié)果。綜上，本研究所得序列是來(lái)自古代樣品的古DNA序列。

3.2 鹿亞科動(dòng)物的系統(tǒng)演化關(guān)系

近年來(lái)，關(guān)于鹿亞科各屬種的系統(tǒng)演化關(guān)系，前人的蛋白質(zhì)電泳、染色體組型、線粒體控制區(qū)序列等研究表明麋鹿和豚鹿與鹿屬動(dòng)物的進(jìn)化關(guān)系較近，認(rèn)為麋鹿和豚鹿應(yīng)該并入鹿屬，歸并后的鹿屬為單系發(fā)生(Wang and Du, 1983; Randi et al., 2001; Liu et al., 2003)。本研究加入古代鹿科動(dòng)物樣品序列后構(gòu)建的系統(tǒng)發(fā)育樹(圖3)顯示，麋鹿與鹿屬的坡鹿互為姐妹支，且自展支持率為100%, 支持Randi et al. (2001)將麋鹿歸為鹿屬的觀點(diǎn)；但是豚鹿在系統(tǒng)發(fā)育樹中的位置與鹿屬各種距離較遠(yuǎn)，且豚鹿與花鹿歸為同一分支，因此，本研究不支持王宗仁、杜若甫(1983)和劉向華等(2003)將豚鹿劃歸鹿屬的觀點(diǎn)。

在鹿屬內(nèi)部，王宗仁、杜若甫(1983)認(rèn)為梅花鹿與馬鹿的關(guān)系最近，現(xiàn)生鹿屬動(dòng)物中水鹿是最原始的，白唇鹿與梅花鹿幾乎同時(shí)從水鹿的祖先種分化而來(lái)，馬鹿則是最后分化出來(lái)的。圖3中，梅花鹿并非直接從水鹿分化而來(lái)，梅花鹿分支與馬鹿阿拉善亞種、天山亞種、甘肅亞種和東北亞種共同聚類形成的分支構(gòu)成姊妹群，而馬鹿塔里木亞種處于二者的根部。該結(jié)果顯示馬鹿塔里木亞種在早期從馬鹿和梅花鹿的共同祖先分化而出，而后依次分化出梅花鹿和馬鹿其他亞種，從線粒體基因組水平印證了王宗仁等的結(jié)論。

3.3 馬鹿種群的基因交流

關(guān)于馬鹿的系統(tǒng)進(jìn)化關(guān)系，大多數(shù)觀點(diǎn)認(rèn)為馬鹿分為東西兩個(gè)進(jìn)化種群：西部群體包括歐洲種群和我國(guó)新疆南部塔里木亞種，東部種群包括北美種群、亞洲及我國(guó)新疆北部的種群(Polziehn and Strobeck, 2002; Muhmut et al., 2002; Ludt et al., 2004)。本研究的系統(tǒng)進(jìn)化樹顯示塔里木馬鹿與其他生活在中國(guó)的4個(gè)馬鹿亞種來(lái)源于共同祖先的兩個(gè)不同分支，支持馬鹿在進(jìn)化中分出東西兩個(gè)種群的觀點(diǎn)。塔里木馬鹿較早地分化出來(lái)，在一個(gè)封閉的地理環(huán)境下獨(dú)立進(jìn)化，因地理隔離而有別于其他馬鹿種群。之后，東北亞種、甘肅亞種、天山亞種和阿拉善亞種相繼分化，這也證明了中國(guó)馬鹿是從中東和歐洲返回中國(guó)的過(guò)程中從西向東逐漸分化的(Sheng, 1992)。

現(xiàn)分布于寧夏和內(nèi)蒙古交界的賀蘭山中段的阿拉善馬鹿種群，是我國(guó)唯一幸存的馬鹿阿拉善亞種的有效種群(Li et al., 1998; Wang et al., 1999)。受賀蘭山周圍被城市、沙漠、河流隔斷形成孤島的地理特征所影響，該種群已成為一個(gè)隔離種(Qiao et al.,2019)。本研究涉及的兩個(gè)古代鹿科動(dòng)物樣品，校正年代分別為(3800±30)和(5100±30)aBP, 采集地黑龍江肇東市和賓縣與賀蘭山相距約1900 km。黑龍江地區(qū)的古代樣品在系統(tǒng)發(fā)育樹上與賀蘭山的現(xiàn)生馬鹿阿拉善亞種聚為一支，而與現(xiàn)生馬鹿東北亞種親緣關(guān)系較遠(yuǎn)，表明在歷史時(shí)期兩地的馬鹿種群存在較密切親緣關(guān)系，可能發(fā)生過(guò)種群消退、遷徙、引入等，也可能曾有過(guò)高于現(xiàn)生種群的基因交流；兩例古代樣品代表的古代黑龍江馬鹿種群，不是現(xiàn)生馬鹿東北亞種的母系直系祖先。本研究后續(xù)工作將對(duì)不同歷史時(shí)期、不同地點(diǎn)的馬鹿古代樣品的線粒體基因組及核基因組進(jìn)行研究，以期為馬鹿各亞種的種群演化歷史提供更為全面的分子證據(jù)。