不同预处理方式下肉桂枝叶内生细菌物种差异分析

曾维星 程 贤 毕良武 卢言菊 陈玉湘 赵振东 张笮晦

（1. 中国林业科学研究院林产化学工业研究所，生物质化学利用国家工程实验室，国家林业和草原局林产化学工程重点开放性实验室，江苏省生物质能源与材料重点实验室，江苏 南京，210042；2. 南京林业大学江苏省林业资源高效加工利用协同创新中心，江苏 南京，210037；3. 广西庚源香料有限责任公司，广西 防城港 538100）

肉桂（Cinnamomum cassia）系樟科（Lauraceae）樟属（Cinnamomum）植物，是我国人工种植的经济树种，主要分布于广西以及广东2 省区[1-2]。人工种植的肉桂一般生长7～10 年后进行砍收，将肉桂枝干剥皮，经过晾晒风干后得到了重要传统的辛香料“桂皮”；而其余的肉桂枝叶则可以通过水蒸气蒸馏法制得另一个重要的工业产品“肉桂油”[3]。根据行业内加工企业长期积累的经验，肉桂枝叶在蒸馏提取之前需要先在通风阴凉处存放45～60 d，肉桂油提取得率大幅增加。前期研究发现，将酵母菌溶液喷洒在采摘后的肉桂枝叶表面，可以在短时间内增加肉桂油的提取得率[4]，初步推断肉桂枝叶采收后，内生菌继续影响肉桂醛等精油成分的内在生物转化。植物内生菌指至少有一段生命周期存活于植物组织中，并对宿主没有负面影响的一类微生物[5-7]。这一类微生物群落大量分布在植物的根、茎、叶、花、果实和种子等器官、组织的细胞或细胞间隙，近些年来，研究学者在许多农作物和经济作物中发现且分离出各种各样的内生菌[8-19]。而且在不同的植物中，内生菌在种类以及物种量上相差甚远，差异性显著，即使是在相同的一株植物中，其内生菌在种类和物种量上亦会存在差异，并且内生菌所处的生活环境以及条件的改变同样会影响甚至改变内生菌的种类和物种量[20]。最近几年，植物内生菌的生物转化作用机制逐渐受到人们的重视。刘斌等对肉桂内生菌转化肉桂醛进行了研究，该研究发现肉桂内生菌Pseudomonassp. RGEB06 具有较好的催化肉桂醛加氢生成肉桂醇的能力[21]。因此，肉桂枝叶的生物预处理是否影响了肉桂枝叶的内生菌，进而增加肉桂油的提取得率，值得进一步探索。本研究通过系统的差异分析，旨在发现生物预处理手段对肉桂枝叶内生菌群落结构的影响，不仅为生物预处理手段干预肉桂枝叶的作用机制研究提供指导，还可以为开发利用微生物资源提供依据。

1 材料与方法

1.1 材料与试剂

采摘放置1 周的肉桂枝叶，取自广西梧州市藤县金鸡镇；实验室自来水、蒸馏水；酵母粉，购自湖北安琪酵母股份有限公司；酵母液，酵母粉和蒸馏水按质量体积比1∶2 配制；无水乙醇、棉球、MP-soil试剂盒，购自南京荣世德贸易有限公司；75%乙醇，由无水乙醇和蒸馏水按体积比3∶1 配制。

1.2 仪器与设备

Milli-Q 纯水仪，美国Millipore 公司；SWCJ-2F 型双人双面净化工作台，购自苏州净化设备有限公司；MISEQ 测序仪Illumina Miseq，美国Illumina 公司。

1.3 生物预处理

在同一天，均各取大小均匀的肉桂枝叶30片，分别用于如下3 种处理：

空白组（C）：叶片不做任何预处理，静置1 周，备用。

水预处理组（W）∶按m（肉桂枝叶）∶m（水）=5∶1 的比例将水均匀的喷洒在枝叶表面，放置在阴凉干燥处，静置1 周。

