基于GBS-SNP 的武夷茶树（Camellia sinensis，Synonym：Thea bohea L.）遗传分析及标记开发

李力，罗盛财，王飞权，黎巷汝，冯花，石玉涛，叶江华，刘菲，赵佳林，李舒莹，张渤*

1. 武夷学院茶与食品学院，福建 南平 354300；2. 武夷山市农业农村局，福建 南平 354300；3. 武夷学院茶叶科学研究所，福建 南平 354300；4. 福建农林大学园艺学院，福建 福州 350002

武夷茶是山茶科（Theaceae）山茶属（Camellia）的一类成员，以“武夷”的英文音译而被命名为Bohea或VarBohea（Thea BoheaLinnaeus. 1762）。武夷茶生长于中国福建省北部的武夷山，武夷山被认为是乌龙茶和红茶的发源地[1-2]。武夷山地区的先民从有性系的武夷茶树群体中筛选分离出品质优良的稀有单丛单株，通过无性繁殖传承至今，成为现在著名的“武夷名丛”[3]。由于武夷茶树生长在岩石土壤中，其鲜叶制作的茶叶具有“岩韵”（风味浓郁，香味持久）的品质特点，被称为武夷岩茶，是中国十大名茶之一，具有巨大的经济和文化价值[4]。

历史上记载的武夷茶树品种有上千种。然而，许多种质资源不断消失，通过对茶树性状的长期观察，目前仅收集到百余个种质资源[5-6]。此外，武夷山当地的育种过程中存在对优良品种过度选择的情况，将导致群体遗传多样性越渐狭窄，最终影响品种性状的综合改良，再加上存在茶树同物异名的现象，进一步阻碍了武夷茶树种质资源的开发。因此，对武夷茶的遗传背景分析和品种鉴定是一项重要的工作。传统的茶树种质资源分类基于形态或农艺性状，容易受到环境条件的限制[7]，这给品种的确定带来了一定的困难[8]。国际植物新品种保护联盟（UPOV）对作物品种建立了“差异性、统一性和稳定性”测试制度[9]，在难以通过表型区分确认新品种的情况下，DNA 指纹图谱可以帮助提高品种的识别。新修改的《中华人民共和国种子法》于2022 年3 月1 日正式施行，首次建立实质性派生品种（Essentially derived variety，EDV）制度[10]，要求对品种的基因或基因型组合进行特异性检测鉴定。建立真实可靠的DNA 鉴定检测技术对新《种子法》修订条例的实施具有重大战略意义。

目前，茶树的特异性DNA 鉴定手段仍相对落后于主要农作物，早期的茶树特异性DNA 鉴定方法多采用第一代或第二代分子标记，包括扩增片段长度多态性（AFLPs）、限制性片段长度多态性（RFLP）、随机扩增多态性DNA（RAPD）、简单序列重复序列（SSRs）等[11-14]，武夷茶树也是如此[15-17]。虽然这些传统方法各有优势，但所使用的标记数据量有限，不适合大样本量的群体分析，不能充分反映物种的遗传多样性。相比之下，第三代分子标记技术单核苷酸多态性（Single nucleotide polymorphism，SNP）[18-20]具有一定的优势，已逐步成为主流的分子标记[21]。Lin 等[22]采用高通量SNP 技术构建了4 个主要乌龙茶产区的DNA 指纹图谱，揭示了4 个产区的遗传关系。Liu 等[23]从茶树的表达序列标签（EST）数据库中筛选出49 个SNP 位点，分析了武夷山及其邻近地区137 份茶树种质的遗传多样性丰富程度及遗传关系。这表明高通量SNP技术在寻找和验证茶树亲本关系方面具有很大的发展空间。然而，茶树的SNP 分子标记开发发展较晚，数量有限，仍需不断完善补充。

传统的SNP 筛选方法需要大量的人力资源和较长的时间，基因分型（Genotyping by sequencing，GBS）测序技术可以快速识别高密度SNP 位点，广泛应用于作物种质鉴定、群体结构分析和全基因组关联分析[24-25]。本研究收集了126 份武夷茶品种/品系，结合来自12 个省份的223 份优异茶树种质资源，共349份茶树样本进行GBS 测序，利用高质量SNPs对武夷茶进行遗传多样性及背景分析，并开发一种简便、快速的茶叶SNP 分子标记鉴定，旨在为茶树研究领域EDV 制度的有效实施提供参考。

