芳樟叶片和枝的转录组分析及生物炭复合肥对芳樟醇合成调控基因表达特性的影响

向 双， 陈德强， 孙维红， 肖 琳， 邹双全，①

(福建农林大学： a. 林学院， b. 自然生物资源保育利用福建省高校工程研究中心， 福建 福州 350002)

芳樟(Cinnamomumcamphoravar.linalooliferaFujita)是典型的芳香植物，其香气的关键成分为芳樟醇。芳樟醇是全球最常用和消耗较多的香料之一[1]，具有多种生物活性[2-4]，广泛应用于食品、医药和化妆品等行业[5]，市场需求量巨大。芳樟醇可分为天然芳樟醇和合成芳樟醇2类，其中，合成芳樟醇为右旋体、具有消光性，其品质远差于左旋体的天然芳樟醇[6]，但由于天然芳樟醇植物源提取物产量较低，加之芳樟自然资源有限，导致天然芳樟醇产量供应不足，使天然芳樟醇的开发利用受到限制。通过基因工程技术对芳樟进行改良，明确其芳樟醇合成调控基因，可为高芳樟醇芳樟品种资源的培育提供新思路和新途径。

近年来，植物的转录组测序在挖掘功能基因、探明活性成分代谢通路和明确次生代谢物合成关键酶基因等方面得到广泛应用[7-8]，通过转录组学研究已经解析了黄酮类[9]、甾体皂苷类[10]和三萜类[11]等成分的生物合成途径。本项目组前期研究结果表明：芳樟精油中的芳樟醇含量受多种因子的影响，不仅因品种、品系(无性系)不同而异[12-14]，也因采收时间不同而变化[15]；并且，对芳樟植株采取不同的栽培措施(例如施肥时间和肥料种类)后，其精油含量和芳樟醇含量也有明显的差异[16-17]。而目前这些影响因子的作用机制尚不明确，特别是对芳樟醇合成过程中分子调控机制缺乏全面而深入的研究，不利于芳樟资源的进一步开发利用。

为了明确芳樟醇合成途径中的关键调控基因，探讨芳樟醇合成相关基因对外在因子的响应机制，本文利用高通量测序平台构建芳樟叶片和枝的转录组数据库，对unigenes进行功能注释，挖掘其芳樟醇合成过程中的关键调控基因；并以“生物炭复合肥对芳樟醇含量有影响效应”的研究结果(另文发表)为基础，比较施肥前后芳樟幼树叶片和枝中芳樟醇含量和相关调控基因表达特性的变化，以期为芳樟中芳樟醇关键基因克隆和功能验证及其品种的遗传改良等奠定研究基础。

1 材料和方法

1.1 材料

供试材料为芳樟优质品系‘MD1’，由厦门牡丹香化实业有限公司培育，种植于福建省泉州市安溪县半林国有林场(东经117°57′、北纬24°55′)。样地面积10 m×10 m，按株距1 m、行距1 m种植60株2年生幼树；幼树平均株高77.11 cm，平均地径15.89 mm，平均冠幅80.95 cm。样地为近年采伐迹地，土质疏松，地力较肥沃，坡度较缓，适合芳樟种植。

根据项目组前期的研究结果确定生物炭复合肥配方，包括3.40 g生物炭(在500 ℃高温厌氧条件下热解玉米秸秆2 h制备)、115.00 g尿素(总N含量大于等于46.0%)、15.00 g过磷酸钙(P2O5含量大于等于16.0%)和15.00 g氯化钾(K2O含量大于等于60%)，其中，生物炭购自四川省久晟农业有限责任公司，尿素购自安阳中盈化肥有限公司，过磷酸钙购自瓮福有限责任公司，氯化钾购自中国农业生产资料集团公司。

1.2 方法