酵母液预处理组（Y）：配制0.1 g/mL 的酵母水溶液，按m（肉桂枝叶）∶m（水）=5∶1 的比例将酵母水溶液均匀的喷洒在枝叶表面，放置在阴凉干燥处，静置1 周。

1.4 提取内生菌

取空白组的30 片肉桂枝叶，将表面用蒸馏水冲洗干净，再用棉球蘸取75% 乙醇擦拭肉桂枝叶，除去枝叶表面的微生物，待叶片表面水分蒸发后，用剪刀将所有肉桂枝叶均匀剪碎，平行分为6 份，置于已杀菌的自封袋中密封，6 份样品命名为C1～C6。用同样的方式平行提取水预处理组（6 份样品命名为W1～W6）与酵母液预处理组（6 份样品命名为Y1～Y6）。最后使用MP-soil试剂盒，按照说明书流程对上述18 份样本进行总DNA 的提取，至此完成了3 组内生菌的提取。

1.5 数据处理

数据的处理借助美吉生物云平台。

1.5.1 数据去杂

通过操作窗口设置1 个平均值，当序列经过窗口后直接截去其尾端碱基，从而达到过滤掉所读取到的序列末端质量值小于20 的碱基的目的；然后对交叠的序列进行拼接；拼接好的序列交叠部分最大错配比率是0.2，高于这个比率的序列直接弃用；序列的首尾两端存在靶向条码，可用于区分各个样本以及辨别序列的方向，条码一般不允许匹配出错，最大引物匹配出错数是2。

1.5.2 得到最优数据

最后读取得到的首尾两端序列数据，两两相互匹配成读取片段拼接为一条序列，并且对这个过程严格进行质检，根据序列首尾两端的靶向条码和引物序列加以区分样本，并且对有效序列进行正序列方向处理以保证序列的质量，得到的数据最优。

1.5.3 数据评估分析

对18 个样本的最优数据序列分别统计所有小组总数以及共有小组数进行曲线分析（Pan/Core分析）评估。为了对所测得序列数据进行注释与评估，本研究采用OTU（operational taxonomic units）标志方法，分析统计测序数据，因为样本较多，所以对每组各个样本进行聚类，这样利于清楚直观的观察其中任一个样本测序结果中物种、物种数等的关键信息。经过聚类的样本，其测序序列相接近的就归类为同一个小组（OTU）。标准是一个OTU 的相似性为97%，统计所有样本的小组总数（Pan OTU）以及所有样本的共有小组数（Core OTU），进行曲线分析。最后通过判断Pan/Core OTU 物种曲线以及稀释曲线是否达到平缓，从而评估这次供测序的样本量是否足够。利用Venn 图进行物种组成分析，统计参试样本中共有以及特有物种和物种数，以及借助层级聚类分析方法分析不同样本间内部的群落组成、结构的相似和差异性，并绘制了树状图。

1.5.4 组间数据显著差异性检验

依据对群落丰度的分析统计，对不同样本间物种群落进行假设检验，对其中的单个或多个物种进行显著性等级评估，得到不同样本间存在差异显著性的物种；分别采用Wilcox 秩和检验（Wilcoxon rank-sum test）和克氏秩和检验（Kruskal-Wallis H test）同时对2 组或3 组样本进行差异显著性检验。多级物种差异判别分析，借助LEfSe 软件（发现高维生物标识和揭示基因组特征），检测存在显著差异的特征物种丰度，并且定向追踪存在显著差异的特征物种族群，该软件使用判别分析（LDA）评估各个物种丰度对显著差异影响的大小进行可视化。

2 结果与分析

2.1 测序数据统计

本研究对18 个待测样本完成测序后，每个样本测序所得序列数、碱基数和平均长度见表1，18 个样本的序列数与碱基数分别在1.5 万～2.5 万条、574 万～915 万条范围内，序列的平均长度在377 bp 左右，最短序列长度为205 bp，最长序列长度为488 bp。上述样本测序信息表明：同组样本的测序结果平行性良好，实验结果真实可靠。

表1 3 组样本内生菌测序数据统计表Table 1 Sequencing data statistics of endophytic bacteria in 3 groups of samples