1 材料与方法

1.1 样本收集与DNA 提取

样品来自武夷山市龟岩茶树种质资源保护基地（福建省优异茶树种质资源保护区闽茶圃004，N 27°43'42.46"，E 118°0'14.40"）和武夷学院茶树种质资源圃（N 27°61'01.34"，E 117°96'63.51"），包括126 个武夷茶品种/品系及223 个来自12 个省份的茶树品种/品系，共349 个（包含同一品种/品系的6 对）。样本信息见附表1（扫文后二维码查看附表1）。

新鲜叶片基因组 DNA 采用植物基因组DNA 提取试剂盒（康威CW0553S）提取，DNA浓度和纯度采用 NanoDrop 紫外分光亮度计（Thermo Scientific，USA）测定。DNA 完整性使用0.8%琼脂糖电泳检测，基因组DNA 保存于-80 ℃。

1.2 基因分型测序及建库

GBS 库的构建参照文献[26]，用限制性内切酶（EcoRⅠ和NiaⅢ）酶解质检合格的基因组DNA，在酶切后的片段两端加上相应的适配接头，利用PCR 扩增两个接头之间的片段，共构建了349 个对应样本的多重库，每个库DNA 样本都有唯一的适配接头。将样品混合后进行电泳回收纯化，纯化产物质检合格后使用Illumina NovaSeq 6000 测序平台（上海凌恩生物科技有限公司）进行双末端测序。

测序获得的原始序列数据（Raw read）去除低质量序列，有效的高质量序列（Clean read）数据经BWA v0.7.10 软件[27]比对茶树基因组[28]，利用SAMTOOLS v1.9 软件将sam 转换为bam 文件。GATK v4.1.2.0 软件用于检测SNP，初始SNP 质量过滤标准为“QD＜2.0 ||MQ ＜40.0 || FS ＞60.0 || SOR ＞3.0 ||MQRankSum ＜-12.5 || ReadPosRankSum ＜-8.0”。使用VCF v0.1.11 软件对SNP 进行严格过滤[29]，经过测序深度4X，缺失率（Miss rate）＜10%、次等位基因频率（MAF）＞0.05，多态信息含量（PIC）＞0.15，SNP 位点两侧150 bp 无变异的过滤，最终共获得12 937 679个SNP。从12 937 679 个SNP 中筛选出遗传多样性较高且均匀分布在15 个染色体上的核心SNP（去除重组率为0 的冗余位点、MAF＞0.15、Miss rate＜5%、PIC＞0.15）用于后续分析。使用PowerMaker v3.25 软件[30]计算观察到的杂合度（Observed heterozygosity，Ho）值和基因多样性（Gene diversity，GD）值。

1.3 茶树群体遗传多样性及结构分析

使用Admixture 软件分析种群结构[31]。根据隶属度系数（Q）的取值对所有品种进行划分，根据交叉验证误差（CV 误差）确定最优K 值。主成分分析（PCA）使用Plink 软件进行，使用 R v3.4 进行可视化绘图[32]。利用MEGA 软件的邻接法（ Neighbor-joining method，NJ）构建系统发育树，软件分析均采用默认参数设置，使用iTOL 在线软件编辑系统发育树[33]。

1.4 基因流分析

根据茶树文献资料[34-35]按茶树原产地（以省/直辖市为单位）对349 个茶树进行分组，为探索武夷山地区茶树的历史遗传基因流，进一步将福建省的地区划分为福建武夷山区（FJW）、闽东区（FJE）、闽南区（FJS）和福建未知区（FJU）。为减少误差，少于5 个品种的地区组被排除分析。最后，共采用11地区组进行分析（附表1）。使用TreeMix v1.01软件进行基因流分析[36]，分别设定假设有1～11 次基因流事件（-m=1～11），以箭头表示不同组之间的基因渗入。

1.5 遗传相似度分析