1.2.1 施肥处理及样品采集 于2019年6月9日进行施肥处理，按株施肥，每株施肥量148.40 g，均一次性施入。分别于6月9日施肥前以及7月9日、9月9日和11月9日(即施肥后第1、第3和第5个月)采集无病虫害的当年生叶片和枝，在样株东、南、西、北4个方位分别采集叶片和枝，各自混合装袋，每袋100 g，每份样品3袋，记为3次重复，用于精油提取和芳樟醇含量测定。转录组分析样品用5 mL冻存管收集，每份样品3管，记为3次重复，置于-80 ℃下保存、备用。

1.2.2 精油提取和芳樟醇相对含量测定 采用水蒸气蒸馏法[18]提取芳樟叶片和枝的精油，并采用HPLC法测定精油中芳樟醇的相对含量。

精密称取43.5 mg芳樟醇标准品(CAS号78-70-6，纯度在97%及以上，购自美国Acros Organics公司)，加入甲醇(色谱纯，德国Merck公司)溶解并定容至10 mL，配制成质量浓度4.35 mg·mL-1的芳樟醇标准品溶液；分别吸取2、4、6、8、10和12 μL标准品溶液，按下述色谱条件进行检测。色谱条件：Agilent 1100高效液相色谱仪(美国Agilent公司)，Dikma Diamonsil C18色谱柱(250 mm×4.6 mm，5 μm)；流动相为V(甲醇)∶V(水)=70∶30混合液，等度洗脱，流速1 mL·min-1；检测波长210 nm，柱温30 ℃，进样量10 μL。记录色谱图，并以标准品峰面积为纵坐标(y)、标准品质量为横坐标(x)绘制标准曲线，回归方程为y=(6×108)x+61 162(R2=0.999 9，线性范围为4.35～52.20 μg)。

称取100 g叶片或枝，于50 ℃烘干至含水率5%～13%，粉碎后置于500 mL挥发油提取器中，加入200 mL超纯水加热提取60～90 min，蒸馏后取上层精油，经无水硫酸钠干燥后称量，按照公式“芳樟精油提取率=芳樟精油质量/样品鲜质量×100%”计算芳樟精油提取率。精密称取各精油样品45 mg，用甲醇定容至10 mL，过滤，滤液即为供试样品溶液；按照上述色谱条件进样测定，每份样品测定3次，记录各样品中芳樟醇的峰面积。

根据芳樟醇标准曲线和样品中芳樟醇的峰面积，按照公式“芳樟醇相对含量=(样品溶液中芳樟醇的峰面积/芳樟醇标准品的峰面积)×(芳樟醇标准品的浓度/样品溶液的浓度)×100%”计算样品中芳樟醇的相对含量。

1.2.3 转录组测序与组装 根据精油提取率和芳樟醇相对含量测定结果，以施肥前及施肥后第3个月的叶片和枝为材料，用RNA试剂盒(美国Omega Bio-Tek公司)提取总RNA；用NanoDrop微量分光光度计(美国Thermo Fisher公司)检测其纯度，并用Agilent 2100生物分析仪(美国Agilent Technology公司)检测其浓度和完整性。

随后，用Oligo(dT)磁珠从总RNA中分离和纯化mRNA，并在富集后进行片段化处理，用随机引物和逆转录酶以mRNA为模版合成cDNA，将cDNA片段末端修复，构建芳樟转录组的cDNA文库。用DNBSEQ高通量测序平台(武汉华大基因科技有限公司)对检测合格的cDNA文库进行测序；过滤低质量、接头污染及未知碱基N含量过高的reads，得到高质量reads；用Trinity软件[19]对高质量reads进行从头组装，并用TGICL软件[20]将组装后的转录本进行聚类去冗余，最终得到unigenes。

1.2.4 Unigenes功能注释 用NCBI BLASTx软件[21]将unigenes与NR(NCBI非冗余蛋白质序列数据库)、NT(核酸序列数据库)、Swiss-Prot(经实验验证的蛋白质序列数据库)、Pfam(蛋白质家族数据库)、GO(基因本体数据库)、KEGG(京都基因与基因组百科全书)和KOG(真核生物直系同源基因数据库)7个数据库分别进行比对，并对unigenes进行功能注释、预测和分类。

