基于SNP,分子标记的波叶杜鹃遗传特征分析

曹毓蓉，马永鹏，张秀姣，刘雄芳，刘德团，张 垚，李正红，马 宏,4**

(1.中国林业科学研究院 高原林业研究所，云南 昆明 650224；
2.南京林业大学 林学院，江苏 南京 210037；
3.云南省极小种群野生植物综合保护重点实验室，中国科学院 昆明植物研究所，云南 昆明 650201；
4.国家林业和草原局 资源昆虫培育与利用重点实验室，云南 昆明 650224)

波叶杜鹃（Rhododendron hemsleyanum）为杜鹃属（Rhododendron）常绿植物，是中国特有物种，残存于四川省峨眉山海拔1 300 m 以下的中山湿性常绿阔叶林陡壁悬崖边.该物种大型白色钟状花冠簇生于枝顶，叶片长椭圆形，先端纯尖，基部耳状心形，边缘波状，故名波叶杜鹃，是极具开发潜力的观叶观花植物.通过对标本记录地点及其它可能分布区的全面调查后发现，现存个体数仅163 株，居群数量少（远低于稳定存活界限的最小生存种群）.其生境严重破碎化，分布地附近大多分布有各种人为建筑，有的甚至开垦为菜地，受外界环境胁迫状况日益严峻.目前，波叶杜鹃在几个历史记录点均已消失，面临着极高的灭绝风险.因此，波叶杜鹃是典型的极小种群野生植物（plant species with extremely small populations,PSESP）[1-2]，已被The Red List of Rhododendrons 和《中国高等植物受威胁物种名录》列为极危物种（Critically Endangered,CR）[3-4]，开展波叶杜鹃拯救性保护研究迫在眉睫.

研究物种遗传多样性和遗传结构，是科学合理地利用种质资源和改进遗传特性的理论基础，也是制订有效保护措施的前提条件[5-6].随着分子生物学技术的快速发展和不断完善，分子标记因其不受发育时期及环境因素的制约，被普遍应用到品种鉴定、群体遗传多样性分析、分子遗传育种以及进化研究等方面[7-9].其中，单核苷酸多态性（single nucleotide polymorphism,SNP）分子标记具有多态性位点数量多、分布范围广、遗传稳定性强、测序通量高等特点[10-11]，已应用于蓝果杜鹃（R.cyanocarpum）[12]、宽杯杜鹃（R.sinofalconeri）[13]等杜鹃属物种的群体遗传学研究中.本研究采用限制性双酶切位点测序（double digest restriction-site associated DNA scquencing,ddRAD-seq）技术开发SNP 位点[14].该技术具有操作简单、建库周期快、通量高且实验成本低等诸多优点[15-17]，可以从任何物种（包括非模式物种）简化的全基因组范围内挖掘大量的SNP[18-19]，可为在DNA 层面上研究波叶杜鹃的遗传背景提供足够的信息位点.以此揭示该物种的遗传多样性水平和遗传结构，探讨种群演化历史，加深对其进化过程的理解，以期为波叶杜鹃种质资源的保护和分子标记辅助育种提供科学的理论依据.

1.1 试验材料对波叶杜鹃自然分布范围内仅存的5 个居群，以及峨眉山植物园自峨眉山引种的1个居群进行采样（表1），摘取成株幼嫩健康叶片，用硅胶干燥带回实验室后在-20 ℃冰箱中保存.样品采集按种群遗传学的原理和方法，遵循代表性、可对照性原则取样（样本间距 1～1.96 m）[20].

表1 波叶杜鹃的采集地信息Tab.1 Collection sites information of R.hemsleyanum

1.2 DNA 提取和ddRAD 文库构建总基因组DNA 提取采用十六烷基三甲基溴化铵（CTAB）法[21].采用1.2 %的琼脂糖凝胶电泳检测DNA 的完整性，QUBIT 3.0 荧光定量仪测定样品DNA 的浓度和纯度.

