張 序 張秀姣 馬永鵬 李正紅 萬(wàn)友名 馬 宏*

（1. 中國(guó)林業(yè)科學(xué)研究院資源昆蟲研究所，昆明 650233；2. 南京林業(yè)大學(xué)林學(xué)院，南京 210037；3. 中國(guó)科學(xué)院昆明植物研究所，云南省極小種群野生植物綜合保護(hù)重點(diǎn)實(shí)驗(yàn)室，昆明 650201）

寬杯杜鵑（Rhododendron sinofalconeri）系杜鵑屬（Rhododendron）常 綠 杜 鵑 亞 屬（Subgen.Hymenanthes）常綠杜鵑組（Sect.Ponticum）杯毛杜鵑亞組（Subsect.Falconera）常綠植物，分布于中國(guó)云南東南部及越南北部［1］。其花色為杜鵑屬植物中少見的黃色，鐘狀花冠簇生于枝頂，綴以褐色斑點(diǎn)狀花藥，具有較高的觀賞價(jià)值，同時(shí)也是難得的育種親本。

近年來(lái)，由于道路建設(shè)、農(nóng)田開墾和修建風(fēng)力電站等人類活動(dòng)，寬杯杜鵑野生種群遭到直接破壞，生境破碎化和人為采挖現(xiàn)象嚴(yán)重。種群年齡結(jié)構(gòu)極度失衡，絕大多數(shù)為老年個(gè)體，偶見幼苗，為典型的衰退性種群；加之寬杯杜鵑對(duì)生境的海拔（＞2 600 m）、濕度等要求高，其面臨著較大的滅絕風(fēng)險(xiǎn)?！抖霹N紅色名錄》將其收錄為易危物種［1］。根據(jù)IUCN（3.1 版）［2］標(biāo)準(zhǔn)和野外調(diào)查結(jié)果，我們將其提升為極危等級(jí)［CR B1b（iv，v）c（iii，iv）］。遺憾的是，有關(guān)寬杯杜鵑的研究?jī)H見花粉儲(chǔ)藏活性方面［3］。因此，開展寬杯杜鵑保護(hù)生物學(xué)研究已刻不容緩。

新一代測(cè)序技術(shù)（next-generation sequencing，NGS）的快速發(fā)展提高了單核苷酸多態(tài)性（singlenucleotide-polymorphism，SNP）的發(fā)現(xiàn)效率，展現(xiàn)出高通量SNP 基因分型在植物分子生物學(xué)研究的潛力?；蚍中停╣enotyping-by-sequencing，GBS）技術(shù)是一種基于NGS 技術(shù)的SNP 基因分型方法。簡(jiǎn)化的GBS 需要限制性內(nèi)切酶在測(cè)序前消化DNA，以生成簡(jiǎn)化的表示庫(kù)。由于參考基因組的可用性并不是實(shí)施GBS 的必要條件，而且成本較低、實(shí)驗(yàn)操作相對(duì)簡(jiǎn)單，高通量測(cè)序獲得的含量SNP 攜帶豐富的全基因組信息，因而已廣泛用于植物遺傳多樣性分析、系統(tǒng)進(jìn)化、分子標(biāo)記定位和遺傳圖譜構(gòu)建等研究［4～6］。

本研究以來(lái)自6 個(gè)亞種群的36 個(gè)寬杯杜鵑個(gè)體為試驗(yàn)材料，通過(guò)GBS 技術(shù)對(duì)中國(guó)現(xiàn)存種群進(jìn)行基因分型，揭示其遺傳多樣性和遺傳結(jié)構(gòu)，為今后原地和遷地保護(hù)策略和開發(fā)利用提供科學(xué)指導(dǎo)。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)材料采自云南省開遠(yuǎn)市大黑山（23°64'N，103°50'E，海拔2 800 m）與文山老君山（23°36'N，103°91'E，海拔2 900 m）寬杯杜鵑種群，隨機(jī)取樣，個(gè)體間距離＞15 m，成株幼嫩葉片保存于硅膠中。將供試寬杯杜鵑材料按照地理來(lái)源區(qū)分為2 個(gè)種群與6 個(gè)亞種群，大黑山種群（包括A，B，C這3 個(gè)亞種群，分別包含樣本D1～D5，D6～D10，D11～D13），老君山種群（包括A1，A2，A3 這3 個(gè)亞種群，分別包含樣本L1～L11，L12～L17，L18～L23）。

1.2 DNA提取和GBS文庫(kù)構(gòu)建

