油茶抗性无性系茶皂素含量对茶籽象幼虫肠道菌群结构的影响*

张守科 李子坤 殷 昊 张 威 舒金平 王浩杰 周旭东 汪阳东2,

(1.浙江农林大学林业与生物技术学院 杭州 311300； 2.浙江省林木育种技术研究重点实验室 杭州 311400； 3.中国林业科学研究院亚热带林业研究所 杭州 311400)

植物化学抗性是林木抗性品种选育与优化研究中的重要目标性状。从提高植物自身抗性入手来解决病虫害问题，是符合绿色防控要求、经济有效的手段(张杰等， 2019； 王晓杰等， 2020)。在经济作物品种抗性选育、优化研究如火如荼的大背景下，我国南方地区重要的木本油料经济作物油茶(Camelliaoleifera)抗虫品种选育工作也已全面开展(李苗苗等， 2017； 张守科等， 2020)。由于前期研究主要专注油茶经济性状筛选，大量纯林推广后，病虫害发生面积也随之大面积增加(张守科等， 2020)。尤其是以蛀果害虫茶籽象(Curculiochinensis)的暴发尤为严重，严重制约了油茶产业可持续健康发展(舒金平等， 2013； 张守科等， 2020)。化学手段防治钻蛀类害虫效果一般且不利于油茶林生态系统及产品的绿色无公害化(舒金平等， 2013)。因此，具备抗虫特性的优良油茶品种的选育以及绿色有效防控技术的研发是当前亟待开展的工作。然而，要达到上述要求则需先行开展茶籽象对油茶抗性响应机制的探究。该研究不仅可为油茶抗虫品种选育提供理论基础，也可为绿色防控技术的集成与应用提供支撑。

茶籽象是山茶属植物果仁内唯一的蛀果害虫，隶属鞘翅目(Coleoptera)象甲科(Curculionidae)象虫属(Curculio)(Zhangetal.， 2018； 张守科等， 2019)。研究者前期针对茶籽象生物学特性、种群动态及防治技术开展了大量探索(舒金平等， 2013； 何立红等， 2014； 李志文等， 2014； 李苗苗等， 2017； Zhangetal.， 2018； 张守科等， 2019)，然而关于茶籽象对油茶化学抗性的适应机制研究鲜有报道。前期开展的山茶属植物抗茶籽象研究工作，基本明确了茶皂素(Tea saponin)是影响山茶属植物茶籽象肠道菌群结构及为害率的主要化学抗性指标(Zhangetal.， 2020)； 且在探究不同油茶无性系间茶皂素含量与5个主推油茶无性系茶籽象为害率关系时也发现茶籽象为害率与茶皂素含量呈现负相关关系(张守科等， 2020)。但是基于4种山茶属植物不同物种对茶籽象肠道菌群结构影响的研究，并未完全考虑物种间遗传背景差异大以及植物化学抗性化合物积累量甚至构型都存在差异等因素，这使得以茶皂素含量作为评价油茶抗虫性强弱的依据存在不确定性。

本研究在前期明确不同山茶属植物茶皂素含量对茶籽象肠道微生物结构产生影响的基础上，进一步探究茶籽象肠道菌群结构对5个茶皂素含量存在差异的油茶无性系的响应。通过将茶籽象回接于普通油茶无性系上，可在排除寄主遗传背景差异的前提下，通过对比油茶无性系茶果内茶籽象老熟幼虫肠道微生物相对丰度、多样性及组成差异筛选差异菌株并初步预测功能，同时初步评价不同菌株在适应寄主化学抗性中的作用。本研究可明确油茶抗性无性系茶皂素含量是否对茶籽象肠道菌群产生影响，从而进一步确定茶皂素含量是否可作为油茶抗虫育种的重要植物化学抗性性状。本研究不仅为油茶抗虫育种提供理论依据，还可为以茶皂素这种植物源毒素为基础的茶籽象生态防控技术研发提供更为直接的理论支持。

1 材料与方法

1.1 试验样地及油茶林概况

选取浙江省青田尖山油茶种植园(28°11′51.61″N，120°23′15.25″E)作为试验样地。该园油茶种植年限超过10年且无人为干扰。试验地均属于亚热带季风气候，为油茶纯林，油茶林基本处于无人为喷施农药防治病、虫、杂草的自然状态，林下分布芒萁(Dicranopterispedata)、杠板归(Polygonumperfoliatum)及禾本科(Poaceae)狗牙根(Cynodondactylon)、白茅(Imperatacylindrica)等杂草。园中多个无性系间混种，均有结实。

1.2 供试植物