使用软件NTSYSpc v2.11[37]进行遗传相似度分析，计算SNP 的分布频率和遗传相似系 数（ Genetic similarity ， GS） ，GS=NS/(NS+ND)。式中，NS为相同SNP 基因型数量，ND为不同SNP 基因型数量[38]。缺失的基因型被视为无效基因型。

根据国际种子联合会[39]的建议，GS 值大于0.9 可作为实质派生关系的有力证据。因此，GS＞0.9 的品种对被认为是近似品种，具有派生关系。349 个茶树中有6 对来自不同生长地的同一品种/品系，作为明确为实质派生关系的阳性对照，将其中最小的GS 值作为具有无可争议的派生关系的阈值。当茶树样本的GS 值大于该阈值，将这对样本明确为具有派生关系。

1.6 简易鉴定SNP 开发

在973 个核心SNP 的基础上，使用Perl方法以349 个样本成对比较的可辨性为过滤条件，选择一组数量最少且可辨性高的SNP标记用于品种简单快速鉴定[40-41]。在SNP 侧翼两个保守区域设计引物，随机进行一代测序验证。

2 结果与分析

2.1 全基因组SNP 筛选

349 份茶树种质资源共获得829.05 G 高质量测序数据，平均每个样品获得2.38 G 数据，Q20（碱基错误率在1%以下）平均占比98%，Q30（碱基错误率在0.1%以下）平均占比93.05%，平均获得高质量序列条数5 737 108 966，平均有5 668 097 099 条序列可以匹配到茶树参考基因组，平均覆盖率为98.83%，表明测序质量较高。经过严格过滤筛选（Miss rate＜20%、MAF＞10%），共获得12 937 679 个SNP。测序获得的原始数据上传 NCBI 数据库（项目号：PRJNA924950）。以去除重组率为0 的冗余位点、MAF＞0.15、Miss rate＜5%、PIC＞0.15 且较均匀分布在茶树基因组15 条染色体上为筛选条件，共获得973 个核心SNPs。

2.2 茶树群体遗传多样性分析

973 个SNP 较均匀地分布在茶树的15 条染色体上，在第1 至第15 条染色体上的个数分别为94、75、62、73、67、63、70、59、69、60、47、44、64、68 和58（图1A）。在349 个茶树中的Ho 值为0.355，53.04%以上Ho 值大于0.3；GD 平均值为0.319，51.28%以上GD 值大于0.3；茶树群体PIC 平均值为0.262。973 个SNP 在武夷茶树品种/品系中的Ho 值为0.334，GD 平均值为0.302，PIC 平均值为0.248（图1B）。

图1 973 个核心SNP 在茶树基因组15 条染色体上分布情况及多样性分析Fig. 1 Distribution and diversity analysis of 973 core SNPs on 15 chromosomes of tea plant genome

基于973 个SNP 位点对349 个茶树进行群体结构分析。基于模型的群体结构分析表明，5 个群体（K=5）为最佳模型（图2A）。通过PCA 对群体结构进行评价，得出了5 个聚类，与K=5 时的结构分析推断结果一致（图2B）。NJ 系统发育树与种群结构分析及PCA结果基本吻合，但少部分个体在分支聚类上有所差异（图2C）。亚群1～5 分别有茶树个体159、72、51、50 个和17 个。亚群1 主要以福建省武夷山地区茶树、浙江地区茶树及安徽地区茶树为代表，亚群2 主要是以福建省闽南地区铁观音和黄旦等茶树为代表，亚群3 主要是以福鼎大白茶及其相关品种类型为代表，亚群4 主要是以湖南地区茶树为代表，亚群5主要是广东地区茶树为代表。武夷山茶树资源大部分集中于亚群1，共108 个，15 个处于亚群2，1 个处于亚群3，2 个处于亚群4。

2.3 基因流分析

在预估的1～11 次基因流事件中（11 个场景），共有5 个场景涉及武夷山地区茶树的基因流动。其中2 个场景发现武夷山地区茶树群体（FJW）的基因主要向浙江地区茶树群体（ZJ）流动（图3B 和3C），没有发现外地省份的茶树群体直接向武夷山地区茶树群体的基因渗入情况，但3 个场景发现存在本省闽南地区（FJS）茶树群体向武夷山地区（FJW）的基因渗入（图3D～F）。

