基于RAD高通量測(cè)序的貴州百里杜鵑保護(hù)區(qū)杜鵑花屬分類*

黃承玲 姚 剛 田曉玲 任永權(quán) 黃家湧 馬永鵬

(1.貴州民族大學(xué)生態(tài)環(huán)境工程學(xué)院 貴陽(yáng) 550025； 2.云南省極小種群野生植物綜合保護(hù)重點(diǎn)實(shí)驗(yàn)室 中國(guó)科學(xué)院昆明植物研究所昆明 650201； 3.貴州民族大學(xué)人文科技學(xué)院 貴陽(yáng) 550025； 4.貴州百里杜鵑科研所 畢節(jié) 551614)

貴州省百里杜鵑保護(hù)區(qū)位于貴州西北部，其自然分布的杜鵑花屬(Rhododendron)植物資源十分豐富，是國(guó)內(nèi)為數(shù)不多的以高山杜鵑為觀賞和保護(hù)對(duì)象的省級(jí)自然保護(hù)區(qū)、省級(jí)風(fēng)景名勝區(qū)和國(guó)家級(jí)森林公園。近年來(lái)，百里杜鵑不僅成為熱門(mén)的觀花旅游目的地，還因其獨(dú)特的杜鵑花資源成為科學(xué)研究的熱點(diǎn)。楊成華等(2006)、Chen等(2010a)、張長(zhǎng)芹等(2015)等都對(duì)百里杜鵑區(qū)域的杜鵑花資源進(jìn)行過(guò)詳細(xì)的調(diào)查研究，但具體到個(gè)別物種，尤其是近年發(fā)表的該區(qū)域的新種的分類地位問(wèn)題仍存在爭(zhēng)議(Marczewskietal.，2016)。以上研究對(duì)百里杜鵑保護(hù)區(qū)的杜鵑花屬僅是依據(jù)形態(tài)學(xué)特性進(jìn)行分類。比如，花芽的著生位置和葉背面是否有鱗片是區(qū)分亞屬的重要形態(tài)特征； 而植株、花冠形狀與顏色，以及葉片、花瓣、雌雄蕊毛的類型(比如：剛毛、絨毛、腺體毛、星狀毛、杯狀毛)是區(qū)分不同亞組和物種的重要特征(胡琳貞等，1991)。然而，形態(tài)學(xué)的一個(gè)突出問(wèn)題是很難把握種間和種內(nèi)形態(tài)變異的幅度(Marczewskietal.，2016)。因而，隨著分子生物學(xué)技術(shù)的發(fā)展，有必要借助分子標(biāo)記技術(shù)對(duì)百里杜鵑保護(hù)區(qū)杜鵑花屬的分類及存疑進(jìn)行印證和解答。目前對(duì)百里杜鵑的杜鵑花屬系統(tǒng)分類的分子水平研究較少。何云松等(2015)采用葉綠體基因組的trnL-F間隔區(qū)測(cè)序，研究了百里杜鵑保護(hù)區(qū)26種杜鵑花屬植物的系統(tǒng)發(fā)育關(guān)系，其結(jié)果顯示trnL-F可以在亞屬水平區(qū)分物種，而在物種水平僅有少數(shù)種被明確區(qū)分，說(shuō)明杜鵑花屬分類的復(fù)雜與困難。近年來(lái)第2代測(cè)序技術(shù)(高通量測(cè)序技術(shù))因其單次運(yùn)行產(chǎn)出的序列數(shù)據(jù)量大、價(jià)格低廉而廣泛用于植物的遺傳轉(zhuǎn)化、分子育種研究中(時(shí)小東等，2016； 胡玉玲等，2014； 鄧敏捷等，2013)。在第2代測(cè)序技術(shù)基礎(chǔ)上，一項(xiàng)基于全基因組酶切位點(diǎn)的簡(jiǎn)化基因組測(cè)序技術(shù)RAD-seq(restriction-site associated DNA-sequencing)迅速發(fā)展起來(lái)，該技術(shù)通過(guò)對(duì)限制性酶切位點(diǎn)的標(biāo)記進(jìn)行測(cè)序，在大多數(shù)生物中獲得海量的單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)位點(diǎn)(Milleretal.，2007； Tasselletal.，2008； 桂柳柳，2017)，因具有多態(tài)性位點(diǎn)數(shù)量多、分布密度高且操作方便等特點(diǎn)逐漸被應(yīng)用到植物群體遺傳學(xué)研究中(魏明科等，2018； 王洋坤等，2014)。李云飛等(2019)采用RAD-seq技術(shù)開(kāi)發(fā)了85種杜鵑花屬植物的高質(zhì)量SNP位點(diǎn)，通過(guò)高通量測(cè)序探討分類，在亞屬水平的分類取得很好的效果。