ddRAD 文库构建依照Peterson 等[14]的方法进行.将DNA原液稀释至10ng/μL，用EcoRⅠ和MseⅠ两种酶在37 ℃下反应5 h，然后在 PCR 仪中升温至65 ℃，保持20 min 使酶失活.产物温度降至12 ℃后加入T4 连接酶和接头（NEB），在16 ℃下连接4 h，再在65 ℃下反应20 min 使酶灭活，最后12 ℃保温.连接产物等量混合，用1.2%的琼脂糖凝胶中进行分离，割胶回收300～500 bp 的区域.水浴溶胶后用gel extraction kit （Omega Inc.,USA）回收片段.将纯化产物进行PCR 扩增至预期浓度，利用Illumina HiSeq X-ten 平台进行150 bp 双端测序（每个样本产生大约0.5 Gb 数据）.DNA 提取和ddRAD-seq 文库制备由捷锐生物科技有限公司完成（中国广州）.

1.3 SNP 挖掘利用Stacks v.2.0 处理ddRAD 测序数据[22-23].首先通过process_radtags 模块进行数据质控，剔除原始数据(raw data)中低质量以及没有rad 的标签序列.使用ustacks 模块对单样品数据进行cluster 聚类以及去重，将每一个体内测序深度达到要求的且错配数符合设置条件的相似序列保留下来，形成loci(参数设定，-min depth of coverage required to create a stack (m)：2；
-max distance allowed between stacks (M)：2；
-max distance allowed to align secondary reads to primary stacks：4；
-max number of stacks allowed perde novolocus：3；
-minimum alignment length：0.8；
-model type：SNP，-alpha significance level for model：0.01).使 用cstacks 模块将所有个体获得的loci 位点合并成catalogs，参数默认.运行sstacks 模块将单样品loci比对至catalog 位点，得到样品的SNP，allele 和tags 信息.使用population 模块进行SNP 的calling(参数设定，-min number of populations a locus must be present in to process a locus(r)：6；
-min percentage of individuals in a population required to process a locus for that population(p)：0.75；
-max observed heterozygosity required to process a nucleotide site at a locus：0.6；
-min minor allele frequency required to process a nucleotide site at a locus：0.01).由 于 同一短序列上不同位置的SNPs 距离较近，存在极强的连锁不平衡现象，因此利用软件中的脚本(--write_random_snp)在每个位点上随机选择一个SNP 用于后续的分析[24-25].

1.4 数据分析

1.4.1 遗传多样性和与遗传结构分析 使用vcftool v.0.1.15 计算Tajima’sD值[26]，以95%的置信区间（-1.082～2.039）对获得的SNP 进行中性检验[15,27]，运算的滑动窗口大小设置为3 000 bp.利用Stacks 软件的populations 模块计算遗传多样性参数[22]，包括观察杂合度（Ho）、期望杂合度（He）、核苷酸多样性（π）、私有等位基因数目（Private）和近交系数（FIS）等.使用PGDSpider v.2.1.1.5 进行格式转换用于后续分析[28].

利用Arlequin v.3.5.2 软件的AMOVA （analysis of molecular variance）模块[29]对3 896 个SNPs 进行分子方差分析，计算居群间和居群内的遗传变异来源，估算居群间成对的遗传分化系数（FST）.通过MEGA-X 软件[30]根据居群间的遗传距离进行NJ （neighbor joining）聚类分析.利用Structure v2.3.5 软件[31]对所有参试个体进行聚类分析.选择依据Alpha 值的祖先混合模型和依据FST值的相关等位基因频率模型，burn-in 次数为50 000，MCMC 迭代次数100 000.设置群体数目（K）值为1～10，重复独立运算10 次.利用Structure Harvester在 线 软 件[32]（http://taylor0.biology.ucla.edu/structu reHarvester/），根据ΔK数值的分布确定最佳K值聚类[33].利用Plink v.1.9 进行主成分分析（PCA）[34]，通过R v.3.6.1的ggplot2 包进行分析结果的可视化，确定参试样本的亚组结构.

