基于马尾松转录组测序的产脂相关SNPs初步分析

吴志铭，冯源恒，杨章旗

（1.广西师范大学生命科学学院，广西桂林 541006；2.广西壮族自治区林业科学研究院广西马尾松工程技术研究中心 广西优良用材林资源培育重点实验室，广西南宁 530002）

马尾松（Pinus massoniana）是松科（Pinaceae）松属常绿针叶树种，分布于秦岭淮河以南和云贵高原以东等17 个省、自治区和直辖市[1]。马尾松材脂兼用，木材可用于建筑、造纸和木纤维工业用材等领域；松脂加工产品是重要的工业原料，可应用于医药、油漆和粘合剂等。中国70%以上的松脂产量来自马尾松[2-4]。据第六次至第八次全国森林资源连续清查数据统计，马尾松的森林面积不断减少，加上不断采割导致采脂树木产脂能力降低及树木死亡，致使马尾松原脂产量出现萎缩。加快高产脂马尾松选育已成为其遗传改良的重要内容。传统选育通过连续几年对处盛产期的马尾松进行采脂测定，对其产脂能力做出评价，选出高产脂单株。其成本较高、耗时长、范围窄且进展慢，开展分子辅助育种，能有效地缩短马尾松育种年限[5]。

单核苷酸多态性（SNPs）是指单核苷酸在基因组水平上的突变引起的DNA 序列间的多态性，其突变包括单碱基的转换、颠换和插入缺失[6]。SNP 因其在基因组中具有分布广、多样性高和易于分型等特点, 成为基因分型最理想的分子标记[7]。表达序列标签（EST）是源于转录表达的特异功能基因的cDNA 片段。利用获得的EST 序列进行SNP 标记开发，对未进行全基因组测序的动植物个体具有重要意义[8]。

本研究对马尾松二代测序转录组数据中的SNP位点进行挖掘，并分析其功能和通路，为马尾松产脂性状关联分析的开展提供可用的SNP标记。

1 材料与方法

1．1 材料

基于连续3年在南宁市林业科学研究所马尾松高产脂种质资源库的采脂试验结果，采集高产脂无性系桂GZ080B（15．96 g/10 cm）和普通产脂无性系桂GZ078B（8．18 g/10 cm）的3 个组织（顶芽、针叶和韧皮部）材料，迅速放入液氮中冷冻，送至北京诺禾致源科技股份有限公司进行高通量测序（Illumina HiSeqTM 2000/MiSeqTM），获得转录组数据。

1．2 RNA提取与检测

通过试剂盒法提取GZ080B 和GZ078B 的顶芽、针叶和韧皮部的RNA，分别使用1%琼脂糖凝胶电泳监测RNA 降解和污染，NanoPhotometer 分光光度计（IMPLEN，CA，USA）检测RNA 的纯度，Qubit 2．0 Flurometer（Life Technologies，CA，USA）中的Qubit RNA 分析试剂盒测量RNA 浓度，Agilent Bioanalyzer 2100 系统（Agilent Technologies，CA，USA）的RNA Nano 6000分析试剂盒评估RNA完整性（图1）。

图1 RNA凝胶电泳图Fig.1 RNA gel electrophoresis

1．3 SNP位点挖掘

从高产脂无性系和普通产脂无性系的顶芽、针叶和韧皮部3 个比较组成的维恩图中得到2 329 个差异表达基因，该图可直观展现各种组合间的差异表达基因数量。为筛选出含SNP位点的Unigene，从转录组数据中找到各比较组中的差异表达基因统计表（表格有7大数据库注释），根据NR注释结果挑选已知功能的Gene ID，通过Novofinder（北京诺禾致源科技股份有限公司测序结果自带）输入Gene ID搜索含SNP 位点的Gene ID 的Unigene，并对其突变类型数量及所在密码子位置进行统计，最后根据含SNP 位点的Gene ID 搜索其Unigene 的GO 功能注释和KEGG 通路注释结果。从中选出3 个比较组中共有的Unigene、只存在于针叶和韧皮部的Unigene 和韧皮部独有的Unigene进行分析。

1．4 SNP位点Unigene的GO功能分析

根据含SNP 位点的Unigene 的GO 注释结果，经过数据处理，通过OmicShare 在线软件的动态GO 富集分析（https://www．omicshare．com/tools/home/report/goenrich．html），把含SNP 位点的Unigene 的GO 注释结果进行功能分类，最后对GO功能进行统计分析。

1．5 SNP位点Unigene的KEGG代谢通路分析

根据含SNP 位点的Unigene 在KEGG 数据库的注释结果，剔除没有K 编号的Unigene，经过数据处理，通过OmicShare 在线软件的动态KEGG 富集分析（https://www．omicshare．com/tools/home/report/koenrich．html），对含编号的Unigene进行统计分析。

2 结果与分析

2．1 马尾松SNP位点的挖掘

