陳祥平 呂 銀 劉 玲 柯皓天 劉凱旋 王茜齡 任艷紅 陳仁芳*

(1.四川省絲綢科學(xué)研究院，成都 610031；2.西南大學(xué)生物技術(shù)學(xué)院，重慶 400716)

?

中國(guó)桑屬15個(gè)種RAD-seq高通量測(cè)序*

陳祥平1呂 銀1劉 玲1柯皓天1劉凱旋1王茜齡2任艷紅2陳仁芳2*

(1.四川省絲綢科學(xué)研究院，成都 610031；2.西南大學(xué)生物技術(shù)學(xué)院，重慶 400716)

利用RAD-seq對(duì)中國(guó)桑屬15個(gè)種進(jìn)行了高通量測(cè)序，總共得到36.72Gb clean data，總Tags數(shù)2 788 927(reads) ，平均每個(gè)種原始數(shù)據(jù)都在10M以上，Tags都在20萬(wàn)條以上，質(zhì)量值Q30都在90%以上。用Stacks軟件對(duì)15個(gè)種進(jìn)行比對(duì)，獲得68904個(gè)SNPs位點(diǎn)。用最大似然法建樹(shù)，分支圖首先將白桑、廣東桑分出，接著是山桑、魯桑、瑞穗桑，再次分出的是雞桑、細(xì)齒桑、蒙桑和鬼桑，最后分出的是黑桑、川桑、華桑、滇桑、長(zhǎng)穗桑、奶桑。分支圖能將栽培種和野生種完全分開(kāi)；可以將蒙桑和鬼桑、雞桑、華桑、川桑、奶桑分開(kāi)。認(rèn)為白桑、廣東桑屬原始類(lèi)型，長(zhǎng)穗桑、奶桑屬進(jìn)化類(lèi)型；山桑、魯桑、瑞穗桑這三個(gè)種被分在一個(gè)分支，自檢支持率99%，黑桑、川桑這兩個(gè)種被分在一個(gè)分支，自檢支持率56%，長(zhǎng)穗桑、奶桑這兩個(gè)種被分在一個(gè)分支，自檢支持率100%，說(shuō)明這些種之間有較近的親緣關(guān)系，桑屬RAD-seq測(cè)序能大規(guī)模篩查SNPs位點(diǎn)，系統(tǒng)發(fā)育分析的準(zhǔn)確性就更加可靠。

中國(guó)桑屬；RAD-seq；SNP標(biāo)記；系統(tǒng)發(fā)育

桑屬(Morus)為桑科的模式屬，由Linne(1753)建立[1]，其后雖有許多學(xué)者進(jìn)行研究[2-12]，但由于學(xué)者們對(duì)物種的認(rèn)識(shí)不同，分出的種數(shù)差異很大。桑屬的分子系統(tǒng)學(xué)也先后有許多學(xué)者進(jìn)行研究[13-21]。但這些研究由于所用的方法、DNA片段，信息位點(diǎn)有限，不能將所有桑種分開(kāi)。

限制性?xún)?nèi)切酶位點(diǎn)相關(guān)DNA(Restriction-site Associated DNA, RAD)由Miller等[22]2007年提出。它是在新一代測(cè)序技術(shù)基礎(chǔ)上發(fā)展起來(lái)的一種DNA新技術(shù)。該技術(shù)利用限制性?xún)?nèi)切酶對(duì)基因組進(jìn)行酶切, 產(chǎn)生一定大小的片段, 構(gòu)建測(cè)序文庫(kù), 對(duì)酶切后產(chǎn)生的RAD標(biāo)記進(jìn)行高通量測(cè)序[23]。由于RAD標(biāo)記代表整個(gè)基因組特異性酶切位點(diǎn)附近的小片段DNA標(biāo)簽,又由于新一代測(cè)序技術(shù)通量高。因此，通過(guò)對(duì)RAD標(biāo)記測(cè)序能夠獲得成千上萬(wàn)的單核苷酸多態(tài)性(Single nuc-leotide polymorphism, SNP)標(biāo)記[24-26]。RAD-seq 已成功應(yīng)用于生物SNP標(biāo)記的開(kāi)發(fā)、超高密度遺傳圖譜的構(gòu)建、動(dòng)植物重要經(jīng)濟(jì)性狀的QTL定位、輔助全基因組重測(cè)序等研究領(lǐng)域[27-33]，在群體遺傳結(jié)構(gòu)、系統(tǒng)演化分析方面也在溫帶竹子、美國(guó)橡樹(shù)、北美瓶草蚊、馬先蒿、蝶屬、三刺魚(yú)上有研究[34-39]，但未見(jiàn)在桑屬上報(bào)道。因此，測(cè)定了中國(guó)桑屬15個(gè)種的RAD-seq序列 ，進(jìn)行SNPs開(kāi)發(fā)和系統(tǒng)進(jìn)化分析?，F(xiàn)將研究結(jié)果報(bào)道如下。

