基于全基因组重测序的泡桐属植物遗传关系分析

赵 阳，冯延芝，杨超伟，王保平，乔 杰，殷世雨，周海江，李芳东

（1.中国林业科学研究院 经济林研究所，河南 郑州 450003；2.经济林种质创新与利用国家林业和草原局重点实验室，河南 郑州 450003；3.信阳市南湾试验林场，河南 信阳 464031）

泡桐Paulownia是玄参科Scrophulariaceae 泡桐属Paulownia树种的统称，是原产我国的重要速生多用途树种，栽培和利用历史悠久[1]。泡桐具有材质优良、生长迅速和繁殖容易等优良特性，在速生丰产林营造、城乡居民区绿化和国民经济建设等领域都发挥了重要作用[2]。我国泡桐种质资源丰富，自然分布范围广，为泡桐遗传改良提供了有利条件[3]。为了高效开发和利用泡桐种质资源，不断提高新一轮次育种群体目标性状的相关基因频率[4]，进一步优化交配设计[5]、实现增效基因的聚合[6]，从而促进种质创新与培育突破性良种，厘清泡桐属植物的分类和系统发育关系是迫切需要解决的问题。

前期已有学者针对泡桐属植物的分类问题开展了很多研究，除了基于形态学分类、数量分类和孢粉学分析等传统方法外，还利用多种分子技术进行分类[7-9]，试图揭示泡桐属植物的系统发育关系。然而，由于泡桐分布广泛且多异花授粉，在漫长的自然迁移、频繁的种间杂交和长期的人工栽培之后，属内种间亲缘关系变得相当复杂，不同研究方法取得的结果也不尽相同。比如，莫文娟等[3]利用ISSR 分子标记分析了泡桐属植物的亲缘关系，并将其分为毛泡桐类（毛泡桐P.tomentosa、兰考泡桐P.elongata、山明泡桐P.lampropylla、白花兰考泡桐P.elongataf.alba、楸叶泡桐P.catalpifolia、圆冠泡桐Paulownia×henanensisC.Y.Zhang et Y.H.Zhao.Hybr.Nov.、宜昌泡桐P.ichangensis和亮叶毛泡桐P.tomentosavar.lucida）、白花泡桐类（白花泡桐P.fortunei、南方泡桐P.australis）和川泡桐类（台湾泡桐P.taiwaniana、川泡桐P.fargesii）3 类；马浩等[10]基于叶绿体DNA 的RFLP 分析，将泡桐属15 个种分为南方泡桐组（南方泡桐和成都泡桐）、毛泡桐组（毛泡桐和兰考泡桐）和白花泡桐组（白花泡桐、山明泡桐、海岛泡桐、台湾泡桐、建始泡桐、宜昌泡桐、兴山泡桐、川泡桐、楸叶泡桐和白花兰考泡桐）3 组；而卢龙斗等[11]基于RAPD 分析认为白花泡桐和川泡桐的亲缘关系在泡桐属的7个桐种中最为接近。可以看出，前人的研究结果不尽相同，仍不能很好地解决泡桐属内部深层次的系统发育关系。

全基因组重测序（Whole genome re-sequencing）是基于DNA 分子水平对已知基因组序列的物种的不同个体进行高通量测序，可以通过将测试品种的序列与已知参考序列进行比对，确定测试品种遗传构成的差异[12]。单核苷酸多态性（Singlenucleotide polymorphisms，SNP）是指基因组中某个特定位置上发生单个核苷酸变异所引起的DNA 序列多态性，与ISSR（Inter-simple sequence repeat）、RFLP（Restriction fragment length polymorphism）、RAPD（Random amplified polymorphic DNA）和SCAR（Sequence characterized amplified region）等分子标记相比，具有数量多、分布密度高、快速且高通量分型等优点[13-14]。随着测序技术的快速发展、拥有参考基因组物种数量的不断增多以及基因分型成本的下降，SNP 标记日渐广泛地被应用到植物群体遗传学研究中[15-18]。白花泡桐全基因组测序的完成推动了SNP 标记在泡桐研究中的应用。李文杨等[19]基于筛选到的SNP 位点对由毛泡桐和白花泡桐为亲本，连同正反交获得的181 个F1代个体构成的作图群体进行基因型分析，构建了高密度的连锁遗传图谱，为泡桐分子育种研究提供了数据参考。但是，尚未见到有关SNP 标记在泡桐系统发育分析和分类地位研究的报道。因此，本研究在对泡桐属11 个桐种的典型株进行全基因组重测序的基础上开发SNP 标记，解析泡桐属内的亲缘关系和遗传结构，以期为泡桐核心育种群体构建、骨干亲本筛选和种质创新提供依据。

1 材料与方法

1.1 试验材料

