The 3D surfaces of the interactive effects of fermentation time (A), fermentation temperature (B), and inoculation volume (C) on the NK are illustrated in Figure 5. Figure 5(A1) shows with the increase in temperature and time, the activity of nattokinase decreased at first and then increased. As shown in Figure 5(B1), with the extension of the fermentation time, the activity of nattokinase showed an increasing trend. This is because the longer the fermentation time, the more time the bacteria will have to grow and ferment, resulting in more nattokinase being produced by the bacteria. The combined effect of inoculum amount and fermentation temperature revealed that nattokinase activity was affected by fermentation temperatures between 31 and 37 ◦C, but the effect was minor (Figure 5(C1)).

Similar to nattokinase activity, Table S3 (Supplementary Materials) shows that the linear terms A-B, cross-product terms AB and AC, and quadratic terms A2 -C2 had a significantly (*p* < 0.05) positive influence on the free amino nitrogen (FAN)content, while the linear terms C, cross-product terms BC and quadratic terms B2 presented a significantly (*p* < 0.05) negative effect (equation Y2).

The level of FAN increased as the temperature and fermentation time increased (Figure 5(A2)). The fermentation period influences the accumulation of products. The free amino nitrogen concentration first declined and subsequently increased as the amount of inoculation and fermentation time increased (Figure 5(B2)).

Table S3 (Supplementary Materials) shows that the linear term C, cross-product term BC, and quadratic term B2 had a significantly (*p* < 0.05) positive influence on the sensory score, while the linear terms A-B, cross-product terms AB and AC, and quadratic terms A2 -C2 presented a significantly (*p* < 0.05) negative effect (equation Y3). The sensory score of natto climbed as fermentation time, fermentation temperature, and inoculum amount increased, reaching a maximum. Then, it began to decline as these three elements continued to rise (Figure 6).

**Figure 5.** Response surface diagram of NK activity and FAN content as a response to the interaction of: (**A1**,**A2**) fermentation temperature and fermentation time (inoculation amount = 6%), (**B1**,**B2**) inoculation amount and fermentation time (temperature = 34 °C), (**C1**,**C2**) inoculation amount and fermentation temperature (time = 30 h). **Figure 5.** Response surface diagram of NK activity and FAN content as a response to the interaction of: (**A1**,**A2**) fermentation temperature and fermentation time (inoculation amount = 6%), (**B1**,**B2**) inoculation amount and fermentation time (temperature = 34 ◦C), (**C1**,**C2**) inoculation amount and fermentation temperature (time = 30 h).

**Figure 6.** Response surface diagram of sensory score as a response to the interaction of: (**A**) fermentation temperature and fermentation time (inoculation amount = 6%), (**B**) inoculation amount and fermentation time (temperature = 34 °C), (**C**) inoculation amount and fermentation temperature (time = 30 h). **Figure 6.** Response surface diagram of sensory score as a response to the interaction of: (**A**) fermentation temperature and fermentation time (inoculation amount = 6%), (**B**) inoculation amount and fermentation time (temperature = 34 ◦C), (**C**) inoculation amount and fermentation temperature (time = 30 h).

#### *3.10. Determination and Verification of Natto Process 3.10. Determination and Verification of Natto Process*

rately predict the real fermentation data for natto.

After response surface optimization, according to the second-order model formula, the theoretical optimal conditions were: sucrose addition was 1%, inoculation amount was 6.96%, fermentation temperature was 35.44 °C, and fermentation time was 24 h. Under these conditions, the predicted value of nattokinase activity was 132.56 FU/g, the sensory score predictive value was 81.21, and the free amino nitrogen content was 4.34 mg/kg. The preceding parameters were slightly modified as follows: sucrose at 1%, inoculation amount at 7%, fermentation temperature at 35.5 °C, and fermentation period at 24 h, tak-After response surface optimization, according to the second-order model formula, the theoretical optimal conditions were: sucrose addition was 1%, inoculation amount was 6.96%, fermentation temperature was 35.44 ◦C, and fermentation time was 24 h. Under these conditions, the predicted value of nattokinase activity was 132.56 FU/g, the sensory score predictive value was 81.21, and the free amino nitrogen content was 4.34 mg/kg. The preceding parameters were slightly modified as follows: sucrose at 1%, inoculation amount at 7%, fermentation temperature at 35.5 ◦C, and fermentation period at 24 h, taking into account the feasibility of the actual operation.

ing into account the feasibility of the actual operation. In order to verify the effectiveness and feasibility of the model, three replicate tests were performed under the optimal culture conditions. The results show that the activity of nattokinase was 144.83 ± 2.66 FU/g, the sensory score was 82.43 ± 5.40, and the free amino nitrogen content was 7.02 ± 0.69 mg/Kg. All of the above values were in line with the expected outcomes. As a result, we came to the conclusion that the model could accu-In order to verify the effectiveness and feasibility of the model, three replicate tests were performed under the optimal culture conditions. The results show that the activity of nattokinase was 144.83 ± 2.66 FU/g, the sensory score was 82.43 ± 5.40, and the free amino nitrogen content was 7.02 ± 0.69 mg/Kg. All of the above values were in line with the expected outcomes. As a result, we came to the conclusion that the model could accurately predict the real fermentation data for natto.

#### *3.11. Thrombolytic Activity, Fibrin Degradation Ability, and Anticoagulant Activity*

*3.11. Thrombolytic Activity, Fibrin Degradation Ability, and Anticoagulant Activity*  The thrombolytic activities of natto were investigated utilizing the agarose–fibrinogen plate method under the abovementioned ideal circumstances. The fibrinolysis action of natto extract is shown in Figure 7A. The diameter of the dissolution circle given in Table 4 was measured under aseptic conditions. The area of the dissolution circle was calculated according to the diameter (π = 3.14). The area of thrombolysis is the area of the dissolution The thrombolytic activities of natto were investigated utilizing the agarose–fibrinogen plate method under the abovementioned ideal circumstances. The fibrinolysis action of natto extract is shown in Figure 7A. The diameter of the dissolution circle given in Table 4 was measured under aseptic conditions. The area of the dissolution circle was calculated according to the diameter (π = 3.14). The area of thrombolysis is the area of the dissolution circle minus the area of the central pore (1.98 ± 0.39 mm 2 ).

circle minus the area of the central pore (1.98 ± 0.39 mm 2).

**Figure 7.** (**A1**,**A2**) Comparison of the effect of single-bacteria and mixed-bacteria fermented natto on fibrin plate dissolution circle; (**B**) changes in nattokinase activity during fermentation; (**C**) prothrombin condensation inhibition ability. **Figure 7.** (**A1**,**A2**) Comparison of the effect of single-bacteria and mixed-bacteria fermented natto on fibrin plate dissolution circle; (**B**) changes in nattokinase activity during fermentation; (**C**) prothrombin condensation inhibition ability.

**4. Discussion Table 4.** Changes in fibrinolytic activity during fermentation.

proved to be an effective thrombolytic in vitro and in vivo [23]. The health benefits of natto have always attracted people's interest [24], and the improvement in its quality has always Note: Data are mean values of three independent experiments ± standard deviation. Mean values displaying different letters within each row are significantly different according to the Duncan test at 95% confidence level.

> been a goal that people pursue, including improving the flavor and taste of natto, and increasing the content of functional ingredients. The combination of *Bacillus subtilis* GUTU09 and *Bifidobacterium* BZ25 (mixed bacteria) was utilized to ferment natto in this As shown in Table 4, as the fermentation time rose, the corresponding diameter of the dissolving circle increased in proportion to the thrombolytic area, and the thrombolysis became more obvious.

> study. Sensory score, pH, nattokinase, and free amino nitrogen concentration were used as indexes to investigate the impact of various variables on the sensory and functional quality of fermented natto. The indexes of natto decreased with the addition of NaCl (Figure 1(b1,c1)). It may be due to the inhibition of some enzyme activities of the two strains under high- and medium-salt conditions. Yue [25] et al. found that histidine decarboxylase was completely inactivated in high- and medium-salt conditions, but remained active when the salt con-Figure 7B depicts the variations in nattokinase activity throughout single-strain and mixed-strain fermentation. The activity of nattokinase gradually increased as the fermentation period was extended. The maximum activity of nattokinase was found in a 36 h fermentation of mixed bacteria, with an enzyme activity of 161 ± 1.96 FU/g recorded. When comparing the two fermentation methods, single-bacteria fermentation had higher natto enzyme activity before 18 h, and equal the enzyme activity at 18 h. After 18 h, the mixed-bacteria fermentation had higher enzyme activity than single-bacteria, which was similar to the fibrinolytic activity indicated before. Finally, the mixed-bacteria fermentation

method produced higher nattokinase activity, more free amino nitrogen, and a superior sensory evaluation. As a result, it is a nutritious food that is readily available to the general public.

The antithrombotic process was recorded when the crude enzyme extract of natto (151.06 ± 2.34 FU/g) was diluted 50 times, 90 times, and 130 times with no added sample as a control (Figure 7C). The control group reached the maximum fibrin concentration (OD405 = 1.107) at 20 min. When the dilution ratio was 130 times, the formation of fibrin was inhibited and the OD405 was 0.494 at 46 min, and fibrin was completely degraded within 100 min. In contrast, the formation of fibrin in a 90-fold diluted sample solution was inhibited at 38 min, and fibrin was almost completely degraded within 60 min. The inhibitory effect, on the other hand, was better and the time was shorter. When the dilution ratio was 50 times, there was no significant change in the curve, indicating that it had strong anticoagulant activity.

#### **4. Discussion**

The formation of a blood clot in a blood vessel is one of the main causes of cardiovascular diseases [5]. Blood clots are formed from fibrinogen via proteolysis by thrombin, and can be hydrolyzed by plasmin to avoid thrombosis in blood vessels [21]. Natural products with anticoagulant activity can hinder the conversion of fibrinogen into fibrin and play a role in the final stage of the anticoagulation process [22]. Natto offers several health benefits, with nattokinase being the most essential functional component. Nattokinase, a serine protease secreted by *Nattobacillus*, has high fibrinolytic activity and has been proved to be an effective thrombolytic in vitro and in vivo [23]. The health benefits of natto have always attracted people's interest [24], and the improvement in its quality has always been a goal that people pursue, including improving the flavor and taste of natto, and increasing the content of functional ingredients. The combination of *Bacillus subtilis* GUTU09 and *Bifidobacterium* BZ25 (mixed bacteria) was utilized to ferment natto in this study. Sensory score, pH, nattokinase, and free amino nitrogen concentration were used as indexes to investigate the impact of various variables on the sensory and functional quality of fermented natto.

The indexes of natto decreased with the addition of NaCl (Figure 1(b1,c1)). It may be due to the inhibition of some enzyme activities of the two strains under high- and medium-salt conditions. Yue [25] et al. found that histidine decarboxylase was completely inactivated in high- and medium-salt conditions, but remained active when the salt concentration decreased. Meanwhile, because of health concerns associated with high salt intakes, the demand for low-salt foods has surged. High salt intake is linked to diseases such as hypertension and renal impairment [26]. As a result, excessive salt should not be used in the fermentation process. Hypertension can be prevented and controlled by eating low-salt foods [27]. It is very beneficial to people with coronary heart disease and hypertension, and it meets the nutritional needs of modern society.

The effect of raising the sucrose content differs from that of raising the salt concentration (Figure 1). Adding about 1–2% sugar will improve the sensory characteristics. The stringiness, flavor, taste, and chewiness of natto are all affected by the sugar concentration. Despite the presence of oligosaccharides and polysaccharides in the steamed soybean, it lacks a carbon source that BZ25 can utilize directly. Increasing the carbon source sugar that the microbes can consume directly results in accelerated growth and metabolism. Our findings are in line with a previous study by Wu et al. [28]. On the one hand, more sucrose provides a better environment for BZ25 to grow, and acid production as a by-product of the microbe increases rapidly, lowering the pH throughout the fermentation process (Figure 1(b2)) and improving the texture and aroma of natto. Adding sucrose for B9, on the other hand, causes rapid fermentation and the production of numerous enzymes to hydrolyze the macromolecular components in soybeans. This not only promotes BZ25 proliferation but also accumulates a significant number of fermentation products.

In the same fermentation environment, a high strain ratio not only prevents the synergy between the strains from working, but also lowers the quality of the fermented product due to nutrient competition between the two strains [29]. A high strain ratio of BZ25 to GUTU09 also produced much acid (Figure 2(b1)), which led to the decline in sensory quality. The reason may be that the internal oxygen content in the solid fermentation was generally low, and increasing the ratio of *Bifidobacterium* BZ25 to *Bacillus subtilis* GUTU09 can enhance the resistance of *Bifidobacterium* to a small amount of oxygen and make it easy to propagate rapidly and produce more acid in the fermentation. Therefore, it is not beneficial to natto fermentation if the ratio of strains GUTU09 to BZ25 is too large or too small. As a result, the suitable strain ratio was approximately 2:1 or 1:1 (Figure 2(a1,b1,c1)).

The number of initial bacteria in the fermentation process is determined by the amount of inoculation, and the appropriate starting bacteria concentration can shorten the strain's growth and proliferation period, allowing it to enter the fermentation stage sooner. An inoculation quantity ranging from 6 to 8% (Figure 2(a2,b2,c2)) is more appropriate; too high or too low inoculation may impair the natto quality. If the initial bacteria concentration is lower, the time needed for reproduction will be longer and the fermentation speed will be slower, which is not conducive to the accumulation of fermentation products. Higher inoculation concentrations, on the other hand, resulted in a strong and speedy bacterial metabolism, requiring the cell to consume the majority of the nutrients in order to continue its rapid growth and hinder the accumulation of fermentation products. When a large amount of metabolic waste is produced, the bacteria's cell senescence may be accelerated, and its functionality may be lowered [30]. Ultimately, the quality of natto could be compromised.

