基于微生物组扩增子ITS的茯砖茶发花特征真菌菌群演替规律分析

姚衡斌，马敬宜，冯瀚瀚，赵仁亮

（河南农业大学园艺学院茶学系 郑州 450046）

茯砖茶属于后发酵茶，微生物大量参与茯砖茶的后发酵过程[1]。茯砖茶的后发酵工艺与其它黑茶稍有不同，主要分为2 个阶段： 一是黑毛茶阶段，刚采摘的鲜叶通过一系列工艺制成黑毛茶。渥堆是该阶段的主要工艺，是湿热作用、微生物及酶共同作用[2]。二是发花阶段，黑毛茶经过气蒸、压制成型、发花等工艺步骤形成茯砖茶。发花也是一个微生物发酵过程，研究表明，发花阶段多种微生物迅速繁殖，茶砖中物质分解，代谢产物积累，微环境变化，菌群的不断演替，最终由散囊菌属（Eurotium）的1 种或几种成为优势菌种，也就是“金花菌”，“金花”即是这些“金花菌”的闭囊壳。

“菌花香”是优质茯砖茶的特征香气，是评价茯砖茶优劣的重要标准。 实践发现有“菌花”的茯砖茶不一定有“菌花香”。 金友兰等[3]分析发花白茶、发花红茶、茯砖茶的香气成分，结果发现不同茶类发花的香气物质存在显著差异。有研究表明[4]不同产地的茯砖茶“菌花香”也不相同。 推测可能由于不同的茶树品种和生长环境导致原料内涵物质的不同，意味着原料香气前体物质的种类及含量不同，且供后发酵微生物生长的碳源、氮源及多种微量元素的种类及含量不同，因此不同地区的茯砖茶发花微生物菌群差异显著[5]，提示微生物可能是产生茯砖茶“菌花香”的关键因素。刘石泉等[6-7]使用DGGE 技术揭示优势微生物是随着发花过程而动态改变，表明“菌花香”可能并不只是冠突散囊菌（Eurotium_cristatum）作用的结果，有其它微生物或微生物菌群进行代谢或提供胞外酶，多种微生物交替发酵共同形成“菌花香”。因此，对发花过程微生物的解析较为关键。

由于技术限制，过去对发花微生物的研究以优势菌株的分离，并通过形态学鉴定为主[8]，而微生物种类繁多，目前成功分离培养的微生物仅占全部的0.1%～1.0%[9]。后来发展出利用变性梯度凝胶电泳使不同的碱基排列的DNA 片段在凝胶上分离获得微生物多样性的PCR-DGGE 方法[10]。 研究人员使用DGGE 技术研究茯砖茶微生物的群落结构，鉴定出多种酵母菌属、曲霉属、及数种不可培养微生物[11]。 DGGE 技术的缺点是通量较低，不能完整反映复杂的环境微生物群落。 宏基因组学技术[12]利用二代测序技术[13]通量高的特点，将环境中的微生物全部DNA 提取分箱[14]，能够获得较为完整的环境微生物的基因组信息[15]。 然而，宏基因组数据量大，分析较为困难，可以仅扩增微生物的特殊片段[16]以达到获得环境微生物多样性及丰度的目的，这就是微生物组扩增子技术[17]。

白茶经鲜叶萎凋干燥而成[18]，在常温下可以长期保存[19]，和茯砖茶一样有紧压茶产品。 由于低级白茶所用原料与茯砖茶原料嫩度一致，因此选用白茶砖作为对照与黑茶一起发花制成发花白茶与茯砖茶。

本研究采用微生物组扩增子ITS，对信阳罗山同一茶叶原料制成的黑茶与白茶进行发花试验。其中，白茶作为对照，围绕发花过程中的特征菌群演替规律，运用多种分析方法，通过两种茶类发花的对比分析，揭示茯砖茶发花特征真菌菌群演替规律，为阐明“菌花香”与微生物菌群的关系，合理控制发花过程，为我国茯砖茶发花工艺改进提供参考。

1 材料与方法

1.1 样品采集与处理