本试验选用的11 份材料均采自中国林科院经济林研究所所属、位于湖北省钟祥市盘石岭林场的泡桐种质资源库，涵盖了泡桐属所有种的典型株。采集叶片后用去离子水清洗，放入液氮中进行速冻，置于冰箱中-80 ℃下保存备用。试验材料详细信息见表1。

表1 供试试验材料Table 1 Materials used in the experiment

1.2 试验方法

1.2.1 DNA 提取、建库和测序

采用改良的CTAB 法提取叶片样品的基因组DNA[20]，经1%琼脂糖凝胶电泳和紫外分光光度计检测质量和浓度后选取合格的DNA 样品进行GBS 文库构建。将DNA 用MseI 限制性内切酶随机打断成短的DNA 片段后进行末端修复，然后在DNA 片段两端连接dA 尾并连接测序接头。对加上接头的DNA 片段经过AMPure XP 磁珠纯化后选择300～400 bp 范围的片段进行PCR 扩增和电泳回收纯化，最后在Illumina HiSeq 4000 平台上机进行双末端PE150 测序。

1.2.2 泡桐基因组SNPs 挖掘

测序获得原始数据（Raw data）后，去除含接头序列的reads、单端测序序列中未检出的碱基超过该序列长度10%的序列和低质量（≤5）碱基数超过该序列长度50%的序列后，得到高质量的测序数据（Clean reads）。

以白花泡桐为参考基因组（Paulowniafortunei,https://www.ncbi.nlm.nih.gov/genome/?term=paulownia），采用BWA（Burrows-wheeler-alignment tool）软件[21]将Clean reads 比对到参考基因组，然后使用picard 对比对结果标记并去除重复序列。采用GATK 软件进行群体SNP 检测，选择贝叶斯模型检测群体中的多态性位点，通过有效过滤得到高质量SNPs 结果用于后续分析[22]。

1.2.3 数据分析

基于检测获得的SNP，采用VCF2Dis 软件计算遗传距离矩阵，采用MEGA11 软件中邻近法（Neighbor-joining method）构建进化树，采用FigTree 软件对进化树进行可视化。采用plink 软件计算特征向量及特征值[23]，然后使用Seaborn 可视化库绘制PCA 分布图，最后通过fastSTRUCTURE的默认参数筛选K值并推断种群结构[24]。以每kb的变异数代表不同基因编码序列（CDS）的突变率，采用四分位法将基因分组，采用clusterProfiler 包对各组基因进行富集分析。

2 结果与分析

2.1 测序质量

采用Illumina 测序平台对泡桐属11 个种的典型株进行全基因组重测序，过滤掉不符合要求（含有接头和低质量的reads）的数据后，共获得高质量测序数据63.08 GB，平均每个样本获得5.87 GB 的数据。经过与参考基因组比对，碱基错误率在1%以下（Q20）的占比95.75%～96.94%，平均为96.54%；错误率低于0.1%（Q30）的占比89.13%～92.00%，平均为90.97%。GC 含量32.98%～38.99%，平均为36.69%。与参考基因组的比对率为72.54%～95.26%，平均为81.87%，测序深度最高为9.69×，最低为6.97×，平均测序深度为8.42×，＞1×碱基覆盖度平均为95.26%。上述结果表明，各样本与参考基因组的相似度满足重测序分析要求，可进行后续分析。各样本数据见表2。

表2 重测序数据概况Table 2 Overview of the whole genome resequencing data

2.2 SNP 类型与分布

采用GATK 进行SNP 检测，过滤后最终获得7 492 966 个高质量SNPs 用于后续分析。进一步统计上述SNP的变异类型，发现在6种变异类型中，A/G 和C/T 变异类型最多，分别占比32.00%和31.96%（图1A），即4 791 902 个SNPs 发生了碱基转换（Ts），2 701 064 个SNPs 发生了碱基颠换（Tv），Ts/Tv=1.78。进一步统计分析他们在基因组的分布情况，发现有57.8%的SNP 位于基因间区，20.9%的SNP 位于内含子区，6.4%的SNP 位于外显子区，基因上游区域和下游区域内的SNP分别占比5.7%和5.0%（图1B）。

图1 泡桐属11 个桐种变异位点及注释信息的统计Fig.1 Statistics of location and annotation for the SNP and InDel variants of 11 Paulownia samples

2.3 SNP 的种间变异

不同品种泡桐个体中检测到的总变异数各不相同，毛泡桐具有最多的变异数（4 900 517）（图1C），和较低的碱基转换/碱基颠换比值（TS/TV=1.77，在11 种泡桐中仅高于台湾泡桐和川泡桐），表明毛泡桐的基因交换事件发生的较少。不同种泡桐非同义突变与同义突变的比值（Nonsyn/Syn）相差不大，白花泡桐最高，达到1.03，其他桐种为0.94～0.95（表3，图1D）。

表3 不同桐种的变异情况Table 3 Summary of variants from all Paulownia samples

