張濤，李學，張競，印軒鵬，賀水蓮

(云南農(nóng)業(yè)大學 園林園藝學院，云南 昆明 650201)

極小種群野生植物(plant species with extremely small populations，PSESP)指分布地域面積狹小或成間斷分布，且長期以來遭受各類不同外界自然環(huán)境因素的脅迫影響,已經(jīng)遠遠小于穩(wěn)定生存范圍界限內(nèi)的最小生存種群，且隨時都有可能發(fā)生瀕臨滅絕的野生植物種類[1]。馨香玉蘭(Magnoliaodoratissima)為木蘭科(Magnoliaceae Juss.)木蘭屬(MagnoliaLinn.)常綠喬木，其花潔白芳香，枝繁葉茂，有較強的凈化環(huán)境和抗污染能力，為優(yōu)質(zhì)的觀賞樹木[2]。馨香玉蘭僅產(chǎn)于云南省，主要分布于文山壯族苗族自治州的廣南、麻栗坡和西疇等地,其所有居群都分布在自然保護區(qū)外，缺乏針對性的保護。馨香玉蘭的野生資源僅剩下5個居群，均處于高大喬木之下，光照不足，生長前景很不樂觀，結(jié)實率低[3]；種子中存在抑制發(fā)芽的物質(zhì),且抑制效果顯著[4]。由此可見其種群生存狀況堪憂。其瀕危級別已屬極危，并列入我國Ⅱ級重點保護野生植物[5]，也是我國云南特有的極小種群物種，急需開展保護。而今，已有學者對其種群生態(tài)學[6-8]、保護生物學[9]、林學[10]、育苗[11-12]等方面開展了研究。然而，對于馨香玉蘭的分子標記開發(fā)、保護遺傳學研究還很薄弱。隨著分子生物學的發(fā)展，開發(fā)分子標記來對瀕危植物進行遺傳評價，從而提出精確保護策略，是保護瀕危植物的有效手段。

SLAF-seq(specific-locus amplified fragment sequencing)技術(shù)是通過構(gòu)建SLAF-seq數(shù)據(jù)庫，將具有特異性長度的片段(SLAF標簽)篩選出來并進行高通量測序，最后將得到的符合質(zhì)量要求的SLAF片段作為目標物種的全基因組信息，進而依據(jù)這些特異性SLAF標簽能夠在全基因范圍內(nèi)迅速鑒定準確性高的SNP變異位點信息，利用特異性位點可以構(gòu)建遺傳連鎖圖譜、系統(tǒng)進化樹、PCA分析、基因流分析等用于遺傳進化和群體結(jié)構(gòu)分析。SLAF測序技術(shù)在保護遺傳學領(lǐng)域已經(jīng)得到了廣泛利用，并取得了一定的成果。如：Zheng等[15]基于SLAF-seq技術(shù)對紅橘×枳(Citrusreticulata×Poncirustrifoliata)雜交群體進行了遺傳多樣性分析、Xia等[16]基于SLAF-seq技術(shù)對油棕(Elaeisguineensis)進行了遺傳多樣性和群體結(jié)構(gòu)的評價、Chang等[17]基于SLAF-seq技術(shù)對側(cè)柏(Platycladusorientalis)進行了遺傳保護研究。此類工作證實了SLAF-seq測序技術(shù)可以提供大量的SNPs以供群體遺傳學分析[18]。本研究以5個野生居群和2個栽培居群共70份馨香玉蘭資源為材料，采用SLAF-seq建立文庫，通過篩選特異性相關(guān)位點擴增片段，在全基因組范圍之內(nèi)開發(fā)大量穩(wěn)定性高、特異性強的基因SNP位點，可為馨香玉蘭的保護遺傳學、譜系地理學及種群生態(tài)學等研究提供分子標記。

1 材料與方法

1.1 實驗材料

選取云南省文山州西疇、麻栗坡及廣南的5個野生居群以及云南省昆明市金殿公園、昆明植物園的2個栽培居群共70份馨香玉蘭資源為實驗材料(表1)。采集其無病斑的新鮮葉片，用硅膠干燥。每個居群采集10個個體，個體間距離不低于10 m。

表1 7個馨香玉蘭居群的地理位置和取樣情況

1.2 DNA的提取

