马尾松生长性状显著关联标记PCZ090所在基因的挖掘与生物信息学分析

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

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

马尾松（Pinus massoniana）广泛分布于我国南方17个省（区），是我国南方特有的多用途重要工业用材树种[1-2]，具有分布面积广、生长迅速、蓄积量大、经济价值高、用途广和适应性强等特点[3]。最新森林资源清查结果显示，我国马尾松分布区域达200万km2，其中林分面积为804．30万hm2，占全国乔木林面积的4．47%，林分总面积居全国乔木树种首位，蓄积量居第六位[4]。马尾松生长迅速，可在较短的经营期内提供大量木材，是许多木材工业、林产工业和造纸工业的支柱，同时也是荒山绿化和生态林建设的先锋树种[5]。通过选育和推广马尾松良种，可提高单位面积木材产量和质量，缓解木材供需矛盾、保障生态安全及增加林业产业效益[6]。因此加快高材性状的选育已成为马尾松遗传改良的重要内容。

从20 世纪80年代开始，广西壮族自治区林业科学研究院松树课题组就开始了马尾松生长性状遗传改良研究工作。至今已选育高材马尾松将近2 000 份，建立了5 个马尾松2 代种子园，3 个马尾松高世代育种园。马尾松的育种周期相对较长，通过传统的遗传育种方法，短时间内无法做到准确快速选育，难以实现良种的快速更新，限制了马尾松的遗传改良进程，制约了森林生产力的提高速率。因此，为更快选育大量马尾松高材良种，满足林业生产的需求，开发研究与生长性状具有较强相关性的遗传标记，结合分子标记辅助选择，是加快马尾松生长性状改良进程的一种可行方法。通过分子标记辅助选择开展早期选择，可有效缩短育种周期，提高选育效率。建立马尾松分子标记辅助选择技术体系对于高效开展马尾松遗传改良意义重大。

课题组前期进行了马尾松主要经济性状候选基因组关联分析研究，已开发获得与生长性状显著关联SSR 标记14 个。本研究以其中一个位点（PCZ090）为目标，基于高通量转录组测序结果对其所在基因序列进行追溯，测序验证获得其基因全长，并通过生物信息学软件分析PCZ090 基因序列，包括同源序列比对、蛋白结构（一级、二级和三级结构）及结构域和功能位点预测等。对生长性状显著关联标记所在基因进行深度挖掘，是发现控制目标性状重要基因的有效途径，对于深层次探究马尾松生长性状遗传调控机理具有重要意义。本研究可为进一步揭示马尾松生长性状相关基因的作用效应、深入研究马尾松主要生长性状遗传机理、分子标记辅助育种技术在其他林木上的应用提供参考。

1 材料与方法

1．1 马尾松生长性状显著关联SSR位点的获得

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

1．2 马尾松PCZ090标记所在基因序列的追溯

胸径和材积均通过3 种线性模型进行分析，分别为EMMMAX、GAPIT 和GLM，计算得到P值，胸径P值分别为0．012 368、0．002 769和0．000 139；材积P值分别为0．006 429、0．002 796 和6．27E-07；树高只通过GLM 模型进行分析，计算得到P值为0．043 690（表1）。3 种生长性状的P值都小于0．05，说明该标记与生长性状显著关联，符合要求。

表1 SSR关联分析结果Tab.1 SSR association analysis results

PCZ090 标记是根据马尾松转录组测序结果开发所得的EST-SSR 标记。从测序结果中可追溯检索到PCZ090标记的原序列Cluster-7694．28796。

1．3 引物来源及PCR反应条件

根据PCZ090标记的原序列Cluster-7694．28796，设计全长扩增引物PCZ090-F（5′-TTTCCTCCCGAAGTAATTTGC-3′）、PCZ090-R（5′-CACCATGATGGATAAAACCGAA-3′），扩增PCZ090标记的cDNA全长序列，最后将目的片段进行克隆测序。

反应采用50 μL 体系进行：TaKaRa LATaq（5 U/μL）0．5 μL，10 × LA PCR Buffer Ⅱ（Mg2+Free）5．0 μL，MgCl2（25 mmol/L）5．0 μL，dNTP Mixture（2．5 mmol/L each）8．0 μL，Forward Primer（10 μmol/L）2．0 μL，Reverse Primer（10 μmol/L）2．0 μL，cDNA template 1．0 μL，Milli-Q Water 26．5 μL。PCR 反应条件为94 ℃3 min；94 ℃30 s，55 ℃30 s，72 ℃2 min，35 cycles；72 ℃延伸10 min，结束PCR 反应。将获得的目的条带送至上海生工生物工程公司进行测序。

1．4 生物信息学分析