1.4.2 Stairway Plot 推测种群动态历史 通过脚本easySFS （https://github.com/isaacovercast/easySFS）转换VCF 生成SNP 频谱（SFS）格式，用Stairway Plot v.2.1[35-36]从所有位点中随机取67% 的位点，推断有效种群大小（Ne）随时间的变化.已有研究表明杜鹃属物种的世代时间普遍在10 a 左右，如彭土杜鹃（R.ponticum）[37]从幼苗至开花结实约需要10 a，以及Liu 等[12]在模拟蓝果杜鹃的居群动态历史时发现，世代设置在10、20、30 a 时呈现出相似的统计模型.根据Yoichi 等[38]的研究，高丽杜鹃（R.Weyrichii）突变速率为1.581×10-9.因此，将波叶杜鹃世代设置为10 a，突变率设置为1.581×10-9.

2.1 ddRAD 测序与数据处理经ddRAD 测序，78 个波叶杜鹃样本共获得测序数据量89.57 Gb.经过process_radtags 模块剔除低质量序列以及无rad标签序列，共保留了87.04 Gb 的清洗数据用于下一步分析.使用ustacks 模块对单样品数据进行cluster 聚类以及去重，共获得11 477 344 个loci.使用cstacks 模块将所有个体获得的loci 合并，参数采取默认，共得到10 603 712 个catalog 位点.运行popolation 模块对共有SNP 位点进行筛选，共获得3 896 个SNP 位点（每个位点随机选择一个SNP）.

2.2 遗传多样性和遗传结构分析经Tajima’D值对获得的SNP 进行中性检验，受选择性的位点仅429 个.物种水平上基于所有位点的Tajima’D均值是0.478 9，说明稀有等位基因频率低(Tajima’D＞ 0)[39]，其原因可能是极危物种分布区域狭窄，环境相似，其受到的选择压力在单一种群或者临近种群均一化所导致.因而后续所有的分析均基于全部SNP 进行.

表2 列出了波叶杜鹃6 个居群的遗传多样性参数.参试居群的期望杂合度（He）和核苷酸多样性（π）平均值分别为0.226 7±0.002 和0.241 0±0.0032.其中，SSG 居群遗传多样性最高（He=0.240 9±0.002 8,π=0.250 4±0.002 9），ZWY 居 群 的遗传多样性最低（He=0.202 8±0.003 2，π=0.228 0±0.003 6）.种群观测杂合度（Ho）均小于期望杂合度（He），平均值分别为0.163 6±0.002 和0.226 7±0.002.各居群中均有特有等位基因，平均值为16.居群内近交系数（FIS）平均值为0.222 0±0.022 2，均为正值，表明其种群杂合子缺失，纯合子过剩，存在一定程度的自交或近交.

表2 波叶杜鹃的种群遗传多样性参数Tab.2 Population genetic statistics of R.hemsleyanum

居群间的遗传距离为0.28～3.29，平均值为1.786 （表3），居群间的遗传分化系数（FST）为0.036 5～0.114 8，平均值为0.062 6（表3），居群间的遗传分化处于中等水平（0.05＜FST＜0.15）.这一点也得到了分子方差分析的支持（表4），有98.42%的遗传变异来源于居群内，居群间变异仅占总变异的1.58%，表明其在遗传上相似，群体分化程度不高.

表3 波叶杜鹃种群间的遗传分化系数(FST 值，在对角线以上)和遗传距离(对角线以下)Tab.3 Genetic distances (FST,above diagonal) and geographic distances (DR,below diagonal)between R.hemsleyanum populations

表4 波叶杜鹃的分子方差(AMOVA)分析Tab.4 Analysis of molecular variance (AMOVA) in the wild populations of R.hemsleyanum