2.2 物种注释与评估

对18 个样本进行Pan/Core 分析和稀释曲线分析得到图1，和图2，图1 中横坐标代表样本量，图1 左纵坐标代表总的小组数，图1 右纵坐标代表共有的小组数，图2 横坐标代表随机抽取的测序数据量，纵坐标代表观测到的物种数量。从图中可以看到，Pan and Core 曲线分析随着样本数量增加到6 个时曲线趋于平缓，而稀释曲线图中，随着随机抽取的测序物种数增加到2 万个时曲线也趋于平缓，所以根据分析结果验证本次测序样本量足够，物种丰度高。

图1 3 组样本物种OTU Pan and Core 分析曲线Fig. 1 OTU Pan and Core analysis curves of 3 groups of sample species

图2 3 组样本物种稀释曲线Fig. 2 Dilution curves of 3 groups of sample species

2.3 物种组成分析

2.3.1 Venn 图分析

借助Venn 图能够很清楚地展现不同环境下样本中物种组成的相似性和交叠部分。通过对物种组成进行Venn 图分析得到图3，图中重叠部分代表组间共有物种数，非重叠部分为各组特有物种数。从图中可以看出测得3 组物种数均超过300 个，其中C 组物种数337 个，W 组物种数345 个，Y 组物种数386 个。3 组共有物种数210 个，其中C 组与水处理组共有物种数40 个，C 组与Y 组共有物种数31 个，W 组与Y 组共有物种数45 个。具有显著差异的一项数据是3 组特有物种数，C 组特有物种数为56 个，W 组特有物种数为50 个，而Y 组特有物种数达到了100 个。从Venn 图可以很直观地看出经过酵母液预处理的肉桂枝叶内生菌物种组成发生了较大变化，推测酵母液预处理这种生物预处理方式可能会影响肉桂枝叶陈化转化的过程。

图3 3 组样本物种组成Venn 分析Fig. 3 Venn analysis chart of species composition of 3 groups of samples

2.3.2 群落组成分析

本研究在门的级别上对3 组18 个群落多级物种分析，借助柱状图和热图直观展现样本中群落的组成差异，同时在根据不同的分类等级对各个样本进行物种丰度差异分析如图4，图中不同颜色代表不同物种，横坐标代表不同分组，纵坐标代表在门的层级上各物种在群落组成上的丰度占比。结果表明3 组样本中的变形菌门（Proteobacteria）、放线菌门（Actinobacteriota）、厚壁菌门（Firmicutes）、拟杆菌门（Bacteroidota）丰度较高，并且这4 种菌门在3 组样本群落中总丰度占比达到了96%；其中3 组样本中的变形菌门经过水和酵母液预处理后丰度降低，而放线菌门、厚壁菌门、拟杆菌门丰度升高，特别是厚壁菌门在经过酵母液预处理后丰度大幅升高。然后对丰度前20 的物种进行热图分析，如图5，图中颜色梯度代表物种丰度，同一颜色梯度下颜色深浅代表该物种丰度变化，对比C 组发现W 和Y 组的酸杆菌门（Acidobacteriota）、装甲菌门（Armatimonadota）、粘球菌门（Myxococcota）、杆菌门（Patescibacteria）等丰度发生显著变化。

图4 群落组成分析Fig. 4 Analysis on community composition

图5 3 组样本群落多级物种组成Heatmap 图Fig. 5 Sunburst diagram of multi-level species composition of 3 groups of sample communities

2.4 参试样本比较分析

本研究对3 组18 个样本进行样本比较分析得到图6，图中横线长度代表各个样本之间的距离，不同的颜色用以区分不同的样本组。从图中可以看出18 个样本间群落组成相似性的有Y4和Y6，C1和W6，C2和C4，C3和C6，Y3和W1，W5，W2和W3，C5和Y1，Y5；而具有显著差异性的是Y2和W42 个样本，初步推测可能是这2 个样本中的物种群落组成更为丰富、多样。

图6 3 组18 个样本OTU 级别的层次聚类树状图Fig. 6 Hierarchical clustering of multi-level species composition of 18 samples in 3 groups

