基于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 个亚种群间的相互引种，以增加其基因交流和遗传多样性，尽可能抵御遗传漂变产生的负面影响，降低野外灭绝的风险。