基于RNA-Seq技术的湿加松特定时期产脂差异基因筛选及表达分析

周晨晨 李义良* 王 哲 赵奋成 钟岁英 林昌明谭志强 李福明 吴惠姗 郭文冰 廖仿炎 王建革

（1. 广东省森林培育与保护利用重点实验室，广州 510520；2. 广东省林业科学研究院，广州 510520；3. 广东省台山市红岭种子园，台山 529200；4. 江西农业大学园林与艺术学院，南昌 330045）

湿加松（Pinus elliottii × P. caribaea）是湿地松（P.elliottii）与加勒比松（P.caribaea）的杂交后代，源产于澳大利亚昆士兰，其在生长量、干形材质、抗病虫害等性状上均优于亲本，具有突出的杂种优势［1］。我国于上世纪90 年代开始湿加松的良种培育工作，培育的湿加松良种干型通直，15年时材积蓄积量可达225 m3·hm-2［2］。目前，湿加松已在我国广西、广东、福建、海南等地广泛种植，并逐渐北扩，成为我国南方培育材脂兼用林和大径材林的有潜力的重要树种。湿加松同其他松树一样可产生具有防御作用的萜类次生代谢物，人们俗称为松脂。松脂不仅是松科（Pinaceae）植物抵御昆虫、病菌、机械损伤的重要物理及化学防御手段［3～4］，又是我国林业经济的主要商品。我国松脂产量占世界松脂总产量70%左右，年出口创汇达1 亿美元。松脂加工后得到的松香、松节油可广泛用于香精香料、合成橡胶、液体生物燃料等方面［5～6］，其关联产业占国民经济的1/10。是我国重要的可再生自然资源。然而，近年来由于我国松林面积大幅减少、采脂技术不规范、乱采滥采现象严重等原因，导致我国松脂资源极度紧张。为改善我国松脂资源短缺现状，湿加松产脂优良品种的培育工作也逐渐开展，虽然已培育的产脂良种5年连续割脂，产脂量可达4.5 万kg·hm-2［2］，但目前对高产脂良种的研究工作仍处于起步阶段，对湿加松个体间出现产脂差异的分子机制不够清晰，导致湿加松产脂良种早期选择困难，短时期内很难补缺我国松脂资源短板。

RNA-Seq 技术已被广泛使用于林木研究中，可以快速全面反映植物组织在不同方式处理下或不同生长发育期内基因的表达情况，通过对获得的转录本信息进行处理可获得一些重要的功能基因，目前已广泛使用在植物次生代谢产物、抗病虫、非生物抗逆等研究中。刘彬等人［7］用接种松材线虫后表现出易感病及抗病的马尾松转录测序，通过数据分析得到一系列差异基因并筛选出参与抗虫的相关基因。Fox H 等人［8］基于RNA-Seq 技术，揭示了地中海松（P. halepensis）应对干旱的响应机制以及从干旱至恢复过程中的动态变化。本研究利用松脂产量受中度至高度的遗传控制这一特征［9～11］，选取生长量一致但产脂量显著差异的湿加松个体为材料，通过对特定时间下的样品进行转录组分析，筛选调控松脂产量差异的相关候选基因，利用qRT-PCR 技术对筛选的基因进行表达量验证分析，为后续分子标记、产脂功能基因克隆及验证奠定基础。

1 材料和方法

1.1 植物材料

试验材料来源于广东省台山市红岭种子园1996 年培育的湿加松测定林。参照《松脂采集技术规程（LY/T 1694-2007）》，对测定林内的155 株湿加松进行连续多次产脂量测定。利用SPSS 软件分析树高、胸径、材积与产脂量的相关性发现，四者之间存在较强的正相关。15 年时树高、胸径与产脂量的相关性分别为0.204、0.389、0.425（见表1）。因此为降低树体生长对松脂产量产生的影响，最终以材积为控制指标，选取材积无显著差异但产脂量连续多年稳定且有显著差异的高产脂（月产脂量＞1 350 g）与低产脂（月产脂量＜615 g）湿加松作为试验材料。

表1 树高、胸径、材积与产脂量的相关性Table 1 Correlation between tree height，DBH and yield of oleoresin

表2 材积、产脂量方差分析Table 2 Analysis of Variance of Volume and yield of oleoresin