百里杜鵑區(qū)域杜鵑花屬中很多種親緣關(guān)系近，自然雜交相當(dāng)普遍(Zhangetal.，2007； Marczewskietal.，2016)。本研究首先剔除前期研究中證實(shí)為天然雜交的種(張長(zhǎng)芹等，2015)，采集該區(qū)域33種杜鵑花屬植物進(jìn)行高通量測(cè)序并進(jìn)行系統(tǒng)分類，對(duì)個(gè)別物種的形態(tài)與系統(tǒng)分類結(jié)果不一致的原因進(jìn)行探討，以期實(shí)現(xiàn)百里杜鵑保護(hù)區(qū)杜鵑花屬的準(zhǔn)確分類，并為該區(qū)域杜鵑花資源的新品種選育提供指導(dǎo)。

1 材料與方法

1.1 試驗(yàn)材料

用于RAD-seq技術(shù)測(cè)序的杜鵑花屬植物樣本共34份，采樣信息見(jiàn)表1，隸屬于6個(gè)亞屬7個(gè)組10個(gè)亞組(表2)，其中大果杜鵑(Rhododendronglanduliferum)采自貴州盤(pán)縣，其余33種按張長(zhǎng)芹等(2015)的描述和記錄，采自貴州省百里杜鵑自然保護(hù)區(qū)。每個(gè)種采集5株成年個(gè)體的健康幼嫩葉片2～3片，采樣的個(gè)體間距10 m以上，采集的葉片保存于4 ℃冰箱備用。

表1 杜鵑花屬34個(gè)種的采樣信息Tab.1 Sample information of 34 species in Rhododendron

1.2 試驗(yàn)流程

每個(gè)個(gè)體摘取新鮮葉片0.5 g，用改良的十六烷基三甲基溴化銨(hexadecyl trimethyl ammonium bromide，CTAB)法分別提取樣品的基因組DNA并進(jìn)行DNA質(zhì)量檢測(cè)(Doyle，1987)，對(duì)質(zhì)量合格的DNA，參考Patterson 等(2012)流程進(jìn)行建庫(kù)； 用Illumina HiSeq 4000 System(San Diego，CA，USA)測(cè)序平臺(tái)進(jìn)行測(cè)序。

1.3 原始數(shù)據(jù)過(guò)濾

原始測(cè)序數(shù)據(jù)(raw data)中一般會(huì)有少部分的reads帶有測(cè)序引物、接頭等人工序列，為保證后續(xù)分析的準(zhǔn)確性，需要先對(duì)raw data進(jìn)行過(guò)濾，去除影響數(shù)據(jù)質(zhì)量和后續(xù)分析的低質(zhì)量區(qū)域(王成賀，2016)。本研究采用Stacks 2.5.2軟件的process_radtags程序(Rochetteetal.，2017； 2019)，按照barcode 序列和index序列對(duì)個(gè)體進(jìn)行區(qū)分并作初步的質(zhì)量控制，去除含有接頭及低質(zhì)量的序列，得到clean data，將所有得到的reads剪切為135 bp，然后通過(guò)FastQC 0.11.8 (Andrews, 2015) 軟件進(jìn)行個(gè)體reads的質(zhì)量評(píng)估。

1.4 基于無(wú)參考基因組的變異檢測(cè)

