基于簡化基因組測序的麋鹿遺傳資源多樣性及種群結構分析

張樹苗，田恒玖，李夷平，陳 頎，李俊芳，李 林，郭青云，王立波，程志斌，白加德

（1.北京麋鹿生態實驗中心，北京，100076；2.北京松山國家級自然保護區，北京，102115；3.江蘇省大豐麋鹿國家級自然保護區管理處，鹽城，224136）

麋鹿（Elaphurus davidianus）是中國特有物種，國家一級重點保護野生動物，IUCN 紅皮書野外滅絕物種［1-2］，是世界物種重引進項目的成功范例。19世紀末，由于棲息地喪失、氣候變化及過度捕獵，麋鹿在中國本土滅絕［3］。1900 年前后，中國僅存的麋鹿被歐洲各大動物園掠奪飼養。英國烏邦寺收集了分布于歐洲各地的18 只麋鹿進行繁殖，逐步建立起完整的小種群［4］。1985—1987 年，我國政府著手開展麋鹿重引入工作［5］，同期，英國烏邦寺分兩批將38只麋鹿重引入到中國，飼養在北京南海子麋鹿苑；1986年引入39 只到江蘇大豐麋鹿國家級自然保護區。經過30 多年的人工繁殖、復壯及野外放歸，中國麋鹿保護工作取得顯著進展。目前人工繁殖的麋鹿種群已全面覆蓋麋鹿原有棲息地，遷地保護地從當初的2 個增至現在的89 個，種群數量12 000 余只，其中有6個野生種群，數量達5 000余只［4，6］。中國麋鹿種群從本土滅絕經重引入，種群逐漸復壯，再到放歸自然，成為保護野外滅絕物種的成功典范。

小種群易受遺傳漂變的影響［7］，特別是重引入的麋鹿種群具有顯著的奠基者效應和瓶頸效應（bottleneck effect）［8-9］。奠基者效應對麋鹿重引入的擴張具有顯著影響［10-11］，而瓶頸效應發生在種群經歷短暫的數量劇減后，種群的遺傳信息完全來自少數存活下來的個體［12］。因此，瓶頸效應與奠基者效應非常類似，這種效應的嚴重后果就是破壞種群遺傳結構，使原有種群動態特征失去平衡，加劇遺傳多樣性的喪失。正因如此，我國麋鹿種群存在遺傳多樣性缺失，近交衰退嚴重，抵抗疾病能力差，壽命降低及雌鹿流產率高等現象［13-14］。人類活動、氣候變化、環境污染和生境破碎化等會導致野生麋鹿種群數量進一步下降，形成離散種群［15-18］。目前，我國麋鹿保護仍面臨嚴峻挑戰［19-20］。

近年來，隨著高通量測序技術的快速發展，通過基因組重測序手段，在對物種鑒定、分類、遺傳多樣性和遺傳結構的研究中發揮越來越重要的作用［21-26］。本研究基于簡化基因組技術，對北京南海子麋鹿苑半散養種群的10 只麋鹿及江蘇大豐麋鹿國家級自然保護區圈養種群和野生種群的8 只麋鹿，共18 只個體進行遺傳多樣性和遺傳結構分析，在不同生境下，分析中國麋鹿的遺傳多樣性背景，并對中國重引入麋鹿工作作出科學而客觀的評價，以期為麋鹿放歸、種群基因交流等保護措施的開展提供科學支撐。

1 研究區概況

北京南海子麋鹿苑（以下簡稱“北京麋鹿苑”）位于華北平原北部，燕山山脈南緣。江蘇大豐麋鹿國家級自然保護區（以下簡稱“江蘇大豐”）地處黃海之濱，江蘇省東部、鹽城市南部，東臨黃海，具體信息見表1。

表1 研究區基本信息Tab.1 The basic information of study areas

2 材料與方法

2.1 材料與基因組樣本收集

共采集18 只麋鹿個體的血液樣品，其中，10 只采集于北京麋鹿苑半散養種群，8只采集于江蘇大豐圈養（2 區：4 只）和野生種群（3 區：4 只）。將血液樣品各取5 mL 至抗凝管中，通過天根試劑盒（QIAamp Fast DNA Stool Mini Kit（50））提取血液樣本DNA，并通過瓊脂糖凝膠和NanoDrop 2000 紫外分光光度計（Thermo Fisher Scientific）對抽提的DNA 質量和濃度進行檢測。DNA產物置于-20 ℃保存備用。

使用BGISEQ-500測序儀測序，隨機打斷檢測合格的DNA 樣品，篩選符合要求、大小適合的DNA 片段，并將純化后的DNA 片段連接測序接頭，滾環擴增制備DNB，然后在陣列化測序芯片上進行雙端測序。

2.2 測序及參考基因組序列比對

原始下機數據（raw data）經過SOAPnuke過濾后，應用BWA 軟件將得到的clean data 比對到參考基因組上（NCBI 數據庫的Project ID：PRJNA391565［25］），然后使用SAMtools 排序，參數為“-q1-C50-g-tDP，SP，DP4-I-d250-L250-m2-p”，并利用Picard 軟件對重復reads 標記，最后使用GATK（https：//www.broadinstitute.org/gatk/）進行BQSR質量值校正，并移除低質量測序數據。