为了减少不同基因型导致与产脂性状无关的差异，本试验共设置高产松脂组、低产松脂组共两组湿加松为材料，每组设置6个重复。高产松脂组各株编号分别为H1、H2、H3、H4、H5、H6；低产松脂组中各株编号为：L1、L2、L3、L4、L5、L6。两组湿加松材积无差异但产脂量差异显著（见表2）。选择树体产脂高峰的9 月中旬进行采样。由于树体在受到伤害后需经过一段时间相关基因才会启动表达，所以在割脂1 h 后再进行取样。在距离地面1 m 位置取剖面为10 cm×10 cm 的正方形，采用与割脂相同的下降式采脂法，采集树干次生木质部组织并立即放于液氮中保存，冻存后置于-80℃冰箱备用。

1.2 试验方法

1.2.1 RNA提取

参照TIANGEN 公司提供的多糖多酚植物总RNA（RNAprep Pure Plant Plus Kit）提取试剂盒说明书提取试验材料的总RNA，用0.5%琼脂糖凝胶电泳检测RNA 是否降解及污染，用超微量蛋白核酸分析仪—柏精（柏精-BioDrop μlite）检测RNA 的纯度（OD260/280比值在1.8～2.0）。

1.2.2 文库构建及高通量测序

将提取的RNA 纯化并检测质量，对符合测序要求的RNA 构建cDNA 文库。为获取更加完整的湿加松次生木质部信息，利用Illumina HiSeq 2500测序平台对12 个样品进行RNA 测序并按照高产脂、低产脂样品划分混合测序结果。

1.2.3 基因功能注释及表达分析

湿加松为无参基因组，为获得更全面的数据信息，将测序所得raw reads 数据质控后得到的clean reads 利用Trinity 软件进行拼接组装。组装后得到的unigenes 经NCBI Blast+与Nr，NT，Pfam，KOG/COG，Swiss-prot，KEGG，GO 共7 个数据库比对，得到相应的功能注释，同时将Trinity 拼接得到的转录本作为参考序列，采用RSEM 软件计算各个样品的基因表达水平。采用DESeq2 进行分析，筛选差异表达基因条件为阈值padj＜0.05 且|log2FoldChange|＞1。若log2FoldChange＞0 则该基因上调，反之该差异基因下调。将得到的差异基因进行GO功能分析及KEGG pathway 通路分析。

1.3 实时荧光定量PCR

用天根FastKing gDNA Dispelling RT Super-Mix FastKing 一步法除基因组cDNA 试剂盒对1 μg的总RNA 进行反转录，转录体系为20 μL。将反转录得到的产物进行10 倍稀释，使用Applied Biosystems 7500 仪器，荧光染料为ChamQTMUniversal SYBR qPCR Master Mix（Vazyme，China），进行实时荧光定量PCR 试验，各试剂加量及扩增程序参照荧光染料说明书，内参基因为TUB 基因，其余基因采用2-ΔΔCt算法计算在各个样品中的相对表达量。

2 结果与分析

2.1 RNA-Seq测序、拼接及功能注释

对高产脂与低产脂植株的总RNA 进行测序，测序所得序列经过过滤后得到594 108 382 条clean reads，总碱基为89.12 Gb。其中Q30 值达到了92%以上（见表3），clean reads 进行数据整合和拼接后，得到270 815 个transcripts 平均长度为934 bp，以及190 514个unigenes，平均长度为1 208 bp，有78 869个unigenes长度大于1 kB，200～500 bp的Genes 数量最多，N50 长度为1 886 bp，为后续分析提供了高质量的数据（见表4～5）。

2.2 基因功能注释

为更加直观地研究湿加松转录本信息，本研究将得到的190 514 条unigenes 在7 大数据库中进行注释，有144 023 条unigenes 至少被1 个数据库注释，占总量的75.59%，有10 265 条unigenes 在所有数据库中均有注释，占总量的5.38%（见表6）。其中有87 608 条unigenes 注释到NR 数据库，43.78%的unigenes 被注释到北美云杉（Picea sitchensis）中，其次是无油樟（Amborella trichopoda）、水芙蓉（Nelumbo nucifera）等，有36.62%的基因未得到明确的功能注释（见图1）。有30 323条unigenes注释到19 条KEGG 代谢通路中，主要富集的通路有“折叠、分类和降解（Folding，sorting and degradation）”，“碳水化合物代谢（Carbohydrate metabolism）”等。其中有770条unigenes富集到萜类物质合成“（Metabolism of terpenoids and poly ketides）”通路上（见图2）其中有304条unigenes直接参与到萜类物质合成过程中，涵盖了萜类物质合成的全过程。

表3 高产脂、低产脂转录组数据信息Table 3 High-yield and Low-yield oleoresin transcriptome data information

表4 拼接长度分布频数Table 4 Distribution frequency of splicing length

表5 湿加松All-unigenes长度分布Table 5 Distribution of All-unigenes length of P.elliottii×P.caribaea