1 材料與方法

1.1 材料

材料取中國(guó)桑屬15個(gè)種，每個(gè)種取1-5位嫩葉后立即裝入冰盒，帶回實(shí)驗(yàn)室-70℃冰箱保存?zhèn)溆?，各個(gè)種采集地、緯度、經(jīng)度、海拔(表1)。

表1 桑屬采樣信息

1.2 實(shí)驗(yàn)流程

(1)提取每個(gè)樣品基因組DNA；

(2)采用限制性?xún)?nèi)切酶EcoRⅠ(GAATTC)酶切： 0.1-1ug基因組DNA在50uL體系中用20u的EcoRI于37℃消化15min，65℃20min使內(nèi)切酶失活；

(3)連接P1接頭：在消化后的DNA片段兩端加P1 Adapter，再次65℃20min使連接酶失活；

(4)片斷化：帶有不同P1 接頭的樣品混合在一起，采用物理方法打斷成300-500bp的片段，1%瓊脂糖凝膠電泳后回收300-700bp的DNA， 末端平化后加A；

(5)連接P2接頭： P2 adapter 是一個(gè)發(fā)散型的“Y” adapter, 它可以防止缺乏P1接頭的基因組片段擴(kuò)增，即只有兩端接頭種類(lèi)不同的片段才能被選擇性擴(kuò)增；

(6)PCR擴(kuò)增、純化、上機(jī)測(cè)序：取5uL，用P1和P2引物PCR擴(kuò)增， 18個(gè)循環(huán)后，跑膠純化回收300-500bpDNA片段，上機(jī)測(cè)序，測(cè)序平臺(tái)Illumina hiseq4000，測(cè)序方法Illumina/Solexa 聚合酶合成測(cè)序，測(cè)序深度3X，測(cè)序在華大基因進(jìn)行。

1.3 序列分析

(1)序列過(guò)濾

