應用簡化基因組技術對富民枳遺傳多樣性檢測1)

張珊珊 陳劍 楊文忠

(國家林業(yè)局云南珍稀瀕特森林植物保護和繁育重點實驗室(云南省林業(yè)科學院)，昆明，650201)

富民枳是典型的極小種群野生植物[1-2]。目前，新建立了1個富民枳野外回歸示范基地。但是，富民枳野生植株的遺傳結(jié)構(gòu)和遺傳多樣性如何？人工建立的回歸種群是否有效？為了有效保護其遺傳資源應采取怎樣的保護策略？這些問題都亟待解決。因此，如能采取科學方法對富民枳遺傳多樣性進行分析，并對回歸種群進行有效性和合理性的綜合檢測，進而提出能有效保護其遺傳資源的保護策略，這不僅能優(yōu)化富民枳的保護管理方案從而保存富民枳這一物種，而且能提升富民枳保護的科學性和完整性，并能通過這一物種的示范效應，輻射帶動云南乃至全國其他極小種群植物的科學保護。

群體遺傳學是研究群體的遺傳結(jié)構(gòu)及其變化規(guī)律的遺傳學分支學科[3]。在群體遺傳學的研究中，遺傳多樣性的高低是一個重要的指標，反映了物種應對棲息地變化的適應性和改變能力，而遺傳多樣性在群體間的分布，即遺傳結(jié)構(gòu)是受許多因素作用的，另外還與物種自身的生物學特征和進化史有著密切關系[4-5]。因此，我們在保護遺傳資源時需要先了解居群的遺傳結(jié)構(gòu)和遺傳背景，進而再制定合理的保護政策[6-8]。

利用分子生物學技術分析遺傳多樣性，可為有效群體的收集、維持和管理提供支持[9-10]。目前最先進的第二代DNA測序技術可通過產(chǎn)生大量的分子標記數(shù)據(jù)進行成株的群體遺傳多樣性分析，建立詳細的遺傳檔案[11]?；诿盖械暮喕蚪M測序(RAD-Seq，Restriction-site Associated DNA Sequence)是對與限制性核酸內(nèi)切酶識別位點相關的DNA進行高通量測序，可大幅降低基因組的復雜度，降低建庫和測序成本，操作簡便，同時不受參考基因組的限制，可快速鑒定出高密度的SNP位點，實現(xiàn)遺傳進化分析及重要性狀候選基因的預測[11-13]。RAD-Seq可以檢測基因組上未知變異點中新的SNP，發(fā)掘新的和稀有的變異，對解決群體遺傳學[14]、遺傳圖譜構(gòu)建[15]、功能基因挖掘[16]、群體進化[17]等問題，具有重大的科研和產(chǎn)業(yè)價值。

因此，本研究采用基于二代測序技術的RAD-seq方法對富民枳開展群體遺傳學研究，為每個成株建立詳細的遺傳檔案，用SNP位點刻畫其遺傳特征，完成遺傳多樣性分析，進而為回歸種群的遺傳保護提供理論基礎。

1 材料與方法

對富民枳現(xiàn)有的天然種群和野外回歸的人工種群分別取樣15株(編號A～編號O)和30株(R1～R30)，其中，建立野外回歸種群的苗木來源為天然種群的實生苗和扦插苗，但是具體來源不清楚；從回歸種群采集的30個樣品為隨機取樣。

采集樣品均為健康幼嫩的葉片，用硅膠快速干燥并拍照，帶回實驗室，用質(zhì)量分數(shù)75%的無水乙醇擦拭清洗，最后用液氮冷凍貯存以備提取DNA。

1.1 DNA提取及質(zhì)量檢測

