RAD-seq技術(shù)研究鵝掌楸屬種源遺傳多樣性和遺傳結(jié)構(gòu)*

潘文婷 孫建軍 原勤勤 張利利 鄧康橋 厲月橋

(中國林業(yè)科學(xué)研究院亞熱帶林業(yè)實驗中心 分宜 336600)

木蘭科(Magnoliaceace)鵝掌楸屬(Liriodendron)為落葉大喬木，葉大形似馬褂，故又稱之為“馬褂木”，歐洲人稱其為“郁金香樹”，是珍貴用材和優(yōu)良觀賞樹種?，F(xiàn)僅存兩個種，即北美鵝掌楸(Liriodendrontulipifera)和鵝掌楸(Liriodendronchinense)，為孑遺樹種，前者生長于美國東部及加拿大南部的闊葉林中，從平原到山區(qū)呈連續(xù)分布，遺傳資源十分豐富； 后者大多零星分布于我國長江流域以南的亞熱帶中、低山區(qū)，?；焐诔＞G或落葉闊葉林中(郝日明等， 1995)，雖分布范圍較大且遺傳多樣性良好，但因種群數(shù)量稀少，種群片段化嚴重，天然更新不良，處于瀕危狀態(tài)，已被列入我國二級珍稀瀕危保護植物(王章榮， 2005)。鵝掌楸屬中的這兩個種是典型的東亞-北美間斷分布“種對”(Vicariad Species Pairs)(Parksetal., 1990)，是植物群體遺傳學(xué)和分子系統(tǒng)發(fā)育地理學(xué)的理想材料(李康琴， 2013)。

隨分子標記研究方法的增多，對鵝掌楸分子育種研究也逐漸完善，王曉陽等(2011)利用SSR技術(shù)研究鵝掌楸苗期生長雜種優(yōu)勢； 李康琴(2013)基于SSR技術(shù)對鵝掌楸屬進行群體遺傳結(jié)構(gòu)及分子系統(tǒng)地理學(xué)研究； 羅群鳳等(2015)采用同源克隆和RACE技術(shù)克隆北美鵝掌楸的查爾酮合成酶基因(LtCHS)，并對其進行生物信息學(xué)及組織表達進行了分析； 祁荔(2017)利用AFLP及SSR分子標記構(gòu)建鵝掌楸遺傳圖譜及重要性狀QTL初步定位。近年來基于全基因組酶切位點相關(guān)的簡化基因組測序基礎(chǔ)上發(fā)展出來了二代測序技術(shù)RAD-seq (Restriction-site associated DNA sequencing)，該技術(shù)由Miller等(2007)提出，主要利用酶切、序列捕獲降低基因組復(fù)雜度從而獲得部分基因組序列信息(Wangetal., 2017; Jérémyetal., 2020)，作為一種快速有效的手段，實現(xiàn)大規(guī)模的分子標記開發(fā)，既節(jié)約成本又節(jié)省時間，已成功應(yīng)用于 SNP 標記的開發(fā)、動植物重要經(jīng)濟性狀的 QTL 定位、生物進化等研究領(lǐng)域(Fengetal., 2020; Lietal., 2020)。陸葉等(2019)利用RRAD-seq 技術(shù)開展了鵝掌楸基因組 SNP 標記開發(fā)，最終獲得3 501個候選SNP標記； Sheng 等(2021)利用RAD-seq技術(shù)分析了鵝掌楸Myb基因家族的表征及其在非生物脅迫應(yīng)答中的作用，可見RAD-seq技術(shù)能提供更多SNP標記，有效補充了鵝掌楸遺傳信息較少這一問題，但有關(guān)鵝掌楸屬群體遺傳結(jié)構(gòu)和遺傳分化的研究較少。

開展群體結(jié)構(gòu)和遺傳分化的研究能精準評估遺傳多樣性和制定科學(xué)有效的種質(zhì)資源保護策略，李云飛等(2019)采用RAD-seq技術(shù)開發(fā)了85種杜鵑花屬植物的高質(zhì)量SNP位點，通過高通量測序探討分類，在亞屬水平的分類取得很好的效果。黃承玲等(2021)采用 RAD-seq技術(shù)對34種杜鵑花屬(Rhododendron)植物進行分類，證實了RAD-seq技術(shù)在復(fù)雜植物類群的物種分類方面比傳統(tǒng)分子標記具有明顯優(yōu)勢。本研究利用RAD-seq測序技術(shù)，以9個中國鵝掌楸和4個北美鵝掌楸種源為研究材料，通過高通量測序，對13個鵝掌楸屬種源143份樣本進行了遺傳多樣性、遺傳變異、遺傳分化以及群體基因交流等分析，揭示鵝掌楸屬各種源親緣關(guān)系，以期為鵝掌楸屬分子鑒定、種質(zhì)創(chuàng)新及種質(zhì)資源收集與保存提供參考。