raw reads包含低質(zhì)量序列、adapter序列等，需要經(jīng)過(guò)一系列數(shù)據(jù)處理來(lái)，得到clean data，詳細(xì)步驟如下：去掉含有adapter的reads；去掉低質(zhì)量的數(shù)據(jù)(過(guò)濾參數(shù)為質(zhì)量值低于12占整條reads的40%或者以上，刪除整條reads；去除含N(表示無(wú)法確定堿基信息)比例大于3%的reads；去除序列5’端未包含樣本barcode(4-8bp)信息的reads(一個(gè)樣品對(duì)應(yīng)一個(gè)barcode)；剪切掉5′端barcode、酶切識(shí)別位點(diǎn)。如果剪切掉barcode的reads 5’端未含酶切識(shí)別位點(diǎn)，則去除該reads。

(2)用Stacks軟件[40]對(duì)reads進(jìn)行整理，通過(guò)ustacks→cstacks→sstacks→population程序處理，獲得SNPs位點(diǎn)。

(3)系統(tǒng)發(fā)育樹(shù)構(gòu)建：population步程序得到的phylip文件上傳到Cipres Science Gateway網(wǎng)站(https://www.phylo.org/)，RAxML-HPC BlackBox程序(Phylogenetic tree inference using maximum likelihood最大似然法)建樹(shù)。

2 結(jié)果

2.1 數(shù)據(jù)產(chǎn)出

經(jīng)過(guò)上述過(guò)濾raw data處理之后，總共得到36.72Gb clean data，總Tags數(shù)2 788 927(reads)，并根據(jù)個(gè)體barcode區(qū)別，將過(guò)濾后的reads細(xì)分至每個(gè)個(gè)體，數(shù)據(jù)產(chǎn)出(表2)。

表2 數(shù)據(jù)產(chǎn)出統(tǒng)計(jì)

從表2可以看出每個(gè)種原始數(shù)據(jù)都在10M以上，通過(guò)過(guò)濾raw data，用Stacks軟件處理后得到序列Tags都在20萬(wàn)條以上，質(zhì)量值Q30都在90%以上，且每個(gè)種得到Tags數(shù)量都較均勻，均可用于后續(xù)SNP標(biāo)記檢測(cè)和系統(tǒng)發(fā)育分析。

2.2 SNP標(biāo)記的檢測(cè)

用Stacks軟件對(duì)reads進(jìn)行整理，通過(guò)population程序，對(duì)15種進(jìn)行比對(duì)，獲得68 904個(gè)SNPs位點(diǎn)，為后續(xù)的系統(tǒng)發(fā)育分析的準(zhǔn)確性創(chuàng)造了條件。

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

用Stacks軟件population程序得到的phylip文件，上傳到Cipres Science Gateway網(wǎng)站(https://www.phylo.org/)，啟動(dòng)RAxML-HPC BlackBox程序，用最大似然法建樹(shù)，得到中國(guó)桑屬15個(gè)種的系統(tǒng)發(fā)育樹(shù)(圖3)。

從圖3可以看出，白桑、廣東桑首先分出，接著是山桑、魯桑、瑞穗桑，再次分出的是雞桑、細(xì)齒桑、蒙桑和鬼桑，最后分出的是黑桑、川桑、華桑、滇桑、長(zhǎng)穗桑、奶桑。山桑、魯桑、瑞穗桑這三個(gè)種被分在一個(gè)分支，自檢支持率99%，黑桑、川桑這兩個(gè)種被分在一個(gè)分支，自檢支持率56%，長(zhǎng)穗桑、奶桑這兩個(gè)種被分在一個(gè)分支，自檢支持率100%。

圖1 基于RAD-seq中國(guó)桑屬15個(gè)種系統(tǒng)發(fā)育(分支間的數(shù)據(jù)為50%以上的自檢支持率)

3 討論

3.1 桑屬RAD-seq測(cè)序要求

RAD-seq屬新一代高通量測(cè)序，對(duì)DNA質(zhì)量要求較高，樣品濃度：100-200ng/μl，樣品純度：OD260/280=1.8-2.0。因此，采樣要用冰盒采集鮮樣，盡量采嫩葉，立即送回實(shí)驗(yàn)室-70℃冷藏保存，DNA提取最好用Plant Genomic DNA Kit。另外，本次研究主要是開(kāi)發(fā)桑屬SNPs，進(jìn)行系統(tǒng)進(jìn)化分析，要求測(cè)序深度至少3X，Q30至少90%以上。

3.2 桑屬RAD-seq測(cè)序能大規(guī)模篩查SNPs位點(diǎn)

本研究采用Illumina hiseq4000測(cè)序平臺(tái)， Illumina/Solexa聚合酶合成測(cè)序，測(cè)序深度3X，對(duì)中國(guó)桑屬15個(gè)種進(jìn)行RAD-seq測(cè)序，過(guò)濾掉一些低質(zhì)量的reads，用Stacks軟件，population程序，對(duì)15種進(jìn)行比對(duì)，獲得68 904個(gè)SNPs位點(diǎn)，相比用常規(guī)方法(ITS序列116個(gè)SNPs位點(diǎn)，ITS、trnL-F、rps16三個(gè)片段合并261個(gè)SNPs位點(diǎn)[18-19])，SNPs位點(diǎn)可成萬(wàn)倍的增加，后續(xù)的系統(tǒng)發(fā)育分析的準(zhǔn)確性就更加可靠。

3.3 基于RAD-seq中國(guó)桑屬15個(gè)種系統(tǒng)發(fā)育樹(shù)的幾個(gè)新觀點(diǎn)

本文根據(jù)篩查到的68 904個(gè)SNPs位點(diǎn)，用最大似然法建立的系統(tǒng)發(fā)育樹(shù)有幾個(gè)新觀點(diǎn)提供給讀者商榷。

(1) 分支圖可以將栽培種和野生種完全分開(kāi)；

(2) 可以將蒙桑和鬼桑、雞桑、華桑、川桑、奶桑分開(kāi)(常規(guī)用ITS、trnL-F、rps16序列不能將蒙桑和鬼桑、雞桑與白桑分開(kāi)，不能將華桑、川桑、奶桑分開(kāi))；

(3) 白桑、廣東桑屬原始類(lèi)型，長(zhǎng)穗桑、奶桑屬進(jìn)化類(lèi)型；

(4) 山桑、魯桑、瑞穗桑親緣關(guān)系近，黑桑、川桑親緣關(guān)系近，長(zhǎng)穗桑、奶桑親緣關(guān)系近，細(xì)齒桑與蒙桑、鬼桑親緣關(guān)系近；

(5) 中國(guó)植物志將山桑定為白桑的變種[41]，理由是山桑無(wú)自然地理分布，而本研究將山桑、魯桑、瑞穗桑分在一個(gè)分支，自檢支持率99%，魯桑、瑞穗桑也無(wú)自然地理分布。因此，認(rèn)為魯桑、瑞穗桑不作種級(jí)單位為宜；

(6) 滇桑與長(zhǎng)穗桑、奶桑親緣關(guān)系近，且滇桑為長(zhǎng)穗類(lèi)，應(yīng)為長(zhǎng)穗類(lèi)里的一個(gè)種，而不作為蒙桑的變種。

3.4 本研究需要進(jìn)一步改進(jìn)之處

(1)由于川桑全基因已測(cè)序，應(yīng)根據(jù)川桑全基因信息篩選合適的內(nèi)切酶進(jìn)行酶切建庫(kù)，有目的的篩查SNPs位點(diǎn)；

(2)stacks軟件應(yīng)選pstacks程序，而不是選ustacks程序，因?yàn)樯僖延写ㄉＨ蚪M信息；

(3)每個(gè)種應(yīng)取更多樣，在population步，統(tǒng)計(jì)分析桑群體遺傳學(xué)Pi、Fis、Fst等相關(guān)數(shù)據(jù)。

[1] LINNE C V.1753.Species plantarum[M]. 2：986.

[2] MORETTI G.Prodromo di una monografia delle specie del centreMorus[M].1842：564.

[3] SERINGE N C.Description,culture et taille des murierrs[M].1855(98)：423-425.

[4] BUREAU E.MORACEAE.Prodromus systematis naturalis regni vegetabilis.Volume 17 XVIL[M].AIph.de Candolle.(ed.).Paris, France.1873：211-288.

[5] SCHNEIDER C K.lantae Wilsonianae[J].Sarg.PI.Wils.1916，3(2)：296-297.

[6] KOIDZUMI G.Taxonomical discussion on Morus plants[J].Bull Imp Sericult Exp Stat，1917(3)：1-62.

[7] 陳嶸.中國(guó)樹(shù)木分類(lèi)學(xué)[M].北京：科學(xué)出版社,1937：228-231.

[8] 胡先骕.植物分類(lèi)學(xué)簡(jiǎn)編[M]．北京：科學(xué)技術(shù)出版社，1957：56.

[9] 張秀實(shí).桑科新分類(lèi)群[J].植物分類(lèi)學(xué)報(bào)，1984，22(1)：64-76.

[10]吳征鎰，張秀實(shí)．中國(guó)?？频囊恍┬路诸?lèi)單位[J]．云南植物研究，1989，11(1)：24-34 .

[11]曹子余.中國(guó)桑屬(桑科)新分類(lèi)群[J]．植物分類(lèi)學(xué)報(bào)，1991，29(3)：264-267.

[12]曹子余.中國(guó)植物志資料(?？? [J]．云南植物研究， 1995，17(2)：153-154，158.

[13]向仲懷,張孝勇,余茂德,等.采用隨機(jī)擴(kuò)增多態(tài)性DNA技術(shù)(RAPD)在桑屬植物系統(tǒng)學(xué)研究的應(yīng)用初報(bào)[J]．蠶業(yè)科學(xué)，1995(4)：203-208.

[14]馮麗春，楊光偉，余茂德，等.利用RAPD對(duì)桑屬植物種間親緣關(guān)系的研究[J]．中國(guó)農(nóng)業(yè)科學(xué)，1997，30(1)：52-56.

[15]楊光偉,馮麗春,張孝勇,等.桑屬種群遺傳結(jié)構(gòu)變異分析[J]．蠶業(yè)科學(xué)，2003，29(4)：323-329.

[16]趙衛(wèi)國(guó), 潘一樂(lè), 張志芳.桑屬植物ITS序列研究與系統(tǒng)發(fā)育分析[J].蠶業(yè)科學(xué)，2004，30(1)：11-14.

[17]汪偉,王興科,朱昱萍.基于trnL內(nèi)含子序列的桑屬植物分子系統(tǒng)學(xué)初報(bào)[J].蠶業(yè)科學(xué)，2008，34(2):89-103.

[18]陳仁芳，余茂德，劉秀群，等.桑種質(zhì)資源ITS序列與系統(tǒng)進(jìn)化分析[J].中國(guó)農(nóng)業(yè)科學(xué)，2010，43(8):34-42.

[19]陳仁芳,張澤，唐洲，等.桑屬I(mǎi)TS、trnL-F、rps16序列與進(jìn)化分析[J].中國(guó)農(nóng)業(yè)科學(xué)，2011，44(8)：1553-1561.

[20]MADHAV P.Systematics and Reproductive Biology of the GenusMorusL.(Moraceae)[D].Kansas：Kansas State University，2008.

[21]MADHAV P NEPAL,CAROLYN J FERGUSON.Phylogenetics of Morus(Moraceae) Inferred from ITS and trnL-trnF Sequence Data[J].Systematic Botany, 2012，37(2)：442-450.

[22]MILLER M R,DUNHAM J P,AMORES A,et al.Ranid and costeffective polvmornhism identification and genotvping using restriction site associated DNA(RAD)markers[J].Genome Res,2007，17(2):240-248.

[23]H.C.ROWE, S.RENAUT，A.GUGGISBERG.RAD in the realm of next-generation sequencing technologies[J].Molecular Ecology,2011，20：3499-3502.

[24]NATHAN A.BAIRD, PAUL D.ETTER, TRESSA S.ATWOOD, et al.Rapid SNP discovery and geneticmapping using sequenced RAD markers[J].PlosOne,2008，3(10):e3376.

[25]VAN TASSELL CP, SMITH TP, MATUKUMALLI LK, et al.SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries[J].Nat Methods,2008，5(3)：247-252.

[26]J.M.PUJOLAR,M.W.JACOBSEN,J.FRYDENBER G, et al.A resource of genome-wide single-nucleotide polymorphisms generated by RAD tag sequencing in the critically endangered European eel[J].Molecular Ecology Resources,2013，13：706-714.

[27]JOHN W.DAVEY，MARKL.BLAXTER.RADSeq：next-generation population Genetics[J].Briefings in Functional Genomics，2011，9(5)：416-423.

[28]ADAM D.L EACHé,BARBARA L.BANBURY,JOSEPH FELSENSTEIN,et al.Short Tree, Long Tree, Right Tree, Wrong Tree：New Acquisition Bias Corrections for Inferring SNP Phylogenies[J].Systematic Biology Advance，2015，00(0)：1-16.

[29]岳桂東, 高強(qiáng), 羅龍海, 等.高通量測(cè)序技術(shù)在動(dòng)植物研究領(lǐng)域中的應(yīng)用[J].中國(guó)科學(xué)：生命科學(xué)，2012，42(2)：107-124.

[30]HOHENLOHE PA, AMISH JS, CATCHEN MJ, et al.Next-generation RAD sequencing identifies thousands of SNPs for assessing hybridization between rainbow and westslope cutthroat trout[J]. Mol Ecol Resour,2011，11(Suppl.1)：117-122.

[31]PFENDER WF, SAHA MC, JOHNSON EA, et al.Mapping with RAD (restriction-site associated DNA) markers to rapidly identify QTL for stem rust resistance in Lolium perenne[J].Theor Appl Genet,2011，122(8)：1467-1480.

[32]POLAND JA, BROWN PJ, SORRELLS ME, et al.Development of high-density genetic maps for barley and wheat using a novel two-enzyme genotyping-by-sequencing approach[J].PlosOne,2012，7(2)：e32253.

[33]劉艷玲.蓮野生居群遺傳多樣性評(píng)價(jià)及高密度遺傳連鎖圖譜的構(gòu)建[D].武漢：華中農(nóng)業(yè)大學(xué)，2013.

[34]ANDREW L.HIPP, DEREN A.R.EATON, JEANNINE CAVENDER-BARES, et al.A Framework Phylogeny of the American Oak Clade Based on Sequenced RAD Data[J].PlosOne，2014,9(4):1-12.

[35]X.Q.WANG,L.ZHAO,D.A.R.EATON,et al.Identification of SNP markers for inferring phylogeny in temperate bamboos (Poaceae：Bambusoideae) using RAD sequencing[J].Molecular Ecology Resources,2013，13：938-945.

[36]JULIAN CATCHEN,SUSAN BASSHAM,TAYLOR WILSON,et al.The population structure and recent colonization history of Oregon threespine stickleback determined using restriction-site associated DNA-sequencing[J].Molecular Ecology，2013，22：2864-2883.

[37]KEVIN J.EMERSON, CLAYTON R.MERZ, JULIAN M.CATCHEN,et al.Resolving postglacial phylogeography using high-throughput sequencing[J].PNAS，2010，107(31)：16196-16200.

[38]DEREN A.R.EATON,RICHARD H.REE.Inferring Phylogeny and Introgression using RADseq Data：An Example from Flowering Plants (Pedicularis：Orobanchaceae)[J].Syst.Biol，2013，62(5):689-706,

[39]NICOLA J.NADEAU, SIMON H.MARTIN, KRZYSZTOF M.KOZAK, et al.Genome-wide patterns of divergence and gene flow across a butterfly radiation[J].Molecular Ecology，2013，22(3)：814-826.

[40]CATCHEN J,HOHENLOHE P A,BASSHAM S,et al.Stacks:an analysis tool set for population genomics[J].Molecular Ecology,2013，22:3124-3140.

[41]中國(guó)植物志編輯委員會(huì).中國(guó)植物志23卷[M].北京：科學(xué)出版社，1998：2-23.

RAD-seq High-Throughput Sequencing of 15 Species ofMorus

CHEN Xiang-ping1Lü Yin1LIU Ling1KE hao-tian1LIU Kai-xuan1WANG Xi-ling2REN Yan-gong2CHEN Ren-fang2*

(1.TheSilkResearchInstituteofSichuanProvince,Chengdu610031,China；2.CollegeofBiotechnology,SouthwestUniversity,Chongqing400716,China)

In this study, RAD-seq (restriction site-associated DNA sequencing) was performed toidentify the SNP loci between 15 species of the genusMorus. A total of 36.72 Gb clean data were sequenced by Illumina hiseq 4000 and a total of 68904 SNP loci were obtained by Stacks software. The data sets were analysed using the maximum-likelihood method. The cultivated and wild species were completely separate,M.albaandM.atropurpureawere original genera of mulberry, however,M.wittiorumandM.macrourawere of the evolutionary type.M.serrata,M.mongolica, andM.mongolicavar.diabolicahad close genetic relationship,whileM.yunnanensis,M.wittiorumandM.macrourahad close genetic relationship. Analyses based on these RAD tags yielded robust phylogenetic inferences, even with data set constructed from surprisingly few loci. The study illustrates the potential for resolving difficult phylogenetic relationships in genusMorus.

Morus; Restriction site-associated DNA; SNP marker; Phylogeny

* 資助項(xiàng)目：四川省科技廳應(yīng)用基礎(chǔ)研究計(jì)劃項(xiàng)目。通訊作者：陳仁芳，博士，副教授，從事桑樹(shù)學(xué)研究。