本研究供试材料为申林茶业开发有限公司的薮北茶树品种[20]。 2020年6月18日采摘茶树鲜叶，采摘标准为1 芽5～6 叶。鲜叶采摘后分别使用黑茶及白茶初制工艺制成相应毛茶。 2020年12月13日将两种毛茶均潮水至含水量28%，然后汽蒸10 s，压制成砖茶。 立即随机取3 个平行样本，其中黑茶砖D1.1、D1.2、D1.3，白茶砖W1.1、W1.2、W1.3，将剩余茶砖置入全新的烘房，烘房温度28℃，湿度80%，分别在发花的第6 天取样（D6.1，D6.2，D6.3，W6.1，W6.2，W6.3），第9 天取样（D9.1，D9.2，D9.3，W9.1，W9.2，W9.3），第12 天取样（D12.1，D12.2，D12.3，W12.1，W12.2，W12.3），第22 天取样（D22.1，D22.2，D22.3，W22.1，W22.2，W22.3），取样后立即置于-80 ℃条件保存。

1.2 主要仪器与试剂

DYY-8C 电泳仪，北京六一生物科技有限公司；Tanon 3500 凝胶成像系统，上海天能生命科学有限公司。

高保真PCR 预混缓冲液和高效高保真酶，New England Biolabs 公司； 胶回收试剂盒，QIAGEN 公司； 建库试剂盒、NovaSeq6000 测序平台，Illumina 公司。

1.3 方法

1.3.1 真菌鉴定 提取全部样品DNA，检测DNA的纯度和浓度。 将DNA 样本于离心管中稀释至1 ng/μL。 以稀释后的基因组DNA 为模板，使用带Barcode 的真菌ITS1 区标准特异引物ITS5-1737F（5′ -GGAAGTAAAAGTCGTAACAAGG -3′ ） 及ITS2-2043R （5′-GCTGCGTTCTTCATCGATGC-3′），进行PCR 扩增。 PCR 产物使用2%琼脂糖凝胶进行电泳检测，合格的PCR 产物进行磁珠纯化，采用PCR-ELISA 定量，根据PCR 产物浓度进行等量混样，混匀后使用2%琼脂糖凝胶电泳检测PCR 产物，对目的条带使用胶回收试剂盒回收产物。使用建库试剂盒进行文库构建。构建好的文库经Qubit 和Q-PCR 定量，文库合格后测序。

1.3.2 数据分析 从下机数据中拆分各样本数据并去除Barcode 和引物序列，然后使用FLASH（v1.2.7）[21]对各样本的reads 拼接，得到的拼接序列为Raw Tags，使用Qiime 软件（v1.9.1）将Raw Tags中连续低质量（质量阈值≤19）碱基数达到3 的低质量部分截去，并进一步过滤掉其中连续高质量碱基长度小于Tags 长度75%的Tags。过滤掉含低质量碱基较多的序列[22]，得到高质量的Clean Tags。 Clean Tags 序列使用VSEARCH（v2.17.1）去除嵌合体序列[23]，得到最终的有效数据（Effective Tags）。

使用R 软件（v 2.15.3）制作稀释曲线，并确定抽平数，使用Uparse 算法[24]（Uparse v7.0.1001）对抽平（Subsample）后的Effective Tags 进行97%一致性聚类聚类，共聚为566 个OTUs，将OTUs 序列数据上传至国家微生物科学数据中心（https：//www.nmdc.cn/），登录号为NMDCX0000111。

筛选出现频数最高的序列作为OTUs 的代表序列，使用Qiime 软件（v 1.9.1）比对Unit（v8.2）[25]数据库进行物种注释，使用MUSCLE（v3.8.31）[26]比对序列，得到OTUs 序列的系统发生关系。 最后以样本数量最少的为标准进行均一化处理，以便后续分析。

使用Qiime 软件 （v1.9.1） 计算Chao1、Shannon、Simpson、PD_whole_tree、Goods_coverage 等 指数，计算Unifrac 距离、构建UPGMA 聚类树，使用SPSS 软件（v22）进行组间差异分析，使用R 软件的pheatmap 包（v1.0.12）绘制热图。