采用ORF Finding（www．ncbi．nlm．nih．gov/orffinder）在线软件进行ORF 确定及氨基酸序列预测；采用在线工具Expasy ProtParam（https://web．expasy．org/protparam/）分析蛋白的分子量、理论等电点和氨基酸组成等；采用InterProScan（https://www．ebi．ac．uk/interpro/）、NCBI CDD（http://www．ncbi．nlm．nih．gov/Structure/cdd/wrpsb．cgi）在线工具软件进行结构域和功能位点预测；采用ProtScale（https://web．expasy．org/protscale/）在线软件绘制蛋白质的亲疏水性系列谱；采用SignalP（https://services．healthtech．dtu．dk/service．php? SignalP-5．0）在线软件预测信号肽；采用TMHMM Server v．2．0（https://services．healthtech．dtu．dk/service．php?TMHMM-2．0）程序进行蛋白序列跨膜区分析；采用SOPMA（https://npsa-prabi．ibcp．fr/cgibin/npsa_automat．pl?page=npsa_sopma．html）在线 软件预测分析蛋白质二级结构[8]；采用SWISS-MODEL（https://swissmodel．expasy．org/interactive）在 线软 件对蛋白质三级结构进行建模[9-11]；采用NCBI BLAST（https://blast．ncbi．nlm．nih．gov/Blast．cgi）中 的Nucleotide BLAST 进行同源比对；采用MEGA 5．0 构建系统进化树，选择邻接法，Bootstrap 检验使用1 000 次重复。

2 结果与分析

2．1 马尾松PCZ090标记相关基因序列分析

基因全长1 636 bp，其中5′非编码区有432 bp，3′非编码区有409 bp，开放阅读框（ORF）序列为795 bp（433 ～1 227 bp），编码264 个氨基酸，起始密码子为ATG（黄色标注），终止密码子为TGA（黄色标注）（图1）。

图1 PmGSTTCHQD1基因ORF结果Fig.1 ORF result of PmGSTTCHQD1 gene

PCZ090 所在基因编码的蛋白质分子式为C1411H2183N377O387S8，编码蛋白的分子量约为30．88 kDa，理论等电点为9．27，属于碱性蛋白，原子总数为4 366。20 种氨基酸组成中，亮氨酸（Leu）数量最多（13．6%）；半胱氨酸（Cys）数量最少（0．4%）。酸性氨基酸（Asp+Glu）总数为29 个，碱性氨基酸（Arg+Lys）总数为36 个，说明该蛋白带弱正电荷。不稳定指数为46．23，该蛋白为不稳定蛋白；脂肪指数为94．28，亲水性总平均值为-0．333。

2．2 编码蛋白产物的功能分析和亲疏水性比较

根据InterProScan 在线软件分析基因编码产物的结构域和功能位点，发现PCZ090 所在基因具有典型的谷胱甘肽S-转移酶的保守结构域：GST-N 末端结构域（从Met1-Gly80）、GST-C 末端结构域（从Thr148-Pro229）（图2），参与谷胱甘肽代谢过程，具有谷胱甘肽转移酶活性。通过NCBI CDD 进一步确定PCZ090 所在基因是谷胱甘肽S-转移酶家族蛋白成员，位于1 ～234 位氨基酸，包含N 端结构域和C 端结构域，可催化还原谷胱甘肽与多种内源性和外源性烷化剂的结合，包括致癌物、治疗药物、环境毒素和氧化应激产物。依据GST 的系统命名体系，将该基因命名为PmGSTTCHQD1。

图2 PmGSTTCHQD1保守结构域Fig.2 Conserved domain of PmGSTTCHQD1

PmGSTTCHQD1 蛋白的亲水性氨基酸残基明显多于疏水性氨基酸残基，属于亲水性蛋白质，亲水指数为-2．467 ～2．211；亲水区域的最高值（Score = -2．467）出现在第153 个氨基酸残基，疏水区域的最高值（Score=2．211）出现在第65 个氨基酸残基（图3）。

图3 PmGSTTCHQD1基因编码蛋白的亲疏水性序列谱Fig.3 Hydropathy profile of protein coded by PmGSTTCHQD1 gene

通过SignalP 进行信号肽预测。结果表明，马尾松PmGSTTCHQD1基因不存在信号肽酶切位点，是非分泌型蛋白。TMHMM 2．0蛋白质跨膜结构分析结果显示，该基因所编码的蛋白没有跨膜结构（图4）。

图4 PmGSTTCHQD1蛋白跨膜结构Fig.4 Transmembrane structure of PmGSTTCHQD1 protein

2．3 编码蛋白质的二级结构预测

SOPMA 分析表明，PmGSTTCHQD1基因编码的蛋白质二级结构含有57．20%（151 个）的α 螺旋（Alpha helix）、5．68%（15 个）的β 转 角（Beta turn）、10．98%（29个）的延伸链（Extended strand）和26．14%（69个）的无规则卷曲（Random coil）（图5）。其中，α螺旋和无规则卷曲是主要成分。