采用Stacks-2.5.2軟件進(jìn)行數(shù)據(jù)預(yù)處理，包括denovo組裝和不連鎖變異位點(diǎn)的過(guò)濾。Stacks分析流程如下：1)使用Linux系統(tǒng)的cat命令合并雙端數(shù)據(jù)。2)采用ustacks根據(jù)序列間相似性原則和最大似然模型對(duì)每個(gè)樣本進(jìn)行denovo組裝，將完全一致的序列聚為一個(gè)stacks，生成的結(jié)果為loci。具體使用的關(guān)鍵參數(shù)為：-M 5 -m 3 -R。3)使用cstacks對(duì)所有樣本的loci進(jìn)行合并，采用kmer的方法來(lái)進(jìn)行比對(duì)，得到catalog.allels.tsv.gz，catalog.snps.tsv.gz和catalog.tags.tsv.gz文件和log文件。具體使用的關(guān)鍵參數(shù)為：-n 3。4)采用sstacks把所有樣本的序列match到cstacks生成的catalog文件中，得到各個(gè)樣本的*.matches.tsv.gz文件和log文件。此步驟采用默認(rèn)參數(shù)。5)使用tsv2bam把sstacks生成的各樣本的matches.tsv.gz文件轉(zhuǎn)為各locus的bam文件，得到各個(gè)樣本的*.matches.bam文件和log文件。此步驟采用默認(rèn)參數(shù)。6)最后使用populations檢測(cè)所有樣本的SNPs，此步驟具體使用的關(guān)鍵參數(shù)為：--write-random-snp --vcf -plink，得到包含SNP的VCF格式文件。

1.5 系統(tǒng)發(fā)生樹(shù)構(gòu)建

首先采用VCFtools(Daneceketal.，2011)對(duì)上一步stacks得到的SNP進(jìn)行過(guò)濾，所使用的關(guān)鍵參數(shù)包括：--min-alleles 2 --max-alleles 2 --remove-indels --maf 0.05 --max-missing 0.2，--minGQ 30 --minDP 3。得到42 083個(gè)SNP位點(diǎn)，對(duì)這些位點(diǎn)進(jìn)行中性檢驗(yàn)，所使用關(guān)鍵參數(shù)為--TajimaD 2500，參照Tajima’sD的95%置信區(qū)間表(Tajima，1989)篩選中性位點(diǎn)。在得到中性位點(diǎn)的VCF文件后，使用Tassel 5.0(Bradburyetal.，2007)將其轉(zhuǎn)化為phy格式。對(duì)phy格式的SNP矩陣采用concatenation的最大似然法(Maximum Likelihood)進(jìn)行系統(tǒng)發(fā)生分析，使用IQTREE(Minhetal.，2020； Kalyaanamoorthyetal.，2017； Hoangetal.，2018) 在北京超級(jí)云計(jì)算中心A分區(qū)上進(jìn)行運(yùn)算，從242個(gè)DNA堿基替換模型中篩選出的最佳堿基替換模型為T(mén)VM+F+ASC，以超快bootstrap法評(píng)估系統(tǒng)發(fā)生樹(shù)分支的支持度，指定次數(shù)為1 000。對(duì)于IQTREE得到的系統(tǒng)發(fā)生樹(shù)文件，使用FigTree(Priceetal.，2010)進(jìn)行查看和美化處理。

1.6 群體遺傳結(jié)構(gòu)fastStructure分析

通過(guò)plink轉(zhuǎn)換為二進(jìn)制文件(得到.ped和.map文件)，然后以轉(zhuǎn)化好格式的SNP文件作為輸入文件，設(shè)置1≤K≤10，對(duì)fastStructure(Aniletal.，2014)使用for-do語(yǔ)句進(jìn)行群體聚類分析，然后使用grep提取最佳K值，即cross-validation error值最低時(shí)所對(duì)應(yīng)的K值。設(shè)置個(gè)體Q≥0.9作為判斷物種的標(biāo)準(zhǔn)。

1.7 主成分分析(PCA)

主成分分析是群體遺傳學(xué)中常用的分析手段，目前用于主成分分析的數(shù)據(jù)主要為高密度的SNP標(biāo)記，常使用的分析軟件為GCTA(Jianetal.，2011)，通過(guò)對(duì)獲得的中性SNP位點(diǎn)的數(shù)據(jù)集進(jìn)行矩陣轉(zhuǎn)換，利用eigen函數(shù)計(jì)算所有樣本遺傳關(guān)系矩陣的特征值和特征向量，再使用R語(yǔ)言的ggplot2包進(jìn)行分析結(jié)果的可視化，并根據(jù)檢索得到的樣品亞屬信息進(jìn)行分類。

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)特性