图3 基因流分析Fig. 3 Gene flow analysis

2.4 遗传相似性分析

基于973 个SNP 位点进行349 个品种之间遗传相似度分析（图4），经过两两比较，GS 在0.524 0～0.996 9，平均值为0.661 8。其中，GS 在0.6～0.7 的最多（90.60%），为55 018 对；GS 在0.8～0.9 的最少，为80 对；GS＞0.9 的共有136 对。126 个武夷茶品种/品系（含6 对已知实质性派生关系的茶树）之间的GS 在0.617 7～0.996 9，平均值为0.686 0。GS＞0.9 的共24 对，在0.8～0.9 的共7 对，在0.7～0.8 的共831 对，在0.6～0.7 的最多（89.05%），共7 013 对。

图4 349 个茶树之间的遗传相似度分析Fig. 4 Genetic similarity analysis among 349 tea resources

349 个茶树中涉及武夷茶品种/品系GS＞0.9 的有26 对（含已知具有派生关系的6 对），被认为是近似品种，具有派生关系，以6 对来自不同生长地的同一品种/品系的两株茶树作为具有明确派生关系的阳性对照，其中GS 最小的为0.972 3，涉及武夷茶品种/品系遗传相似系数大于0.972 3 的共有23 对（含已知具有派生关系的6 对），被明确为具有无争议的派生关系（表1）。其中FJ84（JM051 肉桂）、FJ204（JM046 红海棠）、FJ73（JM037 王母桃）、FJ200（JM035 金鸡母）、FJ191（黄肉桂）5 个品种/品系之间遗传相似系数均大于0.972 3，明确具有派生关系。大红袍2 号株FJ91（JM062-2 大红袍2）与FJ192（奇丹），大红袍4 号株FJ93（JM062-4 大红袍4）与FJ82（JM049 不知春），FJ108（JM081 醉墨）与FJ48（JM005 雀舌），FJ67（JM029 向天梅）与FJ71（JM034 醉贵姬），FJ76（JM042胭脂柳）与FJ77（JM043 醉八仙），FJ13（1113）与FJ12（1114），以及FJ86（JM054 金毛猴）与FJ149（白毛猴）相似系数均大于0.972 3，明确具有派生关系。而FJ52（JM011 石中玉）与FJ161（0205D），FJ63（JM025 石观音）与FJ187（SR），FJ17（0204）与FJ18（0205）相似系数在0.90～0.97，被认为是近似品种/品系，具有一定派生关系。以上品种/品系对之间均可被视为分别属于同一类型。

表1 126 个武夷茶树中遗传相似系数大于0.9 的品种/品系对及派生关系分析Table 1 The information of cultivar/strain pairs with genetic similarity greater than 0.9 and essentially derived relationship analysis in 126 Wuyi tea cultivars

2.5 简易鉴定茶树品种SNP 的筛选

在973 个核心SNPs 基础上，通过349 个样本的成对比较，选择出具有高辨识度的21个SNPs 可用于区分349 个茶树（图5A），其中18 个SNPs 即可用于126 个武夷茶品种/品系建立DNA 指纹图谱（图5B）。进一步在21个SNP 位点两侧保守区域设计引物（表2），使用21 对引物对随机抽取的10 个茶树样本进行PCR 与一代测序验证，分型结果均显示与位点一致，可用于茶树种质资源的快速鉴定（图6）。

图5 349 个茶树（A）及126 个武夷茶（B）的SNP 指纹图谱Fig. 5 SNP fingerprints of 349 tea cultivars (A) and 126 Wuyi tea cultivars (B)

表2 21 对可区分349 个茶树样本的引物信息Table 2 Information of 21 pairs of primers that can distinguish 349 tea resources

3 讨论

3.1 GBS-SNP 及茶树遗传多样性