Fermentation temperature is one of the important parameters affecting fermentation. The suitable growth temperature for *Bacillus subtilis* and *Bifidobacterium* is around 37 ◦C. The activity of nattokinase and the sensory score decreased as the fermentation temperature of natto increased (Figure 3(a1,b1,c1)). The increasing temperature, on the one hand, causes the water in the soybeans to evaporate more quickly. The surface layer of the beans becomes dry, affecting mucus formation and, as a result, the taste. On the other hand, we used a 250 mL flask, which may have led to fermentation, as the heat could not be properly dissipated over time, and an accumulation of heat formed in the core. The high temperature impacts strain development and enzyme reactions, causing uneven natto fermentation and hence affecting nattokinase synthesis [31]. pH increased with the increasing temperature (Figure 3(b1)). The amine bases produced by the breakdown of soybeans and the organic acids produced by BZ25 are the primary sources of the pH changes. The pH of mixed-bacteria fermented natto was practically the same as that of single-bacteria fermented soybean when the temperature was greater than 37 ◦C, indicating that the effect of BZ25 is not noticeable after 37 ◦C. It is possible that the high temperature prevents BZ25 from growing. BZ25's ability to generate acid may be hampered by the higher temperature. The alterations in free amino nitrogen and nattokinase were varied at different temperatures, as shown in Figure 3(c1). The generation of nattokinase was inhibited by high temperatures, while free amino nitrogen levels increased dramatically. The breakdown of soybean protein by excess protease was the predominant source of free amino nitrogen. As a result, high temperatures suppress nattokinase expression while having little to no influence on protease activity.

One of the determining factors of product quality was fermentation time. As the fermentation duration increased (Figure 3(a2,b2,c3)), sensory score decreased, and nattokinase activity and free amino nitrogen increased. It could be because the bacteria in the fermentation reached a steady state of productivity. The catalytic biological reaction did not cease and the fermentation products continued to increase. Nevertheless, longer fermentation times will cause bacteria to enter senescence, resulting in the accumulation of hazardous compounds [32]. As a result, choosing the right fermentation period can help reduce nutrient waste and the formation of toxic metabolites.

After-ripening time can help to minimize natto's unpleasant smell and make it taste richer and fuller. A longer after ripening time will degrade the product's quality and freshness, resulting in a lower sensory score (Figure 4). As a result, selecting the most appropriate after-ripening time is critical.

Through the PB test (Supplementary Materials Table S1 and Table 3) and response surface test (Supplementary materials Tables S2 and S3), it was determined that fermentation time, fermentation temperature, and inoculation amount were significant factors affecting natto fermentation. PB test and response surface test can effectively be used to estimate the effect of fermentation time, fermentation temperature, and inoculation amount, and their interactions on the natto. According to the findings of the study, the optimal fermentation conditions for achieving the finest natto quality were as follows: soybean 50 g/bottle, NaCl was 0%, sucrose 1%, ratio of GUTU09 to BZ25 1:1, inoculum 7%, fermentation temperature 35.5 ◦C, and fermentation time 24 h. The nattokinase activity in the above conditions was 144.83 ± 2.66 FU/g, which was higher than the nattokinase activity in the single-bacteria fermentation procedure. The level of free amino nitrogen was 7.02 ± 0.69 mg/Kg. Furthermore, these optimal conditions resulted in an overwhelming sensory score, with the mixed bacteria fermentation natto scoring much higher on appearance, stringiness, flavor, taste, and chewiness than the non-optimized natto.

Under the optimal circumstances outlined above, the thrombolytic action of natto was investigated (Table 4, Figure 7(A1,A2)). It indicates that as fermentation progressed, strains produced more nattokinase and its thrombolytic activity was improved. The antithrombotic effect includes two aspects: one is the fibrinolysis of the formed thrombus, which is evaluated by the fibrin degradation ability (nattokinase activity); the other is the anticoagulant activity during the formation of fibrin, that is, the inhibition ability of fibrin on the original coagulation [17]. The fibrinolytic activity and anticoagulant activity may reflect the thrombolytic ability to some extent [33].

Single-bacteria fermentation had higher enzyme activity for the first 18 h, but after that, mixed-bacteria fermentation had higher enzyme activity than single-bacteria fermentation (Figure 7B). When both *Bifidobacterium* BZ25 and *Bacillus subtilis* GUTU09 are inoculated at the same time, BZ25 growth is impeded due to the presence of oxygen in the soybean medium at the start of fermentation [34]. According to Hosoi et al. [35], the viability of Lactobacilli in the presence of Bacillus was greatly improved, and they speculated that the production of subtilisin and catalase may play a part in this improvement. Another explanation for the improvement was that the growth of *B. subtilis* GUTU09 in solid-state substrates used dissolved oxygen, allowing *Bifidobacterium* BZ25 to flourish [36]. Simultaneously, the released protease assisted in the hydrolysis of soy protein, providing a nitrogen source for BZ25. Its amylase and other enzymes degraded the polysaccharide in soybeans to give BZ25 a carbon source. These factors promoted the growth of BZ25 [37]. Furthermore, *Bifidobacterium* growth may produce certain compounds that aid Bacillus metabolism. As a result, BZ25 and GUTU09 boosted each other more than single-bacteria fermentation, which had more advantages and was more favorable for bioactive ingredient accumulation. Our experiments revealed that adding a diluent to a natto extract prevented fibrin production, and that the anticoagulant action increased as the enzyme concentration increased (Figure 7C). Thus, mixed-bacteria fermentation methods result in greater nattokinase activity, free amino nitrogen concentration, and sensory assessment. As a result, it is a nutritious food that is readily available to the general public.

#### **5. Conclusions**

Mixed fermentation increased nattokinase activity, free amino nitrogen and the sensory score of natto. Substantial thrombolytic and anticoagulant effects were also observed. A healthy fermented natto with good sensory characteristics was studied and developed.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/foods10112547/s1, Table S1 PB experiment resµlts; Table S2 response surface experimental results; Table S3 analysis of experimental variance with Nattokinase Activity, amino acid nitrogen and sensory scores as response value

**Author Contributions:** Y.Y.: Conceptualization, writing—original draft. G.L.: data curation, formal analysis, writing—original draft. X.T.: data curation, methodology. L.H.: supervising the experiment, writing—review and editing. C.L.: methodology. X.Z.: writing—review and editing. X.W.: writing review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Key Agricultural Project of Guizhou Province (QKHZC— [2021]278, QKHZC—[2021]184, QKHZC—[2021]142, QKHZC—[2019]2382, and QKHZC—[2016]2580), and the Natural Science Foundation of China (31870002 and 31660010), Qiankehe talents project ([2018]5781 and [2017]5788-11).

**Conflicts of Interest:** The authors have no competing interest to disclose.

#### **References**

## *Article* **Metagenomic Study on Chinese Homemade Paocai: The Effects of Raw Materials and Fermentation Periods on the Microbial Ecology and Volatile Components**

**Linjun Jiang † , Shuang Xian † , Xingyan Liu, Guanghui Shen, Zhiqing Zhang, Xiaoyan Hou and Anjun Chen \***

> College of Food Science, Sichuan Agricultural University, Ya'an 625014, China<EMAIL_ADDRESS>(L.J<EMAIL_ADDRESS>(S.X<EMAIL_ADDRESS>(X.L<EMAIL_ADDRESS>(G.S<EMAIL_ADDRESS>(Z.Z<EMAIL_ADDRESS>(X.H.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +86-0835-2882187

† These authors contributed equally to this work.

**Abstract:** "Chinese paocai" is typically made by fermenting red radish or cabbage with aged brine (6–8 *w*/*w*). This study aimed to reveal the effects of paocai raw materials on fermentation microorganisms by metagenomics sequencing technology, and on volatile organic compounds (VOCs) by gas chromatography–mass spectroscopy, using red radish or cabbage fermented for six rounds with aged brine. The results showed that in the same fermentation period, the microbial diversity in cabbage was higher than that in red radish. *Secundilactobacillus paracollinoides* and *Furfurilactobacillus siliginis* were the characteristic bacteria in red radish paocai, whereas 15 species of characteristic microbes were found in cabbage. Thirteen kinds of VOCs were different between the two raw materials and the correlation between the microorganisms and VOCs showed that cabbage paocai had stronger correlations than radish paocai for the most significant relationship between 4-isopropylbenzyl alcohol, α-cadinol, terpinolene and isobutyl phenylacetate. The results of this study provide a theoretical basis for understanding the microbiota and their relation to the characteristic flavors of the fermented paocai.

**Keywords:** red radish; cabbage; fermented foods; microbial ecology; flavor components

#### **1. Introduction**

Fermentation is a widely used technology to preserve foods, improve nutritional value and extend shelf-life, and fermented vegetables are popular and traditional in Asian countries [1]. Among the fermented vegetables, "paocai" is a beloved food with a history of more than 3000 years, which are typically consumed as side dishes or appetizers and are characterized by their tender, crisp texture and rich health benefits owing to the lactic-acid bacteria (LAB) [2–4]. Vegetables for making paocai are usually immersed in 6–8% (*w*/*w*) sodium chloride solution and allowed to undergo spontaneous anaerobic or microaerobic fermentation by the epiphytic microbes (mostly LAB) present on the raw materials for 6–10 days [5].

The microbial composition and flavor in paocai fermentation systems are related to various factors, especially the raw materials. Previous studies have found that the microbial niches in raw ingredients determine the eventual microbial community [6], this conclusion has been confirmed in some studies on microorganisms in Sichuan paocai; for instance, *Rosenbergiella*, *Staphylococcus*, *Hyphopichia* and *Kodamaea* dominated the fermentation of chili pepper [7], while *Lactobacillus*, *Leuconostoc*, *Achromobacter* and *Pediococcus* were the main microbial organisms during the fermentation of cabbage pickle [8]. This is due to the fact that the spontaneous fermentation of vegetable products is highly dependent on the naturally occurring LAB present on the raw materials. A large amount of different LAB can be accumulated in the aged brine in the process of continuous periodic fermentation, and

**Citation:** Jiang, L.; Xian, S.; Liu, X.; Shen, G.; Zhang, Z.; Hou, X.; Chen, A. Metagenomic Study on Chinese Homemade Paocai: The Effects of Raw Materials and Fermentation Periods on the Microbial Ecology and Volatile Components. *Foods* **2022**, *11*, 62. https://doi.org/10.3390/ foods11010062

Academic Editors: Valentina Bernini and Juliano De Dea Lindner

Received: 5 November 2021 Accepted: 22 December 2021 Published: 28 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the total viable count can be up to 8.79 Log CFU/mL [9]. It is worth noting that the microorganisms were shown to rapidly build a stable fermentation system at the initial stages and finally obtain paocai products with different flavors. Due to the complex microbial ecology in the fermentation system, two or three strains of artificially isolated LAB are often used in industry to accelerate the fermentation process and to achieve a uniform quality of the products [10]. However, based on consumer surveys, the quality of the final products is not as good as homemade paocai. Homemade paocai usually uses aged brine as the starter, which is brine fermented for many years and even decades [11,12]. The microorganisms in the fermentation system metabolize the nutrients in the raw materials to produce acids, that makes the fermented vegetables have a unique flavor [13], which is another source of paocai flavor besides that of the raw materials. Therefore, understanding the effects of raw materials and fermentation periods on the microbial ecology and volatile components is of far-reaching significance for understanding how microorganisms construct stable fermentation systems.

Metagenomics sequencing technology can extract all microbial DNA from the fermentation system, construct a metagenomic PE library and use genomics research strategy to study the genetic composition and gene function of all microorganisms contained in environmental samples. This technique can avoid amplification and sequencing bias, which commonly occurred in other sequencing techniques [14]. Metagenomics, although still underexploited, can reveal the succession and function of microorganisms and establish the relationship between microorganisms and substance metabolism [15–17]. However, to the extent of authors' knowledge, little insights on relationship between "raw material– microorganism-VOCs" in the process of paocai fermentation are available, especially based on metagenomics sequencing technology.

Therefore, in this study, red radish and cabbage, the most common paocai raw materials in Sichuan, were used as raw materials for multiround fermentation, and the composition of microbial and volatile components in paocai fermentation systems were analyzed by metagenomics and gas chromatography–mass spectrometry (GC–MS), in order to reveal the effects of vegetable raw materials on the microbial community structure during Chinese homemade paocai multiple fermentation rounds and their relationship with volatile compounds.

#### **2. Materials and Methods**

#### *2.1. Paocai Preparations*

Red radish (*Red raphanus* L.) and cabbage (*Brassica oleracea* L.) were purchased from the local market in Ya'an, Sichuan Province, China. The fresh vegetables were cleaned under the same conditions to ensure the consistency of the initial flora on the surface of the vegetables. Paocai was made according to the local traditional method, specifically: one kilogram of each vegetable was cut into 2–3 cm pieces and set in 4 L pottery jars with garlic (3%), hot red peppers (3.5%), ginger (3%), and Chinese prickly ash (3%) (all percentages were calculated according to the volume of brine in the pottery jars, in *w*/*v*). Two liters of brine was prepared with cooled, boiled water containing 8% salt, 1% Chinese baijiu and 20% aged Sichuan paocai brine (in terms of *w*/*v*, collected from the homemade paocai of a local family used for more than 5 years). The pottery jars were sealed by adding water to the jar edge and stored at 24 ± 2 ◦C for six rounds of fermentation for about 9 d, 7 d, 5 d, 4 d and 3 d, respectively.

#### *2.2. DNA Extraction, Library Construction and Metagenomic Sequencing*

A total of 50 mL of initial aged brine and paocai brine at the end of each fermentation round were enriched on a 0.22 µm sterile microfiltration membrane. According to the manufacturer's recommendation, the microbial DNA from brine samples was extracted with the E.Z.N.A.® DNA Kit (Omega Bio-tek, Norcross, GA, USA). The concentration and purity of the DNA were determined with Turner Biosystems (Sunnyvale, CA, USA) and a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, TBS-380, Shanghai, China), and the integrity of the DNA was detected by 1% *w*/*w* agarose gel electrophoresis.

DNA was fragmented to an average size of about 400 bp using Covaris M220 (Gene Company Limited, Hong Kong, China) for paired-end library construction, and the paired-end library was constructed using NEXTflexTM Rapid DNA-Seq (Bioo Scientific, Austin, TX, USA).

