马尾松树高生长性状显著关联基因挖掘

龚桂芳，冯源恒，罗群凤，杨章旗

（广西壮族自治区林业科学研究院 广西优良用材林资源培育重点实验室 国家林业和草原局马尾松工程技术研究中心 广西马尾松工程技术研究中心，广西南宁 530002）

马尾松（Pinus massoniana）是我国南方主要的用材树种之一，其自然分布广泛，遗传变异丰富，环境适应性强，曾作为荒山绿化的首选树种在南方各省区被大量种植。第八次森林资源清查数据表明，我国马尾松林总面积为1 001万hm2，蓄积量5．91亿m3，其中人工林面积307 万hm2，蓄积量1．72 亿m3[1]。马尾松也是我国最早一批开展遗传改良研究的树种，早在1958年就开始进行马尾松种源试验[2]。1980年后，马尾松遗传育种被正式列入国家重点科研项目[3]。进入21 世纪后，各主要产区的马尾松遗传改良研究先后进入第2 或第3 个轮回育种阶段[4-7]，但较长的育种周期制约了马尾松的遗传改良进程。通过分子标记进行辅助育种成为当前马尾松育种研究的重点之一。

随着马尾松分子遗传学研究的深入，可将分子标记辅助育种分为两个阶段[8]。第一阶段，以利用SSR 分子标记技术分析育种群体的遗传结构、估算亲缘关系、对选自半同胞家系的优树进行亲本分析为主要内容，并将之作为划分育种群体、研究选择方法的重要依据[4,8-10]。第二阶段，以基因关联分析为主要技术手段，开发SSR、SNPs 等标记，与表型性状进行关联分析[11-12]，从而得到高度关联的分子标记，开展早期选育和基因功能分析。

本研究基于马尾松候选基因组关联分析获得的9 个与树高生长性状显著相关的SSR 分子标记，对其所在的基因组序列区域进行分析，以获得与马尾松树高生长性状相关的重要基因。

1 材料与方法

1．1 显著关联SSR位点的获得

研究采用的9个与马尾松树高生长性状显著相关的SSR分子位点来自研究团队前期进行的候选基因组关联分析研究结果。前期研究对两个生长及产脂量性状存在差异的马尾松无性系的顶芽组织、针叶及树干韧皮部进行转录组差异分析，得到差异表达基因。其目的是从顶芽组织中获得与抽梢及纵向生长相关的候选基因，从针叶中获得与光合效率相关的候选基因，从树干韧皮部获得与径向生长及树脂分泌相关的候选基因。在获得的差异表达序列中，设计开发259 对EST-SSR 引物[11]，从中选出65对SSR引物进行候选基因组关联分析[13]。关联分析试验群体由320株1994年造林的马尾松1代种子园自由授粉获得的子代组成，来自106个家系，在同一个随机交配系统下产生，最大母本贡献率为2．5%，不存在单一亲本贡献率过大的问题。2016年12月进行生长量测定，获得树高表型数据，与65 对SSR 引物进行候选基因组关联分析。以P＜0．05 作为标记与树高性状存在连锁不平衡的标准，共计获得9 个与马尾松树高生长性状显著相关的SSR 分子位点，平均表型变异解释率1．42%（表1）。

1．2 显著关联基因的挖掘方法

1．2．1 显著关联位点所在基因序列的挖掘

根据该批次马尾松EST-SSR 引物[11]的设计档案，查找出9 个SSR 位点所在的序列ID 号，获得9 个SSR 位点所在的转录组第2 代测序基因序列。开展第3 代全长转录组测序，通过序列检索、比对，获得上述基因的全长序列。

1．2．2 显著关联基因的功能注释

为获得全面的基因功能信息，对得到的基因序列进行7 个数据库的基因功能注释，包括：Nr（NCBI non-redundant protein sequences）、Nt（NCBI non-redundant nucleotide sequences）、Pfam（Protein family）、KOG/COG（Clusters of Orthologous Groups of proteins）、Swiss-Prot（A manually annotated and reviewed protein sequence database）、KO（KEGG Ortholog database）和GO（Gene Ontology）。

2 结果与分析

2．1 显著关联基因序列挖掘

根据序列ID 号，通过Novofinder 软件在测序数据中检索获得SSR 位点所在的转录组第2 代测序基因序列。9 个基因序列平均长度为1 867 bp，其中PCZ023所在序列最短（910 bp），PCZ129所在序列最长（3 225 bp）（表2）。