表6 湿加松All-unigenes功能注释Table 6 All-unigenes function annotation of P. elliottii×P.caribaea

2.3 差异基因的筛选及分析

为探究湿加松高、低产脂基因差异表达的分子机理，采用DESeq2进行分析，共筛选出414个差异基因，其中183 个基因上调，231 个基因下调。将414 个差异基因进行GO 富集分析。261 条差异表达unigenes 在GO 数据库中具有功能注释，显著性富集的条目群主要是“氧化磷酸化过程”（oxidoreductase activity，acting on）及“双加氧酶活性”（dioxygenase activity）等。

在生物体内，不同基因相互协调行使其生物学功能。为了探究湿加松差异表达基因主要富集的代谢通道。使用KEGG 数据库，对差异基因进行显著性富集分析，得到97 个差异表达unigenes的功能定义，被映射到54 条KEGG Pathways 中。差异基因显著富集的20条代谢通路有内质网蛋白质加工（Protein processing in endoplasmic reticulum）、亚油酸（Linoleic acid metabolism）、苯丙氨酸、酪氨酸和色氨酸的生物合成（Phenylalanine，tyrosine and tryptophan biosynthesis）、剪接体（Spliceosome）、植物病原相互作用（Plant-pathogen interaction）等（见表7）。

2.4 转运与修饰过程对萜类物质的调节作用

根 据 差 异 基 因 在Nr，Nt，Pfam，KOG/COG，Swiss-prot，KEGG，GO 七大数据库中的功能释义发现，与萜类物质合成有关的基因多表现为转运与修饰过程。

表7 差异表达基因KEGG功能注释Table 7 Differential expression genes KEGG function annotation

差异基因中编码ABC 转运蛋白的基因共6 个（Cluster-37595.96613、Cluster-37595.128007、Cluster-37595.77181、 Cluster-37595.89400、 Cluster-37595.20113、Cluster-37595.127923），其中Cluster-37595.96613、 Cluster-37595.77181、 Cluster-37595.89400、Cluster-37595.20113 基 因 在 高 脂 树中的表达量显著高于低脂树中的表达量，Cluster-37595.128007、Cluster-37595.127923 在 低 产 脂 树中的表达量更高。编码细胞色素P450（Cytochrome P450）蛋白共有3 个基因（Cluster-36047.0、Cluster-31181.0、Cluster-14555.0），3 个基因在高产脂树中的表达量均高于低产脂树（见表8）。

2.5 抗逆蛋白在松脂产量中的调节作用

为更加全面探讨湿加松产脂差异的形成原因，进一步从植物防御系统方面对差异基因进行分析。差异基因注释表明，基因中编码植物抗病的TIR-NBS-LRR 蛋白有2 个（Cluster-2245.0、Cluster-37595.30441），均在高产脂植株中呈现高表达。编码双功能醇落叶松醇还原酶（PLR）蛋白1 个（Cluster-37595.73116），其在高产脂单株中表达量更高（见表9）。

表9 TIR-NBS-LRR、PLR的表达Table 9 The expression of TIR-NBS-LRR and PLR

2.6 qRT-PCR验证与产脂相关的候选基因

为了验证转录组数据结果可靠性以及与产脂有关的候选基因，我们选择与产脂有关的4个差异表达基因进行qRT-PCR 分析（见图3）。发现同一基因在高、低产脂不同单株中表现出差异性。CYP450 基因在H1、H4、H5 的表达量高于H6、H2、H3，在低产脂植株中趋于稳定；PLR 基因在低脂树中表达量较为均一，但高脂树中H1、H2、H6 明显高于H3、H4、H5。ABC transporter 在高脂树中H1、H2、H5 的表达量高于H3、H4、H6，其中H5 与H6的表达量差异显著；异戊烯基转移酶在L5 中的表达量最高，其次是L1、L3 的表达量最低，但也显著高于高脂树中的表达。GGPS 基因在高产脂、低产脂树中的表达趋于一致性，高产脂植株中H1、H2、H3 的表达显著高于H4、H5、H6，低脂树中L5、L6的表达量最高，L2 中的表达量最低。总体来看CYP450、ABC transporter、PLR 基因在高脂树中的表达量显著高于低脂树中的表达量，异戊烯基转移酶在低脂树中的表达量更高，GGPS 基因在高脂树与低脂树的表达量无差异，此结果与转录组测序结果一致。

3 讨论

