中國桑屬15個種RAD-seq高通量測序*

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

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

?

中國桑屬15個種RAD-seq高通量測序*

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

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

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

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

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

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

1 材料與方法

1.1 材料

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

表1 桑屬采樣信息

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

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

(2)采用限制性內(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 是一個發(fā)散型的“Y” adapter, 它可以防止缺乏P1接頭的基因組片段擴(kuò)增，即只有兩端接頭種類不同的片段才能被選擇性擴(kuò)增；

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

1.3 序列分析

(1)序列過濾

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

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

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

2 結(jié)果

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

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

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

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

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

用Stacks軟件對reads進(jìn)行整理，通過population程序，對15種進(jìn)行比對，獲得68 904個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/)，啟動RAxML-HPC BlackBox程序，用最大似然法建樹，得到中國桑屬15個種的系統(tǒng)發(fā)育樹(圖3)。

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

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

3 討論

3.1 桑屬RAD-seq測序要求

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

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

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

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

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

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

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

(3) 白桑、廣東桑屬原始類型，長穗桑、奶桑屬進(jìn)化類型；

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

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

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

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

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

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

(3)每個種應(yīng)取更多樣，在population步，統(tǒng)計分析桑群體遺傳學(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] 陳嶸.中國樹木分類學(xué)[M].北京：科學(xué)出版社,1937：228-231.

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

[9] 張秀實(shí).?？菩路诸惾篬J].植物分類學(xué)報，1984，22(1)：64-76.

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

[11]曹子余.中國桑屬(?？?新分類群[J]．植物分類學(xué)報，1991，29(3)：264-267.

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

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

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

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

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

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

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

[19]陳仁芳,張澤，唐洲，等.桑屬ITS、trnL-F、rps16序列與進(jìn)化分析[J].中國農(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), 羅龍海, 等.高通量測序技術(shù)在動植物研究領(lǐng)域中的應(yīng)用[J].中國科學(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]劉艷玲.蓮野生居群遺傳多樣性評價及高密度遺傳連鎖圖譜的構(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]中國植物志編輯委員會.中國植物志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

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