随机选择样地内具有茶皂素含量梯度的5个具备优质农艺性状的高产无性系(长林3号、 长林40号、 长林53号、 长林55号和长林166号)。将选取的不同无性系茶果收集至实验室分别放置并标记，待茶籽象幼虫发育老熟后，将获取的幼虫置于灭菌土中化蛹，供试。

1.3 茶籽象回接、茶果采集及老熟幼虫肠道解剖处理

选取生长状况一致的普通油茶5种无性系各10株。为确保所有的果实都没有被茶籽象产卵而影响试验数据，每棵树都用2 m × 2 m × 2 m的小孔透明塑料网密封，仅留一个拉链开口供后期接种成虫使用。5月下旬，随机选择实验室无菌土中10对同日羽化的茶籽象成虫(♀∶= 1∶1)接于网罩的供试植株上，任其交配产卵。待霜降前，选定的油茶5个无性系进行随机采样，每株分为上、中、下3个层次，东、南、西、北4个方位，每个层次每个方位随机摘取20个茶果，每株共计240个茶果，不足240个的植株则全部采摘。将获取的成熟果实分组编号，带回实验室。室温下，静置于实验室塑料筐内，待老熟幼虫从果实逸出并统计数量后，将茶果果仁取出，并参考张守科等(2020)方法测定茶皂素含量。存活的老熟幼虫每个无性系随机选取36头，于磷酸缓冲液中解剖肠道，2条肠道为1个样本，每个无性系共18个样本，编号后于- 80 ℃保存备用。

1.4 茶籽象肠道菌群总DNA提取及建库测序

茶籽象肠道微生物总DNA的提取，参考QIAGEN DNeasy Blood & Tissue Kit试剂盒细菌DNA提取方法： 在PBS缓冲液中将茶籽象肠道解剖出来，搜集到装有1 mL PBS的无菌PE管中，每管含有3头茶籽象的肠道以保证菌液DNA质量，匀浆； 将匀浆后的菌液7 500 r·min-1离心，取上清液至新的无菌PE管后提取总DNA。收集到的DNA 原液经1% 琼脂糖凝胶电泳，100 V电泳40 min后溴酚蓝到达琼脂4/5处，取出至UV凝胶电泳成像仪进行条带检测，并检测DNA样品纯度和浓度。采用Illumina公司的TruSeq Nano DNA LT Library Prep Kit试剂盒构建测序文库。首先对扩增产物进行序列末端修复，对文库做最终的片段选择与纯化。上机测序前，需要先对文库在Agilent Bioanalyzer上进行质检，采用Agilent High Sensitivity DNA Kit试剂盒。合格后，采用Quant-iT PicoGreen dsDNA Assay Kit试剂盒在Promega QuantiFluor荧光定量系统上对文库进行定量后将合格的样品上机测序； 使用MiSeq测序仪进行2×300 bp的双端测序，试剂为MiSeq Reagent Kit V3(600 cycles)。由于MiSeq测序读长较短的特性，目标片段测序长度定为200～450 bp。

1.5 数据分析

利用质控软件Cutadapt (V1.9.1)将所有16S序列测序质量检验，并根据样本间不同的引物进行分组。去除条形码及引物后，通过UCHIME算法对端序列进行合并和质量过滤。

利用UPARSE软件将序列聚类为分类单元(operational taxonomic units,OTUs)，序列相似度阈值为 >97%，同时删除一部分嵌合体以保留代表性OTUs。在QIIME2软件中将OTU序列与对应数据库序列比对，利用Silva数据库注释获取每个OTUs所对应的分类学信息。将抽平后的OTUs矩阵用来后续分析。依照OTUs划分和分类地位鉴定结果，获取各样本在每个分类水平的物种组成。比较样本在各分类水平的微生物类群数量，使用R软件绘制成分类等级树，比较各样本在同一水平的分类单元数的差异。使用QIIME2软件对每个样本计算Alpha多样性指数的Chao1、Good_coverage、Observed_species、Shannon、Faith_pd、Simpson及Pielou_e指数。Beta多样性分析主要检验不同样本之间群落结构的相似性。主要通过非度量多维标度分析(non-metric multidimensional scaling，NMDS)，对群落数据结构进行自然分解并通过对样本排序，从而观测样本之间的差异及相似性。并根据样本Bray_curtis距离矩阵，对不同分组的样本距离均值进行统计学比较Anosim差异检验(Analysis of similarities test)。使用R软件比较分析高含量茶皂素无性系长林40 + 长林166组与低含量茶皂素无性系长林55 + 长林53 + 长林3组茶籽象肠道菌群物种差异及标志物种并绘制曼哈顿图。

