王宸之，杜沁岭，周思懿，黄　静，兰秋雨，秦　文，杨文钰，刘　江，张　清,＊

（1.四川农业大学食品学院，四川　雅安　625014；2.四川省作物带状复合种植工程技术研究中心，成都　611130）

豆腐是优质植物蛋白来源[1]，经浸泡、研磨、煮沸、过滤、点浆等步骤制成传统填充豆腐[2]。研究表明，在豆浆过滤环节，47%大豆原料过滤后留在豆渣中，豆渣具有较高水分活度，极易腐败变质，过滤豆渣被制成动物饲料或遗弃，浪费原料[3]。为解决传统豆腐原料利用率低问题，全豆填充豆腐研制成为热点。全豆填充豆腐（Packed whole-soybean curd）无豆浆过滤步骤，将大豆原料全部或大部分保留蛋白凝胶。学者针对全豆豆腐凝胶性能，制备工艺、加工条件等开展研究[4-6]。

全豆填充豆腐提高原料利用率，但存在食品安全隐患，大豆生长和种植过程中，表面携带大量微生物，煮沸可杀灭大部分微生物，但仍有部分耐热芽孢残存在豆腐成品中。此外，全豆豆腐含有大量豆渣，初始菌数高于传统去渣豆腐，影响货架期。预测微生物学通过数学模型描述特定环境因素对食物中微生物生长影响[7]，反映微生物生长及衰退趋势，预测产品货架期，在果蔬[8]、海产品[9]生产等领域应用广泛。目前预测微生物学研究主要集中于食品腐败菌及食源性致病菌，Juneja等采用Baranyi模型、Huang模型、modified Gomp⁃ertz模型及三参数线性一级模型，建模预测煮熟大豆中蜡样芽孢杆菌生长动力学特性[10]，利用准确因子（Af）、偏差因子（Bf）、相关系数（R2）、均方根误差（RMSE）评估模型精密度。结果表明，Baranyi模型可精确描述蜡样芽孢杆菌热休克孢子状态萌发、生长、繁殖动力学曲线，Modified Ratkowsky二级模型拟合结果发现贮藏温度与微生物生长速率呈较高相关性（R2=0.98），因此Baranyi一级模型与Modified Ratkowsky二级模型可以用于豆类食品中孢子萌发、繁殖动力学特性研究。预测微生物学模型可用于生长曲线拟合，与理化指标结合分析食品腐败变质原因并预测食品货架期。Matara⁃gas等分析水果酸奶中菌群及理化指标发现，乳酸菌是主要菌群，此外，在限定贮藏温度条件下（5～20℃），酵母数量随时间推移不断上升，但在酵母生长初始阶段出现延滞生长[11]。理化检测结果显示，在冷藏温度下（5℃，10℃），酸奶中可滴定酸和pH均保持恒定，当温度上升时，可滴定酸含量上升，pH下降。预测微生物学模型验证试验表明，17℃温度下，酸奶货架期预测需结合数学模型、感官试验及理化检验。针对豆腐中腐败菌预测微生物学研究较少，Lee等分离和鉴定韩国商业化豆腐中腐败菌，结果表明引发豆腐腐败变质主要微生物是蜡样芽孢杆菌[12]，与Juneja等研究结果一致。在豆腐制作过程中，大豆原料微生物污染状况对豆腐品质有重要影响，后续预测微生物学模型验证试验表明，温度是影响豆腐货架期关键因素。但Rossi等发现豆腐中特异性腐败菌是乳酸菌，即使同一种类食品，腐败菌种类及生长特性不同[13]。预测微生物学现有研究内容包含：①模型在食品中应用、验证、评估。②不同模型精确度对比及常见食物中适用性研究。③预测微生物学模型应结合理化、感官检验技术对食品货架期综合评判。

目前，全豆豆腐腐败菌和货架期研究尚未见报道。因此，本文以全豆填充豆腐为研究对象，分离和鉴定全豆豆腐中腐败菌。在预测微生物学技术基础上，建立腐败菌生长一级模型和二级模型，绘制腐败菌生长动力学曲线，预测和验证不同贮藏温度下全豆豆腐货架期。