2.4 群体进化树与主成分分析

基于SNP 标记，采用VCF2Dis 计算各桐种间的遗传距离，结果见图2。泡桐属植物种间的遗传距离为0.15～0.59，其中，毛泡桐与白花泡桐、台湾泡桐与白花泡桐的遗传距离分别为0.59 和0.58，属于亲缘关系较远的2 组桐种。宜昌泡桐与山明泡桐，鄂川泡桐与山明泡桐的遗传距离均为0.15，属于亲缘关系较近的2 组桐种。根据种间遗传距离，构建泡桐属植物的系统发育树，结果表明11 个桐种可以被分为3 大支（图3）。其中，台湾泡桐和川泡桐为一支（Ⅰ），毛泡桐单独为一支（Ⅱ），白花泡桐、楸叶泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、华东泡桐、建始泡桐和兰考泡桐为一支（Ⅲ）。通过主成分分析方法再次分析11 个桐种间的亲缘关系，根据主成分1 和主成分2 同样将11 份材料分为3 个类群：毛泡桐为一个类群，台湾泡桐和川泡桐为一个类群，其他桐种为一个类群（图4），结果与系统发育分析结果高度一致。这表明，本研究鉴定到的SNPs 在泡桐属植物亲缘关系评价中可靠性较高。

图2 基于SNP 计算的11 个桐种的遗传距离Fig.2 Genetic distances of 11 Paulownia samples based on SNPs

图3 基于SNP 构建的11 个桐种的Neighbor Joining 进化树Fig.3 Neighbor Joining phylogenetic tree of 11 Paulownia samples based on SNPs

图4 基于SNP 构建的11 个桐种的主成分分析Fig.4 Principal component analysis of 11 Paulownia samples based on SNPs

2.5 群体遗传结构分析

群体遗传结构分析是遗传关系分析的一种重要手段[18]。为了进一步解析11 个桐种间的亲缘关系，采用FastSTRUCTURE 软件对他们进行群体遗传结构分析。通过假定不同的祖先来源数目K，可以将不同群体之间是否含有相同祖先成分表示出来。在Structure 群体遗传结构图中，每一列竖格代表一个样本的遗传背景，每种颜色的色块代表估计的一个祖先，每种色块占该竖格的比例代表该祖先对该样品遗传背景的贡献比例[25]。为了确定合适的分群数量，假定祖群值K为1～9，将不同群体之间是否含有相同的祖先成分表示出来。fastSTRUCTURE 中chooseK 程序可根据种群结构边际似然近似值（Marginal likelihood）自动筛选最佳K值，结果显示，3 是边际似然值达到最大时的最小K值，即最合理的祖先种群数为3，表明11 个桐种可能来自于3 个祖先（图5）。该结果与系统发育树和主成分分析结果相互吻合，即来源相同的3 个类群可以独立成群，类群Ⅰ包括毛泡桐；类群Ⅱ包括川泡桐和台湾泡桐；类群Ⅲ包括白花泡桐、楸叶泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、华东泡桐、建始泡桐和兰考泡桐。值得一提的是，类群Ⅲ中的华东泡桐显示了来自类群Ⅱ的混合遗传成分，说明类群Ⅱ和类群Ⅲ可能存在基因交流（图6）。

图5 基于fastSTRUCTURE 计算不同K 值时种群结构的边际似然值Fig.5 Marginal likelihood value of Paulownia samples with different K values calculated by fastSTRUCTURE

图6 基于fastSTRUCTURE 推断的11 个桐种的群体分层Fig.6 Population stratification based on fastSTRUCTURE analysis of 11 Paulownia samples

2.6 SNP 相关基因的功能分析

蛋白质编码序列的突变可能会影响基因功能和不同的生物学过程。采用变异密度法（每kb变异数）计算泡桐属植物蛋白质编码基因的突变率（Mutation rates，MR）[26]，根据突变率分布的四分位距，这些基因被分为高、较高、较低和低突变率4 组。高突变率组包括4 597 个基因（MR ≥6.0），较高突变率组包括4 033 个基因（6.0 ＞MR ≥4.0），较低突变率组包括8 395 个基因（4.0 ＞MR ≥2.0），低突变率组包括9 878个基因（MR ＜2.0）。进一步GO 富集分析表明，进化速度不同的基因，其功能类别也存在差异。高突变率的基因被富集到细胞壁代谢途径，包括木葡聚糖代谢调控、次生细胞壁和木栓素生物合成过程相关基因，如β-己糖胺酶、葡聚糖半乳糖基转移酶、GMA12/MNN10 家族半乳糖基转移酶和UDP-糖基转移酶等33 个基因；花粉代谢途径，包括花粉壁组件、花粉管和花药发育相关基因，如ABORTED 转录因子、ABC 转运蛋白和PHDfinger 等31 个基因；次生代谢物代谢通路，包括苯丙烷生物合成过程、次生代谢物生物合成过程和次生代谢过程相关基因，如苯丙氨酸解氨酶、GST 转运蛋白和AP2 等35 个基因；形态建成相关通路，包括植物形态器官发生调控、幼苗发育调控和侧根形成相关基因，如铵转运蛋白、ABC-2型转运蛋白和MYB 等53 个基因（图7）。