1 材料與方法

1.1 試驗材料

試驗材料共有13個鵝掌楸屬種源，其中9個來自中國7個省，另外4個來自美國，相關(guān)種源地理信息詳見表1。所有種源種植在中國林業(yè)科學(xué)研究院亞熱帶林業(yè)實驗中心年珠實驗林場鵝掌楸種源試驗林中。分別于2016年10月和2017年5月進行采樣，每個種源隨機采取10～14株樹，每株樹取5～6片葉，共采集了143份鵝掌楸屬樣本，其中有97個樣本采自9個中國鵝掌楸種源，還有46個樣本采自4個北美鵝掌楸種源。樣本用自封袋裝好并倒入硅膠干燥、保存。

1.2 基因組DNA的制備

全部143份鵝掌楸屬樣本的基因組DNA采用改良的CTAB法(Liuetal., 2008)提取，用 1.0%瓊脂糖凝膠檢測DNA質(zhì)量。

1.3 酶切建庫

利用RAD-seq技術(shù)對143份鵝掌楸屬樣本的基因組DNA進行簡化測序。采用限制性內(nèi)切酶AvaII與MspI酶切組合對全基因組DNA完全酶切(Petersonetal., 2012)，(插入片段范圍： 500～600 bp，將酶切產(chǎn)物進行5’末端修復(fù)，同時對5’末端進行磷酸化修飾，3’末端加A，使之與接頭5’端T互補，提高接頭連接效率，阻止接頭自連。然后連接測序接頭，將連接產(chǎn)物錨定在flowcell上，進行橋式擴增。用瓊脂糖凝膠電泳進行片段大小選擇，通過PCR擴增，增大文庫量，建好的文庫用Illumina HiSeq X10進行測序。

表1 各參試種源地理信息及參試數(shù)量Tab.1 Geographic information and sample number of 13 provenances in the experiment

1.4 統(tǒng)計分析

對原始數(shù)據(jù)進行識別、過濾、質(zhì)量評估，再進行單樣本測序數(shù)據(jù)的 RAD 標記開發(fā)、樣本群體的RAD標記開發(fā)、多態(tài)性標記SNPs開發(fā)等，最后開展群體多樣性分析，包括進化樹構(gòu)建及群體結(jié)構(gòu)分析等。

主要通過Stacks軟件(http:∥creskolab.uoregon.edu/stacks/)(Catchenetal., 2013)進行RAD標記開發(fā)?；蚪M經(jīng)酶切處理打斷為多個小片段，每個片段相當(dāng)于一個標記位點，同一位點的序列通過相似性聚類，形成一個Stack。每個Stack中存在數(shù)條高深度片段，其余全為低深度片段。高深度片段即為潛在基因型，低深度片段可能由于測序錯誤導(dǎo)致。

采用Mega5.1軟件(http:∥www.atgc-montpellier.fr/phyml/)(Tamuraetal., 2011)中的Upgma算法和Neighbour-joining分別構(gòu)建13個種源的關(guān)系樹和143份樣本的進化樹。通過 Admixture軟件(Lexanderetal., 2009)對143份樣本的群體結(jié)構(gòu)進行聚類分析，該軟件運算速度快，創(chuàng)建Plink(http:∥pngu.mgh.harvard.edu/～purcell/plink/)的輸入文件，用Plink進行SNP過濾，將樣品的分群數(shù)(K值)設(shè)定為1-6進行聚類，根據(jù)CV error最小值確定分群數(shù)。

使用POPGENE 1.32軟件(http:∥www.ualberta.ca/～fyeh/download.htm.)(Yehetal., 2000)計算觀察雜合度(Ho)、期望雜合度(He)、核苷酸多樣性(π)和基因分化系數(shù)(Gst)等主要遺傳結(jié)構(gòu)的度量參數(shù)。

2 結(jié)果與分析

2.1 RAD-seq測序及質(zhì)量評估

本研究通過RAD-seq測序技術(shù)，完成13個鵝掌楸屬種源的143份樣本的測序，總堿基數(shù)為396.95 Gb，樣本平均獲取的堿基數(shù)為2.76 Gb； 測得有效數(shù)據(jù)為1 323.17 Mb，樣本平均獲取的有效數(shù)據(jù)為9.21 Mb； 共開發(fā)RAD標簽 2 611 618個，樣本的平均標簽數(shù)為141 222個，平均測序深度為25.02 x(表2)。所測序列GC含量較低(平均值為26.51%)，Q30數(shù)據(jù)較高(平均值為95.07%)，表明錯誤率低，測序結(jié)果可靠。通過對樣本間進行多態(tài)性分析，共得到高質(zhì)量的SNP位點4 454個。

表2 RAD-seq測序數(shù)據(jù)統(tǒng)計Tab.2 Statistics of RAD-seq sequencing data of 13 provenances

2.2 種源遺傳多樣性統(tǒng)計

對13個鵝掌楸屬種源的遺傳多樣性相關(guān)參數(shù)進行統(tǒng)計(表3)，結(jié)果表明，Ho為0.013 2～0.049 6，均值為0.037 9；He為0.022 2～0.059 7，均值為0.047 7；π為0.023 6～0.064 0，均值為0.050 6。采用隸屬函數(shù)法(彭松等， 2014)對Ho、He和π進行綜合評價(表4)，其中遺傳多樣性排前三位的種源依次為EX，NK和BK，DBS種源的遺傳多樣性最低。

表3 各種源遺傳結(jié)構(gòu)參數(shù)統(tǒng)計Tab.3 Statistics of genetic structure in different provenances

表4 13個鵝掌楸屬種源遺傳多樣性各指標的隸屬函數(shù)值Tab.4 Subordinative function values of genetic diversity parameters among 13 provenances in Liriodendron

2.3 遺傳結(jié)構(gòu)分析

對13個鵝掌楸屬種源的Gst和Nm數(shù)據(jù)(表3)進行分析，結(jié)果表明Gst從0.196 8到0.439 8，鵝掌楸種源之間存在較大的遺傳分化(Gst=0.241 9)以及中等的基因流(Nm=0.805 1，其中種源EX的Nm>1)，北美鵝掌楸種源之間存在很大的的遺傳分化(Gst=0.388 6>0.25)以及低水平的基因流(Nm=0.397 0)，表明北美鵝掌楸遺傳變異主要存在于種源間。

基于篩選的有效SNP，通過Mega5.1 軟件中的Upgma 算法和Neighbour-joining分別構(gòu)建13個種源的關(guān)系樹(圖1)和143份樣本的進化樹(圖2)，結(jié)果顯示，13個鵝掌楸屬種源被分為3個類群，其中XY、YN、EX、SZ和MN 5個種源為一個類群； LY、DBS、LS和WYS 4個種源為一個類群； NK、LYS、MSL和BK 4個種源為一個類群。

圖1 13個種源的關(guān)系樹Fig. 1 Relation tree of 13 provenances in Liriodendron

圖2 143份樣本的進化樹Fig. 2 Neighbour-joining phylogram illustrating genetic relationships among 143 individuals每一個分枝代表一份種質(zhì)材料，樹枝長度代表兩個物種間的進化距離，樹分支上的數(shù)字代表該分支的支持的百分數(shù)。Each branch represents one piece of germplasm, and the length of the branch represents the evolutionary distance between two species, the number on a branch of the tree represents the percentage of support for that branch.

圖3 分群數(shù)為2～4的聚類圖Fig. 3 Bar plot of 2-4 clusters identified with adegenet R package左圖圖中每種顏色代表一個群體，每列代表一個樣品的情況； 圖中展示了143個樣品分群值從2～4的聚類情況。右圖中為每個K值對應(yīng)的CV error值，K為3的時候最小。The picture left shows the clustering of 143 samples with cluster values from 2 to 4 (each color represents a subgroup，Each column represents a sample). The figure on the right shows the CV error value corresponding to each K value (When K is 3, the CV error value is minimal).

通過Admixture 軟件對143份樣本的群體結(jié)構(gòu)進行聚類分析，假設(shè)樣品的分群數(shù)(K值)為1～6，對應(yīng)的CV error值分別為： 0.465 4，0.178 6，0.142 1，0.153 1，0.163 2，0.168 9(圖3)。其群體結(jié)構(gòu)最優(yōu)分群數(shù)為3，分群的結(jié)果與關(guān)系樹和進化樹分析所得結(jié)果基本一致。

3 討論

本研究通過RAD-seq測序技術(shù)，完成13個鵝掌楸屬種源的143份樣本的測序，平均測序深度為25.02 ×，共獲得總堿基數(shù)396.95 Gb，平均每個樣本的堿基數(shù)為2.76 Gb，獲得高質(zhì)量SNPs標記4 454個。比李云飛等(2019年)獲得的3 501個SNP 標記多，這與本研究有13個種源的研究材料，而李云飛等(2019年)為2個種源有關(guān)，但均比以往鵝掌楸的分子標記開發(fā)技術(shù)獲得的候選SNPs數(shù)量有顯著提高。這些標記進一步增加了鵝掌楸屬的基因組資源量，可更準確地界定群體分類，為鵝掌楸屬資源的遺傳保存、育種及改良應(yīng)用提供依據(jù)(李云飛等， 2019； 黃承玲等， 2021)。

以獲得的高質(zhì)量SNP位點為基礎(chǔ)，進一步分析13個鵝掌楸屬遺傳多樣性和遺傳結(jié)構(gòu)。遺傳多樣性排前三位的種源依次為EX，NK和BK，DBS 種源的遺傳多樣性最低。但總體與報道的鵝掌楸遺傳多樣性數(shù)值(李康琴， 2013； 趙亞琦等， 2014)相比較低，可能是由于遺傳多樣性的估算受樣本數(shù)量、采樣策略及基因分型數(shù)據(jù)等因素的影響，本研究樣本數(shù)量較少，且為種源試驗林內(nèi)采樣，因此遺傳多樣性較低。

對13個鵝掌楸屬種源的Gst和Nm數(shù)據(jù)進行分析，結(jié)果表明Gst從0.196 8到0.439 8，鵝掌楸種源之間存在較大的遺傳分化以及中等的基因流(Gst=0.241 9、Nm=0.805 1，其中種源EX的Nm>1)，北美鵝掌楸種源之間存在很大的的遺傳分化以及低水平的基因流(Gst=0.388 6>0.25、Nm=0.397 0)，表明鵝掌楸屬遺傳分化變異主要存在于種源間； EX種源的Nm值大于1，表明基因流可以防止由遺傳漂變引起的群體間的遺傳分化，但各種源間Nm均值只有0.679 6； 且多數(shù)種源Ho小于期望雜合度He。群體間的高度分化的根源在于由生境片段化而引起的遺傳漂變和群體間低水平的基因流(Zaoualietal., 2012)，可見鵝掌楸屬存在由小群體效應(yīng)和片段化影響導(dǎo)致的瀕危現(xiàn)象，該結(jié)論與李康琴(2013)的研究的結(jié)果一致。

通過對13個鵝掌楸屬種源的143份樣本進行聚類分析和群體結(jié)構(gòu)分析，可將13個鵝掌楸屬種源分為3個類群，其中 XY、YN、EX、SZ和MN 5個種源為類群1； LY、DBS、LS和WYS 4個種源為類群2； NK、LYS、MSL和BK 4個種源為類群3。經(jīng)過分析可見類群1和類群2分別處在中國西部和東部，該結(jié)論進一步驗證了朱曉琴等(1995)和李康琴(2013)的研究結(jié)果，可見鵝掌楸屬遺傳結(jié)構(gòu)的形成與其地理隔離和片段化分布有關(guān)。中國東部種源群(類群2)較中國西部種源群(類群1)遺傳多樣性高(Ho： 0.036 1>0.034 9)，但均低于北美鵝掌楸(類群3)遺傳多樣性，該結(jié)論與李康琴(2013)的結(jié)論一致。進化樹和聚類圖均顯示3個類群間出現(xiàn)了少數(shù)的個體混雜現(xiàn)象，表明群體間有基因交流的現(xiàn)象，因此，利用這些材料進行雜交育種時，既要考慮其親緣關(guān)系也要注意他們的遺傳結(jié)構(gòu)。

4 結(jié)論

采用RAD-seq技術(shù)對鵝掌楸屬13個種源143份樣本進行測序，獲得了高質(zhì)量SNP位點，進而分析鵝掌楸屬遺傳多樣性、群體結(jié)構(gòu)分析及進化樹的構(gòu)建，鵝掌楸屬13個種源被分為3個類群，其中北美鵝掌楸和鵝掌楸可顯著區(qū)分、中國東部鵝掌楸和西部鵝掌楸的差異也很顯著。本研究結(jié)果表明高效的RAD-seq技術(shù)獲得的大量SNP位點可應(yīng)用于鵝掌楸屬物種的區(qū)分，進一步證實了鵝掌楸屬遺傳結(jié)構(gòu)的形成與其地理隔離和片段化分布有關(guān)，且存在由于小群體效應(yīng)和片段化影響導(dǎo)致的瀕?，F(xiàn)象； 同時這些SNP標記可進一步用于鵝掌楸屬Q(mào)TL定位、遺傳連鎖圖譜、分子輔助育種等研究，進而為鵝掌楸屬的遺傳演化研究提供理論基礎(chǔ)，為今后鵝掌楸屬種質(zhì)創(chuàng)新及種質(zhì)資源收集與保存提供參考依據(jù)。