2 结果与分析

2.1 茶籽象存活量与油茶抗性无性系茶皂素含量的关联

为进一步明确茶籽象危害率与普通油茶茶皂素含量的关联，在回收5个油茶无性系茶果时，对茶籽象存活数量进行统计，并与无性系茶果果仁茶皂素积累量进行关联分析。通过对霜降后油茶果果仁茶皂素含量测定，茶皂素含量为3.27～9.62 mg·g-1,长林3号与长林53号茶皂素含量差异不显著(P< 0.05)，范围在(3.61± 0.34) mg·g-1～ (3.75 ± 0.23) mg·g-1(图2A)；而长林55号茶皂素含量为(5.52 ± 0.33)mg·g-1，长林40号和长林166号茶皂素的最高含量分别为(7.75 ± 0.16)mg·g-1和(9.47± 0.15)mg·g-1。对不同无性系幼虫存活量统计分析，虫口数量为26～85个(图1)。长林55、53及3号的幼虫存活量差异不明显(P> 0.05)，分别为(77 ± 8)个、(70 ± 10)个及(69 ± 11)个，长林40和长林166号的幼虫存活量相对较少，分别为(43 ± 6)个和(34 ± 8)个，两者差异不显著(P> 0.05)。关联分析表明2组数据的关系为y= 0.252 5x2-10.096x+ 107.47(R2= 0.718 9)，茶籽象存活量呈现随茶皂素含量上升而下降的趋势(图1)。

图1 茶果茶皂素含量与幼虫存活量相关回归分析Fig.1 Correlation regression analysis of tea saponin content and larval survival

图2 普通油茶供试无性系及茶籽象肠道菌群属水平注释结果Fig.2 Clones of Camellia oleifera involved in the experiment and annotated results of weevils’ gut microbiota at the genus level

2.2 油茶抗性无性系对茶籽象肠道菌群结构的影响

对5组茶籽象肠道样本抽平后OTU表属水平分类数量统计发现，长林40号和长林166号茶果内茶籽象肠道菌群相似度较高，其中两者均含有较高水平的不动杆菌属(Acinetobacter)。其中，长林40号中不动杆菌属相对丰度为6.23%，长林166号不动杆菌属相对丰度为10.19%(图2B)。相较于另外3组，长林40和166组内沙雷氏菌属(Serratia)相对丰度较低，分别为5.81%和6.04%； 而长林55、53、3号沙雷氏菌属则为优势菌属，分别为13.41%、16.56%和14.68%(图2B)。另外，长林40号和166号中，假单胞杆菌属(Pseudomonadaceae)相较于另外3组，相对丰度也相对较高，分别为4.94%和3.74%(图2B)。将5组样本进行物种组成最高分类水平为属水平的分类等级树分析发现，前100丰度的OTUs中，丰度前10的OTU的注释结果与前者分析结果一致(图2C)。其中，大部分菌5组占比基本一致，仅不动杆菌属、假单胞杆菌属所包含的OTUs在长林40和166样本中占比较高(图2C)。而沙雷氏菌属、阪岐肠杆菌属(Cronobacter)、茎菌属(Caulobacter)、食酸菌属(Acidovorax)、草酸杆菌属(Aquabacterium)、固氮螺菌属(Azospirillum)以及魏斯氏菌属(Weissella)所包含OTUs多分布于长林55、53及3号茶籽象肠道菌群样本中(图2C)。

2.3 油茶抗性无性系对茶籽象肠道菌群多样性的影响

为明确5种无性系对茶籽象肠道菌群多样性的影响，分别对茶籽象肠道菌群进行了Alpha多样性分析、Beta多样性分析以及基于Bray_curtis距离矩阵的组间Anosim差异分析。其中，Chao1、Good_coverage、Observed_species、Shannon、Faith_pd、Simpson及Pielou_e指数等7个Alpha多样性指数分析中仅Shannon、Faith_pd、Simpson及Pielou_e等4个指数P值小于0.05(图3)； Shannon指数分析中，长林166组肠道样本多样性与长林40号样本差异不大(P> 0.05)，与其他各组均存在极显著性差异(P< 0.001)； 而长林40号样本仅与长林3号样本肠道菌群多样性存在显著性差异(P< 0.01)； 长林55、53及3号样本中Shannon指数差异不大(P> 0.05)(图3)。Faith_pd指数分析中，长林166号样本与长林3，长林53号样本与长林55，长林40号样本与长林3差异极显著(P< 0.001)； 长林55号样本与长林3号样本肠道菌群多样性存在较显著差异(P< 0.05)(图3)。与其他指数不同，Simpson指数中长林166组与其他各组均有显著差异(P< 0.05)，而其他剩余4组之间则多样性差异不大(图3)。而Pielou_e指数分析中，长林40组样本多样性仅与长林166组差异不大，其他各组均差异显著(P< 0.05)，长林166组样本与除长林40组之外的3组样本也存在显著差异(P< 0.05)(图3)。