對(duì)杜鵑花屬34個(gè)種的樣本進(jìn)行RAD-seq技術(shù)測(cè)序，共獲得測(cè)序數(shù)據(jù)量15.1 Gb，平均每個(gè)個(gè)體455.58 Mb。使用fastp軟件對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)量控制，得到10.22 Gb的清洗數(shù)據(jù)，平均每個(gè)個(gè)體307.80 Mb。另外，數(shù)據(jù)的Q20、Q30也分別從過(guò)濾前的94.44%、87.11%提高到過(guò)濾后的97.43%、94.32%。數(shù)據(jù)產(chǎn)出、測(cè)序深度、覆蓋度和多態(tài)性位點(diǎn)的具體信息見(jiàn)表2。

2.2 群體結(jié)構(gòu)分析

混群分析結(jié)果顯示，當(dāng)K=4時(shí)，曲線出現(xiàn)最低值(交叉驗(yàn)證錯(cuò)誤率cross-validation error=0.013 332 7)； 同時(shí)，考慮到樣本所歸屬的分類單元，K=3(cross-validation error=0.030 077 3)和K=5(cross-validation error=0.033 354 2)也做了群體結(jié)構(gòu)柱狀堆疊圖(圖1)。由圖1可知：當(dāng)K=3時(shí)，藍(lán)色、綠色和紅色代表3類不同的遺傳組成，其中藍(lán)色對(duì)應(yīng)映山紅亞屬的全部，綠色對(duì)應(yīng)常綠杜鵑亞屬的大部分，紅色則包含全部的杜鵑亞屬、糙葉杜鵑亞屬、馬銀花亞屬和羊躑躅亞屬，并包含部分常綠杜鵑亞屬個(gè)體； 當(dāng)K=4時(shí)，馬銀花亞屬與羊躑躅亞屬和有鱗類杜鵑(杜鵑亞屬和糙葉杜鵑亞屬)分開(kāi)，而常綠杜鵑亞屬所有種也與有鱗類杜鵑完全分開(kāi)； 當(dāng)K=5時(shí)，分辨率降低，馬銀花亞屬、羊躑躅亞屬和有鱗類杜鵑無(wú)法區(qū)分，而常綠杜鵑亞屬內(nèi)部分為3個(gè)部分，但與映山紅亞屬無(wú)法區(qū)分。綜合來(lái)看，K=4時(shí)可將5個(gè)亞屬中的4個(gè)亞屬很好地區(qū)分，且支持形態(tài)證據(jù)對(duì)這些種的劃分； 而常綠杜鵑亞屬內(nèi)部分為2組反映了該亞屬?gòu)?fù)雜的系統(tǒng)進(jìn)化關(guān)系，映山紅亞屬無(wú)法與部分常綠杜鵑亞屬區(qū)分說(shuō)明二者具有較近的親緣關(guān)系。

2.3 SNP位點(diǎn)主成分分析

PCA 聚類結(jié)果(圖2)與fastStructure結(jié)果(圖1)相似，34個(gè)杜鵑花屬植物樣本聚成5類，即常綠杜鵑亞屬、羊躑躅亞屬、馬銀花亞屬和映山紅亞屬的各自物種聚在一起，與基于形態(tài)性狀的劃分完全吻合； 糙葉杜鵑亞屬與杜鵑亞屬聚在一起，反映了2個(gè)亞屬親緣關(guān)系較近。

表2 杜鵑花屬34個(gè)種RAD測(cè)序的數(shù)據(jù)特性Tab.2 Data characteristics of 34 species in Rhododendron by RAD sequencing

續(xù)表2 Continued

圖1 基于28 983個(gè)SNP位點(diǎn)的杜鵑花屬34個(gè)種的fastStructure聚類Fig. 1 Results from fastStructure analysis of 34 Rhododendron species based on 28 983 SNPs

圖2 基于28 983個(gè)SNP位點(diǎn)的杜鵑花屬34個(gè)種的PCA聚類Fig. 2 Results from PCA analysis of 34 Rhododendron species based on 28 983 SNPs