Adapters containing the full complement of sequencing primer hybridization sites were ligated to the blunt end of fragments. Paired-end sequencing was performed on Illumina NovaSeq (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using NovaSeq 6000 according to the manufacturer's instructions. Sequence data associated with this project have been deposited in the NCBI Short Read Archive database (accession number: PRJNA731323).

#### *2.3. Metagenome Data Integration*

The raw reads from the metagenome sequencing were used to generate clean reads by removing adaptor sequences, trimming and removing low-quality reads (length < 50 bp, or with a quality value < 20 or having N-bases) using the fastp [18] (https://github.com/ OpenGene/fastp, version 0.20.0, accessed on 1 December 2020) on the free online platform of Majorbio's Cloud Platform (https://cloud.majorbio.com/, accessed on 1 December 2020). Megahit (parameters: kmer\_min = 47, kmer\_max = 97, step = 10) (http://www.l3-bioinfo. com/products/megahit.html, version 1.1.2, accessed on 6 December 2020) based on the principle of succinct de Bruijn graphs was used to assemble the optimized sequences [19]. Among the splicing results, contigs ≥ 300 bp were selected as the final assembly result. MetaGene was used to predict the open reading frames (ORFs) of contigs in the stitching results (http://metagene.cb.k.u-tokyo.ac.jp/, accessed on 11 December 2020) [20].

The predicted ORFs with lengths over 100 bp were retrieved and translated to amino acid sequences using the NCBI translation table (https://www.ncbi.nlm.nih.gov/ Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1, accessed on 15 December 2020). CD-HIT (http://www.bioinformatics.org/cd-hit/, version 4.6.1, accessed on 15 December 2020) was used to classify all sequences with 90% sequence identity and 90% coverage as nonredundant gene catalog [21]. Reads after quality control were mapped to the representative genes with 95% identity using SOAPaligner (https://github.com/ ShujiaHuang/SOAPaligner, version 2.21, accessed on 20 December 2020) [22], and gene abundance in each sample was evaluated.

#### *2.4. Species Annotation*

Diamond (https://github.com/bbuchfink/diamond, version 0.8.35, accessed on 21 December 2020) was employed for taxonomic annotations by aligning nonredundant gene catalogs against the NCBI NR database using blastp as implemented in DIAMOND v0.9.19 with an e-value cutoff of 1 × 10−5 [23].

#### *2.5. Determination of VOCs Using GC–MS*

The analysis of the VOCs was performed by HS–SPME combined with GC–MS (7890A GC, Agilent, USA, 5975C MS; Agilent Technologies, Santa Clara, CA, USA). The paocai brine (5 mL) were accurately weighed into 15 mL headspace glass vials containing 1 g of salt, the vials were immediately placed in a heating block to equilibrate for 30 min at 50 ◦C and the VOCs were extracted using a 65 µm 50/30 um DVB/CAR/PDMS optical fiber for 10 min at 50 ◦C. VOCs were desorbed from the SPME fiber at 250 ◦C for 5 min in the injector in splitless mode. The carrier gas was helium, and it was supplied at a flow rate of 1.2 mL/min (constant linear velocity). The column temperature was programmed as follows: maintained at 36 ◦C for 3 min, increased at 5 ◦C/min to 65 ◦C, increased at 3 ◦C/min to 155 ◦C, increased at 10 ◦C/min to 200 ◦C, and maintained at 200 ◦C for 3 min. The MS condition was set as: 280 ◦C for the transfer line, 230 ◦C for the ionization, 70 eV for the ionization energy [10]. The quantification analysis was done by using cyclohexanone (27.424 g/L, *w*/*v*) as an internal standard. The NIST11.L standard mass spectral database

was used to identify VOCs based on retention time and mass spectral similarity match. Each sample was analyzed in triplicates at each sampling time.

#### *2.6. Statistical Analysis*

The nonmetric multidimensional scaling analysis (NMDS), the analysis of similarities (ANOSIM), the heatmap and the bar plot were performed in the Vegan package of R (version 3.3.1). The differences between classes were analyzed using Linear discriminant analysis Effect Size (LEfSe) (http://huttenhower.sph.harvard.edu/galaxy/root?tool\_id=lefse\_upload, Version 1.0, accessed on 1 December 2020). SIMCA-P (Version14.0 (Umetrics, Umeå, Sweden) was used to analyze VOCs by partial least squares discriminant analysis (PLS-DA).

#### **3. Results and Discussion**

#### *3.1. Metagenomic Sequencing Statistics and Quality Control*

A total of 39 samples were sequenced by the Illumina NovaSeq 6000 sequencing systems, and 1708.39 × 106 raw reads, about 257.97 GB, with an average of 6.61 GB of each sample were generated. After the quality processing and eliminating the sample host gene, the sequence utilization of each sample was more than 98.48, including 787.85 × 106 clean reads in the red radish samples and 776.19 × 106 clean reads in the cabbage samples. A total of 966,355 contigs were assembled with Megahit, and the single sample's contigs were between 9431 to 47,483. N50 and N90 were two indices for the distribution of contig lengths within a draft assembly, and the longer the assembly of samples the better. In all samples, the minimum value of N50 reached 1677 bp, the average value was 6680 bp, the N90 minimum was 421 bp, and the average value was 528 bp. Furthermore, a total of 1,778,731 genes were obtained by ORF prediction of the contigs, and the sequence length was 1,150,522,526 bp (Table S1). All predicted genes were clustered using CD-HIT (parameters: identity, 95; coverage, 90), resulting in a total of 122,012 nonredundant gene set. Figure 1A shows the sequence length distribution of nonredundant gene sets, mostly concentrated between 201~600 bp, with an average sequence length of 557.88 bp, and the information was used for subsequent species and functional annotations to reveal microbial community of samples.

#### *3.2. Microbial Structure*

The sequences comprising the total reads corresponded to five domains, eight kingdoms, 47 phyla, 99 classes, 207 orders, 353 families, 706 genera, and 2068 species, of which four domains were known. These four were identified to be bacterial, fungal, archaeal and virus domains, of which the bacterial domain had the highest proportion at 99.31% and the lowest was for the fungal domain. In addition, eight kingdoms in all samples were detected, of which the major phyla were *Firmicutes* and *Ascomycota* with *Firmicutes* the most abundant, contributing to over 98% in all samples, consistent with previous reports of other fermentation vegetables [24,25]. Further, it was interesting to note that with the progress of fermentation, the content of *Ascomycota* decreased gradually in the red radish samples, but there was no obvious change in the cabbage samples. The major classes, order and family were *Bacilli* and *Saccharomycetes*, *Lactobacillales* and *Saccharomycetales* and *Lactobacillaceae* and *Saccharomycetaceae*, consistent with the phyla identified.

**Figure 1.** (**A**) Sequence length distribution of nonredundant gene catalog. (**B**) NMDS analysis at species level. Analysis of similarities between different groups by ANOSIM (**C**,**D**) of red radish paocai and cabbage paocai in the 3rd and 4th round, separately. (**E**) LDA score of the two kinds of red radish paocai and cabbage paocai. (**F**) Wilcoxon rank-sum test of different species in red radish paocai and cabbage paocai. \* 0.01 < *p* ≤ 0.05, \*\* 0.001 < *p* ≤ 0.01, \*\*\* *p* ≤ 0.001.

As for the genus level, a total of 706 microbe genera were detected, of which 319, 521 and 664 genera were detected in the initial red radish and cabbage brines, respectively. Among them, 303 genera occurred in all paocai samples, wherein 38 genera only appeared in the red radish samples, 171 genera only appeared in the cabbage samples and 303 genera appeared in both paocai samples (Figure 2A). The microbial diversity of paocai was found to be relatively simple, and there were only three genera with abundance more than 1% of the diversity, and they belonged to the *Lactobacillus*, *Pediococcus* and an unclassified *Lactobacillaceae* genus. The average abundance of *Lactobacillus* in the initial brine, red radish paocai and cabbage paocai reached 92.35%, 85.10% and 81.84%, respectively. Furthermore,

cai.

*Foods* **2022**, *11*, x FOR PEER REVIEW 6 of 14

brines. Among the detected species 826 of them occurred in all paocai samples (Figure 2C), whereas 129 species appeared only in the red radish samples, and 505 species only in the cabbage samples. A further 564 species commonly appeared in all the red radish and cabbage paocai samples, but only 11 of them were recorded with relative abundances of more than 1%. In addition, a higher number of microbial species were present in the cabbage paocai than in the red radish paocai with a brine of the same age; moreover, with the increase of fermentation time, the number of microorganisms in the red radish paocai increased first and then decreased, while the number of microorganisms in the cabbage paocai did not obviously change, but after six rounds of fermentation, the number of microorganisms in the cabbage paocai was significantly higher than that in the red radish pao-

As for the individual detected species, the *Lentilactobacillus buchneri* (9.9%) and *Lacto-*

Further, *Lactiplantibacillus pentosus,* which is used as the starter in the production of

increased gradually until the third round of fermentation, then decreased continuously

*bacillus acetotolerans* (12.5%) in the initial brine were found to be the dominant species. *Lactiplantibacillus paraplantarum* species were found to be gradually decreasing with the progression of fermentation rounds from 29.1% to 8.6% in the cabbage paocai and from 29.1% to 7.9% in the red radish paocai samples. With regard to the *Pediococcus ethanlidurans*, *Pediococcus damnosus* and *Secundilactobacillus paracollinoides* species, their proportions were found to be gradually increasing with fermentation rounds. Specifically, *P. ethanlidurans* species registered a significant increase from the first to the last round of fermentation from 0.7% to 20.67% in the red radish and 0.94% to 22.39% in the cabbage paocai samples, respectively. This occurrence indicates that *P. ethanlidurans* might have a better growth competitiveness compared to *L. paraplantarum* under the final round of fer-

mentation conditions characterized by a high salinity and acidity.

*Lactobacillus* was the predominant bacteria in the multiround fermentation, and was found to be in agreement with other paocai studies [8]. However, with the progress of fermentation, *Lactobacillus* gradually decreased, accounting for only 67.4% of red radish paocai at the end of the sixth round of fermentation and 69.5% of cabbage paocai. Although previous studies have different conclusions about whether *pediococcus* is the dominant bacterium in fermented vegetables [12,26,27], the proportion of *Pediococcus* genus was found in general to be increasing along with the fermentation. Specifically, its abundance was found to increase from 2.82% to 28.37% in the red radish and 2.58% to 26.28% in the cabbage paocai brines (Figure 2B), respectively. These data indicate that although they belong to the same phylum (Firmicutes), *Pediococcus* may be more stress-resistant or competitive than *Lactobacillus* in the later stage of fermentation. until the last round of paocai fermentation (from 38.8% to 18.0% and from 47.4% to 12.2% in red radish samples and cabbage samples, respectively). This could be due to the preparation process and some limitations as reported in previous studies [28,29]. The microbial community in the later stages of fermentation for both the red radish paocai and cabbage paocai was found to be similar. The proportion of *Lentilactobacillus parafarraginis* was found to be relatively higher in the cabbage paocai (15.2%) than that in the red radish paocai (7.3%), whereas the relative abundance of *S. paracollinoides* was higher in the red radish paocai (17.7%) than that in the cabbage paocai (1.7%) (Figure 2D).

**Figure 2.** Venn plot (**A**) and bar plot of relative abundance (**B**) of microbial communities at the species level in the initial brine and two kinds of paocai and (**C**,**D**) at the genus level in the initial brine and two kinds of paocai. R indicates red radish paocai, B indicates cabbage paocai, IB indicates initial brine. Note: others in B and D means other microbes.

At the species level, a total of 2068 microbial species were detected by sequencing, of which 870, 1525 and 1923 species were detected in the initial, red radish and cabbage brines. Among the detected species 826 of them occurred in all paocai samples (Figure 2C), whereas 129 species appeared only in the red radish samples, and 505 species only in the cabbage samples. A further 564 species commonly appeared in all the red radish and cabbage paocai samples, but only 11 of them were recorded with relative abundances of more than 1%. In addition, a higher number of microbial species were present in the cabbage paocai than in the red radish paocai with a brine of the same age; moreover, with the increase of fermentation time, the number of microorganisms in the red radish paocai increased first and then decreased, while the number of microorganisms in the cabbage paocai did not

obviously change, but after six rounds of fermentation, the number of microorganisms in the cabbage paocai was significantly higher than that in the red radish paocai.

As for the individual detected species, the *Lentilactobacillus buchneri* (9.9%) and *Lactobacillus acetotolerans* (12.5%) in the initial brine were found to be the dominant species. *Lactiplantibacillus paraplantarum* species were found to be gradually decreasing with the progression of fermentation rounds from 29.1% to 8.6% in the cabbage paocai and from 29.1% to 7.9% in the red radish paocai samples. With regard to the *Pediococcus ethanlidurans*, *Pediococcus damnosus* and *Secundilactobacillus paracollinoides* species, their proportions were found to be gradually increasing with fermentation rounds. Specifically, *P. ethanlidurans* species registered a significant increase from the first to the last round of fermentation from 0.7% to 20.67% in the red radish and 0.94% to 22.39% in the cabbage paocai samples, respectively. This occurrence indicates that *P. ethanlidurans* might have a better growth competitiveness compared to *L. paraplantarum* under the final round of fermentation conditions characterized by a high salinity and acidity.

Further, *Lactiplantibacillus pentosus,* which is used as the starter in the production of homemade paocai to improve the quality [5], was observed to be the dominant strain, and increased gradually until the third round of fermentation, then decreased continuously until the last round of paocai fermentation (from 38.8% to 18.0% and from 47.4% to 12.2% in red radish samples and cabbage samples, respectively). This could be due to the preparation process and some limitations as reported in previous studies [28,29]. The microbial community in the later stages of fermentation for both the red radish paocai and cabbage paocai was found to be similar. The proportion of *Lentilactobacillus parafarraginis* was found to be relatively higher in the cabbage paocai (15.2%) than that in the red radish paocai (7.3%), whereas the relative abundance of *S. paracollinoides* was higher in the red radish paocai (17.7%) than that in the cabbage paocai (1.7%) (Figure 2D).

#### *3.3. Microbial Composition Comparison*

The heatmap analysis was utilized to reveal the relative relationship of the top 40 microbial genera detected with the colors representing the species abundance (Figure 3A). It can be clearly observed from the heatmaps that in both evaluated samples, the abundance of *Lactobacillus* and *Pediococcus* in bacteria and *Kazachstania* in eukaryotes was relatively stable during all six rounds of fermentation. Moreover, *Weissella* was found to be relatively rich during the whole fermentation process, but its relative abundance decreased as the fermentation progressed. This decrease in *Weissella*'s abundance can be attributed to the increase in the types of species and the subsequent competition for nutrients, as was also observed in other studies [30,31]. In addition, the abundance of *Geotrichum* and *Kazachstania* decreased gradually during the fermentation in the two fermented vegetables, whereas the abundance of *Enterococcus*, *Streptococcus*, *Clostridium* and *Bacillus* showed a trend toward a gradual increase in the fermentation process, with no change in their relative abundance in both fermentation systems even though they exhibited significantly opposite microorganisms.

In fact, metagenomics sequencing allows the obtained data to be compared with a database to extract the information of microorganisms at the species level [32,33]. Within the top 40 abundant species detected (Figure 3B), 30 species belonged to the genus of LAB. Within the LAB species, a varying behavior of abundance was observed along the fermentation process with species such as *L. paraplantarum*, *Levilactobacillus brevis*, *L. acetotolerans* decreasing in abundance during the fermentation in both samples.

A UPGMA cluster analysis of the samples was performed based on the identified genus and species, demonstrating the higher similarity between the red radish samples.

**Figure 3.** *Cont*.

To further understand the microbiota structure during paocai fermentation, NMDS and LDA Effect Size (LEfSe) was used to analyze the microbial differences. The NMDS analysis (Figure 1B) showed that the microbial community of the red radish and cabbage paocais displayed an obvious spatial structure during the fermentation process, where the composition of microorganisms in different vegetable fermentation systems exhibited greater differences (stress = 0.016 < 0.05). With the progress of the fermentation, the distance between spatial positioning of different fermented vegetable samples of the same round gradually increased, indicating that the influence of raw materials on the microbial composition increased with the increase of fermentation time. These observed results were consistent with the ANOSIM analysis, with the R value closer to 1 confirming the greater differences between groups with a *p*-value less than 0.05 indicative of the high reliability of

the test (Figure 1C,D). Towards the end of the third and fourth rounds of paocai fermentation, the different samples were significantly separated (R = 0.46, *p* = 0.009) with a further increase in distance at the end of the fifth and the sixth fermentation rounds (R = 0.89, *p* = 0.003). The aforementioned data further confirm that the microbial composition of different vegetable paocais presented continuous changes in the multiround fermentation with the differences significantly increasing along the fermentation time.

LEfSe (Figure 1E) was used as an algorithm for high-dimensional biomarker discovery and identifying genomic features characterizing the differences between two or more biological conditions [34]. In this study, LEfSe was performed to evaluate the differences between the bacterial community structures of the different paocais. In the paocai, a total of 31 taxa were found to represent a remarkable difference in their relative abundance, with an LDA (log10) > 2. At the species level, *S. paracollinoides* and *Furfurilactobacillus siliginis* were found to be the characteristic bacteria of the red radish paocai, whereas a total of 15 species of microbes were found characteristic of the cabbage paocai, of which 13 species belonged to the Firmicutes (namely, *L. parafarraginis*, *L. buchneri*, *Levilactobacillus zymae*, *Levilactobacillus namurensis*, *Lentilactobacillus diolivorans*, *Ligilactobacillus acidipiscis*, *Lacticaseibacillus casei*, *Lentilactobacillus parabuchneri*, *Lentilactobacillus hilgardii*, *Lentilactobacillus farraginis*, *Paralactobacillus selangorensis*, *Lentilactobacillus sunkii*, *Lentilactobacillus otakiensis*) and two species to that of Proteobacteria (namely, *Kluyvera ascorbata*, *Pectobacterium carotovorum*). Further, the differential microorganisms were screened by a Wilcoxon ranksum test (Figure 1F), which revealed that only nine species were significantly different in the red radish paocai and cabbage paocai, with *L. parafarraginis*, *S. paracollinoides*, *L. namurensis*, *P. selangorensis* and *F. siliginis* being the significantly different species in the paocais tested. In fact, all the nine species detected were also screened by LEfSe and regarded as biomarkers for the subsequent analysis.

#### *3.4. Profiles of VOCs during Different Paocai Fermentations*

The compositions of VOCs in the paocai samples are shown in Table S2. A total of 62 and 59 volatile substances were identified in the red radish and cabbage paocais. As shown in Figure 4A, the main VOCs detected were alcohols, hydrocarbons, esters, sulfide and ketones during the whole fermentation process. Alcohols were the most abundant compounds among those detected and ranged from 74.1% to 79.3%, with a total of 17 different components. The second largest group was heterocycles which ranged from 9.8% to 13.1% and included 23 component types.

Figure 4B shows 22 major VOCs that accounted for more than one of the two paocai volatile profiles. In general, the types of VOCs identified were similar, both types included compounds like linalool, α-terpineol, eucalyptol, D-carvone and other substances. Among these, linalool registered the highest proportions reaching up to 35% and 41.6% in the red radish and cabbage paocais, respectively, similar to previous studies [8]. Linalool contributes to the unique flavor of paocai and is described as having a strong scent, similar to bergamot oil or French lavender, with a mixture of woody, floral aromas with a touch of spiciness [35]. Other alcohols were also evaluated to have a positive effect on paocai's flavor with α-terpineol contributing a floral typically lilac odor; eucalyptol with a eucalyptus, herbal and camphor odor; and nerol with notes of roses [36,37].

Red radish and cabbage belong to the Cruciferae family of plants, and typically have a strong aroma, bitter taste and pungent taste similar to that of mustard. The compounds responsible for these special flavors were found to be mainly the secondary defense metabolite, glucosinolate and its hydrolysate isothiocyanate [38]. The composition and content of isothiocyanates is variable based on the type of vegetables subsequently imparting characteristic flavor components and compositions in the fermented paocais. Both tested paocais were found to contain butyl isothiocyanate and amyl isothiocyanate, with significantly higher values in the red radish paocai. Among the sulfides, dimethyl disulfide and dimethyl trisulfide were uniquely found in the red radish paocai, whereas 3-butenylisothiocyanate was found to be the unique sulfide in the cabbage paocai.

Among the various potential detected substances, those with a low content (lower than the threshold values) were excluded in the statistical analysis to narrow down the range of key volatiles. In detail, 36 potential volatile substances were screened out, which included 14 alcohols, 3 ketones, 5 esters, 10 hydrocarbons, 3 sulfides, and 1 phenol compounds.

To further screen the different VOCs in two kinds of paocai, a PLS-DA model was established, as shown in Figure 5A,B. Variable projection importance (VIP) was used to evaluate the contribution of variables. Substances with VIP value ≥1.0 were selected and combined with the conditions with a significant difference between the two groups (*p* < 0.05). The point distance between the VOCs of paocai fermented with two different raw materials was obvious, which indicated that the VOCs of paocai fermented with different raw materials were different. The total contribution rate of the first two principal components on variable X was 52.7% (R2X = 0.527), and the contribution rate on variable Y was 95.7% (R2Y = 0.957). At the same time, the prediction ability based on cross-validation was 94.0% (Q2 = 0.94), indicating that the model had a good explanatory degree. The stability and reliability of this model were verified by the cross-validation model with 200 permutation validations. As shown in Figure 5B, all Q2 values on the left were lower than the original point on the right, and the regression line of point Q2 intersected with the vertical axis below zero, and the intercept of −0.306 was less than 0.05, indicating that this PLS-DA model did not exhibit over-fitting and had good representativeness. A total of 13 compounds were screened out as the difference substances and included compounds like (−)-4-terpineol, α-terpineol, nerol, 4-isopropylbenzyl alcohol, carveol, α-cadinol, terpinolene, isobutyl phenylacetate, 1-isothiocyanatobutane, dimethyl disulfide, 1,4-dithiane, dimethyl trisulfide and 2,4-ditert-butylphenol.

Further to explore the correlation between microorganisms and VOCs, a Spearman correlation analysis was carried out with a relative abundance exceeding 5% and the nine biomarkers previously identified. It can be seen from Figure 5C that the microorganisms have different effects on the VOCs of paocai fermented by two kinds of vegetables. Nerol was also significantly correlated with six kinds of microorganisms in the red radish paocai to different degrees. In the cabbage paocai, the correlations between nerol, α-cadinol, isobutyl phenylacetate and the microorganisms were consistent, that is, they have a significant positive correlation with seven species and a significant negative correlation with three species at the same time, indicating that these 13 microorganisms have the same regulatory effect on the changes of these three substances. In addition, 1-isothiocyanatobutane, 1,4 dithiane and dimethyl trisulfide, which account for a high proportion of compounds in the red radish paocai, were related to six kinds of microorganisms, especially *L. paraplantarum* species. However, 2,4-ditert-butylphenol in the cabbage paocai and α-cadinol in the red radish paocai were not related to microorganisms.

From the perspective of the relationship between microorganisms and VOCs, *F. siliginis*, *L. parafarraginis*, *L. namurensis*, *L. zymae* and *S. paracollinoides* were found to be the key microorganisms for distinguishing the two kinds of paocai, and 4-isopropylbenzyl alcohol, α-cadinol, terpinolene and isobutyl phenylacetate were the key VOCs for distinguishing the two kinds of paocai made from red radish and cabbage.

**Figure 5.** (**A**) PLS-DA score scatter plot and (**B**) permutation tests of VOCs. (**C**) Heatmap of correlation between microbial species (relative abundance exceeding 5% and the nine biomarkers previously identified) and different volatile flavor substances in the two kinds of paocai. R indicates red radish paocai, B indicates cabbage paocai. \* 0.01 < *p* ≤ 0.05, \*\* 0.001 < *p* ≤ 0.01, \*\*\* *p* ≤ 0.001. **Figure 5.** (**A**) PLS-DA score scatter plot and (**B**) permutation tests of VOCs. (**C**) Heatmap of correlation between microbial species (relative abundance exceeding 5% and the nine biomarkers previously identified) and different volatile flavor substances in the two kinds of paocai. R indicates red radish paocai, B indicates cabbage paocai. \* 0.01 < *p* ≤ 0.05, \*\* 0.001 < *p* ≤ 0.01, \*\*\* *p* ≤ 0.001.

#### **4. Conclusions**

vor profiles of the studied paocais.

**4. Conclusions**

In this paper, metagenomic sequencing and a HS–SPME method based on a GC–MS method were used to systematically study the changes of microbial flora and VOCs in paocai fermented with different vegetables in six rounds of fermentation, combined with a multivariate statistical method to explore the intrinsic associations of "raw materials microorganism—VOCs". The results revealed that after multiple rounds of fermentation, a higher number of microbial species were present in the cabbage paocai than in the red radish paocai with a brine of the same age. At the species level, *S. paracollinoides* and *F. siliginis* were found to be the characteristic microbial species of the red radish paocai, as well as 13 species belonging to the Firmicutes phylum (namely, *L. parafarraginis*, *L. buchneri, L. zymae*, *L. namurensis*, *L. diolivorans*, *L. acidipiscis*, *L. casei*, *L. parabuchneri*, *L. hilgardii*, In this paper, metagenomic sequencing and a HS–SPME method based on a GC–MS method were used to systematically study the changes of microbial flora and VOCs in paocai fermented with different vegetables in six rounds of fermentation, combined with a multivariate statistical method to explore the intrinsic associations of "raw materials microorganism—VOCs". The results revealed that after multiple rounds of fermentation, a higher number of microbial species were present in the cabbage paocai than in the red radish paocai with a brine of the same age. At the species level, *S. paracollinoides* and *F. siliginis* were found to be the characteristic microbial species of the red radish paocai, as well as 13 species belonging to the Firmicutes phylum (namely, *L. parafarraginis*, *L. buchneri*, *L. zymae*, *L. namurensis*, *L. diolivorans*, *L. acidipiscis*, *L. casei*, *L. parabuchneri*, *L. hilgardii*, *L. farraginis*, *P. selangorensis*, *L. sunkii* and *L. otakiensis*) and two species belonging to Proteobacteria (namely, *K. ascorbata* and *P. carotovorum*) for the cabbage paocai. By constructing a PLS-DA model, it was found that there were 13 kinds of different VOCs in the paocai from two raw materials. In addition, the correlation between the microorganisms and

*L. farraginis*, *P. selangorensis*, *L. sunkii* and *L. otakiensis*) and two species belonging to Pro-

ing a PLS-DA model, it was found that there were 13 kinds of different VOCs in the paocai from two raw materials. In addition, the correlation between the microorganisms and VOCs in the cabbage paocai was stronger than that in the red radish paocai, and the most significant differences were found to be in 4-isopropylbenzyl alcohol, α-cadinol, terpinolene and isobutyl phenylacetate. The metagenomic sequencing in this study allowed the detection, identification and screening of the microbial ecology as well as the characteristic microorganisms of two different fermentation systems of the vegetable paocais. Further, the study also effectively allowed the identification of the VOCs and the establishment of their relationship to the microorganisms, which provided insights into the microbial flaVOCs in the cabbage paocai was stronger than that in the red radish paocai, and the most significant differences were found to be in 4-isopropylbenzyl alcohol, α-cadinol, terpinolene and isobutyl phenylacetate. The metagenomic sequencing in this study allowed the detection, identification and screening of the microbial ecology as well as the characteristic microorganisms of two different fermentation systems of the vegetable paocais. Further, the study also effectively allowed the identification of the VOCs and the establishment of their relationship to the microorganisms, which provided insights into the microbial flavor profiles of the studied paocais.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/foods11010062/s1, Table S1: Statistics profiles of the metagenome, Table S2: Content of volatile organic compounds in the two kinds of paocai.

**Author Contributions:** Conceptualization, L.J. and S.X.; methodology, L.J.; data curation, S.X.; resources, X.L.; writing—original draft preparation, L.J., S.X., G.S. and Z.Z.; visualization, X.H.; supervision, A.C.; project administration, A.C.; funding acquisition, A.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Research on Key Technologies of Food Processing of Characteristic Biological Resources in Arid Valley of Southwest China (No. 2017YFC0505106-3).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in NCBI database at accession PRJNA731323. deposited in https://www.ncbi.nlm.nih.gov/sra/PRJNA731323 (accessed on 21 December 2020).

**Acknowledgments:** We thank Fei Xu and Sichuan Liji Lebao Food Co., Ltd. for their technical support.

**Conflicts of Interest:** The authors declare no competing financial interest.

#### **References**

## *Article* **Bacterial Community of Grana Padano PDO Cheese and Generical Hard Cheeses: DNA Metabarcoding and DNA Metafingerprinting Analysis to Assess Similarities and Differences**

**Miriam Zago 1 , Lia Rossetti 1 , Tommaso Bardelli 2 , Domenico Carminati 1 , Nelson Nazzicari 1 and Giorgio Giraffa 1,\***

**Abstract:** The microbiota of Protected Designation of Origin (PDO) cheeses plays an essential role in defining their quality and typicity and could be applied to protect these products from counterfeiting. To study the possible role of cheese microbiota in distinguishing Grana Padano (GP) cheese from generical hard cheeses (HC), the microbial structure of 119 GP cheese samples was studied by DNA metabarcoding and DNA metafingerprinting and compared with 49 samples of generical hard cheeses taken from retail. DNA metabarcoding highlighted the presence, as dominant taxa, of *Lacticaseibacillus rhamnosus*, *Lactobacillus helveticus*, *Streptococcus thermophilus*, *Limosilactobacillus fermentum*, *Lactobacillus delbrueckii, Lactobacillus* spp., and *Lactococcus* spp. in both GP cheese and HC. Differential multivariate statistical analysis of metataxonomic and metafingerprinting data highlighted significant differences in the Shannon index, bacterial composition, and species abundance within both dominant and subdominant taxa between the two cheese groups. A supervised Neural Network (NN) classification tool, trained by metagenotypic data, was implemented, allowing to correctly classify GP cheese and HC samples. Further implementation and validation to increase the robustness and improve the predictive capacity of the NN classifier will be needed. Nonetheless, the proposed tool opens interesting perspectives in helping protection and valorization of GP and other PDO cheeses.

**Keywords:** Grana Padano cheese; generical hard cheeses; bacterial diversity; DNA metabarcoding; DNA (meta)fingerprinting; predictive models; neural network

#### **1. Introduction**

Protected Designation of Origin (PDO) products represent the excellence of European agricultural food production and are the result of the interplay between environmental (e.g., climate) and human factors (e.g., production techniques handed down over time), which are typical of a given territory. To this regard, many cheeses benefit from the PDO quality label. PDO cheeses are subject to specific production conditions, which producers spontaneously adhere to by joining the Consortia, which establish the transformation criteria through specific shared rules and guidelines. In this way, the consumer is guaranteed in terms of transparency and traceability, while benefiting from high-quality products [1]. The overall quality of cheeses, including PDO cheeses, is the result of many concomitant factors, such as the quality of the raw material, the farming methods, and the processing technology, which, in some cheeses (such as Grana Padano, Parmigiano Reggiano, Silter, Pélardon, Poro), involves the use of undefined microbial cultures. The interaction between these elements contributes to shape the qualitative and quantitative microbiological content

**Citation:** Zago, M.; Rossetti, L.; Bardelli, T.; Carminati, D.; Nazzicari, N.; Giraffa, G. Bacterial Community of Grana Padano PDO Cheese and Generical Hard Cheeses: DNA Metabarcoding and DNA Metafingerprinting Analysis to Assess Similarities and Differences. *Foods* **2021**, *10*, 1826. https://doi.org/ 10.3390/foods10081826

Academic Editors: Valentina Bernini and Juliano De Dea Lindner

Received: 30 June 2021 Accepted: 6 August 2021 Published: 7 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of the ripened product, giving it a specific composition [2–8]. On the other hand, the increasing frequency of imitations and frauds involving PDO foods, including cheeses [9], addressed the search for finer and more sensitive methods of discrimination.

In 2020, 210,000 tons (or approximately 5,200,000 cheeses) of Grana Padano (GP) were produced. With a growing export trend (40% of total production, +3.3% compared to 2019), GP is increasingly consumed both in EU and non-EU countries (https://www. granapadano.it, accessed on 5 July 2021) [10]. Grana Padano is a cooked, long ripened PDO cheese made from raw and partially skimmed milk, which is fermented by lactic bacteria present in natural starters ('sieroinnesto'). The sieroinnesto, which is one of the typical elements of Grana Padano, is prepared by backslopping, i.e., using part of the drained whey from the previous day's cheese making, which is left for 18–24 h at 42–45 ◦C until a final acidity of approx. 60 ◦SH/100 (pH 3.3–3.6) is reached [2,11,12]. On the other hand, GP is also one of the most imitated and counterfeited cheeses. This has stimulated the search for increasingly sensitive analytical approaches to GP mapping and characterization. A recent investigation carried out on a limited number of cheese samples using untargeted metabolomics revealed differences in chemical fingerprints between PDO and non-PDO grana-like cheeses [13]. In other studies, geographic or technological differences were observed in the microbiota of cheese, suggesting it as a possible tool to establish cheese authenticity and diversity [14–17]. In a previous study, the bacterial taxa present in GP were highlighted by 16S rRNA gene sequencing (DNA metabarcoding) [18]. In the present study, the usefulness of the cheese microbiota to distinguish GP cheese from generical hard cheeses, i.e., cheeses whose appearance could be confused with Grana Padano PDO cheese, was investigated. The structure and the genotypic fingerprinting of the bacterial taxa of 119 GP samples were evaluated by 16S rRNA gene sequencing (DNA metabarcoding) and RAPD-PCR from total cheese DNA (RAPD-PCR metafingerprinting), respectively, and compared with 49 samples of generical hard cheeses retrieved from retail.

#### **2. Materials and Methods**

#### *2.1. Sampling of Cheeses*

One hundred nineteen samples (2-kg slice each), representative of all producers of GP and the entire geographical area of production, were collected. Cheeses had been produced in August 2018 and had a ripening time of 6–7 months. Forty-nine samples of generical hard cheeses (from now on 'HC'), labelled as EU and extra-EU products, were retrieved from retail. All cheese samples were stored at −20 ◦C. After thawing at 4 ◦C for 18 h, slices (10 g each) including three different sections of the cheese (i.e., outer, central, and inner) were sampled and put into sterile containers.

#### *2.2. Total DNA Extraction*

The 10-g cheese slices were treated as described by Zago et al. [18]. Briefly, the samples were homogenized twice in a Stomacher 400 Circulator (Seward Laboratory, London, UK) with sterile sodium citrate (NaCt 2% *w*/*v*, pH 7.5). The cheese homogenate was centrifuged (14,400× *g* for 7 min at 4 ◦C) for fat removing. The pellet was then resuspended in Triton X-100 (2.5% *v*/*v*), washed twice with phosphate-buffered saline (pH 7.5), and centrifuged (8700× *g*, 7 min, 4 ◦C). Finally, after the pellet was resuspended in Tris-EDTA (0.1 M, pH 8), total dsDNA was extracted using a QIAcube HT automated station (Qiagen, Milan, Italy) using QIAamp 96 QIAcube HT kit (Qiagen, Milano, Italy). Total dsDNA was quantified fluorometrically (QubitTM, Life Technologies, Monza, Italy).

#### *2.3. DNA Metabarcoding*

Total DNA extracted from the 119 GP and the 49 HC samples was subjected to DNA metabarcoding analysis (IGATech, Udine, Italy) by sequencing of the variable V3–V4 regions of the 16S rRNA gene using an Illumina MiSeq platform, as described previously [18].

#### *2.4. RAPD-PCR Metafingerprinting*

RAPD-PCR from total DNA of each sample was carried out according to Zago et al. [18]. PCR products were separated by QIAxcel electrophoresis using dedicated DNA Screening Gel Cartridges (Qiagen, Milan, Italy). Two QX Alignment Markers (15 bp– 5 kb and 15 bp–600 bp, in a 1:1 ratio) and the QX DNA Size Marker (100 bp–2.5 kb) were included on each run. The repeatability and reproducibility of this method was evaluated by repeated amplification and analysis of four different DNA cheese samples, with both primers for every analysis performed.

#### *2.5. Data Analysis and Bioinformatics Processing*

Reads were de-multiplexed based on the Illumina indexing system. Following the QIIME pipelines, the USEARCH algorithm (version 8.1.1756, 32-bit) allowed the following steps: chimera filtering; grouping of replicate sequences; sorting sequences per decreasing abundance; and OTU identification, with a species-level taxonomic resolution. When the taxonomy assignment did not reach the species level, the genus or family name were reported. After removing OTUs <5 reads [19], alpha (α) diversity (richness and Shannon indexes), beta (β) diversity (principal component analysis-PCA and principal coordinate analysis-PCoA), and the rarefaction curves were estimated, on the resulting OTU table, by means of R (http://www.r-project.org/index.html; accessed on 7 April 2021), using "vegan" [20] and "agricolae" [21]. Statistical differences (*p* value ≤ 0.05) and evaluation of the influence of the two cheese samplings on the microbial indexes were evaluated by ANOVA followed by the Tukey HSD test. Relative abundance for each OTU across all cheeses was calculated, and "subdominant" and "dominant" OTUs were discriminated according to Zago et al. [18]. Taxonomic analysis was carried out through "reshape2" and "ggplot2" packages [22,23]. RAPD-PCR profiles were imported and analyzed by BioNumericsTM (version 7.6, Applied Maths, Sint-Martens-Latem, Belgium), as previously described [18].

A Neural Network (NN), trained to discriminate GP cheese and HC samples, was implemented with the aim to assess the feasibility of automatic classification for cheese samples and to compare the discriminatory power of metabarcoding and metafingerprinting analyses. The NN was implemented from scratch using Keras python library [24]. The input layer received PCoA data and contained 64 fully connected ReLU units; the second layer had ReLU 32 nodes; the third layer had 16. The last layer had only one sigmoid unit and performed the final binary classification. NN were trained on 80% of the available data, picked as a random stratification (thus maintaining the GP cheese/HC ratio).

Adam optimizer was used to minimize "binary cross-entropy" loss function. The following other performance metrics were measured: Area Under Curve (AUC), true and false positive counts (TP, FP), and true and false negatives count (TN, FN). We also derived Precision, True Positive and True Negative Rates (TPR, TNR), and Accuracy and Binary accuracy. All reported performance metrics were averaged over 10 repetitions of the training process, so as to avoid possible biases due to the random selection of validation set samples. Given the unbalance of tally classes, class weights were computed via the 'compute\_class\_weight' function from the sklearn [25] package and passed to the optimizer.

As a post-hoc investigation, we decided to repeat the NN training tuning class weights so as to explore the limits of the detection power of a hypothetical automated NN-based screening system. Specifically, it measured the possibility to raise True Positive Rate (TPR), i.e., the fraction of HC that were correctly identified, without significantly lowering the True Negative Rate (TNR), i.e., the fraction of GP cheese that were correctly classified. Given that, by definition, there is a tradeoff between TPR and TNR, we fixed the target thresholds to have both metrics higher than 0.9.

**3. Results**

#### **3. Results** reads per sample on average (range 80,636–805,340). A total of 477 OTUs, of which 130

*3.1. DNA Metabarcoding Analysis*

#### *3.1. DNA Metabarcoding Analysis* were further split into 48 dominant (≥ 1% total reads) and 82 subdominant (0.1–1% total

*Foods* **2021**, *10*, x FOR PEER REVIEW 4 of 10

the target thresholds to have both metrics higher than 0.9.

lowering the True Negative Rate (TNR), i.e., the fraction of GP cheese that were correctly classified. Given that, by definition, there is a tradeoff between TPR and TNR, we fixed

cheese, or HC) were analyzed. Overall, 35,842,968 reads were sequenced, with 218,554

One hundred and nineteen GP cheese samples and 49 samples of generical hard cheeses (i.e., cheeses whose appearance could be confused with Grana Padano PDO cheese, or HC) were analyzed. Overall, 35,842,968 reads were sequenced, with 218,554 reads per sample on average (range 80,636–805,340). A total of 477 OTUs, of which 130 were further split into 48 dominant (≥1% total reads) and 82 subdominant (0.1–1% total reads) taxa, were identified (data not shown). reads) taxa, were identified (data not shown). 3.1.1. Species Abundance According to the aim of this study, the samples were divided into GP cheese and HC. Considering the relative abundance of the dominant bacterial species (48 taxa; ≥ 1% total reads) found both in GP and in HC, *Lacticaseibacillus rhamnosus* was the prevalent species,

#### 3.1.1. Species Abundance followed by *Lactobacillus helveticus, Streptococcus thermophilus*, *Limosilactobacillus fermen-*

According to the aim of this study, the samples were divided into GP cheese and HC. Considering the relative abundance of the dominant bacterial species (48 taxa; ≥1% total reads) found both in GP and in HC, *Lacticaseibacillus rhamnosus* was the prevalent species, followed by *Lactobacillus helveticus, Streptococcus thermophilus*, *Limosilactobacillus fermentum*, *Lactobacillus delbrueckii, Lactobacillus* spp., and *Lactococcus* spp., with average values between 2% and 46% in GP cheese samples. Notably, the dominant species were less abundant in HC samples compared to GP cheese (Figure 1). In the cheese microbiota, many contaminating bacteria deriving from raw milk or the process environment, such as potentially pathogenic lactic acid bacteria (LAB) (i.e., *Streptococcus uberis* and *Lactococcus garviae*) and other non-LAB taxa (i.e., *Micrococcus*, *Staphylococcus*, *Acinetobacter*, *Pseudomonas,* and many enterobacteria), were also detected (Figure 1). *tum*, *Lactobacillus delbrueckii, Lactobacillus* spp., and *Lactococcus* spp., with average values between 2% and 46% in GP cheese samples. Notably, the dominant species were less abundant in HC samples compared to GP cheese *(*Figure 1). In the cheese microbiota, many contaminating bacteria deriving from raw milk or the process environment, such as potentially pathogenic lactic acid bacteria (LAB) (i.e., *Streptococcus uberis* and *Lactococcus garviae*) and other non-LAB taxa (i.e., *Micrococcus*, *Staphylococcus*, *Acinetobacter*, *Pseudomonas,* and many enterobacteria), were also detected (Figure 1). A different distribution between GP cheese and HC samples of the 82 subdominant species (0.1–1% total reads) was also observed (Figure S1).

**Figure 1.** Average values of relative abundance of the 48 dominant taxa retrieved in Grana Padano (GP) and similar hard cheese (HC) samples. **Figure 1.** Average values of relative abundance of the 48 dominant taxa retrieved in Grana Padano (GP) and similar hard cheese (HC) samples.

3.1.2. Alpha and Beta Diversity A different distribution between GP cheese and HC samples of the 82 subdominant species (0.1–1% total reads) was also observed (Figure S1).

3.1.2. Alpha and Beta Diversity

Richness (i.e., the number of different species in each sample) and Shannon indexes are reported (Figure 2). A greater uniformity of species was retrieved in GP cheese, compared to HC, as stated by the statistical significance of the Shannon diversity index (d.f. = Richness (i.e., the number of different species in each sample) and Shannon indexes are reported (Figure 2). A greater uniformity of species was retrieved in GP cheese, compared to HC, as stated by the statistical significance of the Shannon diversity index (d.f. = 1, F = 16.89, *p* < 0.001), whereas there was not a significant difference in the number of species

1

1

as indicated by Richness (d.f. = 1, F = 2.02, *p* = 0.157). The PCA showed no differences between GP cheese and HC, which were clustered together (Figure 3a). Conversely, the PCoA reported some differences between GP cheese and HC, as the latter were isolated on both sides of the panel (d.f. = 1, F = 4.70, *p* < 0.001; Figure 3b).

**Figure 2.** (**a**) Shannon diversity and (**b**) Richness of Grana Padano cheese (GP) and generical hard cheese (HC) samples. Different letters indicate significant (*p* < 0.05) differences between samplings.

**Figure 3.** (**a**) Principal component analysis (PCA) and (**b**) Principal coordinate analysis (PCoA) based on operational taxonomic unit (OTU) relative abundance of Grana Padano cheese (GP, red) and generical hard cheeses (HC, green). The first component (horizontal) of PCA accounts for 28.14% of variance, and the second component (vertical) of PCA accounts for 23.43% of variance.

#### *3.2. RAPD-PCR Metafingerprinting Analysis*

(**a**) (**b**) After RAPD-PCR amplification of the total microbial DNA extracted from all the cheeses, the resulting patterns were analyzed to detect sample-associated profiles and/or specific bands, with a repeatability and reproducibility of 70%. The presence of a panel of shared bands (labelled I in the left side of Figure S2) was highlighted. It included five series of common bands, present in more than 65% of the whole set of samples (GP cheese and HC). Moreover, a panel of unshared bands (or single bands), marked II in the right side of Figure S2, was observed. The matrix obtained from the data was analyzed by PCA, and statistically significant differences were observed between samples (d.f. = 1, F = 3.49, *p* < 0.001), which allowed to clearly separate the GP cheese and HC at the top and bottom of the panel, respectively (Figure 4a). The same data set was subjected to PCoA, and this trend was confirmed (d.f. = 1, F = 18.43, *p* < 0.001; Figure 4b).

1

(**a**) (**b**)

**Figure 4.** (**a**) Principal component analysis (PCA) and (**b**) Principal coordinate analysis (PCoA) based on band matching and Pearson correlation analysis of the metafingerprinting data of Grana Padano cheese (GP, red) and generical hard cheeses (HC, green). The first component (horizontal) of PCA accounts for 7.99% of variance and the second component (vertical) of PCA accounts for 7.31% of variance.

#### *3.3. Implementation of a Classifier*

DNA metabarcoding and RAPD-PCR metafingerprinting PCoA data were used to implement a classifier that could provide a fast, reliable, and automatic categorization of cheese samples into GP cheese or HC. Results for several performance metrics, measured on a randomly picked 20% of the data acting as validation set, are reported in Table 1. The same metrics measured on the training set are reported in Table S1. In general, the metafingerprinting dataset performed better than the metabarcoding dataset, e.g., for Area Under Curve (0.934 vs. 0.657, respectively) and Accuracy (0.928 vs. 0.637, respectively). Metafingerprinting showed a very low average level of False Positives (0.262), indicating that, on average, on about four out of five runs, no False Positives were detected, and in the fifth run, only one GP cheese was erroneously classified as HC.

**Table 1.** Performance statistics on the validation set for a Neural Network classifier trained on either DNA metabarcoding or RAPD-PCR metafingerprinting PCoA data. The table reports statistics averaged over 10 rounds of training on an 80/20 stratified split of training and validation set. The rightmost column reports metrics when the default class weight for class GP cheese is adjusted by a multiplicative coefficient equal to 0.25, chosen so that both TPR and TNF resulted above a threshold of 0.9.

FN: number of False Negatives. FP: number of False Positives. TN: True Negatives. TP: True Positives. AUC: Area Under Curve. Loss: final value of loss function, optimized during the training phase. Precision: TP/(TP+FP). TPR, True Positive Rate: TP/(TP+FN). TNR, True Negative Rate: TN/(TN+FP). Accuracy: (TP+TN)/(TP+TN+FP+FN). Balanced Accuracy: (TPR+TNR)/2.

As a post-hoc analysis, we decided to repeat the NN training on the metafingerprinting PCoA dataset using class weights more skewed toward the detection of HC samples, with the declared target of having both TPR and TNR metrics above the threshold of 0.9. We tested several configurations of weights (data not reported) and found that using an adjustment coefficient equal to 0.25 on the GP cheese class weight, the desired performance was achieved. On the DNA metabarcoding PCoA dataset, no weight adjustment was found to be able to bring both TPR and TNR above the 0.9 target threshold.

#### **4. Discussion**

Differences in microbial profiling could be useful to establish which microbial components may be responsible for the authentication of PDO cheeses [16,17,26]. On this basis, the microbial diversity of a large sampling that included the entire production area of the GP cheese was investigated and compared with that of 49 generical hard cheeses (HC) taken from retail. Two metataxonomic methods were applied, i.e., 16S rRNA gene sequencing and metagenotyping using RAPD-PCR from total cheese DNA, to identify and implement a robust classifier to differentiate GP cheese from HC. Metataxonomic data revealed a total of about 477 OTUs, including GP cheese and HC, but an overall greater relative abundance in the former ones. This trend occurred among both dominant (among which *Lcb. rhamnosus*, *L. helveticus*, *S. thermophilus*, *Limosilactobacillus fermentum*, *L. delbrueckii, Lactobacillus* spp., and *Lactococcus* spp. prevailed) and subdominant species. While detected LAB taxa within both dominant and subdominant taxa belong to the typical microbiota generally recovered from hard cheeses, the finding of residual psychrotrophic bacteria is not uncommon. Indeed, raw milk for Grana Padano is normally kept under refrigerated conditions before collection and subsequent transport to processing sites [18]. The presence of bacterial DNA of potentially pathogenic LAB and non-LAB taxa is not uncommon, but it has no safety significance, as these bacteria are inactivated by the combined effects of cooking temperature, the curd acidification during the early phases of cheesemaking, and the long ripening times under the harsh conditions (low moisture, low activity water) that characterize the production of hard cheeses [27–29]. The Shannon index, which accounts for both microbial richness and evenness, was significantly influenced by the two samplings, highlighting a greater uniformity of species in GP cheese but a similar richness. Differential analysis based on Bray–Curtis dissimilarity between all samples and calculated using species relative abundancies delineated a significant, although incomplete, separation between the two groups.

Cheese production conditions within a PDO area were selected for specific microbial populations. The microbiota of cheeses, especially those (such as Grana Padano) obtained with artisanal processes and from raw milk, plays a fundamental role in defining the qualitative characteristics and safety parameters of the final products. Its composition, structure, and modulation are the result of different selective pressures, e.g., microbiological quality of milk, technology, type of starter used, and processing environment [30]. Therefore, differences in one or more of these factors can be decisive in shaping specific microbiological profiles in cheese. Grana Padano cheese is obtained from raw milk produced in a large production area, which includes most of the provinces of northern Italy included in the Po Valley. Although animal feed is substantially similar, the microbiological composition of milk aimed at GP cheese production can be influenced by management practices at the farm level and by seasonal, climatic, and environmental variations. For example, the presence of some LAB species, such as *Lactiplantibacillus plantarum, Lentilactobacillus parabuchneri, Lentilactobacillus parafarraginis, Lentilactobacillus hilgardii,* used also as silage starters, can be related to corn silage fed to cows producing milk for GP [31]. The microbial content of raw milk is subsequently modulated by the selective action of the technology and the practice of using undefined whey starter cultures [32,33]. The greater microbial heterogeneity of raw milk for GP, coming from a large production area, and the 'balancing' action exerted by the application of a very similar technology and the addition to raw milk of whey starter cultures, which are usually prepared with comparable methods among the different

dairies, are critical to explain the higher species abundance and uniformity observed in GP cheese. On the other hand, HC not being subject to the constraints of a PDO can be produced from milk outside the production area (with its specific microbial content) and with technologies that, although substantially similar, often rely on the use of selected starter cultures and/or the application of thermization, pasteurization, bactofugation, or microfiltration of milk [34]. These selective pressures could be decisive in explaining the lower OTU abundance and species uniformity in HC. The above trends were corroborated, and the separation between the two groups of samples, which was highlighted by the Bray–Curtis dissimilarity test, was noted.

DNA fingerprinting methods, such as RAPD-PCR, are often applied to evaluate the endemicity, or the prevalence, of a given strain. Furthermore, this technique was found useful in discriminating soil microbial communities and estimating their relatedness [35,36]. In a previous study, RAPD-PCR proved useful to highlight a pattern of bands present in all the samples, as well as more specific bands, which aggregated groups of samples, or distinguished single samples, within the microbial community of GP [11]. In the present work, RAPD-PCR was applied to fingerprint the overall bacterial community of GP cheese and HC samples. The obtained metagenotypes were evaluated as possible tools to differentiate the two sampling groups, assuming, unlike metataxonomic analysis, that this technique was able to identify strain- or group-specific differences within complex microbial communities. Data processing of RAPD-PCR profiles using PCA and Bray–Curtis dissimilarity analysis allowed to obtain a clear and statistically very significant separation between GP cheese and HC samples. Metataxonomic and meta-fingerprinting data were then used as inputs to train and validate a two-class (GP cheese vs. HC) classifier based on a neural network as computational model. While metataxonomic data did not allow for reliable classification, the discriminatory power of metafingerprinting enabled to build an extremely robust model (very high binary accuracy). When trained by metagenotyping data, the model correctly classified GP and HC samples. The origin of this differentiation is not currently known, although it is likely that strain-specific peculiarities in the microbial community, determined by previously outlined ecological, geographic, or technological selective pressures, might explain it. The molecular (meta)fingerprint of the entire microbial community could be promising to assist to authenticate GP cheese and to distinguish it from imitation products as part of an attempt to hinder any counterfeits. Further validation will be needed to increase robustness of the classifier and confirm, with unknown samples, its discriminating ability.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/foods10081826/s1, Figure S1: Average values of relative abundance of the 82 subdominant taxa retrieved in Grana Padano (GP) samples and similar hard cheeses (HC) samples; Figure S2: Band matching cluster analysis showing the relationship between Grana Padano and similar hard cheese samples and M13-/BOXA1R-RAPD-PCR bands; Table S1: Performance statistics on the training set for a Neural Network classifier trained on either metabarcoding or fingerprinting PCoA data.

**Author Contributions:** Conceptualization, M.Z. and G.G.; methodology, M.Z. and L.R.; software, M.Z. and N.N.; validation, N.N.; formal analysis, L.R. and T.B.; investigation, L.R. and T.B.; writing original draft preparation, M.Z.; N.N., and G.G.; writing—review and editing, M.Z.; D.C., and G.G.; supervision, M.Z. and G.G.; funding acquisition, G.G. and D.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Agricultural, Food, and Forestry Policies (MIPAAF) as part of the Project entitled "New Technologies for Cheese Production (NEWTECH)" (D.M. 0021568 of 28 July 2017).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Detection of Enteric Viruses and Core Microbiome Analysis in Artisanal Colonial Salami-Type Dry-Fermented Sausages from Santa Catarina, Brazil**

**Roberto Degenhardt 1,2, Doris Sobral Marques Souza 1,3, Leidiane A. Acordi Menezes 1,4 , Gilberto Vinícius de Melo Pereira 5 , David Rodríguez-Lázaro 6,7 , Gislaine Fongaro 3 and Juliano De Dea Lindner 1,\***

**Abstract:** Microbial fermentation plays an important role in the manufacturing of artisanal sausages and can have major effects on product quality and safety. We used metagenomics and culturedependent methods to study the presence of Hepatitis E virus (HEV) and Rotavirus-A (RV-A), and fungal and bacterial communities, in artisanal Colonial salami-type dry-fermented sausages in Santa Catarina state, Brazil. Lactic acid bacteria (LAB) and yeast dominated the microbiome. *Latilactobacillus sakei* and *Debaryomyces hansenii* were ubiquitous and the most abundant species. The DNA of some foodborne pathogens was found in very low concentrations although viable cells of most of these species were undetectable by cultivation methods. The characteristics of the raw material and hygiene of the artisanal sausage manufacturing process resulted in high loads of beneficial microorganisms and the absence of HEV and RV-A viruses as determined by RT-qPCR assays. In conclusion, high LAB load in sausages was more relevant to preventing pathogen growth than the ripening time and/or physicochemical characteristics. However, the presence of *Clostridium* spp. and other pathogens in some samples must be taken into account for the development of future preservation methods; appropriate LAB starter cultures and health surveillance are required in the production process to prevent foodborne outbreaks.

**Keywords:** swine and pork production chain; Hepatitis E virus; Rotavirus-A; metagenomic analysis; food safety

## **1. Introduction**

Fermented sausages (dry and semi-dry) make up a substantial proportion of pork produced and consumed in western Santa Catarina State. This is in part a consequence of the European origin and cultural characteristics of the population in this region in the south of Brazil. Several small and medium-scale production companies are dedicated to the transformation of pork meat, employing artisanal production techniques including natural fermentation. However, these production practices can allow the survival and mul-

**Citation:** Degenhardt, R.; Sobral Marques Souza, D.; Acordi Menezes, L.A.; de Melo Pereira, G.V.; Rodríguez-Lázaro, D.; Fongaro, G.; De Dea Lindner, J. Detection of Enteric Viruses and Core Microbiome Analysis in Artisanal Colonial Salami-Type Dry-Fermented Sausages from Santa Catarina, Brazil. *Foods* **2021**, *10*, 1957. https:// doi.org/10.3390/foods10081957

Academic Editor: Arun K. Bhunia

Received: 16 July 2021 Accepted: 18 August 2021 Published: 22 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

tiplication of spoilage and pathogenic microorganisms if hygiene and sanitary conditions and process parameters are not strictly controlled.

There is a great diversity of dry-fermented sausages produced in the world. Different raw meat types, formulations, additives, processes, casing, and drying periods are used [1]. Colonial sausages produced in Santa Catarina State are characterized by low acidity and high moisture content when they are offered for sale, but they continue the ripening process outdoors in the market, gradually losing moisture [2]. These uncontrolled conditions can lead to a risk for the multiplication of undesirable microorganisms.

The production of artisanal sausages involve fermentation by the indigenous microbiota of meat (composed mainly by lactic acid bacteria—LAB) producing organic acids, and consequently a decrease in the pH that allows preservation of the product. *Lactobacillus*, *Pediococcus*, *Enterococcus*, and *Leuconostoc* are believed to be the principal genera involved in sausage fermentation [3,4], and diverse LAB are found in food products according to the manufacturing practices used. The preservation ability of LAB is based on competition for nutrients and the production of antimicrobial metabolites, hydrogen peroxide and bacteriocins. The selective influence of diverse factors (e.g., pH, redox potential, water activity (aw), temperature, relative humidity, oxygen availability) in fermented sausage will contribute to determining the microbial consortium present.

The microbial communities in fermented sausage play significant roles in the flavor, texture, quality, and safety of sausages [5,6]. Understanding the ecology of fermented sausages is clearly key to understanding the physical and chemical changes that occur during fermentation and ripening [7]. High-throughput sequencing of the microbiota from salami-type dry-fermented sausages of different origins indicates that LAB and coagulase-negative cocci, including both micrococci and coagulase-negative staphylococci (CNS) are the predominant groups in spontaneously fermented artisanal sausages [8]. *Debaryomyces, Penicillium*, and *Aspergillus* are the yeast and mold genera most commonly found [6,9–11]. Studies involving metagenomics for this type of fermented product are scarce in Brazil, and therefore, the ecology of Brazilian artisanal dry-fermented sausages needs further investigation.

Hepatitis E virus (HEV) and Rotavirus—A (RV-A) belong to the human enteric viruses group. They are non-enveloped and icosahedral viruses, and are very persistent in the environment [12]. HEV has a positive-sense single-stranded RNA and RV-A a doublestranded and segmented RNA (11 segments) that facilitates reassortment events among species [12]. Both are zoonotic viruses that infect many animal species as well as humans. They are mainly transmitted by the fecal-oral route, although HEV is also transmitted parenterally (rare), or by the ingestion of contaminated meat from viremic animals [12,13].

HEV and RV-A infect domestic pigs, which serve as reservoirs, and are emerging risks for human infection [13,14]. These viruses are excreted at high concentrations in animal stools (mainly from swine) so animal waste can contaminate water supplies, and when used in crops may contaminate food. Indeed, they persist as infectious particles in the environment, e.g., water and sewage, for several weeks at lower than body temperatures [13,15–17]. HEV circulation has already been reported in humans and domestic swine in Brazil [17–19].

The present study aimed to detect and quantify HEV and RV-A in artisanal Colonial salami-type dry-fermented sausages. We also analyzed the presence of fermenting, spoiling, and pathogenic microorganisms in the sausages. Both culture-dependent and -independent (metagenomic) approaches were used.

#### **2. Materials and Methods**

#### *2.1. Colonial Salami-Type Dry-Fermented Sausage Production*

Colonial salami-type dry-fermented sausages were collected from 13 production sites (Table 1) under the state inspection system (SIE) in the Vale do Rio do Peixe zone, Santa Catarina state, Brazil (ca. 27◦100 S; 51◦300 W). The sausages were prepared using the same manufacturing process in all sites. Pigs were slaughtered and the carcasses underwent temperature equalization for up to 24 h in cooling chambers to establish rigor mortis. The carcasses were exposed and shanks and palettes were used to prepare the sausages. The meat was comminuted in a meat grinder and seasoned with NaCl, curing agent (NO2), and spices.

**Table 1.** Origin of Colonial salami-type dry-fermented sausages in Santa Catarina State, ripening time and physical-chemical characteristics.

WPS: water phase salt; aw: water activity; N: natural; S: synthetic.

The production of these sausages does not involve addition of commercial starter cultures for pan-curing. The fermentation is due to the natural microbiota of the meat and environment. After homogenization in a mixer, the meat is extruded into a natural bovine or a synthetic cellulose casing. The sausages are placed in a smoker where the drying is started, and the curing is completed under smoking. After smoking, the sausages are ready for sale. Thereafter, drying took place in the market, subject to the environmental conditions of the market stall site. The 13 sausage samples used had drying times of between one and 40 days, with half of the samples having been dried for more than seven days.

#### *2.2. Physical-Chemical Evaluation*

The physical-chemical characteristics of the sausages (Table 1) were determined using the official methods of the Ministry of Agriculture, Livestock and Supply (MAPA) [20]. The pH was determined using a potentiometer (Quimis, model Q400A, Diadema, SP, Brazil). The moisture content was determined by a gravimetric method, placing the samples in a drying oven (SP Labor, model SP-100/150, SP, Brazil) at 105 ◦C for 24 h [21]. The NaCl content was determined by the mercury metric method [20]. The aw was determined by grinding the samples and placing them into a digital hygrometer at 25 ◦C (Aqualab Model-Series 3 TE, Decagon Devices Inc., Pullman, USA). Water phase salt (WPS) was calculated from the values of NaCl and moisture content, according to Koral and Köse [22], using the following equation:

$$\text{WPS\%} = \text{(NaCl\%} / (\text{NaCl\%} + \text{moisture\%})) \times 100$$

#### *2.3. Microbiological Analyses*

The samples were tested for several microbial groups: total coliforms (TC), thermotolerant coliforms (TTC), coagulase-positive staphylococci (CPS), sulfite-reducing clostridia (SRC), LAB, Enterococci, *Salmonella* and *Listeria monocytogenes.* TC were counted by spread plating on Violet Red Bile Lactose agar (VRBL, Oxoid Limited, Hampshire, United Kingdom), incubated aerobically at 37 ◦C for 24 h [23]; TTC by pour plating on double-layered VRBL agar, incubated aerobically at 44 ◦C for 24 h; CPS by spread plating on Baird-Parker agar (BPA, Neogen, Lansing, MI, USA), incubated aerobically at 35 ◦C for 24–48 h [24]; SRC

by pour plating on Iron Sulfite agar (ISA, Neogen), incubated anaerobically at 37 ◦C for 24–48 h [25]; LAB by pour plating on De Man, Rogosa and Sharpe agar (MRS, Neogen), incubated anaerobically at 30 ◦C for 72 h [26]; Enterococci by spread plating on m-Enterococcus agar (Neogen), incubated anaerobically at 35 ± 2 ◦C for 24–48 h; *Salmonella* by incubating a 25 g pre-enrichment sample in buffered peptone water (BPW, Neogen) at 34–38 ◦C for 18 h, selective enrichment by adding 0.1 mL of the BPW to modified semi-solid Rappaport-Vassiliadis medium (MSRV, Neogen) at 41.5 ± 1 ◦C for 24 h and 1 mL of BPW to Muller-Kauffmann Tetrathionate-Novobiocin broth (MKTTn, Neogen) at 37 ± 1 ◦C for 24 h, plating-out in xylose lysine deoxycholate agar (XLD, Neogen) at 37 ◦C for 24 h and brilliant green agar (BGA, Neogen) at 37 ◦C for 24 h, with confirmation by biochemical and serological testing [27]; *Listeria* spp. by incubating a 25 g pre-enrichment sample in demi Fraser broth (DF, Neogen) at 30 ◦C for 24 h, selective enrichment by incubating 0.1 mL of the DF in Fraser broth (Neogen) at 35 ± 1 ◦C for 24–48 h, then plating-out on Listeria agar (ALOA, Laborclin, Pinhais, Paraná, Brazil) at 37 ◦C for 24 h, with confirmation by biochemical testing [28].

Yeasts and molds were counted by spread plating on dichloran glycerol agar (DG18, Neogen) containing chloramphenicol (0.1 g/L) and incubating at 25 ◦C for 5–7 days [29]. Macroscopic observation was used to distinguish between yeasts and molds among the colonies growing on the plates.

#### *2.4. Viral Detection by RT-qPCR*

Virus was detected by mechanical disruption of tissues followed by silica-membrane RNA extraction, as described by García et al. [14] and Rodríguez-Lázaro et al. [30]. Briefly, 25 g samples of sausage (Table 1) were collected, chopped using a sterile scalpel blade and homogenized; aliquots of 0.25 g were taken from the samples, and they were inoculated with ~105 PFU (Plaque Forming Units) of MNV-1 (Murine Norovirus-1) as sample process control virus (SPCV). Each sample (including SPCV) was homogenized in 1 mL of lysis buffer RLT (Qiagen, Hilden, Germany) containing 0.14M β-mercapto-ethanol in a tissue homogenizer, and incubated for 5 min at room temperature. The samples were then centrifuged at 10,000× *g* for 20 min at 4 ◦C and the supernatants (800 µL) were transferred to new tubes for nucleic acid extraction with the RNeasy® Mini kit (Qiagen, Germany), following the manufacturer's protocols. The resulting nucleic acid samples were subjected to RT-qPCR, both undiluted and diluted tenfold to reduce the inhibitory effects of any contaminants in the samples.

The TaqMan technique was used for RT-qPCR, as previously described by Baert et al. [31] for MNV-1 as SPCV, Jothikumar et al. [32] for HEV and [33] for RV-A. All amplifications were performed in a StepOne Plus® Real-Time PCR System (Applied Biosystems, Waltham, MA, USA), and each sample was analyzed in triplicate. Ultrapure water was used as the non-template control for each assay.

#### *2.5. Metagenomic Analyses*

DNA was extracted from samples with the DNeasy Power Soil Pro kit (Qiagen, Germany), following the manufacturer's instructions. Total DNA extracted was used as the template for Next Generation Sequencing (NGS), on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA).

The microbiome of the samples was studied by high throughput sequencing for bacterial taxonomic identification based on conserved and variable regions V3/V4 of the 16S rRNA gene with primers 341F (CCTACGGGRSGCAGCAG) and 806R (GGAC-TACHVGGGTWTCTAAT) [34,35]. Extracted DNA was subjected to PCR using a previously described protocol [36] developed by Neo-prospecta Microbiome Technologies (Florianopolis, Brazil). The first PCR used sequences based on Illumina's TruSeq adapters, which allows a second PCR with primers with indexed sequences. PCRs were carried out in triplicate using Taq Platinum (Invitrogen, Waltham, MA, EUA) with the following conditions: PCR1—95 ◦C/5 min, 25 cycles (95 ◦C/45 s, 55 ◦C/30 s, 72 ◦C/45 s) and final extension

72 ◦C/2 min; PCR2—95 ◦C/5 min, 10 cycles (95 ◦C/45 s, 66 ◦C/30 s, 72 ◦C/45 s) and final extension 72 ◦C/2 min.

The community of filamentous fungi and yeasts was studied by amplifying the ITS1 region using ITS1F (GAACCWGCGGARGGATCA) and ITS2R (GCTGCGTTCTTCATC-GATGC) primers [37], and the same conditions for PCR1. For PCR2, the conditions were: 95 ◦C/5 min, 15 cycles (95 ◦C/45 s, 66 ◦C/30 s, 72 ◦C/45 s) and final extension 72 ◦C/2 min. The final reaction products were purified using a protocol developed by Neoprospecta involving magnetic beads. The samples were grouped in libraries quantified by Real-Time PCR using KAPA Library Quantification (KAPA Biosystems, Woburn, MA, USA). Library pools were sequenced in a MiSeq sequencer system (Illumina Inc., USA) using V2- 300 cycle kit, in a paired-end run, without normalization of libraries. The raw data (DNA sequences in fastq files) were analyzed through bioinformatics workflow considering accumulated error in the sequencing to be at most 1%. To identify the species of microorganisms present in the samples, the DNA sequences obtained were compared with a proprietary database containing well-characterized DNA sequences.

Sequencing data for each sample were processed with the Quantitative Insights into Microbial Ecology (*Qiime*) software package [38]. The sequencing output was analyzed by a read quality filter; reads with an average Phred score < 20 were removed and then 100% identical reads were clustered. Clusters with fewer than five reads were excluded from further analysis, to remove putative chimeric sequences. The remaining good-quality sequences were further clustered at 97% similarity to define operational taxonomic units (OTU). OTUs were classified by comparison with a custom 16S rRNA/ITS1 database (NEORefDB, Neoprospecta). Sequences with at least 99% identity to the reference database were taxonomically assigned. Each OTU was given an arbitrary unit corresponding to the relative abundance (r.a.) calculated as the number of reads pertaining to each species as a percentage of the total number of sequences read for each sample. Species with a r.a. below 1% were grouped in the category "others". The r.a. values were compared between OTUs using heatmaps prepared on *Qiime*.

#### *2.6. Statistical Analyses*

Differences in virus parameters between samples collected from different sites and seasons were analyzed by nonparametric analysis of variance (Kruskal–Wallis) performed with Statistic 7.0. Pearson correlation and linear regression tests, ANOVA, and Student's *t*-tests were performed using GraphPad Prism 5.0 (USA). Differences were considered statistically significant at a *p*-value ≤ 0.05.

#### **3. Results and Discussion**

#### *3.1. Physical-Chemical Characteristics*

Selected physical-chemical characteristics (pH, moisture, NaCl, WPS and aw) of the sausage samples are summarized in Table 1. The sausages can be classified into two subgroups according to the criteria of Brazilian legislation [39]: (i) dry sausages, with moisture ≤ 40% and aw ≤ 0.92, and (ii) semi-dry sausages, for which these parameters are not stipulated. A particular feature of these locally produced sausages is that they are offered for sale as semi-dry sausages and, due to the environmental conditions, they start to present dry sausage characteristics as the ripening continues. The dry sausages (L01, L04, L07, L10, L12) had undergone seven or more days of ripening, whereas the semi-dry sausages (L02, L03, L05, L06, L08, L09, and L11) had undergone up to six, except for samples L02 and L06; these two sausages had been maintained under refrigeration at the site of sale, reducing dehydration.

All the samples presented pH between 5.1 and 6.4 (Table 1), characteristic of the sausages produced in the region. This weak acidity is also characteristic of Italian sausages, produced in the Mediterranean region only with pork meat. Low pH and aw are the two main obstacles to spoilage and the development of pathogenic microorganisms in fermented sausages [40]. If these factors are weak, as is the case in our samples due to the

relatively high pH of the samples (≥5.1), the association between aw and the competitive microbiota (mainly LAB) becomes very important to the safety of the product. The WPS is not frequently used to assess the potential for multiplication of spoilage and/or pathogenic microorganisms in meat sausages in Brazil, but it is a potentially valuable indicator for small producers, due to its relationship with water activity.

The minimum salt concentration in the samples was 3.01 (sample L08) and the maximum was 6.15 (sample L04). The upper limit of sensorial consumer acceptability for salt in cured meats is 3.5% to 5.0% with 1.5% to 2.0% being optimal. The salt concentration was therefore acceptable in the samples, except for sample L04. The higher concentration of salt in this case appears to be a consequence of excessive drying: the moisture content was only 30.2%. In artisanal cured products, salt is still the main resource used for conservation, and this is the main reason for the high salt content in these sausages, mostly between 3.7% and 5.04% (Table 1). These values are compatible with those required for the preservation of sausages and related products [41].

#### *3.2. Viral Detection*

The samples of artisanal Colonial salami-type dry-fermented sausages from different production sites (Table 1) were tested by RT-qPCR for HEV and RV-A. Neither virus was detected in any of the sausage samples, and the mean viral extraction efficiency (measured as SPCV -MNV-1 recovery) was 30%.

HEV is more frequently detected in feces than serum and liver samples from infected animals [42], but there is nevertheless a risk of foodborne hepatitis E infection from consumption of pork products. Boxman et al. [43] found HEV RNA in 14.6% of pork sausages collected from a Dutch market from 2017 to 2019, and HEV RNA was detected in 18.5% of salami (*n* = 92) samples; sequence-based typing was successful for 33 samples, and all were genotype 3c. Berto et al. [44] found HEV RNA in four samples of liver sausage in France, and three of them had infectious particles of HEV as assessed by cell culture assay. Colson et al. [45] found HEV RNA of genotype 3 in seven samples of raw figatelli purchased in French supermarkets. However, HEV RNA was not detected in any of the 22 sausages from different production sites analyzed in this work. This contrasts with Souza et al. [17], who detected HEV in effluent after swine manure digestion by psychrophilic anaerobic biodigesters (PABs) in Concórdia city, in Santa Catarina state. This demonstrates the circulation of HEV in the pig-finishing farms in the zone where the Colonial sausage we studied was produced.

Lowering pH during meat fermentation is believed to inhibit the growth of certain pathogens. Wolff et al. [46] tested the pH resistance of a genotype 3 strain of HEV. Only minimal inactivating effects were found at pH conditions common in sausages during curing (4.5 to 6.5), simulated using D/L-lactic acid. High salt concentrations are also used to preserve meat products and to inactivate foodborne pathogens. Wolff et al. [47] tested the effects of NaCl, sodium nitrite, and sodium nitrate concentrations on the infectivity of HEV. Conditions consistent with those in fermented sausages were simulated. Treatments with up to 20% NaCl for 24 h at 23 ◦C, with and without the addition of 0.015% nitrite or 0.03% nitrate, did not inactivate the virus, demonstrating that HEV is highly stable at salt concentrations used to preserve raw meat products. Therefore, the acid pH and high salinity of the products we tested were unlikely to be responsible for the absence of HEV.

Organs from naturally HEV-infected pigs were tested by García et al. [14] and the HEV virus was not detected in loin samples: meat even from infected animals may therefore be only a low risk for HEV transmission. The main ingredient used in the preparation of Colonial salami is the loin, and this may explain the absence of HEV RNA from the samples tested, despite evidence of the circulation of HEV in pig farms in the region [17]. Other organs would need to be analyzed to assess the HEV contamination of the pigs in the region.

RV-A RNA was not detected in the Colonial sausage samples analyzed, evidencing the good sanitary conditions during slaughter and sausages preparation. RV-A zoonotic transmission to humans usually occurs via direct contact with infected animals or fomites [15]. RV-A causes diarrheic diseases in swine and meat contamination could occur in the slaughterhouse or during sausage preparation by handlers and fomites. RV-A is the most prevalent rotaviral cause of diarrhea in swine [48]. Other RV species that can infect pigs were not investigated in this work.

#### *3.3. Microbiological and Metagenomic Analyses*

The combined culture-dependent and -independent approaches are best for analyses of the microbiota of fermented foods [49]. Culture-dependent analysis showed high levels of LAB and yeasts and the presence of enterococci and molds in all samples. There were low or undetectable counts of undesirable microbiota, in accordance with Brazilian standards (Table 2). In general, the samples showed similar cell count values that were independent of the ripening time. Likewise, there was no correlation between ripening time and microbiota in the metagenomic analyses. The pH in all samples was higher than 5.0, indicating that LAB community was mostly bacteria with low acid production capacity, resulting in low acidity fermented sausages. Although low acidity does not favor the inhibition of undesirable microbiota, it enhances the water-retention capacity of meat proteins, contributing to the maintenance of moisture during ripening [50]. The yeast community is important for the development of flavor and inhibition of toxigenic molds.

\* cfu·log−1 ; \*\* Presence (P) or absence (A) in 25 g of sausage; TC: total coliforms; TTC: thermotolerant coliforms; CPS: coagulase-positive staphylococci; SRC: sulfite-reducing clostridia; LAB: lactic acid bacteria.

> The 16S amplicon target sequencing identified more than 60 species each with a r.a. above 0.1% (Figure 1 shows the 31 species with over 1% r.a., in one or more samples). Species with a very small number of reads (r.a. < 1%) were classified into the category "others".

> LAB and CNS were the most abundant groups in our samples (Figure 1), consistent with previous reports for dry sausages from Spain and Italy [6,8]. LAB was predominant in 11 of the 13 sausage samples: *Latilactobacillus sakei* was the most abundant taxon in samples L01, L03, L05, L06, L08, L10, L12, and L13 (r.a. of 43.56%, 65.85%, 49.57%, 38.38%, 42.70%, 62.75%, 40.30%, and 92.14%, respectively); *Latilactobacillus curvatus* was predominant in sample L02 (r.a. 60.31%); *Lactococcus lactis* in sample L04 (r.a. 33.33%), and *Pediococcus pentosaceus* in sample L11 (r.a. 36.39%).

**Figure 1.** Relative abundance (%) of bacterial species in artisanal Colonial salami-type dry-fermented sausages. Species with a very small number of reads (relative abundance < 1%) were classified into the category "others".

> Most samples (Figure 1) contained three or four species with r.a. above 5%, although in sample L13 *Lb. sakei* was massively predominant and the microbial community was not very diverse. *Lb. curvatus* was found at r.a. of between 5% and 12% in samples L01, L07, L09, and L12, *Lactococcus lactis* at between 9% and 14% in samples L05 and L12, *Lactococcus garvieae* at around 13% in samples L08 and L011, *Leuconostoc mesenteroides* at 17% in L08, *Lactiplantibacillus plantarum* at 7% in L02 and *Weissella hellenica* at 18% in L10. *Lb. sakei, Lb. curvatus* and *Lp. plantarum* are the species of LAB commonly found in sausages with no added starters. Other species of *Weissella, Leuconostoc*, *Lactococcus*, and *Pediococcus* were also found, but as part of the subdominant microbiota in our samples. The fermentation in these sausages is spontaneous, so the microorganisms detected come from the meat and the surrounding environment [6].

> The manufacturing process of artisanal Colonial salami-type dry-fermented sausages from Santa Catarina includes spontaneous fermentation, cure and drying as conservation. Thus, the obstacles to the development of spoilage and pathogenic microbiota are the addition of preservatives (NaCl and sodium nitrite), pH, aw, and the competitive microbiota, mostly LAB. The protective effects of LAB are due in part to their production of bacteriocins with antimicrobial activity against pathogens [51]. Hebert et al. [52] showed that *Lb. curvatus* CRL705 isolated from Argentinean artisanal fermented sausage produces bacteriocin with anti-listeria activity. Palavecino et al. [6] describe coagulase-negative staphylococci (CNS) with antimicrobial activity. The presence of CNS (*Staphylococcus saprophyticus*, *Staphylococcus galinarum*, and *Staphylococcus equorum*) in our samples (Figure 1) may be one of the reasons for the absence of pathogens common in meat products, such as *Clostridium botulinum*, *Staphylococcus aureus*, *Listeria monocytogenes*, and *Escherichia coli* [53].

> In general, the sausage samples (Table 2) that showed the best hygiene indicator results (TC, TTC, CPS, SRC, and *Salmonella*) were drier (e.g., L10) and had longer ripening time (e.g., L02). This is in contrast to the observation for the high LAB semi-dry sausage subgroup, which corroborates the metagenomic data (Figure 1). *Listeria* spp. were detected in 38.4%

of the samples, but *L. monocytogenes* was not detected at a significant r.a. in any sample (Figure 1). *Listeria* spp. are detected in many fermented sausages with low acidity and high aw [54–56]. *Listeria* spp. counts determined by culture-dependent methods are generally low in fermented sausages [57,58], and culture-dependent methods employ enrichment and cell concentration steps. Therefore, the possibility of isolating a microorganism is greater than for the metagenomic analysis when it is present at low r.a.

In addition to LAB, species of the *Micrococcaceae* family, mainly belonging to the genus *Staphylococcus*, are commonly found in artisanal sausages. *S. saprophyticus, Staphylococcus xylosus*, *S. equorum*, and *Staphylococcus carnosus* are often predominant. *Staphylococcus succinus* is also frequently observed in sub-dominant populations [7,59]. For samples L07 and L09 (Figure 1), the predominant microorganism was *S. saprophyticus* (r.a. 42.60% and 61.65%, respectively). This species was also found in samples L01, L06, L08, L10, and L12 (r.a. from 15% to 26%). *S. succinus* was found in sample L09 at an r.a. off 20%.

*Brochothrix thermosphacta, Carnobacterium* spp., clostridia, and *Pseudomonas* spp. are among the bacteria most frequently causing spoilage of refrigerated pork meat [60]. *Pseudomonas fragi* (Figure 1) was found in samples L01, L03, L05, L06, and L011 with r.a. between 10% and 23%. *Carnobacterium divergens* (r.a. 1.87%) and *Clostridium perfringens* (r.a. 27.76%) were found in sample L04, in which the abundance of lactobacilli was low. *B. thermosphacta* was detected in samples (Figure 1) with a low abundance of lactobacilli (L03, L05, and L11). However, for sample L09, *S. saprophyticus* was predominant, and although the r.a. of lactobacilli were low, no sequences belonging to *B. thermosphacta*, *C. divergens* (r.a. 1.87%), *Pseudomonas* spp. or clostridia were detected.

Taxonomical assignment of the yeasts (Figure 2) showed the dominance of *Debaryomyces hansenii* in L02, L04, L06, L07, L09, L10, L11, and L12 (r.a. 63.49%, 44.25%, 82.44%, 80.01%, 65.71%, 35.20%, 28.92% and 10.85%, respectively) followed by *Candida zeylanoides,* predominant in samples L01, L03, and L12 (r.a. 73.91%, 19.13%, and 58.23%, respectively), *Yarrowia galli* predominant in L05 (r.a. 42.15%), and *Kodamaea ohmeri* in L08 (r.a. 53.31%). *D. hansenii*, *C. zeylanoides* and *Yarrowia lipolytica* were found in all the samples, although they were less abundant than the other yeast groups. *Y. galli, Pichia kudriavzevii, K. ohmeri, Candida parapsilosis* and *Candida tropicalis* were also found with relevant r.a. in the subdominant populations of the samples.

*D. hansenii* is the most common yeast in fermented sausages and contributes to two important characteristics: production of volatile compounds during ripening, especially esters, branched alcohols, and aldehydes [61,62], and bioprotective antifungal activity against toxigenic penicillia [63]. Species belonging to the genera *Candida, Yarrowia*, and *Pichia* are also frequently found [6,51,64].

Development of mold on the surface of the casings during ripening of fermented sausages is common, and molds are present in the environment at the production and ripening sites. Molds are important in the development of flavor, but some fungi of the genera *Aspergillus* and *Penicillium* can produce undesirable mycotoxins [65]. For example, a combination of high moisture, low temperature and the predominance of penicillia has been found to be conducive to the production of penicillic acid, a carcinogen [66]. The filamentous fungus found at the highest abundance was *Aspergillus cibarius* in L02, L11 and L12 (3.38%, 42.51%, and 35.30%, respectively) followed by *Penicillium thymicola* (r.a. 8.78% in L12), and *Geotrichum candidum* (r.a. 4.90% in L04). *Aspergillus* and *Penicillium* are the main genera found in fermented sausages, while *Geotrichum* is found with lower frequency [6,51,67–69].