2 结果与分析

2.1 茯砖茶与白茶砖发花情况

两种类型的茶砖发花情况见图1。 茯砖茶D9已有零星“金花菌”出现，D12 出现大量“金花菌”，D22 的“金花菌”茂盛。 白茶砖中W12 有零星“金花菌”，W22 的“金花菌”较为茂盛。

图1 两种类型茶砖发花过程Fig.1 The flowering process of two types of tea bricks

2.2 扩增子测序结果

两类茶砖发花过程样本测序结果见表1。 经质控所得Effective Tags 在Raw PE 中所占比例在59.43%～67.24%之间。 各样本Effective Tags 经97%聚类注释后，获得OTUs 数量在109～247 之间，单个样本的ITS1 序列数均超过60 712 条，因每个样品获得的ITS1 序列数不同，故对数据进行抽平。 从各样本中随机抽取一定测序量的数据（10，6 922，13 834，20 746，27 658，34 570，41 482条序列），统计所代表的OTUs 的数目，绘制稀疏曲线（图2）。当测序序列数达41 482 条时，各样本曲线趋于平缓，表示测序深度充足且合理，因此按41 482 对数据进行抽平。

图2 各样本物种稀疏曲线Fig.2 Rarefaction Curve of samples

表1 所有样本测序及质控结果Table 1 Sequencing and quality control results of all samples

各样品非加权UPGMA 聚类分析见图3。茯砖茶D12 与D22 聚为一类，白茶砖W12 与W22 聚为一类，表明在发花最后10 d 两类茶砖真菌菌群多样性存在差异。 由非加权Unifrac 距离及物种多样性堆积图可知，各样本在不同发花时期真菌多样性差异不大。

图3 基于非加权Unifrac 距离的UMGMA 聚类树Fig.3 Unweighted Unifrac distance-based UMGMA cluster tree

再对各样品进行加权UPGMA 聚类分析，同时分析在门水平上的物种相对多样性及丰度。 如图4 所示，D1 与W1 聚为一类，D6 与W6 聚为一类，D22 与W22 聚为一类，表明发花前期及后期真菌群落丰度在茯砖茶及白茶砖间无显著差异，而发花中期（第9 天，第12 天）存在差异。 由图4堆积图所示，伴随发花过程，子囊菌门相对丰度持续升高，至22 d 发花结束，子囊菌门丰度占比在两类茶中均超过90%，表明发花过程与子囊菌门真菌的相对丰度上升有关。

图4 基于加权Unifrac 距离的UPGMA 聚类树Fig.4 Weighted Unifrac distance-based UPGMA cluster tree

（续表1）

2.3 真菌群落多样性和丰度指数

对两类茶砖在发花过程各样品OTUs 进行Shannon、Simpson、Chao1、PD_whole_tree 等指数的统计分析，结果显示：两茶类Shannon、Simpson 指数均下降，茯砖茶中Chao1 指数和PD_whole_tree下降，而白茶砖中无显著变化，Goods_coverage 显示测序覆盖度均在99.9%以上（表2）。

表2 两茶类发花过程中的真菌群落多样性和丰度指数Table 2 Fungi community diversity and abundance index in the flowering process of the two tea types

2.4 总体真菌群落多样性和丰度

从微生物组扩增子（ITS）中共检测566 个OTUs，分属5 门、19 纲、45 目、99 科、160 属、214 个菌株，其中有131 个菌株鉴定至种水平。门水平上丰度依次是子囊菌门（Ascomycota）（60.4%）、担子菌门 （Basidiomycota）（35.7%）、 球囊菌门（Glomeromycota）（1.3%），壶菌门（Chytridiomycota）及毛霉门（Mucoromycota）丰度低于1%。 以可检测到的最小真菌分类水平为基准，按丰度由高到低将总丰度前100 的真菌排序见表3。Top10 的真菌占总丰度的93.77%，Top50 的真菌占总丰度的97.20%，Top100 的真菌占总丰度的97.35%，其它真菌总丰度小于0.001%，未检出分类的OTUs 占比2.6144%，可见丰度Top10、Top50、Top100 均超过90%，可以较好地代表总体样本。