1　材料与方法

1.1　材料试剂与仪器

大豆（南豆12，蛋白含量为46%，水分含量11%，粗脂肪含量14%），由四川省南充市农业科学院提供；葡萄糖内酯（Glucono-δ-lactone，GDL），食品级，购自安徽省兴省医药食品有限公司；胰酪大豆胨琼脂培养基（TSA），分析纯，购自杭州微生物试剂有限公司；细菌DNA提取试剂盒，DP302，购自天根生化科技（北京）有限公司；高纯度低电渗琼脂糖，Tsingke；DNA凝胶回收试剂盒，Tsingke； DNA Marker， 型 号 DL2000， TaKaRa。测序仪，3730XL，Applied Biosystems；离心机，Micro17，Thermo Legend；PCR 仪，2720 thermal cycler，Applied Biosystems；凝胶成像仪，JY04S-3C，北京君意东方电泳设备有限公司；高速粉碎机，XM-800Y，浙江省永康市铂欧五金制品厂。

1.2　全豆填充豆腐制备

大豆原料经筛选、清洗后于60℃烘干备用。烘干后大豆以高速中药粉碎机粉碎1 min（转速12 000 r·min-1）后制成豆粉，豆粉与饮用水按照1∶7.5质量体积比混合均匀后制成悬浮液，悬浮液升至95℃保温5 min后立即冰浴降温至55℃制成熟豆浆，量取80 mL熟豆浆于100 mL小烧杯内，加入豆浆体积0.4%（m/v）GDL，GDL与豆浆混合均匀后于55℃水浴保温1 h形成全豆填充豆腐。

1.3　全豆填充豆腐中腐败菌分离、纯化与鉴定

豆腐腐败菌分离与鉴定采用Lee等方法[12]，将制备全豆填充豆腐置于37℃条件下放置24 h，取5 g豆腐样品，以无菌生理盐水梯度稀释，选取2～3个适宜稀释度稀释液1 mL倾注TSA培养基，37℃条件下培养24 h，挑取不同形态菌落TSA平板划线纯化形成单菌落。

腐败菌鉴定采用16S rDNA测序法，腐败菌基因组DNA提取参照天根DP302（细菌）试剂盒方法，提取后DNA PCR扩增，PCR长度：1 500 bp；引物设置：27F（AGTTTGATCMTGGCTCAG），1492R（GGTTACCTTGTTACGACTT）；PCR反应体系：T5 Mix，25 μL；PF（10P），1 μL；PR（10P），1 μL；gDNA，1μL；ddH2O，22μL；PCR扩增条件：94℃，5 min；98℃，10 s；55℃，15 s；72℃，15 s；72℃，5 min；4℃，保持。其中，98℃，10 s至72℃，15 s间共循环30次。扩增后产物琼脂糖凝胶电泳：1%琼脂糖，150 V、100 mA，20 min。切割所需DNA目的条带，纯化PCR产物，用引物直接测序，FASTS序列在Reverse complement系统中反向序列互补，互补序列在EZbiocloud系统（https://www.ezbiocloud.net）中序列对比分析，确定腐败菌种类。

1.4　全豆豆腐腐败菌生长一级模型构建

微生物生长一级模型反映恒定环境条件下微生物数量与时间关系，即微生物生长曲线。常见一级模型有Logistic模型、Gompertz模型以及Modified Gompertz模型等。本试验将全豆豆腐分别置于15、20、25和30℃条件下，每间隔一段时间取样，将5 g样品梯度稀释后倾注TSA平板，37℃条件下培养24～48 h，测定豆腐中菌落数，建立微生物数量与时间对应关系。选用Logistic（式1）和Modified Gompertz（式2）两种模型拟合腐败菌生长情况，比较两模型适用性。