1.2.5 基因表达特性分析 用Bowtie2软件[22]分别对芳樟叶片和枝转录组进行序列比对，并用RSEM软件计算基因的表达水平，得到每个unigene的表达水平(FPKM)。用DESeq2软件[23]，根据筛选条件“阈值为|log2FC|>1”和“错误率(FDR)小于0.01”，对与芳樟醇合成有关的调控基因进行筛选，并对各基因的表达水平进行分析。

1.3 数据处理和分析

用EXCEL 2017软件对芳樟转录组unigenes数进行统计并制图，用IBM SPSS Statistics 25软件进行方差分析，采用Dunnett检验进行多重比较。

2 结果和分析

2.1 转录组测序和组装及unigenes功能注释

2.1.1 转录组测序和组装结果 通过DNBSEQ测序平台，从芳樟叶片和枝的转录组中共获得72 560个unigenes，平均长度为1 330 bp，unigenes总长度达到96 543 844 bp，N50值为2 054 bp，GC含量为43.19%。

由unigenes长度分布数据(表1)可见：在72 560个unigenes中，长度小于300 bp的unigenes有12 086个，仅占unigenes总数的16.7%；长度3 000 bp以上的unigenes有6 441个，长度大于1 000 bp的unigenes数占unigenes总数的49%以上，表明芳樟叶片和枝的转录组序列完整性较高，可用于unigenes功能注释。

表1 芳樟叶片和枝转录组unigenes的长度分布1)

2.1.2 Unigenes功能注释结果 对获得的芳樟叶片和枝转录组unigenes进行数据库比对注释，结果显示：在72 560个unigenes中，共有52 613个unigenes被注释到7个数据库中，占unigenes总数的72.51%。其中，在NR数据库中注释到50 149个unigenes，占unigenes总数的69.11%，是注释unigenes数最多的数据库；分别有36 090、39 232、40 698、40 096、39 995、39 186 个unigenes注释到NT、Swiss-Prot、KEGG和KOG、Pfam和GO 6个数据库中，注释的unigenes数分别占unigenes总数的49.74%、54.07%、56.09%、55.26%、55.12%和54.00%。同时注释到7个数据库的unigenes数有21 055个，占unigenes总数的29.02%。

2.1.2.1 序列相似性分析 将芳樟叶片和枝转录组unigenes与NR数据库中一些植物转录组unigenes进行序列相似性比较。结果显示：芳樟与荷花(NelumbonuciferaGaertn.)的同源序列匹配度最高，相似序列比例达到27.09%；其次为博落回〔Macleayacordata(Willd.) R. Br.〕，相似序列比例为14.91%；与葡萄(VitisviniferaLinn.)、油棕(ElaeisguineensisJacq.)和海枣(PhoenixdactyliferaLinn.)的同源序列匹配度也较高，相似序列比例分别为5.85%、5.43%和3.95%；与其他植物的同源序列匹配度均较低，这部分序列比例为42.76%。

2.1.2.2 GO注释分析 在GO数据库对芳樟叶片和枝转录组unigenes功能进行注释，结果见图1。

由图1可见：共有39 186个unigenes注释到生物过程(biological process)、细胞成分(cellular component)和分子功能(molecular function)3大类24亚类中，其中注释到分子功能的unigenes最多。在注释到生物过程类别的21 993个unigenes中，注释到细胞过程(cellular process)的unigenes最多，有9 520个；在注释到细胞成分类别的28 829个unigenes中，注释到膜部分(membrane part)的unigenes最多，有12 323个；而在注释到分子功能类别的44 405个unigenes中，注释到结合(binding)和催化活性(catalytic activity)的unigenes较多，分别有20 456和18 905个。