利用CTAB 方法提取DNA，采用PstⅠ-HF/MspI 對(duì)DNA 進(jìn)行酶切，使用Qubit 測(cè)定PCR 產(chǎn)物濃度，濃度需大于5 ng·μL-1，利用UGbS-Flex 技術(shù)［7］進(jìn)行測(cè)序文庫(kù)構(gòu)建。

1.3 單核苷酸多態(tài)性（SNP）鑒定

依據(jù)建庫(kù)樣品與條形碼接頭對(duì)應(yīng)關(guān)系拆分為單樣品原始數(shù)據(jù)，將混池下機(jī)測(cè)序得到的原始數(shù)據(jù)（raw reads）使用fastqc（v0.11.7）軟件［8］進(jìn)行質(zhì)控。使用stacks（v2.1）軟件包［9］中的process_radtags 程序（主要參數(shù)-r--renz_1--adapter_mm 1），剔除混池下機(jī)Raw Reads 含有接頭序列的Reads，并依據(jù)建庫(kù)樣品與barcode 對(duì)應(yīng)關(guān)系拆分為單樣品Raw Reads。使用fastx toolkit（v0.0.14）軟件包［10］中的fastx_trimmer 程序（主要參數(shù)-f-l），移除酶切位點(diǎn)序列以及3′端fastqc 質(zhì)控質(zhì)量分?jǐn)?shù)＜20 的所有堿基。使用bowtie2（v2.3.4.1）軟件［11］（主要參數(shù)：--maxins1000--no-discordant--no-mixed）將得到的高質(zhì)量數(shù)據(jù)軟件比對(duì)到參考基因組上（ftp：//parrot.genomics. cn/gigadb/pub/10.5524/100001_101000/100331/Genome/Rhododendron_delavayi_genome.fa.gz）。利用GATK（v3.8-1）軟件Unified 程序預(yù)測(cè)樣品中的SNP 和INDEL（Insertion and Deletion）位點(diǎn)，再通過(guò)SelectVariants 程序（主要參數(shù)-restrictAllelesTo BIALLELIC-select"QD＞10.0"）對(duì)預(yù)測(cè)結(jié)果進(jìn)行篩選，得到初步SNP 和INDEL 結(jié)果。為了降低SNP 和INDEL 檢測(cè)的錯(cuò)誤率，使用vcftools（v0.1.13）軟件對(duì)獲得SNP 分型結(jié)果進(jìn)行過(guò)濾（主要參數(shù)--maf 0.01--minDP 4--max-missing 0.8）［12］。

1.4 遺傳數(shù)據(jù)分析

利用treebest［13］（v1.9.2）軟件計(jì)算遺傳矩陣并構(gòu)建進(jìn)化樹，并通過(guò)bootstrap 法進(jìn)行檢驗(yàn)可靠性（重復(fù)1 000 次）。使用plink2（v2.0）軟件進(jìn)行PCA分析。通過(guò)PLINK（v1.9）軟件［14］、ADMIXTURE（v1.3.0）軟件進(jìn)行種群結(jié)構(gòu)分析。采用R 語(yǔ)言擴(kuò)展包genepop（v1.0.5）計(jì)算兩個(gè)寬杯杜鵑種群樣本的哈迪-溫伯格平衡（Hardy-Weinberg equilibrium，HWE）、期望雜合度（expected heterozygosity，He）、觀測(cè)雜合度（observed heterozygosity，Ho）、多態(tài)信息含量（polymorphism information content，PIC）、等位基因數(shù)（number of alleles，Na）、有效等位基因數(shù)（number of effective alleles，Ne）、F 統(tǒng)計(jì)及種群間的遺傳分化系數(shù)（Fst）、基因流（Nm）；運(yùn)用vcftools（v0.1.14）軟件進(jìn)行單位點(diǎn)計(jì)算核苷酸多樣性（Pi）。

2 結(jié)果與分析

2.1 SNP位點(diǎn)挖掘

經(jīng)GBS 測(cè)序，36 個(gè)寬杯杜鵑樣本的測(cè)序總數(shù)據(jù)量為57.82 Gb，清理低質(zhì)量序列后，得到高質(zhì)量序列數(shù)據(jù)共51.49 Gb，平均每樣本為1.43 Gb。測(cè)序質(zhì)量較高（Q20＞=95.43%，Q30＞=88.97%），GC 分布合理，種群樣本與參考基因組的平均比對(duì)率為30.86%。測(cè)序得到所有reads 對(duì)全基因組的覆蓋度為1.8%～3.3%；測(cè)序深度大于1的測(cè)序reads對(duì)全基因組的覆蓋度為1.40%～2.14%，測(cè)序深度大于3的測(cè)序reads對(duì)全基因組的覆蓋度為1.15%～1.58%。過(guò)濾后共獲得103 133個(gè)高質(zhì)量SNP位點(diǎn)。