将获得的基因序列进行检索和比对，以获得完整的基因序列。分析结果表明，9个基因均在第3代全长转录组测序结果中比对到相应序列（表3）。9个基因序列平均长度为1 706 bp，比第2代测序结果略短。其中PCZ023 所在基因序列最短（761 bp），PCZ157 所在基因序列最长（3 831 bp）。9 个基因的第2 代测序结果与第3 代全长转录组测序结果比对一致性平均为99．69%，均具有高度的一致性。

表3 马尾松树高性状显著关联基因序列比对情况Tab.3 Sequence alignment of genes significantly associated with tree height traits of P.massoniana

2．2 显著关联基因的功能注释

为进一步分析与树高性状显著关联基因的功能，将得到的9 个基因在7 个数据库中进行基因功能注释分析。功能注释结果为PCZ002 所在的基因在云杉（Picea asperata）基因组中存在相似序列，其功能属于转录因子类编码基因；PCZ023所在的基因在云杉基因组中存在相似序列，其功能属于60S 核糖体蛋白大亚基编码基因；PCZ090所在的基因在油松（Picea tabuliformis）基因组中存在相似序列，其功能属于TCHQD 类谷胱甘肽S-转移酶编码基因；PCZ099所在的基因在云杉基因组中存在相似序列，其功能属于ATP酶编码基因；PCZ129所在的基因在白云杉（Picea glauca）基因组中存在相似序列，其功能属于泛核蛋白编码基因；PCZ142所在基因在白云杉基因组中存在相似序列，其功能属于60S 核糖体蛋白大亚基编码基因；PCZ157所在基因在白云杉基因组中存在相似序列，其功能属于转录因子类编码基因；PCZ187所在的基因在云杉基因组中存在相似序列，其功能属于氧化还原电子传递链酶类编码基因；PCZ187所在的基因在云杉基因组中存在相似序列，其功能未知。

3 讨论与结论

基因组关联分析是基于基因的连锁不平衡原理，将基因型与观测表型进行群体水平的统计学分析，根据统计量或显著性P值筛选出最有可能影响该性状的分子标记位点，挖掘与性状变异相关基因的一种研究方法。所得的标记位点极有可能与性状变异相关的基因是强度连锁，甚至处于该基因序列上。本研究基于候选基因组关联分析研究结果，采用的分子标记来自候选基因。通过该方法获得的与树高性状显著关联的标记位点极有可能处在控制该性状变异的基因上。通过挖掘标记位点所在的基因，得到控制树高生长主效基因的几率比采用全基因组关联分析与简化基因组关联分析方法更高。

本研究中，挖掘的9 个与树高性状显著关联的基因均在7个数据库中检索到高度同源的基因。其中，5 个在云杉基因组中发现同源基因，3 个在白云杉基因组中发现同源基因，1 个在油松基因组中发现同源基因，说明上述基因可能是松属植物特有的基因家族类型。

对马尾松半双列杂交家系遗传测定试验的分析结果表明，马尾松树高性状加性效应高于显性效应[14]，说明马尾松树高性状是典型的数量性状，由众多的基因甚至基因家族共同控制。树木的光合作用、呼吸作用、激素调控、水分及营养元素吸收和抗逆性等都会对树高生长产生重要影响。本研究中，有8 个基因在数据库中获得功能注释。其中PCZ002 与PCZ157 所在基因属于转录因子类，PCZ023 与PCZ142 所在基因属于60S 核糖体蛋白大亚基编码基因类，PCZ099 与PCZ187 所在基因可能参与了能量代谢及跨膜运输，PCZ090所在基因属于TCHQD 类谷胱甘肽S-转移酶编码基因，广泛参与植物体内解毒及抗逆境胁迫等功能[15-17]，说明参与树高生长过程的基因种类多样。值得关注的是，这些基因多与维持植物细胞基本功能的基因表达与代谢功能相关，而非预想的与细胞分裂、植物激素合成相关。由此推测，植物细胞基本代谢功能旺盛是保证植物高生长的原动力。

本研究采用通过第3代测序技术开展转录组测序获得基因全长的方法。第3代测序技术实现DNA聚合酶内在自身的延续性，一个反应就可以测非常长的序列。第2代测序仅可测上百个碱基，但第3代测序可测几千个碱基，并可对RNA 进行直接测序，大大降低体外逆转录产生的系统误差，且精度非常高，达到99．999 9%，克服了第2 代测序中因拼接过多造成的错误。与传统的通过RACE 技术等克隆获得基因全长的方法相比较，具有一次性获得海量基因的完整序列明显优势。本研究中的9 个基因，其第2代测序结果与第3代测序结果均存在长度差异，这可能是因为第2代测序以其他模式物种为模板进行拼接造成的拼接错误，也可能是因为RNA 序列自身在转录过程中发生剪切、拼接等。