本研究收集了126 个武夷茶品种/品系及223 个来自其他不同地区的茶树种质资源，共349 个茶树样本，利用GBS-SNPs 方法分析了武夷茶的基因遗传多样性与背景，共鉴定出12 937 679 个高质量SNPs，所有样本均通过质量评估，平均准确率达98.8%。高通量和高质量的SNP 表明，GBS 是进行种群遗传多样性分析的有效方法。本研究进一步筛选出多样性较高（缺失率＜5%、MAF＞0.15、PIC＞0.15）并均匀分布在15 个染色体上的973 个SNP 位点进行遗传多样性分析。973 个SNPs 在349个茶树样本中杂合度Ho 平均值（0.355）和GD 平均值（0.319）均较高，表明这973 个SNP信息丰富，具有较高的基因多样性，可用于遗传多样性分析。973 个SNPs 在349 个茶树与126 个武夷茶树中的PIC 平均值分别为0.262与0.248，茶树群体多态性信息含量中等，这暗示着茶树群体中存在一定数量亲缘关系较近的茶树。以12 937 679 个高质量SNPs 和973个SNPs 构建的群体遗传结构、主成分和系统发育关系均表明，349 个茶树样本可划分为5个种群。这5 个茶树种群主要基于茶树之间的亲缘关系进行聚类，而不是树型或叶形等形态特征，这与以往的研究一致[42]。

图6 随机抽取10 个茶树样本进行一代测序的分型结果示例Fig. 6 The sequencing and genotyping result of the random sampling of 10 tea plant samples

3.2 武夷茶历史遗传背景与基因流

中国西南地区是茶树的发源地[28]，同源茶种以不同的传播演化路径发展。早前有学者根据茶树的生态型和地质变化提出茶树五大演化区[43]，即黔区系群、三江区系群、南岭区系群、武夷区系群及江汉区系群，指出茶树在经历自然传播后，因地质和气候变化在第四纪冰期后形成隔离分布和演化。武夷山地处福建西北，在中生代白垩纪时候由海洋逐渐变为陆地，到第四纪末冰期地壳上升，产生褶皱和断裂成山间谷地。武夷山的地质年代，形成土壤的岩性以及存在野生茶等一系列特点，暗示武夷山茶树是由“隔离分布”的同源演化成“武夷类群”[44]。本研究中基因流分析未发现其他省份茶树直接向武夷山茶树群体的基因流动，其结果符合该说法。此外，基因流分析显示武夷茶主要是向浙江地区茶树基因流动，这与种群结构分析中大部分浙江茶树品种与武夷山茶树聚在一起的结果一致，也与Zhang 等[45]研究结果一致。因此，福建北部武夷山与浙江茶区在地质和演化过程可能为同一演化区，这也暗示着浙江地区的部分茶树种与武夷山茶树种存在一定的基因交流，在近代的茶树育种中，人为的跨省之间的茶树杂交育种非常常见，这些品种是否经由武夷山茶树传播发展而来还有待进一步研究。

在基因流分析中，福建省内存在闽南地区茶树向武夷山茶树的基因流动。一些武夷茶品种/品系，例如紫罗兰、正太阳、金毛猴、白毛猴、金丁香和留兰香等在种群结构分析中也与闽南茶树聚在一起。根据相关文献[46-47]记载，明朝年间，部分闽南茶农大举内迁武夷山地区，部分武夷茶品种/品系可能是由闽南茶流入后经过杂交传播发展而来。此外，近期培育的小红袍茶树因在遗传相似度分析中显示与闽南地区的毛蟹茶树遗传相似系数为0.983 6，具有明显的派生关系，未被列入本研究的126个武夷茶树之中。尽管如此，茶树的原产地溯源分析反映的是茶树品种的遗传背景，以便于后期的茶树种质资源利用与开发。

3.3 茶树EDV 及武夷茶亲缘关系

随着育种技术的发展，针对原始品种进行修饰性育种可能出现大量的派生品种。尽管派生品种可提高原始品种的遗传贡献，但也使育种的遗传基础变窄，不利于作物遗传改良，甚至可能造成对粮食安全的威胁。2022 年3 月1日，新《中华人民共和国种子法》第28 条和第90 条首次规定了EDV 制度[10]，对作物品种的遗传基因相似度提出了要求。茶树为异花授粉植物，其遗传组成高度杂合，表现型上千姿百态，对茶树亲缘关系的分析有助于茶树种质资源的合理开发利用，然而茶树品种的EDV 研究未见报道。本研究首次依据遗传相似度对349 个茶树进行了EDV 初步调查。在349 个茶树的遗传相似度分析中，共发现136对茶树遗传相似系数大于0.9，其中26 对涉及武夷茶，这说明在茶树品种资源的大量引种、频繁的种质杂交以及品种选育过程中过度使用优良品种，已出现一定量的茶树EDV，造成“同物异名”的情况[48-49]。