松脂中95%以上是由萜类物质组成，因此人们认为萜类物质合成途径是松脂合成的主要通路。异戊烯基转移酶是萜类物质合成中的一类关键酶。根据产物不同，异戊烯基转移酶可分为香叶基二磷酸合酶（GGPS）、法尼基二磷酸合酶（FPS）、香叶基香叶基二磷酸合酶（GGPPS）［12］。思茅松（Pinus kesiya var.langbia nensis）中GGPS 基因表达量与产脂量呈显著正相关［13］。刘美佳等基于珠子参（Panax japonicus）转录组结果，克隆出一条FPS 的cDNA 序列，并构建珠子参过表达载体导入珠子参细胞中，发现转基因植物的甾醇及三萜含量都有所升高［14］。杨锐等在丹参（Salvia miltiorrhiza Bge.）中过表达GGPPS 基因，发现丹参中丹参酮的含量显著高于对照组［15］。本实验中异戊烯基转移酶在低脂树中的表达量更高，此研究结果与Bai Q S 等人［16］的研究结果相似。当树体遭遇伤害后，高脂树本身储存的松脂大量流出可隔绝伤口与外界的接触，而低脂树树体内的松脂流动性较弱，树体储存的松脂较少，为保护自身免受外界的进一步伤害而增加萜类物质含量以封闭伤口形成物理屏障，因此异戊烯基转移酶在低脂树中的表达量更高。

萜类物质是目前植物体内种类最多的一类次级代谢产物，其种类多样离不开各种修饰酶的作用，其中细胞色素P450 是修饰萜类物质种类多样的一类关键酶，97%的萜类物质的多样性是通过细胞色素P450（CYPs）氧化催化完成［17］。CYP450是一类家族基因，针叶树种中，CYP450 多参与树脂酸的修饰与合成［18～20］，本次测序差异基因中细胞色素P450 上调表达，表明高脂树的CYP450 表达量更高。松脂由树脂道上皮细胞合成并在树脂道中储存。松脂从上皮细胞转移到树脂道的过程中需ABC 转运蛋白的参与。Westbrook 等人［21］发现ABC 转运蛋白与松脂的跨位点流动有关，本研究中转录组与qRT-PCR 数据表明，参与转运的ABC 转运蛋白上调表达，与Liu［22］、陈晓明［23］等人的研究结果一致。这些研究结果表明湿加松松脂产量差异主要表现在萜类物质的修饰与转运过程。

对植物体而言，松脂是其重要的物理及化学防御手段。TIR-NBS-LRR 蛋白属于R 基因的一种，主要参与植物体对病虫的防御过程。李海霞等人［24］从毛白杨（Populus tomentosa）中分离出的一条TIR-NBS-LRR 抗病基因（PtDRG01）并将其导入到拟南芥（Arabidopsis thaliana）中，发现转基因烟草株系的发病程度明显低于非转基因烟草，同时转基因株系中PtDRG01的基因表达水平显著高于非转基因株系。张颖等［25］从野生刺葡萄中克隆出10 条NBS-LRR 类抗病基因片段，同时证实了NB6、NB7的表达受白腐病的调控。王猛［26］也提出TIR-NBS-LRR 的表达在马尾松抗松材线虫中发挥作用。这些结果说明TIR-NBS-LRR 蛋白参与了植物体抗病虫的防御过程。木脂素是植物体重要的防御物质，其合成受双功能松脂醇—落叶松醇还原酶（PLR）和PR基因的调控。本试验中双功能松脂醇—落叶松醇还原酶在高脂树中表达量显著高于低脂树中的含量。野生金荞麦（Fagopyrum cymosum.）的研究中发现，PLR 可催化松脂醇—落叶松脂醇还原开环异构落叶松脂并参与植物体的防御［27］。本研究中编码TIR-NBS-LRR 蛋白的基因及PLR 基因在高产脂树中的表达量显著高于低产脂植株，说明高产脂植株对抵御病虫伤害时表现出更强的抗性。研究表明辐照西洋参中随着ATP 浓度会影响其次生产物皂苷的合成［28］。马尾松中高产树脂的ATP 含量要比低产树脂高［23］。本研究发现GO 项主要富集在氧化磷酸过程及双加氧酶过程中，氧化磷酸化可为植物体生命活动提供必要的ATP，ATP 是细胞进行代谢活动最直接的能量来源，对代谢产物的生成有调控作用，高脂树的代谢活动能力更强。

湿加松松脂合成途径中，部分相关调控基因在高脂树与低脂树中的表达量不同。调控萜类物质转移的细胞色素P450 等基因、ABC 转运体家族基因、参与防御系统的PLR 基因、TIR-NBS-LRR 基因上调表达直接或间接促进了松脂总量的增加，异戊烯基转移酶在低产脂植株中表达量更高。本试验可以作为从转录层面探讨湿加松产脂差异的形成机制，为后续相关分子标记和功能基因克隆奠定重要基础。