表3 丰度前100 的真菌鉴定最小分类水平排序Table 3 The top 100 fungi in abundance were identified by sorting at the minimum taxonomic level

（续表3）

（续表3）

2.5 发花过程真菌丰度动态变化

茯砖茶与白茶砖在发花过程中真菌丰度变化如图5 所示。 热图将Top50 的真菌聚为3 类（Ⅰ、Ⅱ、Ⅲ）。Ⅰ类真菌在两类茶中起始丰度较低，在发花各阶段分别有丰度峰，白茶砖中散囊菌属、假灰绿曲霉、海洋嗜杀酵母、卡利比克迈耶氏酵母等的丰度峰出现在6～9 d，而在茯砖茶中丰度峰出现在9～12 d，22 d 的散囊菌属丰度最高，最终成为优势菌种。 Ⅱ类真菌在白茶砖中的丰度明显高于茯砖茶，且在发花初期及中期丰度较高，茯砖茶的第2个丰度峰在9 d，而白茶的第2 个峰在12 d，在22 d 丰度降至最低。 Ⅲ类真菌在茯砖茶中的丰度显然高于白茶，且茯砖茶丰度峰集中在6 d，白茶砖中没有明显丰度峰。 交叉散尾鬼笔、伞菌目、粉侧担菌在白茶中1～6 d 比茯砖茶的丰度稍低，黑曲霉在茯砖茶中前1～12 d 维持较高丰度，而在白茶砖中丰度维持较低水平。

图5 Top50 真菌的丰度动态热图Fig.5 Heatmap of the abundance of Top50 fungi

2.6 Top10 的真菌动态变化

丰度Top10 的真菌占总体样本丰度的93.77%，分别以堆积图（图6）及热图（图7）表示丰度动态变化。 图6 所示，在发花初期，散囊菌属与假灰绿曲霉在发花初期丰度较低（3%～6%），D9 与W6 中散囊菌属大量繁殖，而D12 与W9 中散囊菌属丰度下降，假灰绿曲霉大量繁殖，22 d 散囊菌属丰度再度上升，在整个发花过程中，此两真菌交替成为优势菌种。 粉侧担菌丰度较高（59.2%，49.2%），在茯砖茶中不断下降，在白茶砖中12 d 有所回升，22 d 丰度降至较低水平（2%左右）；而同为担子菌的交叉散尾鬼笔丰度变化与粉侧担菌保持同步，且丰度平均为粉侧担菌的35%（标准偏差0.12）。黑曲霉在茯砖茶中D1～D12 相对丰度由5.2%缓慢降至3.8%，在D22 骤降至0.08%，而在白茶中由W1 的0.15%缓慢降至W22 的0.07%，是Top10真菌中两类茶差别较大的菌种。 图7 所示热图将Top10 真菌聚为两组，一组为散囊菌属、假灰绿曲霉、黑曲霉3 种菌丰度变化各不相同；另一组7 种真菌在两类茶内丰度变化相似，在茯砖茶中D1 至D9 丰度较高，且呈下降趋势，D12 至D22 丰度骤降至较低水平，白茶砖中W1 至W9 丰度下降，W12 丰度有所上升，W22 骤降至较低水平。

图6 Top10 真菌的相对丰度动态堆积图Fig.6 Dynamic accumulation relative abundance dynamics of Top10 fungi

图7 Top10 真菌的相对丰度动态热图Fig.7 Dynamic heatmap of relative abundance of Top10 fungi

3 结论与讨论

本研究中，两类茶在22d 均发花成功，且测序结果表明不同茶类的发花真菌多样性没有显著差异。有研究表明不同地区原料、不同环境的发花微生物有所不同[27-28]。 本试验控制原料相同，发花环境一致，证明毛茶加工工艺不同并不导致发花真菌多样性的显著差异，而茯砖茶与白茶砖 “菌花香” 的不同可能是由微生物菌群丰度动态变化不同所致。