国际种子联合会（ International seed federation，ISF）[39]将GS 值大于0.9 作为实质派生关系的基本阈值，基于SNP 标记遗传相似度的派生关系阈值的设定建议高于 ISF颁布的标准，且无法对所有作物品种适用统一的EDV 判定标准，必须根据不同作物品种进行确定[50]。茶树品种具体的实质派生关系阈值需要进一步深入研究。实质派生品种除表型与原始品种存在差异外，基因型/基因型组合与原始品种大致相同。考虑到即使同一品种在不同的地方栽培时间较长，也会存在极少量基因突变的积累，若两个表型存在差异的个体的遗传相似系数大于两个不同栽培地的相同品种的遗传相似系数，则可以被确认为具有派生关系。本研究基于遗传相似度分析结果发现，大红袍2 号株（FJ91）与奇丹（FJ192）相似系数为0.988 7，大红袍4 号株（FJ93）与正本不知春（FJ82）相似系数为0.986 7，均具有派生关系，属于同一类型。4 株大红袍茶树与北斗（FJ33）、雀舌（FJ48）的遗传相似系数均较低（0.7 左右），这与先前的研究结果一致[51]。醉墨（FJ108）与雀舌（FJ48）的形态特征与芽叶生长特性十分相似，在叶江华等[17]的ISSR 分子标记鉴定中，相似系数仅为0.82，被认为是不同类型品种。本研究基于SNP 的分子标记显示两者相似系数为0.989 7，应为同一类型的品种，这可能是因为ISSR 分子标记是多等位基因系统，内部存在复杂变异，易产生误差的缘故[52]。留香涧不知春（FJ42）与正本不知春（FJ82）相似系数为0.686 1，为不同类型茶树，这与先前研究一致[17]。此外，本研究还发现，肉桂（FJ84）、红海棠（FJ204）、王母桃（FJ73）、金鸡母（FJ200）和黄肉桂（FJ191）5 个品种/品系之间均具有明确派生关系，属于同一类型。金毛猴（FJ86）与白毛猴（FJ149），向天梅（FJ67）与醉贵姬（FJ71），胭脂柳（FJ76）与醉八仙（FJ77），1113（FJ13）与1114（FJ12）的遗传相似系数均大于0.97，也具有明确派生关系，属于同一类型。石中玉（FJ52）与0205D（FJ161），石观音（FJ63）与SR（FJ187），0204（FJ17）与0205（FJ18）的遗传相似系数在0.90～0.97，均为近似品种，也可认为具有一定派生关系。武夷茶种质资源历史悠久、种类繁多，在引种栽培过程中出现错乱或混杂在所难免。本研究通过SNP 分子标记的遗传相似度分析诠释了武夷茶品种/品系之间的关系，不仅为茶树育种提供了依据，而且对遗传多样性的保护具有重要作用。

3.4 茶树的SNP 分子标记鉴定

SNP 分子标记技术具有基因型信息丰富、检出效率高且准确性高等优势特点，目前已被广泛地应用于作物品种的鉴定[53]。尽管近年来茶树的SNP 分子标记技术进展快速，并有相关的报道[54-56]，但茶树品种资源丰富，基因组庞大且复杂度高，还需要不断地补充完善。本研究筛选出的973 个SNP 位点均匀地覆盖在茶树基因组上，具有较高的基因多样性，可对349 个茶树进行有效区分和指纹图谱构建。数量少且辨别度高的SNP 标记有利于降低构建指纹图谱成本，本研究从973 个位点中筛选21 个SNP 位点作为一组可快速100%区分349个茶树样本的组合，并设计相关引物进行验证，且其中18 个位点即可区分126 个武夷茶品种/品系。这些SNP 标记在未来茶树品种分子鉴定中可起到一定的应用，可作为筛选茶树品种真实性、特异性、纯度鉴定以及类群划分的候选位点组合，为分子标记辅助育种中背景分析提供研究基础。