： 生物过程Biological process； ： 细胞成分Cellular component； ： 分子功能Molecular function.图1 芳樟叶片和枝转录组unigenes的GO功能分类结果Fig. 1 Result of GO function classification of unigenes in transcriptomes from leaf and branch ofCinnamomum camphora var. linaloolifera Fujita

2.1.2.3 KEGG代谢途径分析 通过KEGG数据库对芳樟叶片和枝转录组unigenes功能进行注释，结果见图2。

由图2可见：共有40 698个unigenes注释到细胞过程(cellular process)、环境信息处理(environment information processing)、遗传信息加工(genetic information processing)、代谢(metabolism)和生物系统(organismal system)相关的代谢通路中，共5大类19亚类，且主要集中在代谢和遗传信息处理2大类中。

： 细胞过程Cellular process； ： 环境信息处理Environment information processing； ： 遗传信息加工Genetic information processing； ： 代谢Metabolism； ： 生物系统Organismal system.图2 芳樟叶片和枝转录组unigenes的KEGG功能分类结果Fig. 2 Result of KEGG function classification of unigenes in transcriptomes from leaf and branch ofCinnamomum camphora var. linaloolifera Fujita

注释到细胞过程类别的所有unigenes均注释到运输与分解代谢(transport and catabolism)亚类中，共有1 808个unigenes。注释到环境信息处理类别的unigenes共有2 923个，其中注释到信号转导(signal transduction)亚类的unigenes较多，有2 470个。注释到遗传信息加工类别的unigenes共有8 829个，其中注释到翻译(translation)亚类的unigenes最多，有3 553个。注释到代谢类别的unigenes共有23 363个，其中注释到全局和概述地图(global and overview map)亚类的unigenes最多，有8 972个。注释到生物系统类别的所有unigenes均注释到环境适应性(environmental adaptation)亚类中，共有2 057个。

芳樟醇主要通过MEP(甲基赤藓醇磷酸)途径合成前体物质，因此，根据KEGG数据库的注释结果，筛选出与芳樟醇合成相关的基因及其对应的酶，并推测芳樟芳樟醇的合成途径。结果显示：DXS(1-脱氧-D-木酮糖-5-磷酸合成酶)催化丙酮酸和3-磷酸甘油醛生成DXP(1-脱氧-D-木酮糖-5-磷酸)，是芳樟醇合成的第1个限速酶；随后DXR(1-脱氧-D-木酮糖-5-磷酸还原异构酶)催化DXP生成2-C-甲基-D-赤藓醇-4-磷酸，并在CMS(2-C-甲基-D-赤藓糖醇-4-磷酸胞苷酰转移酶)、CMK〔4-(5′-焦磷酸胞苷)-2-C-甲基-D-赤藓糖醇激酶〕、MDS(2-C-甲基-D-赤藓糖醇-2，4-环二磷酸合成酶)、HDS(4-羟基-3-甲基-2-丁烯基二磷酸合成酶)、HDR(4-羟基-3-甲基-2-丁烯基二磷酸还原酶)和GPPS(牻牛儿基焦磷酸合酶)等一系列酶的作用下合成GPP(牻牛儿基二磷酸)，最后在LIS(芳樟醇合酶)的作用下生成芳樟醇。

2.1.2.4 KOG功能预测结果 通过KOG数据库对芳樟叶片和枝转录组unigenes的功能进行预测，结果见图3。

由图3可见：共有40 096个unigenes注释到24个KOG类别中。其中，注释到一般功能预测(general function prediction only)的unigenes有9 104个，数量最多，说明功能预测是芳樟叶片和枝转录组中最主要的功能分类。另外，注释到信号转导机制(signal transduction mechanism)功能的unigenes有5 352个，注释到翻译后修饰、蛋白质转换和伴侣(posttranslational modification，protein turnover and chaperone)功能的unigenes有3 182个，注释到转录(transcription)功能的unigenes有2 544个，注释到RNA加工和修饰(RNA processing and modification)功能的unigenes有1 522个，注释到次生代谢产物生物合成、运输和分解代谢(secondary metabolites biosynthesis， transport and catabolism)功能的unigenes有1 043个。上述unigenes也是芳樟叶片和枝转录组较为重要的基因。