由于散囊菌属ITS 一区不能鉴定到种水平，因此将“金花菌”鉴定到属水平（散囊菌属）。 有研究表明“金花菌”应归为散囊菌属[29]，相关研究较为明晰。本研究认为鉴定出的散囊菌属真菌即“金花菌”的1 种或几种。

通过图5 将Top50 的真菌聚为3 类，同时对比茯砖茶与白茶砖菌群演替的差异，获得茯砖茶发花特征真菌菌群演替规律，即： 茯砖茶发花前期，第Ⅲ类真菌中的黑曲霉、局限曲霉、萨氏曲霉、岐皱青霉、西兰花青霉、汉逊德巴利酵母在D6 出现丰度峰；之后D9 散囊菌属、卡利比克迈耶氏酵母、海洋嗜杀酵母出现丰度峰，黑曲霉、汉逊德巴利酵母、局限曲霉维持丰度；D12 假灰绿曲霉出现丰度峰，黑曲霉继续维持丰度，D22 散囊菌属最终成为优势菌种。由于茯砖茶的原料是黑毛茶，其独特的工艺“渥堆”过程中黑曲霉、酵母菌、青霉、根霉是主要的微生物[30]，研究表明曲霉属在渥堆过程中能产生单宁酶、糖苷酶、漆酶等胞外酶[31]，与本研究茯砖茶中Ⅲ类真菌D1 时期丰度高于白茶砖结果一致，表明茯砖茶发花特征真菌有一部分是渥堆大量繁殖的真菌，并在发花前期进行代谢与繁殖维持一定的丰度。 图5 的第Ⅱ类真菌在茯砖茶中的丰度明显低于白茶砖，可能是由于茯砖茶经过渥堆，某些微生物生长改变了微生物菌群结构[32]，消耗了部分营养物质；渥堆产生的较高温度使得不适的微生物丰度下降，Ⅱ类真菌丰度因此下降，并在发花过程中丰度持续下降，在D9 有小幅度回升，至D12 丰度降至极低水平，Ⅱ类真菌在先前发花真菌研究中鲜有报道，而在茶叶内生真菌研究中有所报道[33]。 对比白茶砖中的丰度可知，Ⅱ类真菌可能在发花过程中是负贡献微生物，该结论指导白茶砖发花生产前处理应比茯砖茶增加灭菌步骤。 Ⅰ类真菌是茯砖茶与白茶砖共有的发花真菌，大多在先前的研究中有报道，如散囊菌属是“金花菌”的集合，假灰绿曲霉是匍匐散囊菌的无性型，赤曲霉则是赤散囊菌的无性型[29]，这3种金花菌在发花过程中交替成为优势菌种，见图5。 首先是散囊菌属（D9、W6），假灰绿曲霉与赤曲霉丰度峰同时出现（D12、W9～W12），最终散囊菌属成为优势菌种（D22、W22）。 有研究表明假灰绿曲霉与赤曲霉聚在一个进化枝[34]，表明二者遗传距离较近，有较相似的演替规律。发花过程中酵母菌的研究鲜有报道。 本研究中Ⅰ类真菌中卡氏酵母与海洋嗜杀酵母同时出现丰度峰 （D9、W6），热带假丝酵母与隐球酵母属在茯砖茶发花过程的前10 d 为优势菌[6]，而本研究中茯砖茶中的丰度较低，在白茶砖中、后期（W9～W22）的丰度较高，表明它们也可能不是茯砖茶发花的优势真菌。 通过对比Ⅰ类真菌在发花过程中丰度峰出现的时间，发现白茶砖中多数真菌丰度峰比茯砖茶中的丰度峰提前3 d 左右出现，这可能是茯砖茶发花过程中菌群演替规律，可能是茯砖茶渥堆改变了Ⅲ类菌群的丰度，同时消耗了茶砖的营养物质，这仍需进一步研究。

研究茯砖茶发花过程中特征真菌菌落演替是解析茶叶发花过程至关重要的一步，微生物组学的应用，能够将大多数微生物鉴定到种一级，能够深刻认知发花过程，对于后续分析各种微生物在发花过程中的作用奠定基础，对揭开真菌与茯砖茶品质风味形成的关系有重要意义。