采用改良CTAB法從富民枳成年植株和人工培育的幼苗葉片中提取富民枳全基因組DNA。由于簡化基因組測序?qū)NA質(zhì)量要求較高，需要對DNA純度、濃度和完整性進行檢測，以滿足建庫要求，具體操作和要求：用Nano Drop 2000C(Thermo Scientific，USA)檢測所提樣品的DNA質(zhì)量,以保證提取的DNA光吸收值A260/280介于1.7～2.1，OD260/230>2，DNA為無色膠狀物，表明DNA無蛋白、色素等雜志污染，比較純凈；使用Qubit對DNA濃度進行檢測，DNA質(zhì)量濃度>30 μg·L-1、DNA總量大于等于1 μg，即滿足建庫要求；用1%瓊脂糖電泳對DNA完整性進行檢測，DNA主帶明亮，無拖尾表示符合要求。

純度、濃度和完整性均符合要求的基因組DNA置于冰箱中-30 ℃暫時保存，直至送到公司進行建庫。

1.2 文庫構(gòu)建

RAD文庫構(gòu)建及測序委托上海美吉生物醫(yī)藥科技有限公司承擔完成。具體建庫步驟如下：(1)稀釋基因組DNA樣品；(2)根據(jù)選定的酶切方案，對檢測合格的各樣品基因組DNA分別配制酶切反應體系，在恒溫混勻儀上適溫反應一定時間酶切基因組DNA；(3)配制接頭1(300～500 bp)連接反應體系，在恒溫混勻儀中適溫反應一定時間使接頭1與酶切產(chǎn)物連接；(4)將接頭1連接產(chǎn)物按照適當?shù)谋壤M行混樣；(5)配制打斷體系，使用打斷儀對混樣產(chǎn)物進行打斷，再使用試劑盒對打斷產(chǎn)物進行純化；(6)用瓊脂凝膠電泳對產(chǎn)物進行片段選擇；(7)配制末端修復反應體系，在恒溫混勻儀中適溫反應一定時間對片段選擇后的產(chǎn)物進行末端修復反應，反應后用試劑盒進行純化；(8)配制反應體系，在恒溫混勻儀中適溫反應一定時間對片段選擇后的產(chǎn)物進行反應，反應后用試劑盒進行純化；(9)配制接頭2連接反應體系，在恒溫混勻儀中適溫反應一定時間使接頭2與加“A”后產(chǎn)物連接，再使用試劑盒對連接產(chǎn)物進行純化；(10)配制PCR反應體系對上一步的連接產(chǎn)物進行擴增，PCR產(chǎn)物用瓊脂糖凝膠回收純化，并溶解于適量EB solution中，至此文庫制備完成。構(gòu)建好的文庫進行質(zhì)量和產(chǎn)量的檢測，合格后用Illumina HiSeqTM平臺進行測序，測序策略為Illumina PE150，預計測序量為3 Gb，約為富民枳近源種基因組大小的3倍，周期大約為6～8周。

1.3 高通量測序分析

原始數(shù)據(jù)質(zhì)控和過濾：Illumina HiseqTM平臺測序得到的原始圖像數(shù)據(jù)經(jīng)過原始數(shù)據(jù)轉(zhuǎn)化為序列數(shù)據(jù)，得到最原始的測序數(shù)據(jù)文件。在Illumina HiseqTM測序數(shù)據(jù)下機之后，利用Fastp軟件對下機數(shù)據(jù)進行質(zhì)量控制，過濾其中低質(zhì)量的數(shù)據(jù)，獲得高質(zhì)量的質(zhì)控數(shù)據(jù)。

SNP檢測與位點開發(fā)：利用Stacks v1.42軟件進行變異檢測，獲得高質(zhì)量的SNP位點，然后分別進行個體聚類和群體聚類，利用Stacks v1.42軟件鑒別SNPs位點，再對每個個體的SNP進行統(tǒng)計分析，篩選出每個個體特有的SNP組合。Stacks v1.42軟件的分析流程如下。