2.4 系統(tǒng)發(fā)生樹(shù)構(gòu)建

采用IQ-Tree軟件對(duì)杜鵑花屬34個(gè)種的28 983個(gè)單核苷酸多態(tài)性位點(diǎn)(SNPs)構(gòu)建系統(tǒng)發(fā)生樹(shù)，使用篩選好的模型TVM+F+ASC進(jìn)行最大似然親緣關(guān)系分析(maximum-likelihood tree)，結(jié)果(圖3)與PCA主成分分析(圖2)相似，即34個(gè)種的樣本聚成5支，其中4支分別對(duì)應(yīng)形態(tài)上其各自歸屬的常綠杜鵑亞屬、羊躑躅亞屬、馬銀花亞屬和映山紅亞屬，同樣，糙葉杜鵑亞屬與杜鵑亞屬聚成1支，反映了2個(gè)亞屬親緣關(guān)系較近。一些形態(tài)極其相似的種無(wú)法很好地區(qū)分開(kāi)來(lái)，比如九龍山杜鵑和大果杜鵑，前者大部分性狀如幼枝被剛毛狀腺體、葉形狀與大小、葉柄有腺體等與《中國(guó)植物志》中大果杜鵑的形態(tài)描述都有重疊(胡琳貞等，1991)； 而繁花杜鵑與皺葉杜鵑這2個(gè)種的差別僅在于葉背面毛被的顏色。

圖3 基于28 983個(gè)SNP位點(diǎn)的杜鵑花屬34個(gè)種的系統(tǒng)發(fā)生樹(shù)Fig. 3 Results from phylogenomic analysis of 34 Rhododendron species based on 28 983 SNPs

3 討論

通過(guò)對(duì)百里杜鵑保護(hù)區(qū)杜鵑花屬34個(gè)種的簡(jiǎn)化基因組測(cè)序，以獲得的高質(zhì)量SNP位點(diǎn)為基礎(chǔ)，從基因組水平系統(tǒng)分析了34個(gè)種的系統(tǒng)演化關(guān)系。從fastStructure、PCA聚類和系統(tǒng)發(fā)生樹(shù)的結(jié)果來(lái)看，常綠杜鵑亞屬可以與其他5個(gè)亞屬完全分開(kāi)獨(dú)立為一支，支持率92%，這與Kurashige等(2001)、何云松等(2015)、李云飛等(2019)的結(jié)論一致，表明常綠杜鵑亞屬是一個(gè)單系類群。而葉片無(wú)鱗的馬銀花亞屬和映山紅亞屬聚成一個(gè)分支，表明二者親緣關(guān)系較近； 葉片同樣無(wú)鱗的羊躑躅亞屬獨(dú)立成一支。糙葉杜鵑亞屬和杜鵑亞屬不能完全分開(kāi)，這與Kurashige等(2001)、李云飛等(2019)的結(jié)論一致，這2個(gè)亞屬屬于葉片有鱗類群，印證了李云飛等(2019)在亞屬水平可以用葉片被鱗與否來(lái)解釋基因組水平各亞屬的系統(tǒng)分類學(xué)關(guān)系的結(jié)論。

從組、亞組水平看，常綠杜鵑亞屬包含的各個(gè)亞組的多個(gè)種交錯(cuò)組成3個(gè)分支，這個(gè)結(jié)果與Kurashige 等(2001)、何云松等(2015)、李云飛等(2019)的研究一致，表明各亞組之間的界限模糊； 對(duì)該亞屬不同亞組植物的天然雜交現(xiàn)象進(jìn)行研究(Zhangetal.，2007； Zhaetal.，2010； Yanetal.，2013； Maetal.，2016； 2010； Marczewskietal.，2016)，結(jié)果表明分布在同一個(gè)區(qū)域內(nèi)的常綠杜鵑亞屬不同亞組之間的自然雜交極為頻繁，很多種類出現(xiàn)雜交漸滲、基因水平轉(zhuǎn)移以及不完全譜系分化等現(xiàn)象，使得常綠杜鵑亞屬的亞組以下物種水平的種間關(guān)系很難準(zhǔn)確界定(Milneetal.，2010； Yanetal.，2015)。