图3 不同油茶抗性无性系上茶籽象肠道菌群Alpha多样性差异比较Fig.3 Comparison of Alpha diversity of gut microbiota of C. chinensis in different C. oleifera clones

为进一步比较5个油茶无性系果实对茶籽象肠道菌群多样性的影响，对样本间的Bray_curtis距离矩阵和Jaccard 距离矩阵进行了NMDS分析(图4)。基于Bray_curtis距离矩阵的分析发现，高茶皂素含量的长林166和40号上的茶籽象肠道菌群结构相似度较高，在主成分NMDS1上可以将5组样本显著分成2个类群(Stress = 0.052 < 0.1)(图4)。同样的利用Jaccard 距离矩阵进行的分析表明，在主成分NMDS1上也可将5组样本聚类成2个大的类群，高茶皂素含量的长林166和40号样本单独聚类，剩余3组单独聚类，但优劣用胁迫系数(Stress = 0.16)相较于Bray_curtis距离矩阵较高； 虽在此胁迫系数下依然可以将5组样本很好地聚类，但相较于Bray_curtis距离矩阵解释度不够(图4)。为验证5组样本间Beta多样性之间差异的显著性，选取NMDS分析中优劣用胁迫系数(Stress = 0.052)较低的Bray_curtis距离矩阵展开分析。其中，长林40分别与长林55、53和3存在显著差异(P< 0.05)，且R> 0.9 > 0.75，说明样本内肠道菌群差异极大(P< 0.05)，相似度极低(P< 0.05)； 而长林166组同样也与长林55、53和3存在显著差异(P< 0.05)，这同样佐证Alpha多样性分析结果(表1)。

表1 基于Bray_curtis距离的组间Anosim差异分析Tab.1 Analysis of Anosim differences between groups based on Bray_curtis distance

图4 基于2种距离矩阵对不同油茶抗性无性系上茶籽象肠道菌群的NMDS分析Fig.4 NMDS analysis of gut microbiota of C. chinensis in different C. oleifera clones based on two distance matrices

2.4 茶皂素含量差异对茶籽象肠道菌群属水平物种差异及标志物种的影响

为比较长林166和40号高茶皂素含量组肠道菌群相对于长林55、53和3号茶皂素低含量组样本显著富集可能与降解或耐受茶皂素有关的核心菌株，使用曼哈顿图来展现差异OTUs，并通过注释信息获取差异OTUs的属水平注释信息。筛选出51个OTUs在高含量茶皂素组的茶籽象肠道内显著富集，其中，有1个OTU属于黄杆菌属，1个OUT属于固氮螺菌属，5个OTU属于阪岐肠杆菌属(Cronobacter)，31个OTU属于假单胞杆菌属，13个OTU属于不动杆菌属(图5)。

图5 高含量茶皂素无性系与低含量茶皂素无性系茶籽象肠道菌群物种差异及标志物种分析Fig.5 Species difference and marker species analysis of gut microbiota in C. oleifera clones with high and low content of tea saponin坐标系内每个圆点或圈代表1个ASV(分类单元)，大小代表其相对丰度； 虚线分隔了显著差异(P < 0.05)与不显著的ASV，显著差异的点用彩色圆点或圆环标志，不显著的用灰色圆环表示，该分组内显著上调的用彩色实心圆点显示。Each dot or circle in the coordinate system represents one ASV (taxon), and its size represents its relative abundance. The dotted line separates the significant difference (P < 0.05) from the non-significant ASV, the significant difference points are marked by colored dots or rings, the non-significant points are indicated by gray rings, and the significantly upregulation points within the group are shown by colored solid dots.

3 讨论

在植物与植食性昆虫的协同进化关系中，植物为抵抗昆虫为害进化出包括物理和化学的各种抗性表型(Mason， 2020)。植物化学抗性是由植物应对植食性昆虫为害应激反应产生的或者长期稳定积累的有毒、抗营养和难消化的化合物组成，一般具有趋避害虫、损伤或抑制昆虫消化系统对食物消化的能力(Hammeretal.， 2015； Erbetal.， 2019； Yuanetal.， 2021)，起到化学防御的作用(Agrawal， 2019； Masonetal.， 2019)。植食性昆虫也相应地采取一系列对策以适应来自植物寄主抗性的选择压力(Schumanetal.， 2016； Agrawaletal.， 2019； Mason， 2020)，其中就包括对植物寄主产生明显的取食偏好反应以及适应限度差异(Masonetal.， 2019； Liuetal.， 2020)。