2.2 遺傳結(jié)構(gòu)

基于SNP 信息構(gòu)建的系統(tǒng)進(jìn)化樹分析表明，兩個(gè)種群可按地理分布劃分為兩大類群：Group1（老君山種群）與Group2（大黑山種群），其中L21與老君山其他個(gè)體的遺傳距離較遠(yuǎn)（見圖2）。

主成分分析（PCA）結(jié)果表明，處于主導(dǎo)地位的第一主成分（PC1）與第二主成分（PC2）貢獻(xiàn)率分別為52.62%和12.32%（圖3）。二維聚類結(jié)果表明，寬杯杜鵑2 個(gè)種群中的亞種群間隔較小，根據(jù)PC1 可將總樣本分為2 個(gè)類群，分別為Group1（老君山種群）和Group2（大黑山種群）（見圖3）。

對(duì)36 份寬杯杜鵑材料進(jìn)行Structure 分析，根據(jù)CV（Cross validation error）確定最優(yōu)K 值為2，基于此將試驗(yàn)材料劃分為2 個(gè)類群。Group1 包含老君山種群的13 個(gè)個(gè)體，Group2 包含大黑山種群的23個(gè)個(gè)體。

2.3 遺傳多樣性

2.3.1 種群內(nèi)遺傳多樣性

各種群的遺傳多樣性結(jié)果表明，大黑山種群（A，B，C）、老君山種群（A1，B1，C1）中SNP位點(diǎn)的Na值介于1.454 3～1.682 5，Ne值介于1.291 0～1.319 8，Ho值介于0.182 1～0.194 8，He值介于0.170 9～0.196 1，PIC 值 介 于0.136 8～0.160 4，Pi值 介 于0.205 1～0.207 9。以上指標(biāo)表明各亞種群內(nèi)具有較低的遺傳多樣性（PIC＜0.25）（見表1）?？ǚ綑z驗(yàn)計(jì)算出6個(gè)亞種群的HWE 無(wú)偏估計(jì)的概率值（p）均大于0.05，表明所調(diào)查亞種群均達(dá)到遺傳平衡。

2.3.2 種群間遺傳多樣性

結(jié)果顯示，6 個(gè)亞種群間的Fst值為0～0.1822，大黑山種群與老君山種群之間遺傳分化較大（Fst＞0.15），而各自種群內(nèi)的亞種群間分化較低（Fst＜0.05）（見表2）。兩種群內(nèi)部各亞種群間基因流豐富（Nm＞3），兩種群間平均基因流Nm為1.1674（見表3）。

3 討論

豐富的遺傳多樣性是物種應(yīng)對(duì)環(huán)境變化的基礎(chǔ)［15～16］，本研究發(fā)現(xiàn)寬杯杜鵑老君山與大黑山兩個(gè)種群具有較低的遺傳多樣性（He=0.185 6）。近年研究表明多數(shù)杜鵑屬植物具有較高的遺傳多樣性［17～19］，即使種群中個(gè)體數(shù)量稀少的物種，如長(zhǎng)梗杜鵑與大樹杜鵑，其原因可能是許多殘存于遠(yuǎn)距離孤立種群中的瀕危植物，它們可能保留了從前祖先廣泛連續(xù)分布時(shí)的基因信息［20～21］。較低的遺傳多樣性也見于華頂杜鵑（Rh.huadingense）與Rh.ferrugineum，生境破碎化致使華頂杜鵑種群規(guī)模逐步縮小，從而引起近交衰退，致使其遺傳多樣性喪失［22］；Bruni 等分析發(fā)現(xiàn)無(wú)性繁殖和近親繁殖是Rh.ferrugineum 發(fā)生遺傳漂變的誘因［23］。寬杯杜鵑較低的遺傳多樣性與其較小的種群規(guī)模且生境遭到嚴(yán)重人為破壞，最終導(dǎo)致生境破碎化和遺傳漂變有關(guān)。繁育系統(tǒng)通常被認(rèn)為是影響物種遺傳多樣性的重要因素之一，自交與異交類型繁育系統(tǒng)在杜鵑屬植物中均有記錄，且以往的研究表明，杜鵑屬植物缺乏典型的自交不親和系統(tǒng)，表明杜鵑存在復(fù)雜的交配系統(tǒng)［24］。因此，繁育系統(tǒng)可能對(duì)寬杯杜鵑遺傳多樣性具有的重要影響。