基于SNP 信息进行Structure 分析，由Structure Harvest 软件生成的结果显示，当K=2 时，ΔK出现最大值（图1）.表明波叶杜鹃78 个参试个体的最合理聚类组数为2.在K=2 的情况下，根据参试的单株个体分配到不同类群的比值（Q）（图2），GXP居群聚为类群Ⅰ（绿色），HZP、SSG、XXS、ZLP、和ZWY 居群聚为类群Ⅱ（红色）.除ZLP（1 株）居群中的少数个体（占比1.2%）归属于不同类群（Q＜0.6），98.8%的植株个体谱系清晰.

图1 基于Structure 分析的K 与ΔK 的变化趋势Fig.1 Trend of K and ΔK based on Structure analysis

图2 基于SNP 标记的K=2 时波叶杜鹃的群体遗传结构Fig.2 Genetic structure of R.hemsleyanum population at K=2 based on SNP markers

PCA 聚类结果（图3）与Structure 结果一致，前2 个极轴坐标分别解释了整体遗传变异的54.20%和40.44%，78 个参试样本主要分为两个聚类（cluster），GXP 居群的14 个个体聚为一组（clusterⅠ），来自HZP、SSG、XXS、ZLP、ZWY 居群的个体彼此间显示出较近的亲缘关系，倾向于聚为另一个组（clusterⅡ）.聚类分析中6 个居群逐步聚类（图4），很难再区分亚类；
也表明波叶杜鹃具有较低分化的群体结构.

图3 波叶杜鹃78 个个体的主成分分析Fig.3 PCA of 78 R.hemsleyanum individuals

图4 波叶杜鹃 6 个群体的聚类分析Fig.4 UPGMA clustering of six R.hemsleyanum populations

2.3 种群动态历史分析Stairway plot 基于SNP频谱（SFS）推断的结果表明（图5），波叶杜鹃居群的历史时间最早可追溯至214 kaBP （kaBP,kilionanniversary Before Present），初始有效种群大小（Ne）为19 075 个.至100 kaBP，即中更新世（1～0.1 maBP）（maBP,million-anniversary Before Present）期间，波叶杜鹃的有效种群规模急速扩张，达到原来的6.87 倍（Ne=131 × 103），然后逐渐减少，在末次冰期冰盛期（Last Glacial Maximum,LGM）（26.5～19.0 kaBP）以及冰期结束后的全新世（11.7 kaBP）持续下降，直到大约1 000 年前.

图5 波叶杜鹃种群动态历史分析Fig.5 Demographic history of R.hemsleyanum lineages by Stairway plot (95% confidence interval is indicated by thin grey lines)

3.1 居群的遗传多样性及遗传分化研究物种的遗传多样性，可以反映物种适应环境变化的潜力,为物种资源特别是濒危物种现有的生存方式及状态进行有效的评估[6,40].物种的遗传多样性水平与其繁育系统、生态习性、物种分布区的地域特征、气候变化及演化历史等诸多因素有关[41-42].研究表明，由于遗传漂变和近交的原因，分布广泛的物种遗传多样性高于地理分布狭窄的稀有和濒危物种[42-43]，与草本植物相比，木本植物在物种水平具有较高的遗传多样性[44].对于双等位SNP 分子标记来说，核苷酸多样性（π）是衡量群体遗传多样性的整体指标[45].本研究结果表明，木本极危树种波叶杜鹃的遗传多样性（He=0.226 7±0.002；
π=0.241 0±0.003 2）水平与同样使用ddRAD-seq 技术研究遗传多样性的草本植物相比相对较高，如密花石斛（Dendrobium densiflorum） （π=0.105 6）[46]、地宝兰（Geodorum densiflorum） （π=0.035 9）[46]、Viola uliginosa（π=0.044 0）[47]等.与同属物种相比，杜鹃属物种的遗传多样性普遍存在很大的差异，这可能源于不同的有性生殖水平、不同的种群大小、物种和种群的年龄，以及面临的压力和选择不同[47]，如蓝果杜鹃（He=0.067 5）[12]、宽杯杜鹃（He=0.185 6）[13]、长梗杜鹃（R.Longipedicellatum）（He=0.561 0）、大树杜鹃（R.protistumvar.giganteum）（He=0.602 0）[48]等.波叶杜鹃的遗传多样性处于中等水平，其中，ZWY 居群作为栽培群体，在6 个群体中遗传多样性最低，但是群体间遗传一致度变化范围较小，ZWY 居群依然具有丰富的遗传多样性（He=0.202 8±0.003 2,π=0.228 0±0.003 6），说明这些居群间亲缘关系较近.从遗传学的角度分析，丰富的遗传多样性意味着具有较强的适应能力和丰富的遗传改良潜力.波叶杜鹃遗传多样性的维持，多是其繁育系统、生态习性以及自然选择的结果.杜鹃属植物为多年生、长寿命木本树种，频繁的自然杂交且需要传粉昆虫等本身的生物学特性[44,49]，有利于波叶杜鹃和大部分杜鹃属物种维持较高的遗传多样性水平.波叶杜鹃的种群观测杂合度（Ho=0.163 6±0.002 0）均小于期望杂合度（He=0.226 7±0.002 0），说明杂合子缺失（FIS＞ 0），群体间亲缘关系较近，存在一定程度的近交，表明生境片断化的遗传效应正逐渐显现.

