基于全基因组SNP分子标记分析青海云杉遗传结构

张宏斌 吕 东 赵 明 赵 祜 赵兴鹏 李 伟

（1. 甘肃省祁连山水源涵养林研究院，张掖 734000；2. 祁连山特有植物繁育及推广国家地方联合工程研究中心，张掖 734000；3. 北京林业大学生物科学与技术学院，北京 100083）

青海云杉（）是松科（Pinaceae）云杉属（）常绿乔木，为我国特有树种，具有涵养水源、保持水土的作用，是祁连山区主要的水源涵养林和造林更新树种。具有极高的生态价值；木材质量好，可供桥梁、家具及造纸等用材，是重要的用材树种，具有重要的经济价值；四季常青、外观优美，被广泛应用于城市绿化、园林栽植中，具有观赏价值。几十年来，我国林业工作人员对青海云杉的进化演变、空间分布、水分利用、气候响应等进行了一系列研究；并开展了青海云杉母树林与种子园的营建、表型变异、遗传变异及子代测定等方面的遗传改良工作，取得了一些成果。

青海云杉表型变异丰富，在分子水平上也具有较高的遗传变异水平，在群体间和群体内都存在着广泛的变异，具有丰富的遗传多样性，因此，进一步挖掘青海云杉群体变异信息，对青海云杉种质资源保护和遗传改良工作的科学开展具有重要意义。但青海云杉在改良过程中存在许多亟待解决的问题。第一，由于选育周期长，导致良种选育慢良种少，青海云杉遗传改良进程缓慢；已有许多针叶树种相继进入第2 代种子园的建设当中，但由于技术、科技投入不足，青海云杉的遗传改良还停留在初级种子园阶段。第二，由于全基因组较大，青海云杉遗传改良基础薄弱，在分子水平上的研究较少，缺少基因组信息，分子标记多来自挪威云杉（），现有引物通用性差，数量少，导致标记相关分析准确性较低，缺少青海云杉特有的分子标记。近年来，祁连山区生态环境逐步恶化，青海云杉作为主要造林绿化树种，对其良种质量和数量的要求逐渐升高，迫切需求青海云杉良种升级换代，亟待开展青海云杉的高世代亲本选择与遗传评价等基础研究，为青海云杉高遗传改良层次提升奠定基础。

传统的选择优良单株和杂交的育种方式在青海云杉育种过程中所需周期太长，利用目标性状与分子标记的关联分析，筛选出与目标性状紧密关联的标记，为分子标记辅助选择提供依据，可以极大缩减青海云杉遗传改良周期，提升青海云杉遗传改良水平。因此，本研究采用SLAF-seq技术，开发全基因组水平上的青海云杉SNP 分子标记，为青海云杉种质资源的研究、遗传连锁图谱构建以及QTL定位等方面的研究奠定基础。

1 材料与方法

1.1 试验材料

试验材料选取张掖龙渠青海云杉无性系种子园中的106 个青海云杉亲本无性系。种子园建园材料来源于祁连山区的11 个林场（见表1），2019年采集每个无性系单株当年生健康针叶于自封袋中，分别编号，及时存储于液氮中，回到室内，放冰箱里冷藏备用。

表1 试验材料信息Table 1 Experimental material information

1.2 基因组DNA提取与检测

基因组DNA 的提取采用改良CTAB 法，提取的青海云杉基因组DNA 质量浓度不低于20 ng·μL，总量不低于4 μg，保证质量符合建库要求。

1.3 SLAF文库构建和高通量测序

由于目前缺少青海云杉基因组序列的相关信息，选取挪威云杉基因组作为参考基因组（ftp：//plantgenie. org/Data/ConGenIE/Picea_abies/v1.0/），组装出的基因组大小为12 G，GC 含量为37.88%。根据挪威云杉参考基因组，进行电子酶切预测，确定最适酶切方案。

对质量合格的各样品基因组DNA 进行酶切，得到的酶切片段（SLAF 标签）进行3′端加A 处理、连接测序接头后，进行PCR 扩增，通过纯化、混样、切胶等步骤选取目的片段，文库质检合格后用Il⁃lumina HiSeq 进行双端测序。选用水稻日本晴（ssp.）作为试验建库的对照，基因组数据参考来源于：http：//rapdb.dna.affrc.go.jp。由北京百迈克生物科技有限公司完成SLAF文库的构建以及测序。

1.4 SLAF测序数据分析和SNP位点筛选

对获得的Raw Reads 进行质控，得到Clean Reads。将得到的reads利用BWA软件比对到参考基因组上，统计SLAF 标签和多态性SLAF 标签的数量。SNP 分子标记采用GATK 和samtools 两种软件共同开发，以次要基因型频率（MAF）＞0.05、完整性＞0.8 为筛选标准，两个软件开发SNP标记的交集作为最终可靠的SNP标记。