图3 芳樟叶片和枝转录组unigenes的KOG功能分类结果Fig. 3 Result of KOG function classification of unigenes in transcriptomes from leaf and branch of Cinnamomum camphora var. linaloolifera Fujita

2.2 施肥后芳樟醇相对含量及其合成调控基因表达水平的变化

2.2.1 芳樟醇相对含量的变化 施肥前后芳樟叶片和枝中芳樟醇相对含量的变化见表2。比较结果显示：不论是施肥前还是施肥后，叶片中芳樟醇相对含量均高于枝，但差异幅度不大。但施肥后不同时间叶片和枝中芳樟醇相对含量均高于施肥前，差异均达显著(P<0.05)水平；并在施肥后第3个月叶片和枝中芳樟醇相对含量均达到最高，分别为89.89%和84.83%，但在施肥后第5个月叶片和枝中芳樟醇相对含量则显著下降。

表2 施肥前后芳樟叶片和枝芳樟醇相对含量变化

2.2.2 芳樟醇合成调控基因及其表达水平变化 结合芳樟叶片和枝中芳樟醇的代谢途径，进一步筛选关键调控基因，并对施肥前后各基因的表达水平(FPKM)进行分析，结果见表3。

由表3可见：从施肥前芳樟叶片和枝的转录组数据中共鉴定出23个与芳樟醇合成相关的基因，其中，15个基因在枝中的表达水平高于叶片，8个基因在叶片中的表达水平高于枝，说明在芳樟不同部位这些基因的表达具有组织特异性。

由表3还可见：与施肥前相比，施肥后第3个月叶片中有18个基因表达水平上调，枝中有15个基因的表达水平上调，叶片和枝中均上调表达的基因有13个，包括3个DXS基因(DXS1、DXS3和DXS6)、1个DXR基因、2个HDS基因(HDS1和HDS2)、1个HDR基因(HDR2)、2个IDI基因(IDI1和IDI2)、1个GPPS基因(GPPS1)和3个LIS基因(LIS1、LIS2和LIS4)；另外，虽然施肥后叶片中CMS和GPPS2基因以及枝中CMK、MDS和IDI3基因的表达水平下调，但这些基因的总体表达水平较高，表明这些基因都可能参与芳樟醇合成的调控。

表3 施肥前和施肥后第3个月芳樟叶片和枝中芳樟醇合成调控基因表达水平(FPKM)的变化

在施肥前后，叶片中MDS和IDI2基因以及枝中GPPS1基因的表达水平差异显著(P<0.05)，表明施用生物炭复合肥可能对这些基因的表达有一定的影响。且这些基因表达水平均上调，表明施用生物炭复合肥可能对基因表达有促进作用。而其他基因尽管在施肥前后表达水平差异不显著，但也出现不同程度的上调或下调，表明施用生物炭复合肥对这些基因的影响较小。

3 讨论和结论

上述高通量测序结果表明：从芳樟叶片和枝的转录组中共获得72 560个unigenes，平均长度1 330 bp，N50值为2 054 bp，转录组数据较为完整，基因信息较为丰富；而江香梅等[24]对5种化学型樟树叶片转录组进行测序，共获得55 955个unigenes，平均长度584 bp， N50值为1 023 bp，造成这一差异的原因可能与转录组的来源器官、株龄以及栽培措施和生境条件等因子有关。本研究结果表明：通过KEGG数据库比对，芳樟叶片和枝的转录组中有40 698个unigenes被注释到5大类19亚类代谢途径中，涵盖了多糖、甾体皂苷、萜类和苯丙素类等生物合成相关基因。其中，注释到萜类和聚酮化合物代谢的有603个unigenes，明显少于花生(ArachishypogaeaLinn.)转录组(1 033个)[25]，可能与花生的萜类和聚酮化合物比芳樟更丰富有关。