采用改進的CTAB法[21]提取馨香玉蘭樣品的基因組DNA，用NanoDrop 1000分光光度計評估所得DNA的濃度和質(zhì)量，并使用2%SDS-PAGE顯示結(jié)果；然后將DNA樣品稀釋至100 ng/μL，用于隨后的SLAF-seq技術(shù)分析。

1.3 酶切建庫

由于馨香玉蘭的基因組數(shù)據(jù)尚未有研究，因此無法確定其基因組結(jié)構(gòu)大小及GC濃度等相關(guān)信息。本研究對70個馨香玉蘭個體分別進行建庫處理。首先選取近緣種鵝掌楸(Liriodendronchinense)基因組(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/003/013/855/GCA_003013855.1_ASM301385v1/GCA_003013855.1_ASM301385v1_genomic.fna.gz)作為參考基因組進行酶切預測(鵝掌楸基因組大小為1.56 Gb、GC含量為39.06%)，然后利用RsaI+EcoRV-HF?限制性內(nèi)切酶組合進行酶切，對獲得的酶切片段(SLAF標簽)進行3′端加A處理、連接Dual-index[22]測序接頭、PCR擴增、純化、混樣、切膠選取目的片段，經(jīng)文庫質(zhì)檢合格后用Illumina平臺完成檢測。

1.4 Illumina HiSeqTM2500測序及產(chǎn)出數(shù)據(jù)的質(zhì)量分析

經(jīng)文庫質(zhì)檢合格后用Illumina HiSeqTM2500進行測序。測序后獲得的序列通過去接頭、去低質(zhì)量閱讀框和去污處理方法,最后獲得到干凈序列。為充分保證評估后本次酶切實驗結(jié)果的正確性，采用粳稻(Oryzasativassp.japonica)作為酶切實驗結(jié)果對照，進行了相同的數(shù)據(jù)處理參與建庫和測序。為有效保證每個樣本序列測序的數(shù)據(jù)質(zhì)量，對測序后所得的樣品序列數(shù)據(jù)進行GC比例和Q30數(shù)據(jù)分析。

1.5 SLAF標簽的獲得和SNP 標記的開發(fā)

利用BWA軟件[23]將序列測序reads比對在不同參考基因組上，統(tǒng)計不同序列染色體上的SLAF標簽和多態(tài)性SLAF標簽。按照不同序列之間的相識程度，對不同馨香玉蘭樣本的序列進行比對聚類,最后聚集在一起的reads來源于同一個SLAF標簽，不同SLAF標簽間的相似度相對于同一SLAF標簽在不同樣品間的序列相似度較低。一個SLAF標簽的不同樣品間序列有差異(即有多態(tài)性)，被稱之為多態(tài)性SLAF標簽。通常使用GATK[24]和SAMTOOLS[25]兩種篩選方法可以組合開發(fā)出SNP，將兩種篩選方法結(jié)合得出的SNP標記交集作為最終正確的SNP標記數(shù)據(jù)集，篩選的標準值為MAF>0.05。

2 結(jié)果與分析

2.1 建庫評估

電子酶切預測以鵝掌楸為參考基因組，采用組合限制性內(nèi)切酶RsaI+EcoRV-HF?進行酶切，將酶切片段為414～464 bp的序列定義為SLAF標簽，共得到127 393個SLAF酶切標簽(表2)。通過對粳稻檢測數(shù)據(jù)的評估監(jiān)控，從而判斷酶切方法的有效性。粳稻基因組大小為374.30 Mb，并通過SOAP[26]軟件將粳稻的測序reads與其參考基因組進行比對，比對后的結(jié)果見表3。經(jīng)比對，雙端比對效率在90.39%(一條序列兩端在參考基因組上的比對跨度介于500～1 000 bp的reads占總數(shù)的比例)，比對效率基本正常。樣品測序堿基分布情況見圖1。

表2 酶切預測確定的酶切方案信息統(tǒng)計

表3 粳稻測序reads比對結(jié)果統(tǒng)計

圖1 馨香玉蘭測序堿基分布