在专食性象虫这类昆虫的进化进程中，适应限度差异现象尤为突出(Morera-Margaritetal.， 2021)。前期研究发现，茶籽象对山茶属植物的适应限度差异现象表现在发育状况及存活率上(Zhangetal.， 2020)。由于以山茶属属下阶元观测这种适应限度差异现象并未完全排除物种遗传背景差异大的干扰，因此，为探究这一现象并进一步评判茶皂素是否可以作为油茶抗茶籽象的重要植物抗性成分，本研究观测了5个不同油茶无性系对茶籽象老熟幼虫存活量及老熟幼虫肠道菌群结构的影响。通过对茶籽象幼虫个数与茶皂素含量之间的关联分析，发现油茶无性系茶皂素含量对茶籽象幼虫存活量有显著影响。此结果与前期山茶属之间为害率不同的结果基本一致，幼虫在低茶皂素含量的茶果内相较于高茶皂素含量的茶果内发育更为良好(Zhangetal.， 2020)。这进一步论证了“茶皂素是山茶属植物乃至油茶无性系抑制茶籽象发育的重要植物化学抗性成分”的结论。

适应限度差异的产生，除植食性昆虫消化道内自身消化酶、解毒酶发挥作用外(Grbietal.， 2011； Xuetal.， 2019； Guietal.， 2020； Baietal.， 2021)，肠道菌群解毒功能在昆虫抵消植物化学抗性方面的作用也不可忽视(Agrawaletal.， 2019； Masonetal.， 2019； Guptaetal., 2020； Mason， 2020； Shuklaetal.， 2020)。研究发现，昆虫肠道内微生物不仅可以帮助宿主消化营养物质，还可帮助植食性昆虫快速适应一定量的植物有毒次生代谢物(Douglas， 2013； Chengetal.， 2018； Morera-Margaritetal.， 2021)，有助于昆虫在一定程度下有效分解或规避寄主植物有毒次生代谢物，健康地完成营养积累过程(Hansenetal.， 2014； Masonetal.， 2019； Francoeuretal.， 2020； Mason， 2020； Zhangetal.， 2020)。肠道微生物菌群能够促进宿主昆虫适应寄主，并介导了昆虫代谢化学防御增强取食消化植物寄主的能力(Nikohetal.， 2011； Adamsetal.， 2013； Chenetal.， 2018)。前期研究也发现，不同山茶属植物上茶籽象肠道营养代谢相关的微生物在消化植物物质的过程中会受到山茶属植物的影响，可能在无法正常消化营养物质时，虫体无法获取足够的营养物质而营养不良，导致为害率差异的现象(Zhangetal.， 2020)。

研究通过将同一来源的茶籽象回接于茶皂素含量存在差异的油茶无性系上，直接探讨茶皂素对茶籽象肠道菌群结构的影响。不同寄主回接后获取的老熟幼虫茶籽象肠道菌群结构变化差异巨大(图2、3，表1)，从分析菌群Alpha多样性、Beta多样性分析结果来看，茶籽象肠道微生物菌群结构及数量对无性系茶皂素含量反应明显，高含量茶皂素的长林166和长林40回收的茶籽象肠道内与有毒物质代谢相关的不动杆菌属、假单胞杆菌属占据较高权重； 与营养消化功能相关菌属则被抑制； 结合前期研究(Zhangetal.， 2020； 张守科等， 2020)，本研究中5个油茶无性系对茶籽象幼虫肠道菌群结构的影响以及老熟幼虫存活量的差异，进一步明确普通油茶茶果茶皂素含量是影响肠道菌群的化学抗性组分，同时不动杆菌属则是与茶籽象耐受茶皂素毒性密切相关的优势肠道菌属。

4 结论

本研究利用物种茶皂素含量存在差异的油茶无性系，探究了茶皂素含量对茶籽象幼虫肠道菌群结构的影响，发现高含量茶皂素的无性系可以明显影响茶籽象肠道菌群多样性并造成肠道菌群组成结构差异。对茶果果仁茶皂素含量与幼虫存活量的相关性分析表明，油茶无性系茶皂素含量可以作为油茶抗虫品种选育的一个重要指标。这一结论也佐证了笔者前期的研究结果。本研究的结论进一步明确了茶皂素含量在抗茶籽象上的作用，为油茶抗虫品种选育工作以及开发新的防治茶籽象的绿色防控技术提供了理论和技术支持。