遗传分化系数（FST）是衡量等位基因频率在不同群体间分化程度的指标[50].在本研究中，6 个波叶杜鹃居群间的遗传分化系数（FST）平均值为0.062 6（表3），在中度分化范围之内（0.05～0.15）[51]，分化程度不高.基因流是影响种群间遗传分化的主要因素[52]，植物种群间的基因流动主要是通过花粉和种子传播来实现[51,53].波叶杜鹃现有野生居群均位于峨眉山且地理距离较近（0.28～3.29 km），交配系统以异交为主，有效传粉昆虫雄蜂（Bombussp.）等具有较长距离传粉的能力.此外，密集的大型总状伞形花序及花香和腺体，能增加其显著度，吸引传粉昆虫来访，且其种子小、轻、具翅，易被风吹散，均能促进不同种群间基因交流，可能会抵消部分生境严重破碎化造成的不利影响，有助于维持群体稳定，在一定程度上增强了群体对环境变化的适应性[53-54]，这也是波叶杜鹃遗传分化水平不高的主要原因.将波叶杜鹃的6 个居群视为一个类群进行AMOVA 分析发现，居群内和居群间均存在显著的遗传变异，居群内变异远大于居群间，与Liu等[12]利用ddRAD-seq 技术分析蓝果杜鹃4 个居群间遗传分化（FST=0.031 4～0.045 2）的结果一致.这进一步验证了多年生木本植物由于异花传粉等特性，居群内遗传变异水平较高，而居群间遗传变异水平较低的观点[55].因此，在选择较高遗传多样性群体的基础上，应该侧重于群体内个体的选择.聚类分析中，Structure 分析及PCA 分析均支持波叶杜鹃群体并不是遗传组分均质的单一群体，而是具有2 个遗传组分异质的亚居群.鉴于两亚居群间相对较小的遗传分化（FST=0.068 5），以及生境特征和两者间存在的山谷、人为建筑、农田等屏障间隔，推测亚居群可能是由生境的彼此隔离以及微环境的选择进化造成的[40,54].