原始數(shù)據(jù)處理：此步驟運行STACKS軟件數(shù)據(jù)進行質(zhì)量過濾，去除含有模糊堿基(N)和低質(zhì)量堿基的序列，并將序列截成相同的長度，形成質(zhì)控數(shù)據(jù)。在此步驟中，為了保證堿基序列的質(zhì)量，將所有片段截取至140 bp。c參數(shù)用于去除含有模糊堿基的任何堿基序列，產(chǎn)生質(zhì)控數(shù)據(jù)；q參數(shù)用于去除低質(zhì)量值的片段。

個體分析：此步驟運行STACKS中的個體分析程序，分別對每個個體進行分析，生成一系列的位點。利用SNP模型最大似然法對每個位點進行SNP發(fā)掘，生成的數(shù)據(jù)用于群體遺傳學分析，本研究使用SNP模型。在此步驟中，使用d參數(shù)，實行去杠桿化算法，用來解決過多產(chǎn)生的標簽；使用r參數(shù)，實行去除算法，用于去除算法中過度重復的片段和鄰近的錯誤；并設置創(chuàng)建一個片段所需的最小覆蓋深度(m參數(shù))為5；每個片段之間允許堿基間的最大距離(M參數(shù))為2。

數(shù)據(jù)整合：此步驟運行STACKS中的數(shù)據(jù)整合程序，將所有個體的所有位點數(shù)據(jù)進行整合，形成一個目錄，并形成一組一致性位點。在此步驟中，設置產(chǎn)生目錄時樣本標簽之間允許的錯配數(shù)目(n參數(shù))為5。

數(shù)據(jù)匹配：此步驟運行STACKS中的數(shù)據(jù)匹配程序，將個體分析程序中構(gòu)建的每個個體的一系列位點匹配到數(shù)據(jù)整合程序形成的目錄中，以確定每個個體的每個位點所含有的等位基因類型(位點結(jié)合形式)。在這個過程中，如果一個樣本的一個位點被匹配到目錄中的多個位點上，則這個位點將被刪除。

遺傳參數(shù)分析：此步驟運行STACKS中的遺傳參數(shù)分析程序，用于計算群體遺傳學參數(shù)，并能夠?qū)?shù)據(jù)輸出為一系列的標準輸出格式。在每個位點中保留一個隨機的SNP，從而保證每個SNP位點之間是獨立的，避免產(chǎn)生連鎖的SNPs位點。

SNPs分析群體的遺傳分析：利用Stacks v1.42軟件鑒定群體SNP，按照樣品最低測序深度>2，樣品缺失率<0.5，MAF>0.05過濾得到有效的SNP位點，分析富民枳天然居群和回歸種群中的遺傳多樣性參數(shù)，主要包括私人等位基因數(shù)目、平均觀測雜合度(Ho)、平均期望雜合度(He)、核苷酸多樣性(π)和平均近交系數(shù)(FIS)等遺傳多樣性指數(shù)。

2 結(jié)果與分析

2.1 RAD-seq測序數(shù)據(jù)質(zhì)控和過濾

本試驗中富民枳樣品A的測序堿基含量分布如圖1所示。序列的起始位置與測序的引物接頭相連，因此A、C、G、T在起始端會有所波動，后面會趨于穩(wěn)定。模糊堿基N所占比例越低，說明未知堿基數(shù)越少，測序樣本受系統(tǒng)AT偏好影響越小。虛線左側(cè)為測序片段1的統(tǒng)計，虛線右側(cè)為測序片段2的統(tǒng)計結(jié)果。結(jié)果表明，該樣品的文庫構(gòu)建質(zhì)量和測序質(zhì)量均可滿足后續(xù)分析。

橫坐標是測序片段堿基坐標，坐標表示測序片段上從5’到3’端依次堿基的排列；縱坐標是所有測序片段在該測序位置A、C、G、T、N堿基分別占的百分比，不同堿基用不同顏色表示。

圖1 樣品A的堿基組成分布圖