1.5 遗传结构分析

根据开发的青海云杉高质量SNP 标记，采用MEGA X软件，基于邻接法（neighbor-joining），构建各青海云杉样品的系统发育树。利用admixture软件，研究青海云杉群体遗传结构，预先设定亚群数目（）为1～10，进行聚类，并对聚类结果计算交叉验证错误率，依据其谷值确定最优分群数。并且利用EIGENSOFT 软件，进行主成分分析，得到样品的聚类情况。

2 结果与分析

2.1 试验建库评估

根据参考基因组的电子酶切预测，最终确定限制性内切酶为HinCII，长度在314～394 bp的酶切片段序列定义为SLAF标签，预测将获得506 448个SLAF 标签。本研究中水稻日本晴数据的酶切效率高达98.12%，表明酶切反应正常。通过SOAP软件将测序的水稻日本晴reads 与其参考基因组进行比对，显示青海云杉基因组DNA构建的SLAF文库双端比对效率为95.33%，说明青海云杉基因组DNA构建的SLAF建库正常。

2.2 测序数据评估

通过SLAF-seq 测序后，对照组水稻日本晴测序获得1.48 Mb reads 的数据量，106 个青海云杉样品共获得1 375.57 Mb reads 数据。对照水稻日本晴reads 数为1 475 338，GC 含量为40.42%，Q30 为92.69%。各样品的测序结果（见表2），所获得的reads 数目范围为4 411 192～30 427 534，其中，7 号所获得的数据量最多，36 号获得的数据量最少。GC 含量范围在39.19%～43.62%，平均GC 含量为40.27%，测序质量值Q30 范围在90.97%～97.82%，平均Q30 为95.51%。本研究中所测序列的Q30 较高，表明碱基测序错误率低，说明测序质量高，可以用作下一步的数据分析。

表2 青海云杉测序结果Table 2 Sequencing result of P.crassifolia

2.3 SLAF标签与SNP标记开发

通过对测序数据的处理分析，本研究共开发4 058 883个SLAF标签，平均每个青海云杉样本开发了389 459 个SLAF 标签，其中多态性SLAF 标签有1 573 899 个，占SLAF 标签总量的38.78%，标签的平均测序深度为21.21×（见表3）。共开发12 275 765个青海云杉SNP标记，各青海云杉样本的SNP数量为1 890 934～4 487 841，各样品检测到的SNP完整度为15.40%～36.56%，SNP的杂合率为5.41%～10.99%（见表4）。

表3 青海云杉SLAF 标签统计结果Table 3 SLAF tag statistics of P.crassifolia

表4 青海云杉SNP 信息统计结果Table 4 SNP information statistics of P.crassifolia

2.4 群体遗传结构

利用已开发的高质量青海云杉SNP 标记，构建了106 个青海云杉的系统发育树，结果如图1所示，大致可分为6 小组，分别具有5、19、15、25、21、21个青海云杉无性系，来自多数种源的青海云杉在各组中分布比较均匀，不同种源的青海云杉多聚为一类，但每一类在不同种源的占比差别较大。

图1 106个青海云杉材料的系统进化树Fig.1 The phylogenetic tree of 106 P.crassifolia materials

基于SNP 标记，利用admixture 软件分析了106 个青海云杉无性系群体结构，结果表明当=1时，交叉验证错误率最低（图2：A～B），说明这些无性系来源于同一个祖先的可能性较大。

通过主成分分析，绝大多数来自不同种源的青海云杉形成一簇，只有小部分的青海云杉无性系分散分布，并未出现显著的分化（图2：C～E），表明无性系间亲缘关系相近。

图2 青海云杉遗传结构分析A.青海云杉样品聚类结果；B.K值的交叉验证错误率；C～E.PCA聚类Fig.2 Analysis of genetic structure of P.crassifolia A.Clustering results of P.crassifolia；B.Cross validation error rates of K-values；C−E.Principal components analysis

3 讨论

针叶树种的基因组庞大复杂，目前仅有挪威云杉、白云杉、火炬松进行了基因组测序，利用传统的方法进行分子标记的开发，工作量大且复杂，花费时间长，得到标记的数量不能满足试验的要求。袁行栋从已发表的挪威云杉60 对SSR 引物中筛选出了10 对具有多态性的SSR 引物，用作青海云杉遗传多样性分析，但分子标记的数量少，引物之间的通用性较差，导致利用分子标记的相关分析准确性较低。

3.1 青海云杉全基因组SNP分子标记开发

