基于RAD高通量測序的貴州百里杜鵑保護區杜鵑花屬分類*

黃承玲 姚 剛 田曉玲 任永權 黃家湧 馬永鵬

(1.貴州民族大學生態環境工程學院 貴陽 550025； 2.云南省極小種群野生植物綜合保護重點實驗室 中國科學院昆明植物研究所昆明 650201； 3.貴州民族大學人文科技學院 貴陽 550025； 4.貴州百里杜鵑科研所 畢節 551614)

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

百里杜鵑區域杜鵑花屬中很多種親緣關系近，自然雜交相當普遍(Zhangetal.，2007； Marczewskietal.，2016)。本研究首先剔除前期研究中證實為天然雜交的種(張長芹等，2015)，采集該區域33種杜鵑花屬植物進行高通量測序并進行系統分類，對個別物種的形態與系統分類結果不一致的原因進行探討，以期實現百里杜鵑保護區杜鵑花屬的準確分類，并為該區域杜鵑花資源的新品種選育提供指導。

1 材料與方法

1.1 試驗材料

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

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

1.2 試驗流程

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

1.3 原始數據過濾

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

1.4 基于無參考基因組的變異檢測

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

1.5 系統發生樹構建

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

1.6 群體遺傳結構fastStructure分析

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

1.7 主成分分析(PCA)

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

2 結果與分析

2.1 測序數據特性

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

2.2 群體結構分析

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

2.3 SNP位點主成分分析

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

表2 杜鵑花屬34個種RAD測序的數據特性Tab.2 Data characteristics of 34 species in Rhododendron by RAD sequencing

續表2 Continued

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

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

2.4 系統發生樹構建

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

圖3 基于28 983個SNP位點的杜鵑花屬34個種的系統發生樹Fig. 3 Results from phylogenomic analysis of 34 Rhododendron species based on 28 983 SNPs

3 討論

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

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

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

4 結論

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