酶切效率是簡化基因組測序是否成功的關(guān)鍵性指標?；蜉d體上的結(jié)構(gòu)區(qū)域(如環(huán)狀結(jié)構(gòu)域、連續(xù)酶切位點等)復雜、基因組DNA樣品純度較低、酶切綜合利用持續(xù)時間短和質(zhì)量水平不足等這些客觀因素均有可能直接地影響限制性內(nèi)切酶的活力，從而造成酶切位點識別不準確而未能完全切開。根據(jù)reads插入片段中剩余酶切位點的具體檢測結(jié)果比例進行測算，其比例越高，酶切效率越好[27]。粳稻數(shù)據(jù)的酶切效率統(tǒng)計結(jié)果見表4。粳稻對照組的酶切效率為91.46%，表明酶切反應正常。90.39%的雙端比對效率和91.46%的酶切效率說明SLAF建庫正常。

表4 粳稻數(shù)據(jù)酶切效率評估統(tǒng)計

通過對照雙端比對序列在基因組中的定位，可以估算SLAF標簽的實際長度，并繪制對照讀長插入片段的長度分布圖(圖2)，從而估計實際的片段選擇范圍。序列的插入片段長度都在預期的選擇范圍之內(nèi)(圖2)，說明序列測序數(shù)據(jù)質(zhì)量正常，采用的這種測序分析方法可信度高[28]。

圖2 對照序列插入片段分布

2.2 測序數(shù)據(jù)統(tǒng)計與評估

為提高綜合分析的結(jié)果質(zhì)量，以讀長100 bp×2 的數(shù)據(jù)作為后續(xù)的實驗數(shù)據(jù)進行綜合評估，以粳稻的數(shù)據(jù)為參考評估馨香玉蘭資料庫的數(shù)據(jù)準確性。通過Illumina HiSeqTM2500平臺測序，總共得到了229.54 Mb reads的數(shù)據(jù)(表5)，在各個實驗樣本中所統(tǒng)計得到的數(shù)據(jù)讀長長度值的范圍為1 499 485～7 274 846個。其中：來自昆明金殿的JD4樣本所獲得數(shù)據(jù)量最大，為7 274 846個讀長；西疇縣對門山的DM7樣本數(shù)據(jù)量最小，為1 499 485個讀長。樣本測序所獲得的GC比例為41.5%～44.38%，其中GC比例的最大值44.38%為來自廣南縣的GN7樣本、GC比例的最小值41.5%為來自西疇縣竜樹山的LS7樣本、GC比例的平均值為43.01%。由此可見，GC比例普遍較低，說明已達到測序條件。測序質(zhì)量值 Q30 的范圍為78.54%～93.22%，平均值為86.96%。其中：Q30測序的最小值來自于麻栗坡縣小箐的XQ2樣本，其值僅為78.54%；來自于昆明植物園的KIB9樣本的Q30測序值最大，為93.22%，每個樣本的Q30值都在70%以上。表明測序技術(shù)的堿基錯誤率非常低，所獲數(shù)值均合格。

2.3 SLAF 標簽與 SNP 標記的開發(fā)

通過序列分析，從70份馨香玉蘭樣本的基因組中共獲得843 082個SLAF標簽。標簽的平均測序深度約是11.78 ×(表5)，其中多態(tài)性SLAF的標簽有287 843個。對獲得的多態(tài)性SLAF標簽進行分析，共得到2 998 167個馨香玉蘭群體SNP標記，各樣品中SNP標記數(shù)量為476 062～1 392 468個，各樣品檢測到的SNP完整度為15.87%～46.44%、雜合率為3.42%～10.69%。篩選去除完整度≤50%以及次要基因頻率≤5%的SNP，總共獲得了180 650個高度一致性的群體SNP位點，占SNP總量的6.03%。

表5 馨香玉蘭SLAF標簽統(tǒng)計和SNP 信息統(tǒng)計及對照測序結(jié)果

續(xù)表5

2.4 遺傳多樣性分析