2.5 物种差异分析

2.5.1 组间差异显著性检验

图中P值为假阳性概率值，通常情况下P＜0.05 时差异显著。在门的级别上对3 组样本的群落物种进行克氏秩和检验，得到图7，横坐标代表物种不同分组中平均相对丰度，不同颜色的柱子表示不同分组，纵坐标代表不同物种。从图中可以看出变形菌门、厚壁菌门、蓝藻门（Cyanobacteria）的P值小于0.05，为显著差异，可以推断这3 个群落物种在3 组样本中的丰度、物种组成、物种数量可能存在很大差异。再单独进行两组互为比较，得到物种丰度前20 的P值，如表2，从表中可以看出W 组和Y 组对比C 组，其中变形菌门、厚壁菌门的P＜0.05，为显著差异；而W组与Y 组比较发现厚壁菌门的P＜0.05，为显著差异。

表2 3 组样本丰度前20 的群落物种Wilcox 秩和检验结果Table 2 The Wilcox rank sum test results of the top 20 community species in 3 groups of sample abundance

图7 3 组样本的群落物种克氏秩和检验结果Fig. 7 The results of Kruskal-Wallis H test of community species of 3 groups of samples

2.5.2 多级物种差异判别分析

对W 组和Y 组进行多级物种差异判别分析得到图8，图8 中不同字母代表不同物种，不同的颜色节点对应不同分组，节点距离圆心越近表示影响越大。从图中可以看出，在W 组中，e—环脂酸芽孢杆菌科（Alicyclobacillaceae）、n—鞘脂杆菌目（Sphingobacteriales）、y—立克次氏体门（Rickettsiales）3 个群落物种对组间差异存在最显著影响，而甲基杆菌门（Methylobacterium）的丰度对差异效果影响最大；在Y 组中，a—厚壁菌门对组间差异存在最显著影响，同时其丰度对差异效果影响也是最大。

图8 2 组样本LEfSe 多级物种差异分析Fig. 8 LEfSe multi-level species difference analysis of 2 sets of samples

3 结论与讨论

1）本研究采用基因测序手段，对经过特殊处理的3 组肉桂枝叶样本的内生菌进行基因测序，对数据进行OTU/Pan and Core OUT/稀释曲线分析、评估，得到并统计了充足的有效数据。

通过物种组成分析发现具有显著差异3 组特有物种数，C 组特有物种数为56，W 组特有物种数为50，而Y 组特有物种数达到了100。3 组样本中的变形菌门、放线菌门、厚壁菌门、拟杆菌门丰度较高，并且这4 种菌门在3 组样本群落中总丰度占比达到了96%。初步推断这4 种菌门可能是对肉桂枝叶自然陈化转化的过程影响最大；同时经过酵母液预处理以后特有菌种数增加明显，证明这种生物预处理方式能够很大的改变肉桂枝叶的内生菌数量组成，亦或可能间接地影响甚至改变肉桂枝叶的出油率。

2）通过对3 组18 个样本进行比较分析，发现具有显著差异性的是Y2和W42 个样本，推测可能是这2 个样本中的物种群落组成更为丰富、多样，内生菌提取的更加充分完整。

3）通过物种差异分析发现3 组样本中变形菌门、厚壁菌门、蓝藻门为显著差异。W 组和Y 组对比C 组，变形菌门、厚壁菌门为显著差异，可以得出肉桂枝叶经过水和酵母液预处理后的内部的这两种菌类发生了显著性的变化，倘若经过这2 种生物预处理手段的肉桂枝叶的出油率发生变化，可以推断这2 种菌类可能是影响变化的因素。而W 组与Y 组比较发现厚壁菌门为显著差异。最后对WY 两组进行LEfSe 多级物种差异判别分析发现W 组中甲基杆菌门的丰度对差异效果影响最大；Y 组中，厚壁菌门对组间差异存在最显著影响，同时其丰度对差异效果影响也是最大。最终推断厚壁菌门可能是影响肉桂枝叶内部生物转化的最大因素。