将高产脂无性系与普通产脂无性系的韧皮部、针叶与顶芽的转录组数据进行两两对比（图2）。在发现的2 329 个差异表达基因中进行两次筛选，第1次根据NR 数据库的注释结果筛选出366 条Unigene，第2 次筛选出含SNP 位点的Unigene 125 条，共656 个SNP 位点。对656 个SNP 位点的突变类型进行统计，发现转换类型有4 种，颠换类型有8 种。发生转换突变频率较高，其中T/C 和C/T 转换占总SNP位点的33．54%，A/G 和G/A 占28．05%；发生颠换类型的各种突变频率较低，分别为11．59%（C、G）、9．60%（G、T）、8．69%（A、T）和8．53%（A、C）。对SNP位点所在密码子位置进行统计时，发现只有47．99%的SNP 位点在密码子上，在第一位置、第二位置和第三位置发生的突变比例分别为2∶1∶2。为进一步了解SNP 位点的信息，分析每个个体该位点的基因型和突变后的基因型，根据支持该位点的reads个数和GATK3 软件得到的该位点的基因型，若/两边碱基相同，则为纯合位点，若不同，则为杂合位点；基因型在不同部位中相同但在产脂能力不同的马尾松中不同。根据统计结果发现，纯合突变32 个，杂合突变121个。

2．2 马尾松SNP位点所在Unigene的GO分类

图2 差异表达基因Venn图Fig.2 Venn map of differentially expressed genes

为了解筛选出的含有SNP 的Unigene 的功能，对GO 注释结果进行进一步分类。这些Unigene 被注释到3大类41个功能区（图3）。其中93条Unigene 参与生物过程（Biological process）的18 个功能区，49条Unigene 参与细胞成分（Cellular component）的14 个功能区，92 条参与分子功能（Molecular function）的9 个功能区。生物过程中参与代谢过程（metabolic process，73 条）、细胞过程（cellular process，68 条）和单一生物过程（single-organism process,64条）的基因最多，参与生物过程的正、负调节(positive regulation and negative regulation of biological process)和免疫系统过程(immune system process)的基因最少，均只有1 条；细胞成分中参与细胞（cell，25 条）、细胞组分（cell part，25 条）和膜（membrane，24 条）的基因最多，其次是大分子复合物（macromolecular complex，18条）；分子功能中参与催化活性（catalytic activity，71 条）和结合活性（binding，62 条）的基因最多，参与抗氧化活性（antioxidant activity，1 条）和分子功能调节器（molecular function regulator，1 条）的基因最少。含有SNP 的Unigene 主要与马尾松的代谢过程相关。

图3 马尾松SNP位点所在Unigene的GO分类Fig.3 GO classification of Unigene SNP loci of P.massoniana

2．3 马尾松SNP 位点Unigene 的KEGG 代谢通路分析

对转录组数据中125 条含SNP 的Unigene 的KEGG注释结果进行处理，发现57条含K编号，已知KEGG 功能的基因有29 条被注释到5 大类13 个通路中（图4）。其中新陈代谢（Metabolism）有21 条Unigene，环境信息处理（Environmental information processing）有4 条，组织系统（Organismal systems）有2 条，遗传信息处理（Genetic information processing）和细胞过程（Cellular processes）各1条。新陈代谢的Unigene 可分为9 个亚类，氨基酸代谢类（Amino acid metabolism）最多，有8 条；其次是碳水化合物代谢类（Carbohydrate metabolism）和其他次生代谢产物合成（Biosynthesis of other secondary metabolites），分别有7 条和6 条；涉及萜类和聚酮化合物代谢（Metabolism of terpenoids and polyketides）的有4 条。结果表明，KEGG 代谢通路分析与GO 分类得出的结果均与代谢相关。

图4 马尾松SNP位点所在Unigene的KEGG的代谢通路Fig.4 Metabolic pathway of KEGG in Unigene at P.massoniana SNP loci

3 讨 论

本研究根据马尾松高产脂无性系桂GZ080B和普通产脂无性系桂GZ078B的转录组测序结果，在差异表达分析结果中根据NR数据库注释和含有SNP位点两大特点，从2 329个Unigene中筛选出374条Unigene，共2 192个SNP位点，根据试验要求从中选出3个比较组中共有的基因、只存在于针叶和韧皮部的基因和韧皮部独有的基因。在对突变的类型进行统计时，发现C转换为T的频率最高（33．54%），原因是CG中的C常为甲基化状态，自发脱氨后成为胸腺嘧啶。

综合GO分类和KEGG代谢通路分析，从马尾松转录组中筛选出的含SNP 的基因主要是与生物体代谢和分子功能相关，与萜类化合物和聚酮类化合物合成相关的基因较少。SNP标记位点源自编码序列，通过EST 数据库可以直接开发出与功能基因相关的SNP 标记，为进一步的功能基因研究提供依据[9]。本研究可为马尾松产脂相关基因与标记开发等研究提供参考。