单核苷酸多态性（SNP）分子标记是指DNA 序列发生多态性的变化，是由单个核苷酸发生缺失、插入、转换等变异导致的，由于其在基因组中数量众多，且存在范围广阔，在遗传图谱的构建、数量性状定位分析和全基因组关联分析中被广泛应用。SLAF-seq 技术是一种快速开发大量SNP 标记的技术，对于物种的要求较低，没有其基因组信息也可以顺利开发出大量SNP 标记，对于未知基因组物种全基因组分子标记的开发是非常有效的方法，在动植物中已经得到了大量的使用。在针叶树中，段红静通过SLAF-seq技术，首次在杉木群体中开发出166 646 个SNP 标记。Bai 等采用SLAF-seq 技术对广东马尾松种质资源进行分析，在599 164个SLAF多态性标签中，共鉴定出471 660个SNP标记。董明亮利用SLAF-seq技术在华北落叶松全基因组范围内挖掘SNP 标记，共获得6 323 943 个SLAF 位 点，检 测 到324 352 个SNP 标记，其中有122 785 个呈现多态性。王飞等利用SLAF-seq 技术开发华山松SNP 标记，共获得SLAF标签12 952 676 个，多态性SLAF 标签为1 456 486个，开发了3 469 074个群体SNP标记。

在本研究中，首次将SLAF-seq 技术应用于青海云杉群体，共获得4 058 883 个SLAF 标签，开发了12 275 765个青海云杉SNP标记，获得了数量大且分布广泛的分子标记，并基于开发的分子标记，分析了青海云杉无性系群体结构，结果表明不同种源的青海云杉来源于同一个亚群的可能性较大。青海云杉SNP 标记的开发为后续分子育种工作奠定了良好的基础，可为遗传评价、遗传图谱的构建以及目标性状的关联分析提供数据支撑。

3.2 青海云杉种子园遗传结构

种子园是按设计要求营建的，以优树的无性系或家系为材料，以生产优良遗传品质的种子为目的营造的特种人工林。青海云杉种子园的研究包括林木生长，繁殖特性等。对于遗传生物学形状也有一些研究，吕东等以与本研究相同的材料进行无性系结实性状遗传变异研究，表明青海云杉无性系间的结实能力差异明显，其无性系重复力较高，球果主要集中在树冠中上部的冠层分枝上。袁行栋2016年对其进行了遗传多样性研究，表明青海云杉无性系较高的遗传多样性水平。本研究在前人的基础上进行，与之前遗传多样性水平较高的结论一致。

家系来源的11个林区基本涵盖了青海云杉的自然分布区，结合地理分布分析。连城林区（LC）位于所有种源地的最东南方向，其15 个无性系分布在第Ⅱ、Ⅲ、Ⅴ、Ⅵ分组中，而仅Ⅴ组和Ⅵ组又多达11 个。大河口林区位于各林区的中间位置，其9 个无性系在5 个组分里面都有分布。隆畅河林区在最西北角上，10 个无性系中有7 个属于第Ⅳ组分，其余Ⅲ、Ⅴ、Ⅵ组分各有1 个。总体来说，靠中间的种源地青海云杉在各组中分布相对比较均匀一般能聚为一类，靠边缘的种源地虽然仍有多个组分的遗传物质来源，但都不能涵盖所有组分且都有特定组分占大比例。青海云杉无性系群体结构未出现分层情况，这些无性系来源于同一个祖先的可能性较大。

张掖龙渠青海云杉无性系种子园具有较高的遗传多样性水平，可为高世代种子园的建设提供丰富的原材料。该种子园建园初期包括了大量的无性系，拥有丰富的育种资源，但之后没有连续补充新材料，良种基地发展的物质基础不够充足。没有充分重视遗传测定工作，对收集的资源所做的作遗传测定连续性不足，基地的更新换代和发展的理论依据薄弱。建议今后应不断补充新的育种资源，对资源进行有计划的遗传测定，同时可以采用多种管理技术措施；以充实基地发展的物质基础，为基地的更新换代提供足够的理论依据。

4 结论

通过SLAF-seq 技术，共开发了4 058 883 个SLAF标签，其中多态性SLAF标签有1 573 899个，占SLAF标签总量的38.78%，标签的平均测序深度为21.21×，共开发了12 275 765 个SNP 标记，为今后遗传多样性的分析、遗传图谱的构建等提供了基础数据。

对青海云杉无性系遗传结构进行分析，根据系统发育树的结果，青海云杉可分为6 小组，来自不同种源的青海云杉在各组中分布整体比较均匀，但相距较远的种源之间存在一定差别。群体结构分析中，当=1 时，交叉验证错误率最低，说明这些无性系来源于同一个祖先的可能性较大，主成分分析显示大部分来自不同种源的青海云杉无性系聚为一类，说明不同种源的青海云杉无性系间的亲缘关系较近，青海云杉无性系群体结构未出现分层情况，这些无性系来源于同一个祖先的可能性较大。