3.2 种群动态历史波叶杜鹃仅残存于四川省峨眉山（29°31′～29°38′N，103°15′～103°28′E），峨眉山地处川西高原边缘地带，地形地貌特殊，垂直气候带复杂，物种组成丰富[56].已有研究表明，自更新世以来青藏高原的快速隆升以及第四纪冰期和间冰期的反复波动对青藏高原及邻近地区的动植物种类迁移、分化和物种形成等事件产生了重要的影响[57-61].动态历史研究表明，波叶杜鹃最早可追溯至214 kaBP，其有效种群大小（Ne）大约在中更新世之后开始下降，在末次冰期冰盛期（LGM）（26.5～19.0 kaBP）以及冰期结束后的全新世（11.7 kaBP）持续下降（图4）.冰期全球气温的显著下降、冰盖的形成和海陆形势的变化、气候带的相应移动，以及大气环流和洋流的变化，都会直接影响植物的生长、演化和分布[38,62-63].冰期以及全新世以来气候环境的一系列变化，有限的遗传多样性水平可能削弱了波叶杜鹃居群在迁徙过程中的适应潜力，导致其有效种群规模降低，再加上人为活动干扰（发展旅游、基础设施建设等）等不利的环境条件使得其生存压力增大，种群规模进一步减小，只能以极小种群的形式生活在原栖息地，无法恢复到正常规模[60,64].根据种群动态推断的结果以及实地调查发现，目前波叶杜鹃的有效种群大小远小于可以防止适应性丧失而推荐的阈值（Ne≥ 100），尤其小于可以保留进化潜力而建议的有效种群大小的最低值（Ne≥ 1 000）[65].因此，亟待重视和加强对波叶杜鹃这些残存居群的保护.

3.3 保护建议保护物种的遗传多样性，是防止该物种特异种质流失的关键，特别是针对波叶杜鹃这样天然群体数量极小的濒危特有植物，对物种多样性的维持以及长期生存均具有重要的意义.为了有效维持波叶杜鹃的遗传多样性并实现种群规模恢复，在今后的工作中，有必要进一步研究波叶杜鹃的无性繁殖技术，通过组织培养、扦插等方式实现波叶杜鹃的快速繁殖.根据本研究，SSG 居群和ZLP 居群的遗传多样性最为丰富，在种质资源收集和回归引种时应重点考虑.其次，建议开展就地保护与迁地保护、野外回归相结合的策略.在原生地区进行区域划分，严格限制旅游设施建设，防止人为干扰，做好抚育管理，有效保护其自然栖息地.在就地保护的同时，通过测定波叶杜鹃当前生境的生态因子状况，寻找合适的适生区域，实行实生苗的近地迁地保护，以扩大波叶杜鹃的种群分布，实现种群的复壮和人工种群的自然回归，最终避免这一物种的灭绝，并得以有效保护和利用.

本研究首次利用ddRAD-seq 测序技术，系统地分析了波叶杜鹃的遗传背景.分析结果表明，波叶杜鹃在物种和居群水平上，具有丰富的遗传多样性，进一步选择和育种利用的潜力较好.遗传分化水平中等，遗传变异主要来源于居群内.由于近期冰期气候环境的变化，再加上人为活动（发展旅游、基础设施建设等），导致波叶杜鹃的生存环境逐步恶化，使得该物种面临极高的灭绝风险.从对该濒危植物保护的角度考虑，应尽快建立自然保护小区予以有效保护，同时加强人工繁育技术研究，在扩大种群规模的同时维持其遗传多样性，这对波叶杜鹃的种群恢复和开发利用都具有重要的意义.

猜你喜欢 位点杜鹃分化 Pd改性多活性位点催化剂NH3-SCR脱硝反应机理研究分子催化(2022年1期)2022-11-02DNA脱碱基位点的检测方法及其生物学研究进展分析测试学报(2022年9期)2022-09-21多环境下玉米保绿相关性状遗传位点的挖掘中国农业科学(2022年16期)2022-09-19两次中美货币政策分化的比较及启示清华金融评论(2022年4期)2022-04-13云上杜鹃家教世界·创新阅读(2022年3期)2022-04-12一种改进的多聚腺苷酸化位点提取方法电脑知识与技术(2018年19期)2018-11-01红杜鹃幸福·悦读(2016年11期)2016-12-21作为春花的杜鹃含笑花(2015年4期)2015-02-27日本中世末期五山禅僧的“儒·释”论争与其内部分化古代文明(2014年1期)2014-02-23《诗经·大车》正音与谈部中古分化探文学教育·中旬版(2012年4期)2013-02-01 相关热词搜索：杜鹃，遗传，标记，