萜类化合物合成途径分为MEP和MVA(甲羟戊酸)两类途经[26]，其中MEP途径主要合成单萜和二萜及类胡萝卜素和叶绿素等多萜类化合物，MVA途径主要合成倍半萜和甾醇等三萜类化合物。芳樟醇是单萜醇，主要通过MEP途径生成前体物质GPP，然后在LIS的作用下生成芳樟醇；芳樟醇是桂花〔Osmanthusfragrans(Thunb.) Loureiro〕[27]、夜来香〔Telosmacordata(Burm. F.) Merr.〕[28]和小苍兰(FreesiarefractaKlatt)[29]等芳香植物的特征香气成分之一，也是芳樟叶片和枝精油的重要成分之一。本研究鉴定出23个与芳樟醇合成相关的基因，其中，15个基因在枝中表达水平高于叶片，8个基因在叶片中表达水平高于枝，揭示出这些基因在芳樟中的表达具有组织特异性。枝中有15个基因的表达水平高于叶片，但枝中的芳樟醇含量反而低于叶片，推测这15个基因对芳樟醇合成调控的影响作用小于在叶片中表达较高的8个基因。DXS是芳樟萜类合成中第1个限速酶[30]，在拟南芥〔Arabidopsisthaliana(Linn.) Heynh.〕[31]、银杏(GinkgobilobaLinn.)[32]和番茄(LycopersiconesculentumMiller)[33]等植物中，DXS基因过量表达可使各种单萜类化合物产量显著提高。而在芳樟转录组数据中，2个DXS基因表达水平较低，并且施肥后芳樟醇相对含量提高，3个DXS基因表达水平也上调，因此，可以通过调控DXS基因的表达水平改变芳樟醇含量。在芳樟叶片和枝的转录组中，尽管施肥后LIS1、LIS3和LIS4这3个基因的表达水平有所上调，但4个LIS基因的表达水平总体上仍较低，说明LIS基因也是芳樟醇合成的限速酶。值得注意的是，杨艳萍[34]认为DXR是萜类合成中的第2个限速酶，但本研究结果则显示：在施肥前后芳樟叶片与芳樟醇合成相关的调控基因中，DXR基因的表达水平均最高，且在施肥后枝中该基因的表达水平也升至最高，推测DXR基因参与芳樟芳樟醇合成过程的调控，但并不起限速作用，这种现象可能与DXR基因的表达具有种类特异性有关[30]。

芳樟醇的形成、积累和转化过程非常复杂，受到遗传和环境因子的综合影响。黄秋良等[35]的研究结果表明：施用微生物菌剂能够促进芳樟生长和精油积累；而在本研究中，施用生物炭复合肥也可显著提高芳樟叶片和枝中芳樟醇的相对含量，说明施用合适的肥料对植物次生代谢产物合成有一定的促进作用。芳樟叶片和枝中的芳樟醇相对含量在施肥后第3个月达到最高，在施肥后第5个月不同程度下降，这一现象可能与气温变化及芳樟生长期有关。施肥后第3个月正好是9月中旬，气温开始下降，植株即将进入一个生长周期的末期，生长和生理代谢速率均减缓，使芳樟醇的合成同步减缓。并且，施肥后芳樟叶片和枝中的芳樟醇相对含量提高，多数基因的表达水平也不同程度上调，尤其是叶片中的MDS和IDI2基因以及枝中的GPPS1基因的表达水平显著上调，说明施用生物炭复合肥可能对这些基因的表达有一定的影响。

综上所述，施用生物炭复合肥后，芳樟叶片和枝中芳樟醇的相对含量均显著提高，多数基因表达水平也不同程度上调，说明生物碳复合肥可能影响芳樟醇合成调控基因的表达特性，从而影响芳樟醇含量，但其中关键基因的功能还有待进一步验证。