图7 高突变率基因集中所包含细胞壁合成、花粉、次生代谢和形态建成相关基因的GO 富集分析Fig.7 GO enrichment analysis for the genes involved in the pathways of cell way,pollen,secondary metabolite and morphogenesis with high mutation rates (MR)

3 讨 论

我国泡桐种质资源丰富，但很多材料缺乏可追溯的谱系，遗传背景不明，给泡桐属植物核心育种群体的建立、亲本选配和种质创新等带来诸多不便。李芳东等[2]在对全国泡桐属种质资源调查收集过程中，从形态特征、生态习性、分布状况和性状稳定性等方面对泡桐属植物进行了全面比较分析，提出泡桐属主要包括11 个桐种的认识（表1），推动了泡桐的遗传演化及分类学研究。

本研究共鉴定到7 492 966 个高质量SNPs，其中碱基转换与颠换比为1.78，说明碱基转换是泡桐属植物主要的SNP 变异类型，该结果与茶树Camelliasinensis(L.) Kuntze[27]、葡萄Vitis[28]和毛竹Phyllostachysedulis[29]中SNP 的鉴定结果基本一致。研究表明，碱基转换的SNP 对蛋白质氨基酸序列的改变小于碱基颠换的SNP[30]。本研究鉴定的大多数SNP 分布在基因组的基因间区和非编码区，该结果与咖啡和茶树中SNP 的注释结果相似[31-32]。虽然这些SNP 可能不影响基因功能，但却能够体现植物群体间的遗传和进化关系[18]，作为遗传标记应用于泡桐属植物遗传关系分析。

系统进化树、主成分分析以及遗传结构分析的结果高度一致，均支持3 个亚群的分类模式。其中，毛泡桐为类群Ⅰ；川泡桐和台湾泡桐为类群Ⅱ；白花泡桐、楸叶泡桐、山明泡桐、鄂川泡桐、宜昌泡桐、华东泡桐、建始泡桐和兰考泡桐为类群Ⅲ，该结果与范永明[33]基于花序形状、花蕾形状和花色等22 个性状采用数量分类学方法取得的结果基本一致。本研究还发现，地理位置相近的桐种并没有因此聚到同一类群，说明3 个类群可能的基因来源与地理位置没有明确的相关性[34]。

本研究中，高突变率基因（MR ≥6.0）被富集到细胞壁代谢、花粉代谢、次生代谢物代谢和形态建成等生物学过程，共同作用导致了不同桐种在遗传、基因功能和表型上的多样性。比如，木葡聚糖代谢调控、次生细胞壁和木栓素生物合成过程等相关基因影响细胞壁的代谢过程。树干形成层向内分化形成具有生活力的木质部细胞，失活后残留的细胞壁是影响木材品质的重要成分[35]。因此，这些SNP 可能在木材形成过程中相关基因的功能调控中发挥作用，需要后续进一步研究。再如，花粉壁组件、花粉管和花药发育相关基因影响花粉的代谢过程，这些SNP 可能是泡桐属在长期进化过程中出现楸叶泡桐和山明泡桐2 个雄性不育种[1]的诱因之一。楸叶泡桐是泡桐属11 个种中唯一兼具适宜密植和材质优良特性的桐种，是培育装饰材的首选桐种[36]，雄性不育致使其只能用作母本进行杂交育种，限制了杂种优势的利用。对此类关键SNP 位点的筛选和鉴定可为通过分子育种手段恢复楸叶楸叶泡桐的育性提供参考。

本研究以泡桐属11 个桐种的典型单株为对象，利用WGRS 技术挖掘覆盖全基因组的SNP 标记，对解析泡桐属的种间亲缘关系及遗传结构具有重要意义，但也存在一些局限，下一步还需要增加种内不同基因型或不同种源的试验材料参与分析，开展更深入的理论研究，更有效地支撑泡桐种质创新工作。

4 结 论

本研究首次报道了泡桐属11 个桐种间的全基因组遗传变异。通过全基因组重测序鉴定出7 492 966 个高质量SNPs。这些SNP 不仅有助于进行泡桐种质资源遗传变异分析和起源分化探索，同时还对鉴定具有表型效应的等位基因变异具有重要意义。本研究基于这些SNP 对泡桐属植物的亲缘关系和遗传结构进行了系统分析，将为泡桐种质创新、挖掘优良性状相关基因进而培育突破性良种奠定基础。