式中，Nt—贮藏时间为t h时菌数，log（cfu·g-1）；μmax—微生物生长最大比生长速率（h-1）；λ—微生物生长延滞时间（h）。A—相对最大菌浓度,即logNmax/N0，log（cfu·g-1）；N0—样品初始菌数，log（cfu·g-1）；Nmax—样品最大菌数，log（cfu·g-1）；

式（2）中，Nt—贮藏时间为t h时菌数，log（cfu·g-1）；N0—样品初始菌数，log（cfu·g-1）；C—微生物最大菌数（Nmax）与初始菌数（N0）差值，即C=Nmax-N0，log（cfu·g-1）；μmax—微生物生长最大比生长速率（h-1）；λ—微生物生长延滞时间（h）。

1.5　全豆豆腐腐败菌生长二级模型构建

微生物生长一级模型反映恒定环境条件下微生物生长情况，二级模型可反映环境因素变化对微生物生长特性影响，即环境因子影响模型。由于温度对豆腐中微生物生长有显著影响，选用延滞期模型（Lag Time,LT模型，式3）和平方根模型（式4）拟合一级模型所得数据。

式（3）中，LT—延滞期（h）；T—贮藏温度（℃）；aLT，bLT回归常数。

式（4）中，μmax—微生物生长最大比生长速率（h-1）；T—贮藏温度（℃）；Tmin—微生物生长最低温度（℃）；b，回归常数。

1.6　全豆豆腐货架期预测模型构建

全豆豆腐货架期预测基于一级模型，当样品菌数达6 log（cfu·g-1）时，认为豆腐到达货架期终点，采用迭代法求出对应时间，即为全豆豆腐货架期。

1.7　统计与分析

采用Origin 9.0软件对数据非线性回归拟合及误差分析，采用Microsoft Excel 2014对数据单变量求解。

2　结果与分析

2.1　全豆豆腐中腐败菌分离与鉴定结果

将制备全豆豆腐于37℃条件下贮藏24 h后，取样5 g，以无菌生理盐水梯度后倾注TSA培养基，37℃恒温培养24 h，因样品存在芽孢杆菌，若培养24 h后仍无菌落生成，则延长培养时间至48 h。经观察，TSA平板中共生成2种形态不同菌落，以无菌接种环挑取菌落划线于TSA平板上，纯化后获得单菌落T-1及T-2。全豆豆腐中腐败菌在TSA平板上纯化单菌落形态如图1所示。菌落形态学分析如下：T-1经接种环挑起后有粘性，单菌落呈圆形，无光泽，较干燥，边缘不整齐且有明显褶皱。T-2呈灰白色、不透明，表面粗糙，菌落有明显光泽，呈毛玻璃或熔融蜡状，边缘呈扩展状，经接种环挑起后有粘性。菌株分析参考文献[14]，结果表明T-1及T-2含芽孢杆菌，其中T-2为蜡样芽孢杆菌。

为更精确分析腐败菌类别，本试验采用16S rDNA测序法分析T-1及T-2，首先采用DP302（细菌）试剂盒提取细菌基因组，特异性引物扩增16S rDNA序列后，采用琼脂糖凝胶电泳对扩增产物纯化鉴定。琼脂糖凝胶电泳图谱显示（见图2），T-1及T-2扩增后DNA序列约1 500 bp，与16S rDNA序列长度一致，扩增成功。

图1　全豆豆腐中腐败菌在TSA平板上菌落形态Fig.1　Colony morphology of spoilage bacterium isolated from whole-soybean curd on TSA medium

图2　全豆豆腐腐败菌16S rDNA琼脂糖凝胶电泳结果Fig.2　Agarose gel electrophoresis of spoilage bacterium 16S rDNA isolated from whole-soybean curd

切割1 500 bpDNA目的条带，将纯化后PCR产物用引物直接测序，测序结果与EZbiocloud数据库对比。如图3所示，菌株T-1 16S rDNA反向互补序列长度1 381 bp，覆盖范围为47～1 427 bp，覆盖率93.8%。菌株T-2 16S rDNA反向互补序列长度为1 377 bp，覆盖范围50～1 426 bp，覆盖率93.4%。

