基于宏基因测序及分箱技术的复合乳酸菌发酵乳菌种鉴定及定量

程坤，赵婷，刘蕊，周立光，冯会粉，于学健，李宏，姚粟*

1(中国食品发酵工业研究院有限公司，中国工业微生物菌种保藏管理中心，北京，100015) 2(中国合格评定国家认可中心，北京，100062)

发酵乳是以生牛(羊)乳或乳粉为原料，经杀菌、发酵后制成的pH值降低的产品(GB 19302—2010《食品安全国家标准发酵乳》)。发酵乳行业的规模和质量飞速发展，符合“健康中国”、“国民营养计划”等国家战略，是关系国计民生的重要产业。嗜热链球菌(Streptococcusthermophilus)和德氏乳杆菌保加利亚亚种(Lactobacillusdelbrueckiisubsp.bulgaricus)是发酵乳常用发酵剂菌种，二者之间存在共生关系，德氏乳杆菌保加利亚亚种代谢生成的氨基酸可刺激嗜热链球菌的生长，而嗜热链球菌在代谢过程中产生的甲酸可以促进德氏乳杆菌保加利亚亚种的生长[1]。发酵乳通常添加对人体有益的益生菌，如：动物双歧杆菌乳亚种(Bifidobacteriumanimalissubsp.lactis)和干酪乳酪杆菌(Lacticaseibacilluscasei)，有助于促进肠道微生物的定殖，治疗和预防腹泻、缓解便秘、增强人体免疫力等[2]，某一株益生菌在产品中需要达到一定数量才能起到保健作用，因此发酵乳中菌种的鉴定和定量尤其重要。

目前我国对于发酵乳中乳酸菌的使用和监管，标准主要为GB 4789.34—2016《食品安全国家标准 食品微生物学检验 双歧杆菌检验》和GB 4789.35—2016《食品安全国家标准 食品微生物学检验 乳酸菌检验》，采用传统培养方法和生理生化鉴定方法进行发酵乳中乳酸菌的鉴定和定量，具有成本低，易于操作等优点，但需要获得菌种纯培养物，且存在检测结果不准确等问题。建立快速、准确的复合乳酸菌发酵乳菌种鉴定和定量检测方法，成为发酵乳行业亟待解决的热点、难点问题之一[3]。近几年来，高通量测序技术迅速发展，采用宏基因组测序分箱分析可将包含不同生物体的序列进行组装，重现复合样品中的微生物种类，解析产品组成。该方法以复合样品总基因组DNA为分析对象，有效避免了菌种分离培养步骤，目前已被应用于复合益生菌产品的物种组成分析[4]。平均核苷酸一致性(average nucleotide identity，ANI)是通过全基因组序列鉴定细菌物种的常用方法之一，其域值95%～96%与DNA-DNA杂交值70%相对应，可作为细菌种水平界定的标准[5]。ANI分析具有方便、快捷、分辨率高的优点[5]，其与宏基因组分箱分析相结合，为复合乳酸菌发酵乳菌种的种水平精准鉴定提供了一种可靠方法。

复合乳酸菌发酵乳产品的菌种组成与定量分析是产品质量的关键。本文以3款市售复合乳酸菌发酵乳产品为研究对象，利用宏基因组测序和分箱技术、ANI分析等技术，研究复合乳酸菌发酵乳产品中菌种种水平鉴定及相对定量。对复合菌发酵乳产品的质量控制和行业监管具有重要参考意义。

1 材料与方法

1.1 实验材料

1.1.1 复合乳酸菌发酵乳产品

本研究选取3款不同品牌的复合乳酸菌发酵乳产品A、B、C，均购自北京大型超市。产品标签标识菌株类别及活菌数如下：

样品1：产品A(菌种名称：动物双歧杆菌BB-12、嗜热链球菌、德氏乳杆菌保加利亚亚种；保加利亚乳杆菌和嗜热链球菌总数1×106CFU/g，双歧杆菌1×107CFU/g)；

样品2：产品B(菌种名称：保加利亚乳杆菌、嗜热链球菌、干酪乳杆菌，乳酸菌活菌数≥1×106CFU/g)；

样品3：产品C(菌种名称：乳双歧杆菌BL-99、嗜热链球菌、保加利亚乳杆菌；乳双歧杆菌总数≥1×106CFU/g)。

1.1.2 试剂与仪器

DNA提取试剂盒：E.Z.N.A.®Bacterial DNA Kit，美国OMEGA公司；文库构建试剂盒NEXTflex Rapid DNA-Seq Kit，美国Bioo Scientific公司。

Eppendorf微量可调节移液器、5424R小型台式冷冻离心机，德国Eppendorf公司；HG-50高压灭菌锅，日本HIRAYAMA公司；AC2-4S1生物安全柜，新加坡Esco公司；BioDrop核酸蛋白分析仪，豪沃生物科技(上海)有限公司；Qubit 3.0 核酸荧光定量仪，美国ThermoFisher Scientific公司；Illumina HiSeq 3000基因测序平台，美国Illumina公司。