從物種水平看，張長(zhǎng)芹等(2015)認(rèn)為Chen等(2010b)發(fā)表的新種九龍山杜鵑與大果杜鵑相同，從系統(tǒng)發(fā)生樹(shù)看，大果杜鵑與九龍山杜鵑聚成一小支，支持率達(dá)100%，支持大果杜鵑與九龍山杜鵑為同一個(gè)種的推測(cè)。張長(zhǎng)芹等(2015)認(rèn)為百里杜鵑的一些新種的分類地位有待商榷，從黃坪杜鵑、金波杜鵑、枇杷葉杜鵑、匙葉杜鵑等形態(tài)特征推測(cè)這些新種可能是馬纓杜鵑與露珠杜鵑、大白杜鵑或皺葉杜鵑的天然雜交后代； 從系統(tǒng)發(fā)生樹(shù)可以看出，金波杜鵑與馬纓杜鵑聚成一小支，支持率100%，說(shuō)明金波杜鵑與馬纓杜鵑有很近的親緣關(guān)系，普底杜鵑與皺葉杜鵑聚成一個(gè)分支，同樣表明二者有很近的親緣關(guān)系。雖然不能由此明確金波杜鵑和普底杜鵑是否為雜交后代，但至少印證了二者分別與馬纓杜鵑和皺葉杜鵑有親緣關(guān)系的推測(cè)。另外，用RAD序列得到的SNP仍然無(wú)法區(qū)分某些更小分類單元的物種，比如，有鱗大花杜鵑亞組的樹(shù)楓杜鵑和百合花杜鵑，以及銀葉杜鵑亞組的銀葉杜鵑和云錦杜鵑亞組的美容杜鵑等。目前，比較認(rèn)可的觀點(diǎn)是通過(guò)分子標(biāo)記區(qū)分物種間的親緣關(guān)系本質(zhì)上與形態(tài)判斷同等重要，不能顧此失彼。對(duì)于一些形態(tài)差異明顯但系統(tǒng)進(jìn)化分析無(wú)差別的種類，本文不做定論。期待后續(xù)通過(guò)更多、更準(zhǔn)確的分子標(biāo)記(比如全基因組測(cè)序技術(shù))，尋找這些近緣種在整個(gè)基因組水平的位點(diǎn)與結(jié)構(gòu)變異； 用更科學(xué)的系統(tǒng)進(jìn)化樹(shù)構(gòu)建方法進(jìn)一步區(qū)分和驗(yàn)證。

4 結(jié)論

采用RAD-seq技術(shù)對(duì)百里杜鵑保護(hù)區(qū)杜鵑花屬34個(gè)種進(jìn)行測(cè)序，通過(guò)denovo組裝并獲得高質(zhì)量SNP位點(diǎn)，基于這些SNP位點(diǎn)對(duì)杜鵑花屬植物進(jìn)行群體結(jié)構(gòu)分析、主成分分析與系統(tǒng)樹(shù)構(gòu)建，結(jié)果表明RAD-seq技術(shù)產(chǎn)生的大量SNP位點(diǎn)可在亞屬水平和物種水平區(qū)分百里杜鵑的大部分杜鵑花屬植物，同時(shí)證實(shí)RAD-seq技術(shù)在復(fù)雜植物類群的物種分類方面相比傳統(tǒng)分子標(biāo)記具有明顯優(yōu)勢(shì)。杜鵑花屬6個(gè)亞屬中的常綠杜鵑亞屬、馬銀花亞屬、羊躑躅亞屬和映山紅亞屬能明顯分開(kāi)，而糙葉杜鵑亞屬與杜鵑亞屬不能區(qū)分開(kāi)，表明二者親緣關(guān)系較近。對(duì)百里杜鵑保護(hù)區(qū)新發(fā)表的金波杜鵑和普底杜鵑的分類定位進(jìn)行探討，表明金波杜鵑與馬纓杜鵑、普底杜鵑與皺葉杜鵑有較近的親緣關(guān)系，但還需進(jìn)一步證實(shí)。本研究結(jié)果可為復(fù)雜植物類群的物種界定在方法上提供參考，也可為百里杜鵑保護(hù)區(qū)杜鵑花的資源開(kāi)發(fā)與新品種培育提供指導(dǎo)。