表1 不同種群寬杯杜鵑遺傳多樣性度量指標(biāo)Table 1 Genetic diversity indexes of R.sinofalconeri of different groups

表2 種群間的遺傳分化系數(shù)和遺傳距離表Table 2 Genetic differentiation coefficient（Fst：above diagonal）and genetic distance（DR：below diagonal）between populations

表3 種群間的基因流（Nm）表Table 3 Gene flow between populations（Nm）

遺傳結(jié)構(gòu)受多種因素的影響，如繁育系統(tǒng)、遺傳漂變、種群大小、種子傳播、基因流、進(jìn)化史和自然選擇等［25］，種群間的遺傳分化系數(shù)（Fst）分析表明，36 個(gè)個(gè)體在大黑山與老君山兩種群間分化程度較高（Fst＞0.15），該結(jié)果與進(jìn)化樹、主成分分析與種群遺傳結(jié)構(gòu)分析結(jié)果一致，均表明大黑山與老君山兩寬杯杜鵑種群的遺傳距離較遠(yuǎn)。Hamrick等認(rèn)為基因流（Nm）＞1可一定程度抵消遺傳漂變帶來(lái)的負(fù)面影響［26］。雖然兩種群間平均基因流（Nm）為1.167 4，但兩種群間呈高度分化（Fst＞0.15），這可能是由于寬杯杜鵑是多年生喬木，現(xiàn)有的基因流可能源于其祖先種群間的遺傳交流或者共享祖先種群的某些基因型所致。以往的研究發(fā)現(xiàn)，杜鵑屬植物種子通過(guò)風(fēng)的散布距離為30～80 m，而通過(guò)昆蟲與鳥類傳播的杜鵑花粉移動(dòng)距離通常介于3～10 km［27］，但兩個(gè)寬杯杜鵑種群間較遠(yuǎn)的距離（＞80 km）使得花粉與種子的傳播交流難以實(shí)現(xiàn)。長(zhǎng)距離地理隔離和特殊的生境在將來(lái)會(huì)進(jìn)一步阻礙這些高度隔離的殘存種群之間的基因交換，其潛在的負(fù)面后果是減少剩余的小種群之間的基因流動(dòng)和增加遺傳漂變的風(fēng)險(xiǎn)［28］。

就地保護(hù)是保護(hù)瀕危物種的有效措施［29］。現(xiàn)存寬杯杜鵑種群的生境在近年來(lái)皆受到道路建設(shè)等人為破壞的影響，而生物多樣性對(duì)道路建設(shè)等的滯后反應(yīng)，使得在短期內(nèi)難以完全認(rèn)知該干擾的全部影響［30～31］，已有研究發(fā)現(xiàn)杜鵑花根際真菌多樣性可能由于人類活動(dòng)而發(fā)生變化［32］。大黑山種群因盜挖致使種群中僅剩老樹與少量幼苗。王書珍等通過(guò)SSR 分析發(fā)現(xiàn)映山紅（Rh.simsii）老樹與幼苗種群遺傳多樣性最低，以小樹種群多樣性最為豐富［33］。因此，保護(hù)寬杯杜鵑的自然生境，將其并入保護(hù)區(qū)或者建立保護(hù)小區(qū)應(yīng)作為優(yōu)先考慮的策略?；蛄魇怯绊懳锓N遺傳多樣性與遺傳分化的關(guān)鍵因素［34］，我們發(fā)現(xiàn)寬杯杜鵑2個(gè)種群間遺傳分化程度較高，加之地理隔離使其面臨極高的遺傳漂變風(fēng)險(xiǎn)。遷地保護(hù)是保護(hù)野生物種的重要途徑，遷地引種時(shí)應(yīng)盡量涵蓋不同單倍型（Haplotype）和進(jìn)化顯著單元（evolutinary significant units，ESUs）的個(gè)體，同時(shí)，在回歸引種時(shí)注意6 個(gè)亞種群間的相互引種，以增加其基因交流和遺傳多樣性，盡可能抵御遺傳漂變產(chǎn)生的負(fù)面影響，降低野外滅絕的風(fēng)險(xiǎn)。