1.2 实验方法

1.2.1 基因组DNA的提取

取20 mL发酵乳样品，4 000 r/min离心5 min，吸取上清液，12 000 r/min离心10 min，取沉淀按照试剂盒E.Z.N.A.®Bacterial DNA Kit的说明书提取基因组DNA。使用Qubit 3.0核酸荧光定量仪、BioDrop核酸蛋白分析仪分别检测基因组DNA的浓度和纯度，使用1%琼脂糖凝胶电泳检测基因组DNA的完整度。

1.2.2 宏基因组文库构建与测序

通过物理方法将质检合格的基因组DNA随机打断成300～400 bp的片段，使用NEXTflex Rapid DNA-SeqKit按照说明书进行文库制备。使用Illumina HiSeq 3000测序平台(上海美吉生物医药科技有限公司)进行宏基因组测序，宏基因组序列NCBI登录号分别为：SAMN22013808，SAMN22013809，SAMN22013810。

1.2.3 下机数据质控与组装

采用MetaWRAP流程[6]对宏基因组下机数据进行质控和组装，采用Trim Galore(http://www.bioinformat- ics.babraham.ac.uk/projects/trim_galore/)去除reads 的adapter、过滤质量值低于20的reads，参数采用默认设置。使用Megahit[7]软件进行序列组装，采用succinct de Bruijn graph方法进行拼接，k-mer参数从小至大迭代，选择长度≥1 000 bp的contigs用于后续分析。

1.2.4 分箱分析

采用MetaWRAP流程对拼接后的contig开展分箱分析，获得各样品的分箱结果(bin)。利用checkM[8]评估各bin的完整度和污染度，选取完整度>70%和污染度<10%的bin[9]，利用blobology模块绘制各样品中bin的丰度散点图。通过Circos[10](版本0.69-8)软件对各样品bins 进行可视化分析，每个bin的contigs重新基于模式菌株有参组装，展示最终的基因组圈图。

1.2.5 ANI分析

从NCBI genome 数据库下载相关种模式菌株的基因组序列，使用fastANI[11](版本：1.32)计算各bin与相关物种模式株基因组序列之间的ANI值。以ANI值为95%～96%作为细菌物种界定的标准[5]，对各样品的分箱结果进行物种鉴定。

1.2.6 物种定量分析

采用Quant_bins模块分析各bin的相对丰度。

2 结果与分析

2.1 测序数据质控和组装

针对样品的宏基因组测序下机数据，采用Fastp[12]软件，碱基质量值≥20，对raw data进行质控。质控后A、B、C 3个产品的clean data数分别为13 667 335 968、12 529 436 956和12 225 941 409 bp，碱基质量值Q30均大于90%，测序数据质量能够满足后续分析(表1)。

表1 样品宏基因组测序数据统计Table 1 Data statistics for metagenomic sequencing data

对clean data进行组装，各产品获得的contigs数目及相关参数见表2。

表2 样品宏基因组数据组装结果Table 2 Assembly result of the metagenomic sequencing data

2.2 复合乳酸菌产品中物种基因组重建

使用MetaWRAP宏基因组分析流程对组装后的contigs进行分箱分析，获得分箱结果(bin)，即重建的基因组。产品A、B、C的contigs 均分别分箱成3个bin(表3)，与各产品标识所声称添加的菌种数量一致。通过checkM来评估分箱后bin的完整度和污染度，结果表明，各样品分箱产生的bin的完整度均高于85%，污染度均低于5%(表3)。各Bin丰度散点图见图1。

表3 宏基因组分箱后各bin的参数Table 3 The results of metagenomics binning

a-产品A；b-产品B；c-产品C图1 Bin丰度散点图Fig.1 Scatter plots of bins abundance注：一个点代表一个contig，相同颜色的contig来自同一个bin

2.3 复合乳酸菌产品中物种的鉴定

将分箱获得的各bin通过NCBI-Blast进行初步比对，确定其近缘物种，利用fastANI软件计算各bin与近缘物种模式株之间的ANI值(表4)。每个样品中的bin与最近缘物种模式菌株基因组ANI分析的值均>95%，表明每个bin 代表的菌种与模式菌株为同一个种(表4)。

表4 分箱结果与参考模式菌基因组ANI 分析Table 4 The analysis results of ANI between bins and genomes of type strains

产品A和产品C中3个bin通过ANI分析鉴定为嗜热链球菌、动物双歧杆菌和德氏乳杆菌，与该产品声称添加的菌种一致。产品B中3个bin通过ANI分析鉴定为嗜热链球菌、干酪乳酪杆菌和乳酸乳球菌(Lactococcuslactis)，其中，嗜热链球菌和干酪乳酪杆菌与该产品声称添加的菌种一致，未检测到该产品声称添加的德氏乳杆菌，检测到的乳酸乳球菌未出现在该产品声称中。

Contig利用模式基因组进行有参组装后，通过Circos v 0.69-8软件对各样品bin进行可视化，基因组圈图见图2。

图2 产品宏基因组菌种基因组重建Fig.2 Re-constructed genomes from samples by means of the metagenomics binning注：每一个产品的bins 都由横向排列的遗传图谱来表示，相同物种用同一种颜色显示，绿色“√”表示分箱得到的物种与产品声称菌种一致

动物双歧杆菌和干酪乳杆菌在我国《可用于食品的菌种名单(卫办监督发〔2010〕65号)》之列，是

发酵乳常用的益生菌，在A、B和C 3个产品中也均有添加。其中，动物双歧杆菌乳亚种曾用名为乳双歧杆菌(Bifidobacteriumlactis)[13]；德氏乳杆菌保加利亚亚种的曾用名为保加利亚乳杆菌(Lactobacillusbulgaricus)[14]。干酪乳杆菌(Lactobacillucasei)在2020年其分类学地位变更为干酪乳酪杆菌[15]。

德氏乳杆菌在产品B中未被检测到，进一步分析发现，在利用完整度>70%、污染度<10%参数对产品B宏基因组分箱初步结果的过滤前，检测到1个完整度仅为6.61%的bin，通过ANI分析显示为德氏乳杆菌，定量结果显示相对丰度为0.09%，由于bin完整度较低，可能会引起该定量检测结果的偏差，该菌种按照宏基因组数据过滤参数未形成完整bin，说明产品中菌种的相对丰度较低时，可能引起分箱bin的完整度不够而没有被检测到。产品B中检测到的乳酸乳球菌，推测可能为原料、环境、生产过程等环节引入，或未经菌种鉴定误按嗜热链球菌添加进产品，尚需进一步的验证确认。

2.4 复合乳酸菌产品中各物种的定量分析

利用Quant_bins模块计算各样品中bin的相对丰度，结果如表5所示。

表5 各分箱结果的定量分析Table 5 Quantitative analysis of bins

对于产品A，其分箱分析所检测到的嗜热链球菌、德氏乳杆菌和动物双歧杆菌3个物种的相对丰度分别为96.85%、0.12%和3.02%，嗜热链球菌为绝对优势物种。该产品标识声称嗜热链球菌和德氏乳杆菌总量为1×106CFU/g，定量结果表明，产品中二者的相对比例约为807∶1。产品B标识声称活菌总数为1×106CFU/g，通过定量分析显示，分箱分析所检测到的嗜热链球菌、干酪乳酪杆菌和乳酸乳球菌3个物种的相对丰度分别为99.43%、0.27%和0.21%，嗜热链球菌为绝对优势物种。产品C声称添加的物种类别与产品A相同，且其中动物双歧杆菌的添加量为1×106CFU/g。定量结果表明，嗜热链球菌的相对丰度为99.55%，德氏乳杆菌和动物双歧杆菌的相对丰度分别为0.18%和0.14%。复合乳酸菌发酵乳产品中对所添加各菌种定量主要是以乳酸菌活菌数或某一种/株益生菌活菌数为指标，本研究中基于宏基因组的定量分析为相对定量，以活菌数和死菌数总和的百分比表示，随着后生元(postbiotics)[16]概念逐渐被大众所熟悉及相应产品的开发，相较于传统定量检测方法，本文提供了新的思路。

3 结论

本研究以3款市售复合乳酸菌发酵乳产品为研究对象，利用宏基因组测序和分箱分析、ANI分析等技术，研究复合乳酸菌发酵乳产品中菌种种水平鉴定及相对定量。产品A和产品C分箱获得的3个bin通过ANI分析鉴定为嗜热链球菌、德氏乳杆菌和动物双歧杆菌，与产品声称添加菌种一致，产品A 3个物种的相对丰度分别为96.85%、0.12%和3.02%，产品C 3个物种的相对丰度分别为99.55%、0.18%和0.14%。产品B分箱获得的3个bin通过ANI分析鉴定为嗜热链球菌、干酪乳酪杆菌和乳酸乳球菌，与产品声称添加菌种部分一致，3个物种的相对丰度分别为99.43%、0.27%和0.21%。本研究为复合乳酸菌发酵乳产品的菌种鉴定与定量提供了一种可行的分析方法，与传统基于可培养的分析方法相比，具有准确、可靠等优点，且不受菌种死活状态的影响。随着测序技术的不断革新，高通量测序在成本与周期方面必将进一步下降和缩短，宏基因测序及分箱作为一项新兴技术在食品检测领域将进一步推广应用，有望发展成为一种常规分析方法。本研究对复合乳酸菌发酵乳产品的质量控制和行业监管具有重要参考意义。