EZbiocloud系统将覆盖序列与数据库中标准菌株序列对比，排除无效菌株后根据16S rDNA序列覆盖度及序列相似度排序（只列出排序前3菌株），如表1所示，腐败菌菌株T-1序列相似度由高到低依次为：Bacillus subtilis subsp.Inaquosorum（KCTC 13429（T），序列相似度99.93%）、Bacillus tequilen⁃sis（KCTC 13622（T），序列相似度99.93%）、Bacillus subtilis subsp.Subtilis（NCIB 3610（T），序列相似度2/1381）。腐败菌菌株T-2序列相似度由高到低依次为：Bacillus cereus（ATCC 14579（T），序列相似度100.00%）、Bacillus wiedmannii（FSL W8-0169（T），序列相似度99.93%）、Bacillus albus（N35-10-2（T），序列相似度99.93%）。由上述数据可得：菌株T-1为枯草芽孢杆菌（Bacillus subtilis）或特基拉芽孢杆菌（Bacillus tequilensis），T-2为蜡样芽孢杆菌（Bacillus cereus）。

图3　EZbiocloud系统中腐败菌16S rDNA反向序列信息Fig.3　Reverse sequence details of 16S rDNA isolated from spoilage bacterium in EZbiocloud system

表1　全豆豆腐腐败菌16S rDNA测序鉴定结果Table 1　Taxonomic details of spoilage bacterium in whole-soybean curd as identified on the basis of 16 S rDNA gene partial sequence

大豆原料大部分来源于田间，大豆表面附泥土，携带大量土壤微生物，豆浆加热过程使大豆蛋白发生热变性，杀灭主要细菌营养体细胞，但部分芽孢杆菌在逆环境中能生成耐热芽孢残留，环境适宜芽孢重新萌发，产生毒素，危害食品安全。枯草芽孢杆菌可产生淀粉酶及枯草杆菌素，引起食品腐败变质，蜡样芽孢杆菌可引起食物中毒[15]。

Lee等研究报道，成品豆腐额外热处理利于豆腐保存[12]，全豆豆腐因芽孢对逆环境耐受力极强，普通热处理难以杀灭豆浆中芽孢，因此二次杀菌应集中辐照处理，抑制芽孢数量，延长全豆豆腐货架期。

2.2　全豆豆腐中腐败菌生长动力学模型构建

如图4所示，采用Logistic（A）和Modified Gom⁃pertz（B）模型对15、20、25和30 ℃温度下腐败菌生长趋势非线性回归拟合，其中Modified Gompertz（B）模型所拟合曲线呈典型延滞期、对数期和稳定期，而Logistic（A）模型20、25、30℃温度下呈典型对数期和稳定期，15℃下Logistic（A）模型仅呈对数期，与模型拟合精度有关，拟合精度越高，腐败菌实际生长趋势还原度越高。

图4　不同温度下全豆豆腐中腐败菌生长动力学Fig.4　Kinetic curves of spoilage bacterium on whole-soybean curd stored at different temperatures

表2为拟合曲线精度统计学分析结果，15～30℃温度下，Logistic（A）模型拟合精度（R2）依次为0.93、0.92、0.93、0.97，对应 Modified Gompertz（B）模型拟合精度（R2）依次为 0.99、0.98、0.98、0.99。各温度下Modified Gompertz拟合精度均高于Logistic（A）模型，可见Modified Gompertz（B）模型更适用于全豆豆腐中腐败菌生长趋势拟合。此外，Logistic（A）或Gompertz（B）模型拟合结果均显示，随豆腐贮藏温度由30℃降至15℃，腐败菌最大比生长速率（μmax）降低，延滞期（Lag time）增大，温度与腐败菌生长呈负相关趋势。Logistic（A）模型相对最大菌浓度（A值）恒定在8 log（cfu·g-1），但15 ℃下相对最大菌浓度为（9.18±0.87）log（cfu·g-1），因15℃下曲线拟合未收敛所致。Gompertz（B）模型最大菌数（Nmax）与初始菌数（N0）差值较为恒定，在5 log（cfu·g-1）。全豆豆腐初始菌数较为恒定，在3 log（cfu·g-1），因此在15～30 ℃条件下腐败菌最大生长量（Nmax）较为恒定，15～30℃时，温度仅改变腐败菌达到最大生长量时间，不改变最大生长量。