對質(zhì)量剪切后的質(zhì)控數(shù)據(jù)分別進行測序片段數(shù)、總堿基數(shù)、GC含量和Q30比例的統(tǒng)計，共獲得高質(zhì)量的片段數(shù)200 465 643，過濾后得到質(zhì)控數(shù)據(jù)59.07 G，平均每個樣品1.31 G，平均GC質(zhì)量分數(shù)為38.06%，堿基質(zhì)量Q30比例達到92.79%。由于所測序列的Q30數(shù)據(jù)較高，表明堿基出錯率很低，GC分布正常，數(shù)據(jù)量達到分析要求，建庫測序成功。詳細結(jié)果見表1。

表1 富民枳基因組DNA測序數(shù)據(jù)

2.2 SNP檢測與位點開發(fā)

SNP檢測結(jié)果發(fā)現(xiàn)，所有樣品酶切片段數(shù)量為141 728～234 584，平均值為165 034；酶切片段平均深度為19.45～35.33X，平均值為26.00X(表2)，表明SNP位點檢測成功。

表2 富民枳的SNP檢測結(jié)果

依據(jù)材料方法部分給出的SNP檢測標準，我們對樣本進行了過濾，在經(jīng)過填補和質(zhì)控后，共計得到1 589 744個SNP位點，每個樣品檢測到29 351～47 810個SNP位點，平均獲得35 327個SNP位點(表3)。根據(jù)樣品最低測序深度>2，樣品缺失率<0.5、次要基因型頻率(MAF)>0.05的選擇標準，從這些SNP中篩選出41 077個有效SNP，用于群體遺傳分析。樣本中SNP位點的雜合率為80.78%～97.41%，平均雜合率為91.56%；樣本中SNP位點的純合率為2.59%～19.22%，平均純合率為8.44%。另外，每個樣品的轉(zhuǎn)換與顛換的比值(Ti/Tv)為1.48～1.62，平均值為1.52。富民枳轉(zhuǎn)換與顛換的比值均大于1，說明富民枳的變異更多的是發(fā)生在嘌呤與嘌呤或者嘧啶與嘧啶之間的轉(zhuǎn)換。

2.3 遺傳多樣性分析

基于過濾得到的41 077個SNP標記，對富民枳天然種群和回歸種群分別完成了遺傳多樣性分析。結(jié)果表明，富民枳天然種群的私人等位基因數(shù)目為22 782，Ho值為0.458 7，He值為0.371 7，π值為0.455 6，F(xiàn)IS值為0.088；而回歸種群的私人等位基因數(shù)量、Ho值、He值、π值分別為6 699、0.398 3、0.291 1和0.330 6，F(xiàn)IS值為-0.244 4，見表4。

表3 富民枳的SNP信息統(tǒng)計

續(xù)(表3)

表4 富民枳天然種群和回歸種群的遺傳多樣性

3 結(jié)論與討論

基于二代測序的簡化基因組測序技術目前已作為主流方法應用于很多物種SNP標記的開發(fā)[17-20]。富民枳作為云南特有的極小種群植物，曾采用SSR標記研究其遺傳多樣性水平[1，21]，但是這些標記存在通量小、準確性較低、耗時耗力、成本高等局限性[22-23]，而且難以篩選出多態(tài)性足夠高的遺傳標記。而基于RAD-seq技術開發(fā)的SNP標記作為現(xiàn)在分子遺傳學中最重要的分子標記，是反映生物內(nèi)DNA序列變異程度的重要參數(shù)，并經(jīng)統(tǒng)計分析可獲得核苷酸多樣性信息[24]，具有數(shù)量多且分布豐富代表性高、遺傳穩(wěn)定性好、檢測快速、不受基因組序列的限制等特點[18，25-26]。目前，富民枳尚未獲得足夠的核酸序列數(shù)據(jù)，甚至也沒有蕓香科內(nèi)近緣物種的可參考基因組序列，RAD-seq是一種比較理想的研究技術。從本試驗所測得15個個體的有效片斷數(shù)看，數(shù)據(jù)量基本均勻，未出現(xiàn)個體數(shù)據(jù)量差異極大的現(xiàn)象，基本滿足后續(xù)分析的數(shù)據(jù)量要求。經(jīng)過預試驗，我們選擇EcoR Ⅰ進行限制性酶切。對富民枳進行基因組測序分析，共得到17.63 G的質(zhì)控數(shù)據(jù)，Q30比例達到92.71%。由于所測序列的Q30數(shù)據(jù)較高，表明堿基出錯率很低，GC分布正常，表明數(shù)據(jù)量達到預期目標，文庫構(gòu)建質(zhì)量較好，測序成功，可以進行后續(xù)分析。本次研究中我們共獲得SNP位點491 189個，通過樣品最低測序深度>2，樣品缺失率<0.5、次要基因型頻率(MAF)>0.05篩選以后，得到有效SNP位點41 077個。因此，本研究嘗試應用簡化基因組技術RAD-seq對富民枳開發(fā)SNP位點，取得了成功，為富民枳的群體遺傳分析奠定了基礎。

