基于GBS簡化基因組技術的寬杯杜鵑遺傳多樣性分析

張 序 張秀姣 馬永鵬 李正紅 萬友名 馬 宏*

（1. 中國林業科學研究院資源昆蟲研究所，昆明 650233；2. 南京林業大學林學院，南京 210037；3. 中國科學院昆明植物研究所，云南省極小種群野生植物綜合保護重點實驗室，昆明 650201）

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

近年來，由于道路建設、農田開墾和修建風力電站等人類活動，寬杯杜鵑野生種群遭到直接破壞，生境破碎化和人為采挖現象嚴重。種群年齡結構極度失衡，絕大多數為老年個體，偶見幼苗，為典型的衰退性種群；加之寬杯杜鵑對生境的海拔（＞2 600 m）、濕度等要求高，其面臨著較大的滅絕風險。《杜鵑紅色名錄》將其收錄為易危物種［1］。根據IUCN（3.1 版）［2］標準和野外調查結果，我們將其提升為極危等級［CR B1b（iv，v）c（iii，iv）］。遺憾的是，有關寬杯杜鵑的研究僅見花粉儲藏活性方面［3］。因此，開展寬杯杜鵑保護生物學研究已刻不容緩。

新一代測序技術（next-generation sequencing，NGS）的快速發展提高了單核苷酸多態性（singlenucleotide-polymorphism，SNP）的發現效率，展現出高通量SNP 基因分型在植物分子生物學研究的潛力。基因分型（genotyping-by-sequencing，GBS）技術是一種基于NGS 技術的SNP 基因分型方法。簡化的GBS 需要限制性內切酶在測序前消化DNA，以生成簡化的表示庫。由于參考基因組的可用性并不是實施GBS 的必要條件，而且成本較低、實驗操作相對簡單，高通量測序獲得的含量SNP 攜帶豐富的全基因組信息，因而已廣泛用于植物遺傳多樣性分析、系統進化、分子標記定位和遺傳圖譜構建等研究［4～6］。

本研究以來自6 個亞種群的36 個寬杯杜鵑個體為試驗材料，通過GBS 技術對中國現存種群進行基因分型，揭示其遺傳多樣性和遺傳結構，為今后原地和遷地保護策略和開發利用提供科學指導。

1 材料與方法

1.1 試驗材料

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

1.2 DNA提取和GBS文庫構建

利用CTAB 方法提取DNA，采用PstⅠ-HF/MspI 對DNA 進行酶切，使用Qubit 測定PCR 產物濃度，濃度需大于5 ng·μL-1，利用UGbS-Flex 技術［7］進行測序文庫構建。

1.3 單核苷酸多態性（SNP）鑒定

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

1.4 遺傳數據分析

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

2 結果與分析

2.1 SNP位點挖掘

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

2.2 遺傳結構

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

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

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

2.3 遺傳多樣性

2.3.1 種群內遺傳多樣性

各種群的遺傳多樣性結果表明，大黑山種群（A，B，C）、老君山種群（A1，B1，C1）中SNP位點的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。以上指標表明各亞種群內具有較低的遺傳多樣性（PIC＜0.25）（見表1）。卡方檢驗計算出6個亞種群的HWE 無偏估計的概率值（p）均大于0.05，表明所調查亞種群均達到遺傳平衡。

2.3.2 種群間遺傳多樣性

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

3 討論

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

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

表2 種群間的遺傳分化系數和遺傳距離表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）

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

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