表2　不同温度下全豆豆腐腐败菌生长曲线拟合参数Table 2　Fitting parameters of spoilage bacterium growth on whole-soybean curd stored at different temperatures

2.3　温度对全豆豆腐中腐败菌生长特性影响

微生物生长一级模型虽可拟合腐败菌生长趋势，但无法描述温度对腐败菌生长参数影响[16]。因此，采用延滞期模型和平方根模型两个二级模型探究温度对腐败菌生长延滞期（Lag time,LT）和最大比生长速率（μmax）影响。

如图5所示，分别采用延滞期模型和平方根模型对一级模型所得延滞期（Lag time）和最大比生长速率（μmax）非线性回归拟合。将曲线拟合范围范围延长至5～40℃，计算10和35℃对应延滞期和最大比生长速率预测值，构建10和35℃一级模型。

图5A为Logistic方程对应二级模型（平方根模型），平方根模型显示当温度在10℃时，模型出现收敛趋势，经计算得出10℃对应μmax值为0.0005。此外，当温度为5～10℃时，平方根模型函数曲线回升，不符合微生物生长规律，该范围内μmax值不可信，延滞期和最大比生长速率预测值和实际值线性相关系数R2分别为0.98573和0.95249。图5B为Modified Gompertz方程对应二级模型（平方根模型），平方根模型显示当温度在10℃时，模型未出现收敛趋势，延滞期和最大比生长速率预测值和实际值线性相关系数R2分别为0.98491和0.99032，Modified Gompertz最大比生长速率实际值与预测值相关性高于Logistic方程，延滞期实际值与预测值相关性与Logistic方程相近。

表3显示二级模型拟合度、微生物生长参数实际值、预测值，其中Logistic方程对应延滞期模型和平方根模型拟合度（R2）依次为0.98和0.95，Mod⁃ified Gompertz方程对应延滞期模型和平方根模型拟合度（R2）依次为0.98和0.99，Modified Gompertz方程二级模型拟合精度高于Logistic方程。腐败菌生长参数预测值指对应自变量（本试验中为温度）带入二级模型方程后计算出因变量值，而预测值和实际值之间相对误差反映预测值与实际值之间差距，如表3所示，除去15℃下延滞期预测，其余情况下Logistic方程对应延滞期和最大比生长速率相对误差均大于Modified Gompertz方程，因此Modi⁃fied Gompertz方程预测精度更高，更适用于全豆豆腐腐败菌生长趋势预测。

图5　不同温度下腐败菌延滞期及最大比生长速率回归曲线与相关性曲线Fig.5　Regression and correlation curves of lag time（LT）and maximum growth rate（μmax）of spoilage bacterium at different temperature

表3　全豆豆腐腐败菌延滞期及最大比生长速率预测值与实际值Table 3　Observed and predicted lag time(LT)and maximum growth rate(μmax)of spoilage bacteria in whole-soybean curd

2.4　全豆豆腐货架期预测

全豆豆腐初始菌数较高，货架期计算尤为重要。当豆腐中菌数达6 log（cfu·g-1）时，豆腐货架期达到终点[13]，利用Microsoft Excel 2014对一级模型单变量求解，目标值设置为6，迭代后计算自变量即为豆腐货架期。此外，本试验将2.3中计算10及35℃下腐败菌生长参数（延滞期，最大比生长速率）带入一级模型，给出该温度下全豆豆腐货架期预测值。