基于180 650個高質(zhì)量SNP位點，對馨香玉蘭5個野生居群及2個栽培居群進行遺傳多樣性檢測(表6)。結(jié)果顯示：7個馨香玉蘭群體的觀測雜合度(Ho)范圍在0.136 3～0.355 1之間，其中最高的是XQ群體，最低的是JD群體；期望雜合度(He)的范圍在0.296 1～0.357 4之間，其中最高的是KIB群體，最低的是GN群體；Shannon指數(shù)(I)的范圍在0.455 3～0.535 9之間，其中最大值是KIB群體為0.535 9，最小值是GN群體為0.455 3；Nei多樣性指數(shù)的范圍在0.312 6～0.386 3之間，其中最大值是KIB群體為0.386 3，最小值為GN群體為0.312 6。7個馨香玉蘭群體在群體水平上都具有較為豐富的遺傳多樣性，其中KIB群體的遺傳多樣性水平相對較高，GN群體的遺傳多樣性水平相對較低。KIB為栽培居群，其遺傳多樣性指數(shù)相對處于較高水平，推測其個體可能來源于不同的野生居群。

表6 馨香玉蘭的遺傳多樣性指數(shù)

3 討論與結(jié)論

目前，我國在極小種群野生植物保育方面已做了大量工作，云南省在極小種群野生植物保育方面的成績尤為突出。隨著記錄了云南101種極小種群野生植物的主要識別特征、分布現(xiàn)狀、受威脅因素和主要保護建議等信息的圖書——《云南省極小種群野生植物保護名錄(2021版)》[29]的出版發(fā)行，將更進一步地促進云南極小種群野生植物的保育工作。隨著測序技術(shù)的發(fā)展，傳統(tǒng)的分子標記已無法滿足研究需求，越來越多的研究者開始采用高通量測序技術(shù)來進行極小種群的保護遺傳學研究。馬永鵬等[30]為弄清漾濞槭(Aceryangbiense)極小種群的形成與維持機制，更好地指導保護工作，對漾濞槭105個個體開展了重測序，從遺傳分化、基因流、遺傳多樣性、種群歷史、有害突變和近交等方面對漾濞槭進行了深入研究。楊豐懋等[31]對顯脈木蘭(Magnoliafistulosa)進行dd-RAD簡化基因組測序，對顯脈木蘭的遺傳多樣性，種群歷史及第四紀冰期的有效居群大小等進行了研究，對顯脈木蘭提出了精確保護策略。由此可見，使用高通量測序深層次挖掘物種的遺傳多樣性、種群結(jié)構(gòu)及歷史、有害突變等信息具有明顯的優(yōu)勢。

馨香玉蘭基因組龐大而復雜，目前尚無全基因組序列信息公開，開發(fā)傳統(tǒng)分子標記耗時耗力，不能滿足馨香玉蘭保護遺傳學研究需求。因此，本實驗運用SLAF-seq技術(shù)對該物種的特異性分子標記進行開發(fā)，并對其遺傳多樣性參數(shù)進行分析，結(jié)果顯示這些SNP標記在馨香玉蘭不同群體中表現(xiàn)出較為豐富的多態(tài)性，與金蕊等[9]利用RAPD技術(shù)對馨香玉蘭進行遺傳多樣性測定的結(jié)果(Shannon指數(shù)為0.568 6、Nei’s指數(shù)為0.390 3)相比，本研究得到的遺傳多樣性指數(shù)稍低，因為RAPD技術(shù)的穩(wěn)定性較差，是顯性遺傳，不能識別雜合子位點，這使得遺傳分析相對復雜，在遺傳多樣性計算時會因顯性遮蓋作用而使計算位點間遺傳距離的準確性下降，而使用SLAF-seq,穩(wěn)定性和準確性較高。

馨香玉蘭被列為第一批極小種群物種，但跟其它木蘭科植物相比，如大果木蓮(M.grandis)[32](I=0.365 1,H=0.243 3)、樂昌含笑(Micheliachapensis)[33](I=0.475 1,H=0.325 5)等，其遺傳多樣性水平卻相對較高(I=0.5069，H=0.355 3)。這說明馨香玉蘭雖然片段化明顯，僅剩的種群小，但仍具有較強的生存力，導致其瀕危的原因可能是外界環(huán)境的因素，如生存環(huán)境惡劣，上層喬木搶奪生存資源[33]，而不是由物種本身遺傳多樣性單一導致。本研究探索利用SLAF-seq開發(fā)分子標記，既可為馨香玉蘭未來的種群結(jié)構(gòu)、遺傳進化與種群演化歷史等提供高質(zhì)量的分子標記，同時也可為馨香玉蘭分子輔助選育、高密度遺傳連鎖圖譜的建立等方面提供參考。