图5 PmGSTTCHQD1蛋白二级结构预测Fig.5 Prediction of secondary structure of PmGSTTCHQD1 protein

2．4 编码蛋白质的三级结构预测

蛋白质三级结构预测结果显示，PmGSTTCHQD1基因编码的蛋白质以SWISS-MODEL 上来自毛果杨（Populus trichocarpa）的谷胱甘肽转移酶F2 的晶体结构（编号5ey6．1．B）模型作为模板，与模板有较高相似性的是1 ～231位氨基酸，模板覆盖率为81%，两者的序列具有33%的一致性（图6）。

图6 PmGSTTCHQD1蛋白三级结构模型Fig.6 Predicted tertiary structure model of PmGSTTCHQD1 protein

2．5 同源性与系统进化分析

由NCBI 核酸数据库Blast 比对结果可知，PmGSTTCHQD1基因与油松（Pinus tabuliformis）TCHQD类谷胱甘肽S-转移酶（JX962828．1）和日本落叶松（Larix kaempferi）TCHQD 类谷胱甘肽S-转移酶（KF438050．1）的完整mRNA 同源性分别为99．75%和93．09%，覆盖率均为48%；与白云杉（Picea glauca）GQ02824_L21（BT105585．1）的mRNA 同源性为91．01%，覆盖率为49%；与火炬松（Pinus taeda）两个分 离 株5639 和5637（FJ066498．1 和FJ066514．1）的基因组序列同源性分别为99．31%和99．14%，覆盖率分别为26%和21%。

MEGA5．0 系统进化结果表明，PmGSTTCHQD1蛋白与油松和落叶松TCHQD 类谷胱甘肽S-转移酶蛋白关系最近，推测PmGSTTCHQD1 蛋白与油松和落叶松TCHQD 类谷胱甘肽S-转移酶蛋白功能类似（图7）。

图7 PmGSTTCHQD1蛋白系统进化分析Fig.7 Phylogenetic analysis of PmGSTTCHQD1 protein

3 结论与讨论

本文基于高通量转录组测序结果对PCZ090 标记所在基因序列进行追溯，获得了马尾松PCZ090基因的全长。并通过生物信息学软件分析SSRPCZ090 标记所在的mRNA 序列，包括同源序列比对、蛋白质一级、二级和三级结构及功能基因的预测。同源序列比对结果表明，该基因是马尾松谷胱甘肽S-转移酶家族蛋白编码基因成员。

植物谷胱甘肽S-转移酶（GSTs）是一类种类丰富且普遍存在的酶[12]。其参与多种细胞功能，如初生代谢、次生代谢、解毒作用、抗低温、抗氧化损伤及异源物质隔离等[13]。该类酶还能作为配体蛋白参与植物激素代谢途径[14]。在高等植物中，GST 被分为八种亚类：phi、tau、theta、zeta、lambda、EF1Bγ（elongation factor 1gamma）、DHAR（Dehydroascorbate reductase，谷胱甘肽依赖的抗坏血酸还原酶）和TCHQD（Tetrachlorohydroquinone dehalogenase，四氯代氢醌脱卤素酶）[14-15]。本研究挖掘得到马尾松PmGSTTCHQD1基因属于TCHQD 亚家族中的一员。该类谷胱甘肽S-转移酶的研究尚不全面，在模式植物拟南芥（Arabidopsis thaliana）中仅鉴定到1 个TCHQD 类GST，因其与原核酶序列同源，推测其可以催化氯代异源物质的代谢[16]。

GST 作为一种重要的解毒酶类，其催化反应主要分3 步进行：（1）催化特异性底物进行脱卤水反应；（2）已被活化的疏水底物与谷胱甘肽过氧化物酶（GSH）结合，增加底物的亲水性；（3）复合物被转运出细胞质。这一过程在植物中是由ATP 驱动泵驱动复合物进入液泡或者非共质体中[17]。在本研究团队通过候选基因组关联分析得到的与生长性状显著关联的其他位点中，PCZ099所在基因即为ATP酶基因。 由此推测PCZ090 所在基因PmGSTTCHQD1可能参与有毒物质的脱卤水反应，而PCZ099 所在的ATP 酶基因则参与复合物被转运出细胞质的过程，两者共同参与细胞解毒功能。细胞解毒能力很有可能影响马尾松的生长。

林木的生长性状是典型的数量性状，由多基因控制。本研究得到的PmGSTTCHQD1基因是通过关联分析研究发现的对生长性状表型变异解释率较高的基因之一，其编码蛋白所属的GST 酶类主要通过参与细胞解毒机制以提高抗性。其对生长性状的作用还需深入研究。可对关联分析群体中生长性状存在显著差异的个体进行PmGSTTCHQD1基因表达量分析、编码蛋白活性分析等。并开展PCZ090位点不同基因型间表型效应分析，研究其在不同林龄、不同地点的表型差异稳定，为利用该标记进行分子标记辅助选择提供参考。