如表4所示，随温度降低全豆豆腐实际货架期均延长，15℃条件下Logistic方程和Modified Gomp⁃ertz方程对应实际货架期依次为41.38 h（1.72 d）和40.73（1.70 d），预测货架期依次为164.60 h（6.86 d）和40.25 h（1.68 d），由此可见Logistic方程预测值相对误差明显大于Modified Gompertz方程，预测值可信度较低，此外，10℃时Logistic方程预测货架期达12 374.03 h（515.58 d），与实际情况偏差极大，原因为Logistic方程二级模型计算μmax值不准确，35℃下预测货架期为5.22 h。

由图6可知，10℃下全豆豆腐预测货架期为143.96 h（5.99 d），35℃下预测货架期为3.36 h，即10℃下，全豆豆腐货架期可接受。

表4　不同温度下全豆豆腐货架期预测值和实测值Table 4　Predicted and observed shelf life of whole-soybean curd stored at different temperature

图6　10及35℃下全豆豆腐货架期预测模型Fig.6　Shelf life prediction of whole-soybean curd at 10 and 35℃

3　讨论与结论

国内关于豆制品风险评估质量控制体系尚未建立，大豆生长和收获过程中，主要微生物风险因素包括大豆原料与土壤接触（如土壤中常见B.ce⁃reus对原料污染）和灌溉水质量。与豆腐生产相关危险因素包括受污染原材料（大豆），生产过程中水质，生产线污染等。豆腐生产加热过程可杀死存在蜡样芽胞杆菌营养细胞，但蜡样芽孢杆菌孢子和催吐剂毒素耐热，条件合适（环境温暖且潮湿）孢子可萌发，营养细胞在食物中再次生长（可发生在豆腐制造加热过程后）。B.cereus可能是引发豆腐腐败变质特异性腐败菌，且B.cereus可能来源大豆原料本身。为延长豆腐货架期，针对大豆原料杀菌和质量控制尤为关键。

豆腐二次杀菌工艺和冷链运输、贮藏处理可最大限度抑制孢子萌发和萌发后营养细胞生长。Agata等研究表明，毒素产生以蜡样芽孢杆菌大量繁殖为前提，其中呕吐毒素产生需要106～108cfu·g-1B.cereus[17-18]。因此，如蜡样芽孢杆菌被认定为豆腐变质特异性腐败菌（Specific spoilage organism），其在食物中含量可作为评价食品品质标准之一。

预测微生物学模型在实际生产中应用广泛，但食品成分复杂，基于培养基计数机械性模型已不适用。此外，虽然不同恒定温度下微生物生长曲线拟合情况较好，但食物实际储藏环境多变，恒定温度下设计模型难以反映变动环境实际情况，当冷链断裂或贮藏温度变化时，微生物初始生理学状态参数发生改变。综上所述，预测微生物学模型在豆腐中应用应关注波动温度条件下腐败微生物生长特性改变，可结合感官、理化检验，为豆腐货架期和食用品质提供准确判断。

综上所述，本试验结果表明，全豆豆腐中腐败菌经16S rDNA测序后鉴定为T-1:Bacillus tequi⁃lensis或 Bacillus subtilis subsp.Inaquosorum，T-2:Bacillus tequilensis。15～30℃温度下 Modified Gomp⁃ertz模型拟合精度（R2）依次为 0.99、0.98、0.98、0.99，均高于Logistic模型。随豆腐贮藏温度由30℃降至15℃，腐败菌最大比生长速率（μmax）降低，延滞期（Lag time）增大。Modified Gompertz方程对应延滞期模型和平方根模型拟合度（R2）依次为0.98和0.99，拟合精度高于Logistic方程，除去15℃下延滞期预测，Logistic方程对应延滞期和最大比生长速率相对误差均大于Modified Gompertz方程。15℃条件下Logistic方程和Modified Gompertz方程对应实际货架期依次为41.38 h（1.72 d）和40.73（1.70 d），在10℃时，Logistic方程货架期预测值无效，Modified Gompertz方程预测出全豆豆腐货架期为143.96 h（5.99 d）。