2.3 基因型鑒定

SNP檢測：使用Genome Analysis Toolkit（GATK）［27］對比對得到的BAM 文件進行單核苷酸多態（single nucleotide polymorphisms，SNPs）檢測和過濾。為獲得高質量SNP，選用參數QD＜2，MQ＜30，MQRankSum＜-12.5，FS＞200，ReadPosRankSum＜20，QUAL＜30.0和AN＜40 對低質量的SNP 進行過濾。過濾后，利用HaplotypeCaller 軟件，以GVCF 文件生成并統計每只個體形成的SNP基因型。為進一步精準比較北京麋鹿苑和江蘇大豐麋鹿SNP 位點的多態性，使用VCFtools 對種群間未覆蓋位點、個體間多態性低位點及重疊缺失位點進行過濾，參數為MAF＞0.05。

2.4 遺傳多樣性分析

利用R 軟件包統計SNP 位點數目和觀測雜合度（Ho）。

2.5 遺傳結構和分化

使用完整的SNP 數據集，利用VCFtools 計算2 個種群間的群體遺傳分化指數（FST），分析群體間遺傳關系遠近，并計算Tajima’sD值，檢測選擇效應。

使用fastStructure［28］對種群祖先遺傳成分進行分析，首先利用structure Harvester 做K值評估，設置K=2～8，迭代數設為100，每個K值運行20 次，統計K值的最大似然值。Structure 計算結果用Structure Harvester網頁工具分析并進行ΔK計算，公式為：ΔK=mean（|L″（K）|）/sd（L（K）），根據結果繪制ΔK折線圖，找出最佳K值。

使用Plink軟件對2個群體的SNP基因型進行主成分（PCA）分析，利用R 軟件繪制第一主成分與第二主成分的二維圖。

使用FastTree［29］中的鄰接法（neighbor-joining，NJ）進行1 000 次Bootstrap 檢驗，構建系統發育樹（phylogenetic tree），用以研究個體之間的遺傳聚類分析。

3 結果

3.1 測序數據質量及與參考基因組比對

18 只麋鹿個體簡化基因組測序reads 數量為448.98 Mb～711.08 Mb。過濾后reads有效數據比例為98.91%～99.78%。Raw bases為44.90 Gb～71.11Gb，Clean bases 為44.50 Gb～70.94 Gb。CleanQ30bases rate 為92.86%～95.14%，說明樣本建庫測序的質量達到重測序標準，測序數據質量信息見表2。

表2 測序數據質量信息匯總Tab.2 Summary of sequencing data quality information

本研究采用的參考基因組為麋鹿參考基因組，大小約為2.6 Gb（https：//www.ncbi.nlm.nih.gov/assembly/GCA_002443075.1/），包 含194 812 條scaffold，GC 含量為40.5%［25］。18 份麋鹿樣品簡化基因組測序數據與參考基因組比對情況見表3。樣本比對到參考基因組上的比例為98.18%～99.67%，比對到參考基因組唯一位置的比例為94.03%～94.68%，平均測序深度為15.70～24.37 倍，覆蓋率為99.28%～99.54%，說明麋鹿簡化基因組測序基本覆蓋參考基因組，可用于后續分型和群體結構分析。

表3 麋鹿樣品簡化基因組測序數據與參考基因組比對Tab.3 Comparison of simplified genome sequencing data and reference genome of Père David’s deer

3.2 種群遺傳多樣性

經與參考基因組比較，樣本共檢測到SNP 位點37 581 752 個，遺傳多樣性統計見表4。北京麋鹿苑10 只麋鹿的SNP 平均值為2 166 903 個，江蘇大豐8 只麋鹿的SNP 平均值為1 989 090 個。兩地麋鹿SNP 位點的觀測雜合度分別為0.59～0.68 和0.63～0.71。總體來說，北京麋鹿苑的麋鹿遺傳多樣性略低于江蘇大豐。整合北京地區和江蘇地區的兩個麋鹿種群的SNP 標記，最終篩選到均勻分布于每條染色體、多態性好的共有SNP 位點1 456 457 個，用于后續遺傳結構分析。

表4 麋鹿樣本的SNP遺傳多樣性統計Tab.4 Statistics of SNP genetic diversity of Père David’s deer

3.3 種群遺傳結構和分化

使用VCFtools 分析種群間FST，結果顯示江蘇大豐和北京麋鹿苑2 個種群間平均FST為0.085，表明2 個種群間存在中等程度的遺傳分化。18 只麋鹿平均Tajima’sD為1.29，北京麋鹿苑麋鹿種群的平均Tajima’sD為0.90，江蘇大豐麋鹿種群的平均Tajima’sD為0.89，各種群Tajima’sD分布如圖1所示。

圖1 不同種群Tajima’s D分布Fig.1 Distribution of Tajima’s D values in different populations