評估一個物種的遺傳變異程度是遺傳資源保護的基礎和重要內(nèi)容[27]。基于SNP位點的核苷酸多樣性(π)、觀測雜合度(Ho)、期望雜合度(He)和近交系數(shù)等參數(shù)常被用來表征群體的遺傳多樣性大小[20]。核苷酸多樣性(π)是反映生物體內(nèi)基因組DNA序列變異程度的重要參數(shù)，通過對基因組DNA的測序開發(fā)SNP位點，經(jīng)過計算獲得核苷酸多樣性信息[24]。尤其對于雙等位基因SNP位點來說，π值可以全面地測定一個群體的遺傳多樣性[28-30]。由于進化速率的不同，同一個物種中不同的DNA片段也可能具有不同的π值[31]。目前通過對SNP的分析發(fā)現(xiàn)，很多瀕危植物的π值都低于0.1[20，24，32]。本研究基于開發(fā)出的SNP位點，測得富民枳天然種群的平均觀測雜合度(Ho)、期望雜合度(He)和核苷酸多樣性(π)分別為0.458 7、0.371 7和0.455 6，而回歸種群的Ho值、He值、π值分別為0.398 3、0.291和0.330 6，是天然種群遺傳多樣性的87%、78%和73%，這顯然并未達到前人提出的“在100年內(nèi)保持珍稀瀕危植物90%以上的遺傳多樣性才能稱之為“成功保護”的標準。綜上所述，盡管富民枳天然種群和回歸種群的π、Ho和He值均較高，意味著富民枳現(xiàn)存植株在SNP水平上的遺傳多樣性較豐富，保護其遺傳資源具有重要意義。但是，回歸種群的遺傳多樣性仍未達到我們的保護目標，基于之前的研究結(jié)果，本項目計劃將回歸種群的樣本量配置調(diào)整至“65株”有效種質(zhì)資源的29株(或復本)以內(nèi)[21]，使得遺傳資源保護科學化。

富民枳的遺傳多樣性分析可為評價和指導富民枳拯救保護工程中遺傳資源的保護提供理論依據(jù)。為了保證新建種群個體間的遺傳距離保持最大，遺傳資源保護最豐富，就應構(gòu)建一個包括所有有效樣本的樣本集。綜合考慮環(huán)境壓力和資金投入因素，可最終確定每個有效樣本的備份數(shù)及恢復重建種群的規(guī)模。在野外恢復或重建種群時，確保相鄰或相近的植株屬于不同有效樣本，并結(jié)合現(xiàn)存種群密度(特別是不同有效樣本間的距離)，使新建種群個體間的遺傳距離保持最大。結(jié)合恢復或重建種群規(guī)模及其個體間的空間布局，通過ArcGIS在衛(wèi)星影像圖上模擬，確定保護小區(qū)面積、重建種群所需生境面積。