種群結構分析顯示，當K=2 時，為最優K值，推測北京麋鹿苑和江蘇大豐麋鹿種群有了初步分化（圖2），并且北京麋鹿苑中，6 個個體顯示出純合性，表明相對于江蘇大豐麋鹿種群來說，北京麋鹿苑個體的血統可能更為單一。

圖2 麋鹿種群遺傳結構分析Fig.2 Population structure analysis of Père David’s deer

基于SNP 標記的主成分分析結果顯示，北京麋鹿苑和江蘇大豐2 個種群之間具有一定的遺傳分化，推測2 個種群地理引入個體可能存在一定的遺傳差異（圖3）。系統發育樹結構也支持主成分分析結果，2 個種群分別形成單獨的遺傳聚類（圖4）。總體來說，主成分分析與進化樹聚類結果均表明北京麋鹿苑和江蘇大豐2 個種群的遺傳多樣性已出現一定遺傳的差異。

圖3 麋鹿種群的主成分分析Fig.3 Principal component analysis diagram of Père David’s deer

圖4 麋鹿種群系統發育樹Fig.4 Phylogenetic tree NJ-tree of Père David’s deer

4 分析與討論

種群結構、群體分層及基因滲入是全面理解種群結構變異的重要指標。本研究從群體遺傳結構、群體主成分分析和系統發育樹，全面了解北京麋鹿苑和江蘇大豐麋鹿種群間的分類關系及親緣關系遠近，顯示出北京麋鹿苑半散養麋鹿與江蘇大豐麋鹿種群存在中等程度的遺傳分化。在群集最佳值為2（K=2）時，麋鹿種群分化較為顯著。基于成對的SNP 差異NJ系統發育樹也揭示2個重引入種群之間的獨立遺傳簇，這種遺傳分化暗示了北京麋鹿苑與江蘇大豐重引入群體之間存在一定的遺傳差異。

使用遺傳數據推斷種群歷史動態在群體遺傳學研究中具有重要作用。種群中存在許多低頻率的等位基因（稀有等位基因），表明群體存在定向選擇或種群擴張，而種群中高等/中等頻率的等位基因較多，則表明群體存在一定的平衡選擇和瓶頸效應［30］。FST和Tajima’sD是衡量種群結構差異和動態變化的關鍵參數，FST適用于亞種群間的多樣性比較，用于衡量種群分化程度，分化指數越大，差異越大。本研究結果顯示江蘇大豐和北京麋鹿苑2 個種群間平均FST為0.085，表明2個種群間存在中等程度的遺傳分化。同時，本研究分析得出北京麋鹿苑與江蘇大豐麋鹿種群的Tajima’sD＞0，表明2個種群在種群重引進和遷地保護中可能都存在一定的瓶頸效應。

麋鹿群體盡管經歷過嚴重的種群瓶頸，但仍存留較高的遺傳多樣性。本研究發現北京麋鹿苑與江蘇大豐麋鹿種群的觀測雜合度均大于0.50，表明2個種群的遺傳多樣性較豐富。目前，關于野生動物群體遺傳多樣性的研究已經進入基因組水平。Zhang等［25］通過第二代高通量基因組測序技術，構建了麋鹿第一個基因組草圖。針對麋鹿［26］、朱鹮（Nipponia nippon）［31］、大熊貓（Ailuropoda melanoleuca）［32］和 北極熊（Ursus maritimus）［33］，研究人員利用種群基因組測序和比較基因組進行遺傳多樣性的比較分析，結果顯示朱鹮的近交程度最高，麋鹿次之，并認為麋鹿并沒有存在很高的近交程度，這與多年來麋鹿群體的幼體存活率較高的結果［34］一致。

相比于單獨繁育一個大群，多個小群分開繁育所面臨的風險要大得多。中國在復壯麋鹿繁育種群時，不斷對麋鹿實施多地遷地保護，根據長期的棲息地調查和麋鹿種群數量普查工作，目前已確認建立了89 個麋鹿繁育場所，為保護麋鹿儲備了豐富的種源，但各種群間由于地理隔離而缺乏基因交流。本研究結果與這種現象一致，通過基因組測序技術，筆者發現北京麋鹿苑和江蘇大豐種群之間產生了遺傳分化，缺乏基因交流。因此，為更好地保護中國麋鹿種群的健康繁衍，加強各繁育場所的基因交流非常重要，應該增進不同種群間的基因交流。

本研究通過簡化基因組測序技術，獲得海量的第三代分子標記，利用高質量SNP 標記對18 只麋鹿的遺傳資源進行遺傳多樣性和種群結構分析，發現北京麋鹿苑麋鹿的遺傳多樣性略低于江蘇大豐。遺傳結構與進化樹結果均提示北京麋鹿苑與江蘇大豐的麋鹿種群有初步分化，該分化可能由地理隔離因素、遺傳漂變和瓶頸效應共同作用導致。由于本研究選用的麋鹿個體數量較少，SNP 標記組對麋鹿不同分布區的鑒定有效性還需進一步驗證。