# A retrospective multivariate case-control approach to investigate factors influencing the sanitary quality of raw milk cheeses

Alisia Chiado Rana, Géraldine Cazeau, Stéphanie Bougeard, Kevin Bethune, Vincent Chochois, Corinne Teyssier, Charlotte Valat

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Alisia Chiado Rana, Géraldine Cazeau, Stéphanie Bougeard, Kevin Bethune, Vincent Chochois, et al. A retrospective multivariate case-control approach to investigate factors influencing the sanitary quality of raw milk cheeses. Applied Food Research, 2024, 4 (2), pp.100617. 10.1016/j.afres.2024.100617. anses-04816878

HAL Id: anses-04816878

https://anses.hal.science/anses-04816878v1

Submitted on 3 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# A retrospective multivariate case-control approach to investigate factors influencing the sanitary quality of raw milk cheeses

Alisia Chiado Rana<sup>a</sup>, Géraldine Cazeau<sup>b</sup>, Stéphanie Bougeard<sup>c</sup>, Kevin Bethune<sup>d</sup>, Vincent Chochois<sup>d</sup>, Corinne Teyssier<sup>a</sup>, Charlotte Valat<sup>b,*</sup>

## ARTICLE INFO

Keywords:

16S rRNA metabarcoding

Short reads

Multi-block analysis

Raw milk cheese

Production environment

Risk factors

Food-borne zoonotic pathogens

## ABSTRACT

A retrospective case control study conducted in France between 2022 and 2023, involving 20 farms (10 case and 10 control farms), investigated the factors influencing the sanitary quality of raw milk cheeses, focusing on the agricultural aspect and the bacterial microbiota characterised using 16S rRNA short reads metabarcoding from cheese production and the environment during manufacturing. The microbiological quality of raw milk cheeses is defined in Europe by Regulation (EC) No. 2073/2005. In this study, a set of food safety criteria related to the presence of food-borne pathogens such as *Salmonella*, *Listeria monocytogenes* and virulence factors of STEC in cheese for a long-time period (2017 - 2022) prior to sampling was used to select the case farms. Using multi-variate analysis, this study identified multiple bacterial species and critical sampling points such as raw milk, milking machines and wooden shelves that contribute to food-borne contamination of cheeses and highlights the importance of farm management practices such as additional vaccination or transhumance, hygiene and manufacturing practices to ensure the safety of raw milk cheeses. The bacterial species from raw milk cheeses or contact materials along the production chain contributing to the case or control farms were diversified in winter and summer. This study confirms the value of multi-block analysis for managing complex datasets and suggests that 16S rRNA markers from metabarcoding can help to identify critical control points and risk factors.

## 1. Introduction

In 2021, France was the second-largest cheese producer in the European Union (Eurostat, 2022). Part of this production involves cheeses produced from unpasteurised milk that hasn't been cooked over 40 °C, known as raw milk cheeses (FAO, 2022). In France, approximately 6000 artisanal dairy farms produce milk exclusively for the manufacturing of raw milk cheese (CNIEL, 2022) and raw milk cheeses account for 77.5% of PDO (Protected Designation of Origin) production (INAO, 2022). The production of raw milk cheeses is therefore one of the pillars of the French culinary heritage (Donnelly, 2013). The PDO label guarantees harmonised production and farming practices as well as the geographical localisation, covered by the Regulation (EU) No. 1151/2012 (European Commission, 2012). Cow's milk has been found to contain more than 100 genera and 400 species of microorganisms (Montel et al.,

2014). This complex and heterogeneous nature of the microbiota of cow milk and raw milk cheeses is crucial to the art of cheesemaking because it affects the development of flavours and organoleptic properties (Anastasiou et al., 2022; Quigley et al., 2013). Cheese making involves microorganisms naturally present in raw milk (non-starter culture) and additional starter or adjunct cultures produced at an industrial or artisanal level. This microbiota can play a number of roles; the most important of which is to facilitate milk fermentation and contribute to flavour and good texture, mainly involving lactic acid bacteria (LAB). Among LABs, those that are considered safe for humans mainly belong to the genera *Lactococcus*, *Oenococcus*, *Lactobacillus*, *Leuconostoc*, *Pediococcus*, and some *Streptococcus* (Balciunas et al., 2013). Some LABs also play a role in controlling the growth of bacteria (Mariani et al., 2011). Numerous factors can affect the microbiological composition of milk, and by extension, cheese. These factors include the cow's diet (Ponzoni

<footnote><sup>a</sup> Qualisud, Univ Montpellier, CIRAD, Avignon Université, Institut Agro, IRD, Univ de La Réunion, Montpellier, France</footnote>

<footnote><sup>b</sup> Unité Epidémiologie et Appui à la Surveillance, ANSES Laboratoire de Lyon - Université de Lyon, Lyon, France</footnote>

<footnote><sup>c</sup> Unité EPidémiologie, SAnté et Bien-Etre (EPISABE), ANSES Laboratoire Ploufragan-Plouzané-Niort, Ploufragan, France</footnote>

<footnote><sup>d</sup> CIRAD, UMR Qualisud, F-34398 Montpellier, France. Qualisud, Univ Montpellier, Avignon Université, CIRAD, Institut Agro, IRD, Université de La Réunion, Montpellier, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: charlotte.valat@anses.fr (C. Valat).</footnote>
et al., 2009; Pozzi et al., 2023), the mountain pasture, the animal microbiota (e.g. udder surface, teat canal, faeces), the livestock housing conditions or the cheese production environment such as water, contact materials (Cruciata et al., 2018; Sun & D'Amico, 2021), ripening shelves, milking machine (Cruciata et al., 2018; Didienne et al., 2012), air in the barn (Hohmann et al., 2020) and the season (Li et al., 2018). Previous study has analysed the combination of farming practices using regression trees, and breed and diet have been identified as the most influential factors for sensory and health dimensions in cheese assessment (Rey-Cadilhac et al., 2023). In France, we observe two main production systems depending on the season: the traditional mountain pasture system in summer and the indoor system with indoor feeding in winter. In addition, the composition and the structure of the raw milk cheese ecosystem depend on the cheese production. The microbiota of PDO uncooked pressed raw milk cheeses, mostly produced in France, includes LABs such as *Streptococcus thermophilus* and *Lactobacillus delbrueckii* ssp. *bulgaricus*, *Lactobacillus helveticus* and *Lactococcus lactis*, yeasts such as *Geotrichum candidum* and *Debaryomyces hansenii* (Ceugniez et al., 2017; Mariani, 2007), and ripening bacteria such as *Brevibacterium aurantiacum*, *Arthrobacter arilaitensis* and *Corynebacterium variabile* (Cogan et al., 2014; Monnet et al., 2016). The microbiota of raw milk cheeses can promote health (e.g. *Lactobacillus* and *Bifidobacterium*), cause spoilage (e.g. *Pseudomonas*, *Clostridium*, *Bacillus* and other spore-forming or thermoduric microorganisms) or cause human diseases (Oliver & Murinda, 2011; Yoon et al., 2016). Over the last decade in France, contaminated raw milk cheeses have been linked to 34 % of salmonellosis outbreaks, 37 % of listeriosis outbreaks, and 60 % of enterohemorrhagic *Escherichia coli* (EHEC) infections (ANSES, 2022). In fact, dairy cows are a major reservoir of bacterial pathogens, and the major zoonotic bacterial pathogens that can be transferred to raw milk include *Salmonella* spp., *Listeria monocytogenes*, *Staphylococcus aureus*, *Campylobacter* spp. and pathogenic *E. coli* (Farrokh et al., 2013). Although the growth of these bacterial pathogens can be impacted by antagonistic microorganisms resulting in a very low prevalence of these bacteria in raw milk cheese, some pathogens can persist in the cheese-making environment and in some rare cases contaminate cheese during production (Oliver et al., 2005). Therefore, the production of raw milk cheese requires careful attention to bacterial control and advanced hygiene practices to reduce the risk of bacterial pathogen contamination. The regulations governing milk quality and herd health in dairy farming are several, encompassing sanitary measures, production practices, and legislative frameworks. The Codex General Standard for Cheese (CODEX-STAN 283-1978) emphasises the importance of safety and hygiene in the production of cheese, including raw milk cheese. The microbiological quality of raw milk cheeses is defined in Europe by Regulation (EC) No. 2073/2005 (European Commission, 2005). The risk of contamination can be reduced by controlling for conditions that favour pathogen introduction and growth in farms and in the environment of production. Raw-milk cheese safety assessment is complicated by the fact that each farm has a distinct combination of environmental risk factors (e.g., composition and properties of raw milk, technology, seasonal variations, production and farm practices, environment). In the face of this public health challenge, a better understanding of the relationship between the microbial ecosystems present in the raw milk cheese and from the production environment is essential.

The microbial risk assessment of raw milk cheeses has been addressed previously (FSANZ, 2009). The main risk factors are the transmission of bacterial pathogens through fecal contamination, which could be influenced by farm management practices, including milking processes, cleaning and disinfection procedures udder care, together with milking machine cleanliness (temperature of washing detergent and the use of a water softener) (Elmoslemany et al., 2010; O'Mahony & Austin, 1991) and animal health (Ricci et al., 2013). Some other preventive and corrective measures are recommended such as pre-harvest interventions, including milking hygiene practices (Perrin et al., 2015) or probiotic use (El Jeni et al., 2024). Several contact surfaces during

Uncooked Pressed Cheese (UPC) manufacturing such as the milk filtration system (Bücher et al., 2024) or the wooden shelves support significant growth of bacterial species (Busetta et al., 2024; Settanni et al., 2021). In contrast, biological indicators from cheese that have a positive impact on human health have been identified (e.g. β-glucuronidase-positive *E. coli* and β-carotene) (Rey-Cadilhac et al., 2021). However, a number of gaps in the risk assessment have been identified. In particular, the knowledge of bacterial strains circulation at the farm level and during the cheese production, and the identification of the optimal hygiene and manufacturing practices (Farrokh et al., 2013). To fill these gaps, it is important to further investigate the bacterial composition of cheeses throughout the process, depending on both the farm and the cheesemaking practices. Researchers should focus on filling these gaps by further investigating the bacterial composition of cheeses throughout the process, depending on both the farm and the cheese-making practices. In particular, they should focus on non-starter lactic acid bacteria (NSLAB) (Betterta et al., 2023) and the bacterial composition of environmental surfaces (Mendonça et al., 2024).

The aim of this paper was to evaluate, through statistical analysis, potential risk or protective factors and critical control points in UPC food safety. For this purpose, a retrospective case-control study was carried out on 10 control and 10 case PDO small-scale farms from alpine mountain (France). This article focused on PDO cheeses sharing standard production methods and guarantees harmonised practices. To take into account the low prevalence of zoonotic bacterial pathogens in this raw milk cheese production, a set of multiple food safety criteria as been laid down in Regulation (EC) No. 2073/2005, related to the detection of *Salmonella*, *Listeria monocytogenes* and virulence factors of STEC in cheeses prior to sampling for a period of five and a half years (January 2017 - September 2022), was used to select the case farms. The selected case farms were associated with repeated contamination events, with a maximum number of detection events. In contrast, the control farms had no reported cases of contamination during the same period. Two campaigns were conducted in winter and summer to assess the impact of the seasonal variations on factors potentially influencing cheese safety. Multi-block statistical analysis (Bougeard et al., 2011) usually applied to multiple sets of variables, but rarely to bacterial variables, was adapted in our study in order to include quantitative bacterial variables characterised by 16S rRNA sequences metabarcoding (Castellanos-Rozo et al., 2020) and data on both farming and cheesemaking practices using a questionnaire. The blocks have been defined in a farm-to-fork approach in accordance with the production stage. The importance of the correlation of each farming, production and microbiological variable and its respective block with the case or control farms, was evaluated. An extensive on-farm sampling was carried out in this work to investigate the critical control points from the farm to the final raw milk cheese product.

## 2. Material and methods

### 2.1. Case-control study: farms selection

A retrospective case-control study was carried out in artisanal UPC with PDO producing farms, in the French mountains. The same retrospective case-control approach has previously been used to assess risk factors associated with *Salmonella* on the pig farms (Smith et al., 2018) or on broiler farms (Yamaguchi et al., 2024). In this study, 109 PDO producing farms were identified in 2022. Ten control farms and 10 case farms were subsequently included in this study. The control farms were randomly selected from among 23 farms meeting the cheese quality standard required by French legislation based on European Commission (EC) Regulation No. 2073/2005 by analysing the confidential background data on the microbiological criteria set by Regulation (EC) No. 2073/2005 from 2017 to 2021, prior to sampling. All of which were subsequently contacted by telephone and e-mail; 10 agreed to participate. Specific food safety criteria for cheese related with Regulation No.
2073/2005 (European Commission, 2005), such as the detection of the main serotypes of *E. coli* (O157:H7, O26:H11, O45:H2, O103:H2, O111:H8, O121:H19 and O145:H28), the major virulence of STEC (*stx* or *eae* genes), *Listeria monocytogenes* and *Salmonella* spp., were checked during the same five year period to select the case farms. The analytical methods are described in the following part. The 10 case farms with the highest number of detection events related to these food safety criteria and who agreed to participate were selected. In each case farm, the number of events ranged from three to six from 2017 to 2022. The expected low incidence of contamination events, required a long period of five years of retrospective analysis to attribute the sanitary status of farms.

As farm management practices vary according to the season, two different sampling campaigns were carried out, the first one during the winter of 2022 (mid-November to mid-December), when the cows were indoors and the second during the summer of 2023 (mid-June to mid-July), when the cows were outdoors in the pastures. For each campaign, a total of 40 visits were made, with two visits per farm. During the first visit to each farm, a questionnaire on farming and cheesemaking practices was completed with the farmer.

## 2.2. Sample collection

Each farm, in each sampling campaign, was visited twice. The first visit corresponded to the day of milking (day 0) during which samples and measurements of physico-chemical parameters (temperature and pH) were collected. These parameters were measured from raw milk with starters and rennet before coagulation throughout the entire cheese-making process, including the final cheese product. The second visit took place on day 7 for cheese curd and rind sampling. A total of 789 samples were collected along the cheesemaking process and from the environment (Table A1).

At day 0, fifteen sampling points were collected: liquid samples (raw milk, raw milk with starter, whey, brine, cheese smear), swabs (milking machine swab teat cup before washing, milking machine swab teat cup wipes after washing, milking machine teat filter swab of junction before washing, milking machine swab filter of junction after washing, moulding cloth, cheese washing brush), filters (milking machine filter, wooden shelves in the dryer, wooden shelves in the cellar) and solid pre-cheese sample (raw milk coagulum) (Table A1). On day 7, during the second visit, two solid cheese samples (cheese curd and rind) were collected. The liquid samples (RM BS, RM AS, RMC, W, B and CS) were collected and 40 mL of each was centrifuged, the whole resulting pellet was extracted (see DNA extraction section). Swabs and filters were preserved in 2 mL and in 90 mL of Ringer's solution, respectively (Castro et al., 2017; Rubiola et al., 2020). Sampling of the ripening wooden shelves (10 cm2), both in the dryer and in the cellar was carried out using a brushing method (Gaglio et al., 2022; Mariani, 2007). During the second campaign, in summer, an additional two litres of water from the cheese brushing machine was sampled on each farm and filtered on day 0; the resulting 0.22 µm cellulose filter was stored in 30 mL of Ringer's solution. Two 10 g aliquots of solid samples were aseptically collected and stored at -20 °C for subsequent analysis. All samples collected were cooled with ice packs and then stored at -20 °C (Castellote et al., 2015; Wolfe et al., 2014).

## 2.3. Analysis of bacterial pathogens

*Listeria monocytogenes*, *Listeria* spp. in milking machine filters (MMF) and raw milk coagulum (RMC), *Salmonella* spp. in raw milk after the addition of starter (RM_AS) and cheese (CC) were detected using the VIDAS® system (bioMérieux, France) according to the manufacturer's instructions. GENE-UP® STEC (bioMérieux, France) was used for the detection of STEC-encoding virulence genes in RMC. Analyses were performed after enrichment from 25 mL of raw milk or 25 g of RMC, respectively.

## 2.4. Metabarcoding

### 2.4.1. DNA extraction and purification

DNA was extracted from 2 mL of centrifuged pellet from each individual sample, except for swabs which have been pooled from four or two samples (Table A1). Environmental samples, including milking machine filters and dryer/cellar's ripening shelves (WSD, WSC), were mixed in a stomach for 4 min at 230 rpm, and centrifuged at 14,000 g for 10 min (Rubiola et al., 2022; Wang et al., 2019). The pellet was diluted in 500 µL UltraPure water and vortexed. Ten grams of cheese samples (CC or CR) or 20 g of coagulum (RMC) were homogenised in sterile PBS and mixed in a stomach for 4 min at 180 rpm (Ceugniez et al., 2017; Dreier et al., 2021; Zago et al., 2021). Liquid samples were centrifuged at 14,000 g for 15 min at 4 °C. The fat layer and the supernatant were discarded, leaving only the pellet, except for (RMF BS), where the supernatant was analysed separately as suggested in a previous protocol (Bagel & Sergentet, 2022; Lima et al., 2018; Stinson et al., 2021). Pellets were washed four times with PBS, vortexed and centrifuged. DNA from the whole pellet was extracted using the DNeasy® powerFood® Microbial Kit (Qiagen, France) according to the manufacturer's instructions, with an additional enzymatic lysis step (Siebert et al., 2021). The pellet of crude DNA was resuspended in 100 µL water. The nucleic acids were precipitated by adding 1 vol of ice-isopropanol and concentrated into DNA pellets by centrifugation (13,000 rpm, 5 min, 4 °C) and drying (Chemidlin Prévost-Bouré et al., 2021). Finally, DNA pellets were washed and dried as described by Chemidlin Prévost-Bouré et al. (2021) and the dried pellet was dissolved in 15 µL of UltraPure Water.

### 2.4.2. 16S rRNA sequences PCR library preparation

DNA concentration was measured using a Qubit (dsDNA Quantitation, High Sensitivity, Q32851 Thermo Fisher scientific, France) and adjusted to 5 ng/µL. Amplicon library preparation was performed according to the Illumina manufacturer's instructions. An initial PCR was performed using the 16S-specific primers 341F and 785R (Klindworth et al., 2013), complemented by 5' nucleotide overhangs used to anchor indexing primers during the subsequent indexing PCR. For each individual sample (n = 789), a 16S rRNA gene locus-specific PCR was performed in duplicate and pooled. The 20 µL PCR mix contained 10 µL of 2X Taq polymerase Master Phusion Flash High-Fidelity (Thermoscientific, Courtaboeuf, France), 0.25 µM of each primer, 20 ng template DNA and 1.25 µM of each peptide nucleic acid (mPNA and pPNA, PNA Bio Inc, Newbury Park, CA, USA) as described previously (Lundberg et al., 2013). A second short PCR (10 cycles) was performed to add Illumina barcodes after a magnetic bead purification (Clean PCR, Proteigene, Saint-Marcel, France). The final indexed PCR products were then magnetic beads purified and pooled. The resulting library was paired-end sequenced on an Illumina MiSeq sequencer using the MiSeq Reagent Kit v3 (600 cycle; Illumina, San Diego, CA, USA).

### 2.4.3. Sequencing data analysis

The sequencing data were trimmed using Cutadapt (v. 4.0) with the following parameters: "-n 3 -m 150:150 -quality-base=33 -g CCTACGGNGGCWGCGAG -a GGATTAGATACCBDGTAGTC -G GAC-TACHVGGGTATCTAATCC -A CTGCWGCCNCCCGTAGG". Subsequent steps were performed in R (v. 4.3.2). The read quality was checked and filtered using the DADA2 (v1.30.0) functions 'plotQualityProfile' and 'filterAndTrim', respectively. The following parameters were used for 'filterAndTrim': truncQ=2, maxEE=c(2,2), trimleft=c(0, 20), truncLen=c(240, 200), minLen=20, maxN=0, rm.phix=TRUE. For taxa identification, amplicon sequence variants (ASVs) were obtained using the 'dada' function from the DADA2, paired-end reads were merged, and chimeras were removed prior to taxonomy assignment, which was performed against the DairyDB (v. 3.0) database (Meola et al., 2019) using the 'assignTaxonomy' function with default parameters. After taxonomic assignment, the lowest taxonomic level identified for each ASV was
reported. Reports at the family (or higher) level meant ‘one or more
unidentified ASV belonging to this family’. Sample depth was then
rarefied to a uniform 10,165 reads per sample. In addition, a common
Shannon diversity index was calculated for samples from milking ma-
chine and slot swab filters (MMF; MMFS_BW_Sw), ripening shelves in the
dryer and cellar shelves (WSD; WSC) and cheese curd with rind (CC; CR)
by regrouping the ASVs from each sample.

To compare diversity between samples, alpha diversity was assessed.
The Shannon index (Shannon, 1948) was calculated on the rarefied data
using the phyloseq package (v. 1.46.0) (McMurdie and Holmes, 2014).

## 2.5. Statistical analyses

### 2.5.1. Differential bacterial abundance analysis

Differential bacterial abundance analysis was calculated from the
16S rRNA metabarcoding data using the DESeq2 R's package (v. 1.38.3),
with default parameters (Love et al., 2014) for all samples. Prior to the
DESeq2 analysis, the absolute number of reads across all samples was
filtered by retaining all species present in at least four samples, with at
least 20 reads in at least one sample. To perform the DESeq2 analysis, a
pseudo-count of one was added to all the read counts prior to rarefac-
tion. The outputs of the analysis were the bacterial variables, i.e. unique
ASVs associated with each sample type.

### 2.5.2. Variable selection with univariate analysis

Univariate statistical analysis was used to select variables for the
multi-block analysis, integrating both metabarcoding data and ques-
tionnaire responses. This approach was implemented to mitigate the
inclusion of extraneous variables, thereby reducing potential con-
founding factors and increasing the robustness of the multivariate
model.

Farm questionnaire data were compared between case (n = 10) and control (n = 10) farms using the Wilcoxon test for numerical variables and the Chi-squared Fisher exact test for categorical variables, with the null hypothesis was rejected at 20 % (p-value < 0.2). The Shannon index data of each sample were compared between case and control farms using a non-parametric Kruskal-Wallis test, due to the reduced number of observations (n = 20), which allowed the selection of significant variables (p-value < 0.2). This threshold was chosen to deal with the reduced number of farms and to increase the number of variables to be included in the multivariate analysis from 10 to 40.

For the bacterial variables, ASVs with an absolute log2 fold change greater than 2.5 and a Benjamini-Hochberg adjusted p-value (p-adj) less than 0.05 after DESeq analysis were selected. A lower p-value was set for these bacterial variables compared to the other variables in order to manage the large number of taxonomic data, control false positives, and identify differentially expressed ASVs with higher confidence.

### 2.5.3. Dealing with missing values

The MissMDA package (v 1.19) (Josse & Husson, 2016) was used to
deal with missing values. After the imputation process, a post-filtration
event was applied to reduce multicollinearity as follows: the same ASVs
from different samples with a significant correlation (Kendall's rank test;
p-value < 0.05) were excluded. An additional univariate analysis was
performed to ensure that the final set of variables remained statistically
significant and non-redundant.

### 2.5.4. Multiblock analysis

Multi-block analysis was used to: (i) quantify the relative contribu-
tion of selected explanatory variables to the sanitary quality of the
cheese represented by the case and control farms (outcome) and (ii)
examine the contributions of their associated blocks.

All significant variables in the univariate analysis, including both the
questionnaire data and the 16S rRNA sequences from metabarcoding,
were transformed where necessary and filtered prior to the multi-block
analysis. The absolute number of unrarefied reads for each ASVs

associated with the sample, was used as a numerical variable. In order to
integrate numerical and categorical variables in the same analysis,
categorical variables were transformed into dummy variables (coded as
0 or 1) to represent categorical variables that have two or more levels.
Variables with two categories were coded as 0 or 1. For variables with
more than two categories, multiple dummy variables were created to
represent each level of the variable, and only one dummy variable took
on a value of 1 for each observation. Values were centred at zero and
scaled to remove bias due to differences in units or amplitudes between
measures.

The selected variables were divided into four blocks, as described in Fig. 1. The block-1 (B1) included physicochemical (temperature and pH) and bacterial variables (Shannon's index, taxa abundance) from cheese and production environment samples. Block-2 (B2) included the same variables as in block-1 that were significant after univariate analysis, but from samples after the addition of starters. Block-3 (B3) was related to farm practices such as cattle health, udder hygiene, milking machine washing, farm structure, number of cattle up to raw milk collection. Block-4 (B4) was related to the cheesemaking process and including procedures for cheese production, starting from the production room until cheese ripening. The variable Y was the variable to be explained (case versus control farms), while the potential risk factors were distributed across the explanatory blocks (B1...B4). The multi-block partial least square (mbPLS) method was applied using the ade4 R package (v. 1.7) (Bougeard & Dray, 2018) to the specific case of a single response dataset Y (Wold, 1966) to deal with multi-correlations between explanatory variables within the blocks. The mbPLS method calculates components, as linear combinations of all explanatory variables for each block (t1...t4). According to the criterion to be optimised, these components have the highest covariance with the outcome variable Y, representing the sanitary status of farms (case or control). The global component tT sums up the partial components (t1...t4). Equal weights were assigned to each variable in each block (variable centering and scaling) and to the four blocks (block scaling with block inertia) (Bougeard et al., 2011).

The outputs of the multi-block method were the block importance (BlockImp) and the regression coefficient estimate (CE) to quantify the effect of each block and each explanatory variable on the outcome (sanitary status). A bootstrap simulation with 500 replicates was performed to calculate 95 % confidence intervals (CI) of these results. Blocks whose estimates 95 % CI did not include the standard threshold value of 1/K were considered statistically significant in explaining Y. Variables, whose estimates 95 % CI not including the standard threshold of zero, were considered statistically significant in explaining Y (Bougeard et al., 2011).

## 3. Results

### 3.1. Selection of variables of campaign 1 (winter) using univariate analysis

In total, 10,748,080 raw paired-end reads were generated from the 16S rRNA sequences of the whole samples of campaign-1. After quality filtering steps, the reads were reduced to 7491,448 and 69.70 % were further assigned and classified into 44,150 total ASVs (10,039 unique ASVs). Taxonomic analysis revealed a total of 30 bacterial phyla, 295 families and 9503 unique bacterial taxa across all samples from the first campaign. The same most abundant taxa were found in both control and case farms (Fig. 2). At the genus level, most of the sequences were associated with LABs (Streptococcus, Lactobacillus, Lactococcus, Aero-coccus and Leuconostoc) (Burgain et al., 2014) or additional taxa such as Corynebacterium, Intestinibacter, Acinetobacter, Brevibacterium, Kocuria, Staphylococcus, Facklamia, Rothia, Pseudomonas and Tetragenococcus (Fig. 2).

Univariate statistical analysis identified a total of 69 bacterial vari-
ables including 64 unique ASVs associated with bacterial taxa, for which
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Definition of the four blocks used for the multi-block analysis within the production flowchart of uncooked pressed cheese and distribution of the 213 variables selected from the univariate analysis. C1 and C2 are the number of variables selected by univariate analysis for each campaign in relation to the samples within the blocks. Blocks composition: Block 1 (B1): microbiological and physicochemical variables before starter addition; Block 2 (B2): microbiological and physicochemical variables after the addition of the starter cultures to cheese ripening; Block 3 (B3): farming practices; Block 4 (B4): cheesemaking practices. t = time; RH = relative humidity, T = temperature and pH recommended; * temperature and pH measured in product.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Boxplots of relative abundance of the top 15 most abundant taxa (at genus taxonomic rank) in samples from control versus case farms during campaign 1, in winter. Bars are coloured according to the taxonomic rank of the bacterial genus.</figcaption></figure>

the number of assigned reads differed significantly between samples from case and control farms. Of these 64 unique ASVs, 21 belonged to the most common genera (Fig. 2), in particular *Corynebacterium* and *Staphylococcus* followed by *Brevibacterium*, *Leuconostoc*, *Streptococcus* and *Pseudomonas*. Among the significant bacterial variables, 36 variables were related to samples before the addition of starter defined as block-1, and 33 ASVs were associated with samples after the addition of starter (cheesemaking process) defined as block-2 (Table 1). The alpha diversity of the milking machine filter slot taken before washing (MMF_MMFS_BW_Sw) was higher in control farms than in case farms (Table A2). The bacterial community of the ripening shelves (WSD_WSC) and the cheese curd and rind (CC_CR) was less diverse in the case farms compared to the control farms (Table A2). During campaign-1, a number of reads mapping to *Salmonella enterica* greater than 100 was found in 4 samples from control farms, 3 of which belonged to the RM BS sample type. In addition, *S. enterica* ASV was detected in two distinct types of samples from case farms. *Shigella* spp. (> 100 reads) was detected in one sample from case farms and in two distinct types of samples from control farms. No ASV corresponding to *Listeria* spp. was detected. However the presence of STEC, *Listeria* spp. and *Salmonella* spp. was not confirmed in the cheese by the conventional VIDAS and GENE UP methods. In campaign-1, of the 128 variables from the questionnaire, eight differed significantly between case and control farms (Tables 2 and 3). A total of 79 significant variables were included in the multi-block analysis (Table 1).

## 3.2. Selection of variables of campaign 2 (summer) using univariate analysis

Campaign-2 was carried out in summer when the dairy cows were on the mountain pastures. In total, 12,471,211 raw paired-end reads were generated from 16S rRNA sequencing for this second campaign. After quality filtering steps, the reads were reduced to 8954,861 and 71.8 % were further assigned and classified into 63,365 total ASVs (21,088 unique ASVs). Taxonomic analysis revealed the same most abundant taxa as were found in campaign-1 (Fig. 3). However, the relative abundance of genera and the alpha diversity were higher during the campaign-2.

Univariate statistical analysis comparing control and case farms in campaign-2 identified a total of 110 bacterial variables including 101 unique ASVs, 34 in block-1 and 76 in block-2 (Table 1). Seven of the 29 genera taxa associated with block-1 belonged to the top 15 genera taxa (Fig. 3). Of these 101 unique ASVs, 26 belonged to the most common genera (Fig. 2), in particular *Staphylococcus*, *Acinetobacter*, *Brevibacterium* followed by *Aerococcus*, *Corynebacterium*, *Kocuria*, *Lactococcus*, *Pseudomonas*, *Streptococcus*, *Intestinibacter* and *Tetragenococcus*. The alpha diversity of the ripening shelves (WSD_WSC) was significantly higher in the control farms compared to the case farms (Table A2) confirming the observation of the first campaign. On the contrary, the alpha diversity of the brine (B) was lower in the control farms compared to the case farms in this second campaign (Table A2). During campaign-2 was carried out in summer when the dairy cows were on the mountain pastures. In total, 12,471,211 raw paired-end reads were generated from 16S rRNA sequencing for this second campaign. After quality filtering steps, the reads were reduced to 8954,861 and 71.8 % were further assigned and classified into 63,365 total ASVs (21,088 unique ASVs). Taxonomic analysis revealed the same most abundant taxa as were found in campaign-1 (Fig. 3). However, the relative abundance of genera and the alpha diversity were higher during the campaign-2.

<table><caption>**Table 1**
Number of variables included in the multi-block analysis compared to the total number of variables from campaigns 1 and 2.</caption><thead><tr><th rowspan="2">Block (B) number - description</th><th colspan="2">N0. / n.<br>(total variables)</th><th colspan="2">n. / N.<br>(numerical variables from<br>questionnaire)</th><th colspan="2">n. / N.<br>(categorial variables from<br>questionnaire)</th><th colspan="2">n. / N.<br>(Shannon index)</th><th colspan="2">n. / N.<br>(ASVs)</th></tr><tr><th>C1 winter</th><th>C2 summer</th><th>C1 winter</th><th>C2 summer</th><th>C1 winter</th><th>C2 summer</th><th>C1 winter</th><th>C2 summer</th><th>C1 winter</th><th>C2 summer</th></tr></thead><tbody><tr><td>B1 - production before starter</td><td>28/37</td><td>22/34</td><td>0/0</td><td>0/0</td><td>0/0</td><td>0/0</td><td>1/2</td><td>0/2</td><td>36/6409</td><td>34/11,750</td></tr><tr><td>B2 - production after starter</td><td>23/35</td><td>63/84</td><td>0/18</td><td>5/18</td><td>0/0</td><td>0/0</td><td>2/8</td><td>3/8</td><td>33/5396</td><td>76/9895</td></tr><tr><td>B3 - farming practices</td><td>1/4</td><td>9/14</td><td>0/30</td><td>9/28</td><td>4/36</td><td>5/31</td><td>0/0</td><td>0/0</td><td>0/0</td><td>0/0</td></tr><tr><td>B4 - cheesemaking practices</td><td>0/3</td><td>0/2</td><td>3/13</td><td>2/14</td><td>0/31</td><td>0/38</td><td>0/0</td><td>0/0</td><td>0/0</td><td>0/0</td></tr><tr><td><strong>Total</strong></td><td><strong>52/79</strong></td><td><strong>93/134</strong></td><td><strong>3/61</strong></td><td><strong>19/60</strong></td><td><strong>4/67</strong></td><td><strong>5/69</strong></td><td><strong>4/10</strong></td><td><strong>3/10</strong></td><td><strong>69/9503</strong></td><td><strong>110/18,329</strong></td></tr></tbody></table>

N. is the total number of variables, number of unique ASVs, before univariate analysis; n. is the number of variables included into the multi-block analysis after univariate analysis; N0. is the number of variables significant after the multi-block analysis.
<caption>**Table 2**</caption>

Descriptive statistics by case and control farms for significant numerical variables from univariate analysis of campaign 1 (p-value < 0.2, Wilcoxon test).

<table><thead><tr><th rowspan="2">Numerical variable-C1</th><th colspan="4">Case</th><th colspan="4">Control</th><th rowspan="2">p-value</th></tr><tr><th>max</th><th>mean</th><th>min</th><th>N</th><th>max</th><th>mean</th><th>min</th><th>N</th></tr></thead><tbody><tr><td>cleaning_productionEquipmentFreq(n/day)</td><td>2.00</td><td>2.00</td><td>2.00</td><td>10</td><td>2.00</td><td>1.80</td><td>1.00</td><td>10</td><td>0.167</td></tr><tr><td>washingCheeseDuration(day)</td><td>6.00</td><td>4.50</td><td>3.00</td><td>10</td><td>7.00</td><td>5.40</td><td>4.00</td><td>10</td><td>0.037</td></tr><tr><td>ripeningDuration(day)</td><td>17.00</td><td>8.70</td><td>5.00</td><td>10</td><td>21.00</td><td>12.30</td><td>7.00</td><td>10</td><td>0.081</td></tr></tbody></table>

N = number of farms.

<caption>**Table 3**</caption>

Descriptive statistics by case and control farms for significant categorical variables from univariate analysis of campaign 1 (p-value < 0.2, Chi square Fisher exact test).

<table><thead><tr><th>Categorical variables-C1</th><th>N.</th><th>Category</th><th>N. case</th><th>N. control</th><th>p-value</th></tr></thead><tbody><tr><td rowspan="2">feeding_during_milking</td><td rowspan="2">20</td><td>no</td><td>9</td><td>5</td><td rowspan="2">0.141</td></tr><tr><td>yes</td><td>1</td><td>5</td></tr><tr><td rowspan="4">udder_cleaning_disinfecting</td><td rowspan="4">20</td><td>other<sup>a</sup></td><td>3</td><td>3</td><td rowspan="4">0.082</td></tr><tr><td>post</td><td>0</td><td>4</td></tr><tr><td>pre</td><td>3</td><td>0</td></tr><tr><td>pre_post</td><td>4</td><td>3</td></tr><tr><td rowspan="2">udder_cleaning_paper_wipe</td><td rowspan="2">20</td><td>no</td><td>1</td><td>6</td><td rowspan="2">0.057</td></tr><tr><td>yes</td><td>9</td><td>4</td></tr><tr><td rowspan="2">udder_cleaning_pre.soak_disinfecting</td><td rowspan="2">20</td><td>no</td><td>3</td><td>7</td><td rowspan="2">0.178</td></tr><tr><td>yes</td><td>7</td><td>3</td></tr><tr><td rowspan="2">vaccination_option</td><td rowspan="2">20</td><td>no</td><td>8</td><td>3</td><td rowspan="2">0.069</td></tr><tr><td>yes</td><td>2</td><td>7</td></tr></tbody></table>

N = number of farms.

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Boxplots of relative abundance of the top 15 most abundant taxa (at genus taxonomic rank) in all samples from control and case farms during campaign 2, in pasture in summer. Bars are coloured according to the taxonomic rank of the bacterial genus.</figcaption></figure>

2, a number of reads mapping to *Salmonella enterica* greater than 100 was found in three samples from control farms and 6 from case farms, two of which belonged to the MC_Sw samples. In addition, *S. enterica* was detected in two different types of samples from case farms. *Shigella* spp. (> 100 reads) was detected in two MC_Sw samples from case farms. No ASV corresponding to *Listeria* spp. was detected. However the presence of food-borne pathogens in the cheeses was not confirmed by the conventional VIDAS and GENE UP methods. In campaign-2, of the 129 variables from the questionnaire, 25 were selected by the univariate analysis (Tables 4 and 5). In contrast to the first campaign, the selected

variables mainly concerned feeding practices and cheese production in campaign-2 (Table 4). A total of 134 significant variables were included in the multi-block analysis (Table 1).

### 3.3. Multi-block analysis of campaign 1 (winter)

The multi-block analysis of the data collected during the winter of campaign-1 revealed 51 variables significantly associated with case or control farms (Table 1). Of the four blocks analysed, block-3, including variables related to farming practices, was the block with the highest

<footnote>a Direct branching, paper or wood wool or mop wipe cleaning operations not including pre or post soaking. Information related from the "udder_cleaning_pre-soak_disinfecting" was summarized in "udder_cleaning_disinfecting".</footnote>
<caption>**Table 4**</caption>

Descriptive statistics by case and control farms for significant numerical variables from univariate analysis of campaign 1 (p-value < 0.2, Wilcoxon test).

<table><thead><tr><th rowspan="2">Numerical variable-C2</th><th colspan="4">Case</th><th colspan="4">Control</th><th rowspan="2">p-value</th></tr><tr><th>max</th><th>mean</th><th>min</th><th>N</th><th>max</th><th>mean</th><th>min</th><th>N</th></tr></thead><tbody><tr><td>feed_number</td><td>3.00</td><td>1.80</td><td>1.00</td><td>10</td><td>2.00</td><td>1.10</td><td>0.00</td><td>10</td><td>0.048</td></tr><tr><td>feed_total_quantity<br>(kg/day per cow)</td><td>5.00</td><td>4.20</td><td>3.00</td><td>10</td><td>5.00</td><td>2.70</td><td>0.00</td><td>10</td><td>0.018</td></tr><tr><td>protein_total_quantity<br>(kg/day per cow)</td><td>1.65</td><td>1.00</td><td>0.60</td><td>10</td><td>1.01</td><td>0.50</td><td>0.00</td><td>10</td><td>0.045</td></tr><tr><td>mountain_pasture_elevation (m)<sup>a</sup></td><td>1443.00</td><td>1364.33</td><td>1300.00</td><td>3</td><td>1 648</td><td>1465.43</td><td>1400.00</td><td>7</td><td>0.171</td></tr><tr><td>mountain_pasture_duration (month)<sup>a</sup></td><td>7.00</td><td>6.00</td><td>5.00</td><td>3</td><td>5.00</td><td>4.57</td><td>4.00</td><td>7</td><td>0.049</td></tr><tr><td>breed2_dry_cows</td><td>20.00</td><td>5.50</td><td>0.00</td><td>10</td><td>8.00</td><td>1.30</td><td>0.00</td><td>10</td><td>0.147</td></tr><tr><td>calves_female</td><td>12.00</td><td>6.50</td><td>0.00</td><td>10</td><td>9.00</td><td>2.80</td><td>0.00</td><td>10</td><td>0.072</td></tr><tr><td>totCalves</td><td>15.00</td><td>7.30</td><td>1.00</td><td>10</td><td>11.00</td><td>3.50</td><td>0.00</td><td>10</td><td>0.048</td></tr><tr><td>temp_dryer (°C)</td><td>23.00</td><td>19.40</td><td>14.00</td><td>10</td><td>19.00</td><td>18.00</td><td>17.00</td><td>10</td><td>0.014</td></tr><tr><td>ripening_duration</td><td>21.00</td><td>10.20</td><td>6.00</td><td>10</td><td>21.00</td><td>13.10</td><td>7.00</td><td>10</td><td>0.188</td></tr><tr><td>cows_staff</td><td>61.00</td><td>28.40</td><td>11.33</td><td>10</td><td>29.00</td><td>21.30</td><td>13.25</td><td>10</td><td>0.173</td></tr><tr><td>dry_cows_staff</td><td>10.00</td><td>6.20</td><td>2.00</td><td>10</td><td>9.00</td><td>4.20</td><td>1.25</td><td>10</td><td>0.150</td></tr><tr><td>calves_staff</td><td>15.00</td><td>4.60</td><td>0.33</td><td>10</td><td>5.50</td><td>1.50</td><td>0.00</td><td>10</td><td>0.028</td></tr><tr><td>RMC_pH_before_coag</td><td>6.89</td><td>6.70</td><td>6.60</td><td>10</td><td>6.70</td><td>6.60</td><td>6.46</td><td>10</td><td>0.085</td></tr><tr><td>RMC_pHzanmoulding</td><td>5.66</td><td>5.40</td><td>5.17</td><td>10</td><td>5.52</td><td>5.30</td><td>5.00</td><td>10</td><td>0.197</td></tr><tr><td>RMC_temp_before_coag (°C)</td><td>34.10</td><td>33.30</td><td>32.50</td><td>10</td><td>35.40</td><td>33.90</td><td>32.30</td><td>10</td><td>0.094</td></tr><tr><td>RMC_temp_cutting (°C)</td><td>33.10</td><td>32.50</td><td>31.00</td><td>10</td><td>38.40</td><td>33.50</td><td>31.00</td><td>10</td><td>0.176</td></tr><tr><td>CC_temp_day1 (°C)</td><td>20.90</td><td>18.60</td><td>12.40</td><td>10</td><td>19.00</td><td>18.10</td><td>17.00</td><td>10</td><td>0.071</td></tr><tr><td>CC_temp_day2 (°C)</td><td>20.90</td><td>18.60</td><td>12.40</td><td>10</td><td>18.80</td><td>18.00</td><td>17.00</td><td>10</td><td>0.071</td></tr></tbody></table>

N = number of farms; day1 corresponding to the first day after coagulation, thus first day on a wooden shelf.

a Correlated questions to "moutain_pasture". The variables relating to pasture should be treated with caution, as there has been some confusion between lowland and mountain pasture in the questionnaire.

<caption>**Table 5**</caption>

Descriptive statistics by case and control farms for significant categorical variables from univariate analysis of campaign 2 (p-value < 0.2, Chi square Fisher exact test).

<table><thead><tr><th>Categorical variable-C2</th><th>N</th><th>Category</th><th>N. case</th><th>N. control</th><th>p-value</th></tr></thead><tbody><tr><td>mountain_pasture</td><td>20</td><td>No</td><td>7</td><td>3</td><td>0.178</td></tr><tr><td>mountain_pasture_start**</td><td>10</td><td>yes**</td><td>3</td><td>7</td><td></td></tr><tr><td>mountain_pasture_start**</td><td>10</td><td>mar_avr</td><td>2</td><td>0</td><td>0.067</td></tr><tr><td>mountain_pasture_start**</td><td>10</td><td>may_jun</td><td>1</td><td>7</td><td></td></tr><tr><td>udder_cleaning_disinfecting</td><td>20</td><td>other<sup>a</sup></td><td>1</td><td>2</td><td>0.177</td></tr><tr><td>udder_cleaning_paper_wipe</td><td>20</td><td>post milking</td><td>2</td><td>5</td><td></td></tr><tr><td>udder_cleaning_paper_wipe</td><td>20</td><td>pre milking</td><td>4</td><td>0</td><td></td></tr><tr><td>udder_cleaning_paper_wipe</td><td>20</td><td>pre_post</td><td>3</td><td>3</td><td></td></tr><tr><td>udder_cleaning_paper_wipe</td><td>20</td><td>no</td><td>4</td><td>8</td><td>0.170</td></tr><tr><td>udder_cleaning_wiping</td><td>20</td><td>yes</td><td>6</td><td>2</td><td></td></tr><tr><td>udder_cleaning_wiping</td><td>20</td><td>cloth</td><td>4</td><td>2</td><td>0.021</td></tr><tr><td>udder_cleaning_wiping</td><td>20</td><td>other<sup>b</sup></td><td>0</td><td>6</td><td></td></tr><tr><td>udder_cleaning_wiping</td><td>20</td><td>paper</td><td>6</td><td>2</td><td></td></tr><tr><td>litter_type</td><td>20</td><td>no litter</td><td>1</td><td>2</td><td>0.117</td></tr><tr><td>litter_type</td><td>20</td><td>sawdust or</td><td>0</td><td>1</td><td></td></tr><tr><td>litter_type</td><td>20</td><td>shavings</td><td>4</td><td>0</td><td></td></tr><tr><td>litter_type</td><td>20</td><td>straw/slope</td><td>4</td><td>0</td><td></td></tr><tr><td>litter_type</td><td>20</td><td>bedding</td><td>5</td><td>7</td><td></td></tr><tr><td>litter_type</td><td>20</td><td>other<sup>c</sup></td><td>5</td><td>7</td><td></td></tr></tbody></table>

N = number of farms;

a Direct branching, paper or wood wool or mop wipe cleaning operations not including pre or post soaking;

b Direct branching, wood wool, pre or post soaking not including paper or cloth wiping; c = rubber carpet;

** Correlated question. The variables relating to pasture should be treated with caution, as there has been some confusion between lowland and mountain pasture in the questionnaire.

contribution to the distribution of the case-study (block importance index (BI) = 33.6°% [22.1, 49.9]), followed by block-1 (BI = 23.2°% [12.2, 32.0]), block-4 (BI = 22.8°% [14.0, 33.0]) and block-2 (BI = 20.0°% [10.0, 25.3]) related to bacterial and physicochemical variables of product and environmental samples along the cheese production, but no block had a significant impact on control/case status because of a large tolerance interval (Figure B3). However, the number of significant

variables was higher in block-1 = 28/37; block-2 = 23/35 compared to block-3 = 1/4 and block-4 = 0/3. There was no statistical difference in any other variables linked to the cheesemaking process practices (block-4) in campaign-1 (Figure B1). A total of 18 significant variables distributed in block-1 and 12 with block-2 were associated with case farms.

Regarding the bacterial variables from the campaign-1, of the 69 variables associated with block-1 and -2 by univariate analysis, 45 unique ASVs were significant after multi-block analysis (Figs. 4A and 4B). Most ASVs, in block-1 (28 unique ASVs), contributing to case farms were associated with raw milk samples (7/28) and with the milking machine (11/28), in particular the milking machine filter slot swabs before washing (5/28). While, in block-2, ASVs from brine (6/21), followed by washing cheese brushes (4/21) were significantly associated with case farms. The significant ASVs were distributed across 35 genera and seven phyla, including Actinomycetota, Pseudomonadota, Bacillota, Bacteroidota, Fusobacteriota, Gracilibacteria and Thermotogota. Among these, Fusobacteriota, Gracilibacteria and Thermotogota were specifically associated with case status in this campaign. More than 62.2 % of these unique ASVs belonged to Actinomycetota and Pseudomonadota. While the Actinomycetota and Pseudomonadota were among the top 15 phyla observed during the winter campaign, Fusobacteriota and Gracilibacteria were not part of the dominant bacteriological population. Some significant bacterial genera such as *Corynebacterium*, *Staphylococcus*, *Brevibacterium*, *Streptococcus*, *Leuconostoc*, *Rothia* and *Pseudomonas* (Figure B5) belonged to the dominant microbiota (Figure 2). The bacterial taxa contributing the most to control farms were *Lactiplantibacillus pentosus* CE = -0.82 [-1.40, -0.87] and *Corynebacterium lubricantis* CE = -0.58 [-1.00, -0.55] both from MMFS_BW_Sw and *Marinomonas arctica* from brine CE = -0.60 [-0.54, 1.20]. On the contrary, *Fusobacterium necrophorum* from MMTC_BW_Sw, CE = 0.52 [0.49, 0.88], *Pseudoclavibacter* spp. CE = 0.54 [0.53, 0.90] and *Leuconostoc mesenteroides* from RM_BS, CE = 0.54 [0.52, 0.94] or *Ochrobactrum rhizosphaerae* from BWC_Sw, CE = 0.49 [0.44, 0.86], were mostly contributed to case farms (Figs. 4A and 4B). A significantly increased diversity in CC_CR and WSD_WSC was confirmed by multi-block analysis (Figure 4B). Regarding the data from questionnaire associated with farming practices (block-3), only the association between voluntary vaccination (optional to the farmer) and control farms was confirmed by
A

<figure><img src="image_7.png" /><figcaption>**Fig. 4A.** Association between the 37 explanatory variables of block-1 in campaign 1 and the outcome block Y, representing the sanitary status of the farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome Y. Significant positive and negative associations with case farms are lined respectively in blue and red. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero). MMF_MMFS_BW_Sw: Shannon diversity index regrouping the ASVs from milking machine and slot swabs filters.</figcaption></figure>

multi-block analysis in campaign-1 (Fig. 4C). These vaccinations are recommended but not required by the authorities.

### 3.4. Multi-block analysis of campaign 2 (summer)

The multi-block analysis of data collected in summer illustrated 91 variables significantly associated with case or control farms (Table 1) (block-1 = 22/34; block-2 = 63/84; block-3 = 8/14; 4 = 0/2). As was
<figure><img src="image_8.png" /><figcaption>**Fig. 4B.** Association between the 35 explicative variables of block-2 in campaign-1 and the outcome block Y, representing the case farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome variable Y. Significant positive and negative associations with case farms are lined in blue and red, respectively. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero). CC_CR: Shannon diversity index regrouping the ASVs from cheese rind and curd; WSC_WSD: Shannon diversity index regrouping the ASVs from cellar and dryer's wooden shelves.</figcaption></figure>

observed in the first campaign, the block-3 (BI = 35.1°% [23.3, 53.3]) mostly contributed to discriminating between case and control farms (Figure B4). However, the contribution of block-1 (BI = 24.6°% [15.3, 32.0]) and block-2 (BI = 21.2°% [10.0, 27.6]) was higher in summer than in winter (Figure B4). A total of 11 significant variables were associated with case farms in block-1 and 19 with block-2.

A total of 47 unique ASVs were associated with control farms and 29 with case farms. Case farms in block-1 (22 unique ASVs), were associated with raw milk samples (4/22) and with the milking machine (7/22), in particular with the washing of milking machine filters (4/22). In contrast, in block 2 (56 unique ASVs), ASVs from case farms were associated with the washing of cheese brush samples (9/56). All 75 significant ASVs were distributed in 64 genera in six phyla such as Pseudomonadota, Bacillota, Actinomycetota, Bacteroidota, Acidobacteriota and Verrucomicrobiota, and the last two phyla being specifically associated with control farms. More than 57.3% of these

unique ASVs belonged to Pseudomonadota and Bacillota. The top-5 bacterial taxa related to control farms were *Corynebacterium felinum*, CE = -0.54 [-0.95, -0.51] from raw milk (MMTC_BW_Sw), *Akkermansia muciniphila* (CE = -0.46 [-0.76, 0.33]) from MMTC_BW_Sw, *Aminicella lysinilytica* (CE = -0.46 [-0.38, 0.79]) from MMTC_BW_Sw, followed by Lachnospiraceae from MMFS_BW and *Latilactobacillus sakei* from MMF (Fig. 5A). The higher association with case farms was *Arca-nobacterium hippocoleae* in MMTC_BW_Sw (CE = 0.67 [0.52, 1.23]) followed by *Rathayibacter festucae* from RM_BS, *Kocuria salsicia* and *Exiguobacterium sibiricum* both from MMF and *Staphylococcus* spp. from BWC_Sw (Figs. 5A and 5B). Only *Staphylococcus*, *Kocuria* and *Acinetobacter* genera were common in both block-1 and -2 in campaign-2 (Figs. 5A and 5B). Regarding physico-chemical variables in block-1 and 2, control farms maintained a significantly higher temperature during both rennet addition (RMC_Temp_before_coag) and curd cutting phases (RMC_Temp_cutting), resulting in higher temperature in cheese
<figure><img src="image_9.png" /><figcaption>Fig. 4C. Association between the four explicative variables of block-3 in campaign-1 and the outcome block Y, representing the sanitary status of the farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome variable Y. Significant positive and negative associations with case farms are lined in blue and red, respectively. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero).</figcaption></figure>

at day 1 and 2 after coagulation (CC_Temp_day). The pH after rennet addition (RMC_pH_before_coag) was significantly lower in control farms (Fig. 5B).

Most variables related to farming practices (block-3) were associated with case farms (Fig. 5C). The higher association with case farms was found for the type of breed during the dry period (breed2_dry_cow), as well as a high number of mixed feed types (feed_number), daily kilograms per cow of mixed feed (feed_total_quantity), an higher daily intake of protein per cow (protein_total_quantity), litter on slope covered with straw (litter_type_straw-bedding_slope), use of udder disinfection with pre-soaking during the milking phase (udder_cleaning_disinfecting_pre) and use of paper for udder cleaning (udder_cleaning_wiping_paper). Of the ten farms practising transhumance (7 control farms and 3 case farms), the control farms started the mountain grazing period (pasture_start) significantly later (May-June) than the case farms (March-April) (Fig. 5C). However, the results for pasture need to be treated with caution and further research is needed, as there has been some confusion between lowland pasture and mountain pasture in the questionnaire.

## 4. Discussion

Classifying the main protective and risk factors for the sanitary quality of raw milk cheeses requires a relevant selection of indicators covering microbiological data, farming and cheesemaking practices and a statistical analysis managing a large number of variables in combination. To our knowledge, the association between a retrospective control-case study including bacteriological and questionnaire data analysed by multivariate analysis has never been carried out in cheese production. Compared to the previous study that used a regression tree to estimate the relative contribution of farm management practices in dairy farms (Rey-Cadilhac et al., 2023), the multiblock analysis used in this work manages a large set of variables with a relatively small sample size (number of farms) and provides the relative contribution of block of variables to explain case or control status, what the regression tree does not allow.

The results highlighted critical sampling points such as raw milk, milking machine, wooden shelves, associated with specific bacterial species related to case or control farms, while the bacterial ecosystem of cheese washing production water contributed specifically to control farms and brushes to case farms. These results complement previous work showing the influence of bacterial diversity of milking machines on bacterial contamination (Cremonesi et al., 2018; Sundberg et al., 2011; Weber et al., 2019) by extending the study to the entire raw milk cheese production process and combining practices and bacterial populations.

The bacterial taxa associated with the case or control farms depended on the stage of the production process and the campaigns conducted in winter and summer respectively. At both phylum and genus level, most of the significant bacterial variables associated with cases or controls in campaigns -1 and -2 did not overlap, suggesting that the ecosystem of critical points differed significantly between winter and summer production conditions. In addition, the number of bacterial variables correlated with case farms was higher in winter than in summer. Focusing on the milking machine samples in winter, most of the bacteria significantly found in case farms (e.g. *Rhodococcus qingshengii* or *Brachybacterium ginsengisoli*) were typical soil or water inhabitants (Iminova et al., 2022) or pathogen opportunists (*Pseudomonas* spp.). In contrast, in summer, most of the bacteria correlated with control farms (e.g. *Aminocella lysinolytica*, *Lachnospiraceae*, etc.) have the animal gut as their natural habitat (Thomas et al., 2017; Ueki et al., 2015). At the phylum level, Fusobacteriota, Gracilibacteria and Thermotogota were specifically associated with case status in winter, whereas Acidobacteriota and Verrucomicrobiota contributed to control status in summer. Some bacteria persisted among seasons such as *Klebsiella* spp., and *Enterococcus* spp., which contributed to the case farms in both campaigns and *Corynebacterium* spp., which belonged to the dominant microbiota, and was most often associated with the control farms. These results confirm the presence of *Corynebacterium* spp., known to enhance the flavour of raw milk cheeses (Anastasiou et al., 2022; Bertuzzi et al., 2018), among dominant microbiota of cheese curd (Castellanos-Rozo et al., 2020) and rind (Beresford & Williams, 2004) and from specific
<figure><img src="image_10.png" /><figcaption>**Fig. 5A.** Association between the 34 explicative variables of block-1 in campaign-2 and the outcome block Y, representing the sanitary status of the farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome variable Y. Significant positive and negative associations with case farms are lined in blue and red, respectively. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero).</figcaption></figure>

environmental niches on dairy farms (Anand et al., 2014; Flint et al., 2022), particularly on wooden shelves (Guzzon et al., 2017; Mariani, 2007) or in brine (Innocente et al., 2023). These findings highlight the importance of understanding the dynamics and impact of Corynebacterium spp. on milk production in dairy farms. As expected, most of the LAB species identified in this study, including *Enterococcus faecium*, *Enterococcus* spp., or *Tetragenococcus muriaticus* were associated with control farms. Despite the high relative abundance of *Streptococcus* and *Lactobacillus* in the whole samples of the two campaigns, none of the significant taxa corresponded to the starter cultures commonly used in the production of the PDO cheeses studied, such as *Streptococcus thermophilus*, *Lactobacillus delbrueckii* or *Lactococcus lactis* (Mariani, 2007). Interestingly, a higher diversity in product and environmental samples was observed in the control farms, especially in wooden shelves, confirming previous studies (Gaglio et al., 2022; Settanni et al., 2021). Indeed, maintaining a balanced and controlled bacterial diversity in dairy production is essential to ensure food safety and minimise the risk of growth of bacterial pathogens (Breitenwieser et al., 2020; Kable et al., 2016; Skeie et al., 2019), especially on wooden surfaces (Licitra et al., 2017; Mariani et al., 2011).

With regard to farm management practices, voluntary vaccination of dairy cows in winter was associated with the control farms, whereas most of the other management practices analysed contributed to cases of farms that could be potential risk factors of bacterial contamination. In

France, the use of these vaccines (e.g. Salmo-Past, Qfever, Respiratory syncytial virus (RSV) vaccines, etc.) is based on farm's risk assessment, disease prevalence in the region, or specific veterinary advice. In summer, most of the significant variables were related to feeding practices. The case farms gave more feed mixtures and also a higher proportion of protein, due probably to a less varied diet. Both case and control farms allowed their herds to graze in summer. However, the control farms were more likely to practise transhumance, moving their herds to mountain pastures, which probably contributed to the increase in LAB, bifidobacteria and propionibacteria, as previously shown (Secchi et al., 2023), and reduced the use of feed mixtures.

Contrary to the microbiota of the production environment, most of the variables analysed related to the cheesemaking practices were not significantly correlated with control or case farms. Only changes in temperature and pH during the first phase of production were significantly associated with contamination events detected in a past period. Higher temperatures in the RMC before coagulation and during coagulum cutting were measured in the control farms, and pH at both renneting and demoulding was less acid in the case farms, confirming the crucial role of temperature (Panthi et al., 2019) and pH during cheesemaking in cheese quality and safety (Falardeau et al., 2021; Schvartzman et al., 2014). This work confirms the considerable importance of hygiene and udder care, together with the cleanliness of milking machines, for raw milk quality (Elmoslemany et al., 2009; O'Mahony
B

<figure><img src="image_11.png" /></figure>

(caption on next page)
**Fig. 5B.** Association between the 84 explicative variables of block-2 in campaign-2 and the outcome block Y, representing the sanitary status of the farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome variable Y. Significant positive and negative associations with case farms are lined in blue and red, respectively. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero). Variables: temperature during rennet addition (RMC_Temp_before_coag), temperature during coagulum cutting phase (RMC_Temp_cutting), pH at rennet addition and after cheese's demoulding (RMC_pH_before_coag; RMC_pHDemoulding) and temperature in cheese at day 1 and 2 after coagulation (CC_Temp_day1; CC_Temp_day2).

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**Fig. 5C.** Association between the 14 explicative variables of block-3 in campaign-2 and the outcome block Y, representing the sanitary status of the farms (case versus control). The line lengths represent the 95 % confidence intervals (CI) and the point is the coefficient estimating the strength of the association between the predictor variable and the outcome variable Y. Significant positive and negative associations with case farms are lined in blue and red, respectively. Grey lines indicate that there is no significant association between the variable and the outcome (i.e., the confidence interval includes zero). The variables relating to pasture should be treated with caution, as there has been some confusion between lowland and mountain pasture in the questionnaire.

and Austin, 1991). Complementing these previous findings, this study highlights that in summer, case farms were more likely to use pre-soaking disinfection to improve udder cleaning and reduce the risk of clinical mastitis. However, this suggests that pre-soaking may not be as effective as intended in preventing contamination, or may have unintended consequences. The efficacy of pre-soaking could be affected by factors such as improper application, environmental conditions, or potential selection of resistant bacteria. In summer, farms also differed in the use of paper wiping in combination with pre-soaking suggesting that drying may be a critical step in effective hygiene. Previous work suggested that pre-soaking with a foaming product could increase the somatic cell count (SCC) (Barnouin et al., 2004). In contrast, control farms were more likely to use post-soaking disinfection or wood wool wiping. These practices are recommended by professional networks and stakeholders and their effectiveness has been supported by scientific studies. Gleeson et al. (2009) found that post-milking teat disinfection significantly reduced new intramammary infections compared to no disinfection and Williamson and Lacy-Hulbert (2013) showed that wood wool wiping was effective in reducing bacterial counts on teats (Gleeson et al., 2009; Williamson & Lacy-Hulbert, 2013).

Despite the limited number of farms that could be studied in a specific artisanal UPC production and the absence of events during the limited period of sampling time, this work has highlighted the main bacterial variables correlated with past bacterial contamination of cheese and potential critical control points. The multi-block analysis provided new insights into the correlations of variables, including the bacterial ecosystem of both the cheese and the environment, related to past contamination events of raw milk cheeses. However, based on the significant factors identified here, further larger or longitudinal studies dealing with rare potentially recurring events (Hocine & Chavance, 2010) are needed to establish causal relationships between these factors

and cheese contamination. We conclude that ASV markers from metabarcoding analysis are useful in combination with other diagnostics (such as VIDAS, PCR) and multivariate statistical analysis to help identify critical control points or potential bacterial risk factors. However, ASV quantification may also be limited due to the fact that multiple 16S rRNA sequences from the same genome may alter relative abundance results and interpretation (Johnson et al., 2019). In addition the number of ASV is not related to the growth of viable bacteria. Metabarcoding (Illumina 16S V3-V4 region) was not used here to detect bacterial pathogens. Long-read metabarcoding offers greater taxonomic resolution allowing for the potential identification of bacterial pathogens in a single food, ingredient, or environmental sample (Billington et al., 2022), but this technology would have significantly increased the cost of sample analysis. This work confirms that multi-block analysis allows the management of several types of variables including here animal health management, cheese-making procedures, environmental conditions during production and storage, and bacterial taxa. To date, no simulation studies have been carried out to assess the sensitivity of the mbPLS approach in relation to the blocks specified and the number of variables within these blocks. However, the sensitivity of the approach is influenced by the composition of the blocks and the assignment of weights to the variables. Unlike other methods, such as PLS regression, which focus on single explanatory blocks, the mbPLS approach is designed to handle multiple blocks. The mbPLS method is promising for the analysis of complex datasets and offers new possibilities for data analysis, visualisation and interpretation at a more comprehensive level (Mötus et al., 2024).

## 5. Conclusion

Food-borne pathogens cause a wide range of diseases with a
significant impact on human health and the economy. This work demonstrates the importance of characterising the ecosystems of the production environment in order to assess both the risk of contamination by foodborne pathogens and, conversely, the bacteriological factors reducing the risk of contamination. In addition, this work has identified a number of critical control points throughout the cheese production process. Finally, a sampling campaign in a control case study, combined with a suitable multivariate methodology, allowed a conjoint evaluation of bacteriological, physico-chemical parameters as well as farming and production practices. This methodology could be extended to other food-borne viral or bacteriological pathogens and to other types of cheese production.

## Ethics approval and consent to participate

Farmers agreed to participate in the study signing a written form stating that all given information is confidential and will only be published in an anonymised fashion. Questionnaires which do not ask for personal data are analysed anonymously and do not require clearance by the ethics committee.

## Data availability: sequencing data and codes

After acceptance of the manuscript, the raw data and codes supporting the findings of this study will be openly available at https://doi.org/10.18167/DVN1/6ZOPLZ

## Fundings

This work was supported by the University of Montpellier, the Centre de coopération internationale en recherche agronomique pour le développement (CIRAD) and the Agence nationale de sécurité sanitaire de l'alimentation, de l'environnement et du travail (Anses). For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## CRediT authorship contribution statement

**Alisia Chiado Rana:** Writing – original draft, Visualization, Software, Investigation, Formal analysis, Data curation, Conceptualization. **Géraldine Cazeau:** Writing – review & editing, Software, Methodology, Formal analysis, Conceptualization. **Stéphanie Bougeard:** Writing – review & editing, Software, Formal analysis. **Kevin Bethune:** Investigation. **Vincent Chochois:** Writing – review & editing, Supervision, Software, Methodology, Formal analysis, Data curation, Conceptualization. **Corinne Teyssier:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. **Charlotte Valat:** Writing – review & editing, Writing – original draft, Validation, Supervision, Methodology, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The data used in this paper were partly generated by the genotyping and sequencing facilities of the ISEM (Institut des Sciences de l'Evolution-Montpellier) and the Labex Centre Méditerranéen Environnement Biodiversité. We thank all the farmers involved in these campaigns. We thank Didier MONTET for his contribution to the initial conception of this project.

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.afres.2024.100617.

## Data availability

I have share the link to my data/code
Dataverse Lactic4safe (Original data)

## References

Anand, S., Singh, D., Avadhanula, M., & Marka, S. (2014). Development and control of bacterial biofilms on dairy processing membranes: Biofilms on dairy processing membranes.... *Comprehensive Reviews in Food Science and Food Safety*, *13*(1), 18–33.
https://doi.org/10.1111/1541-4337.12048

Anastasiou, R., Kazou, M., Georgalaki, M., Aktypis, A., Zoumpopoulou, G., & Tsakalidou, E. (2022). Omics approaches to assess flavor development in cheese. *Foods* (*Basel*, Switzerland), *11*(2). https://doi.org/10.3390/foods11020188. Article 2.

ANSES. (2022). Raw-milk cheeses: What are the associated health risks and what preventive measures can be taken? *Anses - Agence Nationale de Sécurité Sanitaire de l'alimentation, de l'environnement et Du Travail*. https://www.anses.fr/en/content/ra/w-milk-cheeses-what-are-associated-health-risks-and-what-preventive-measures-can-be-taken.

Bagel, A., & Sergentet, D (2022). Shiga toxin-producing escherichia coli and milk fat globules. *Microorganisms*, *10*(3), 496. https://doi.org/10.3390/microorganisms10030496

Balcıunas, E. M., Castillo Martinez, F. A., Todorov, S. D., Franco, B. D. G. de M., Converti, A., & Oliveira, R. P. de S. (2013). Novel biotechnological applications of bacteriocins: A review. *Food Control*, *32*(1), 134–142. https://doi.org/10.1016/j foodcont.2012.11.025

Barnouin, J., Chassagne, M., Bazin, S., & Boichard, D (2004). Management practices from questionnaire surveys in herds with very low somatic cell score through a national mastitis program in France. *Journal of Dairy Science*, *87*(11), 3989–3999. https://doi.org/10.3168/jds.S0022-0302(04)73539-0

Beresford, T., & Williams, A. (2004). The microbiology of cheese ripening. In P. F. Fox, P. L. H. McSweeney, T. M. Cogan, & T. P. Guinee (Eds.), *Cheese: Chemistry, physics and microbiology* (pp. 287–317). Academic Press. https://doi.org/10.1016/S1874-558X(04)80071-X. Vol. 1.

Bertuzzi, A. S., Walsh, A. M., Sheehan, J. J., Cotter, P. D., Crispie, F., McSweeney, P. L. H., Kilcawley, K. N., & Rea, M. C. (2018). Omics-Based Insights into Flavor Development and Microbial Succession within Surface-Ripened Cheese. *mSystems*, *3*(1). https://doi.org/10.1128/msystems.00211-17. 10.1128/msystems.00211-17.

Betterla, L., Levante, A., Bancalari, E., Bottari, B., & Gatti, M. (2023). Lactic acid bacteria in cow raw milk for cheese production: Which and how many? *Frontiers in Microbiology*, *13*. https://doi.org/10.3389/fmicb.2022.1092224

Billington, C., Kingsbury, J. M., & Rivas, L. (2022). Metagenomics approaches for improving food safety: A review. *Journal of Food Protection*, *85*(3), 448–464. https://doi.org/10.4315/JFP-21-301

Bougeard, S., & Dray, S. (2018). Supervised Multiblock Analysis in R with the ade4 Package. *Journal of Statistical Software*, *86*, 1–17. https://doi.org/10.18637/jss.v086.i01

Bougeard, S., Qannari, E. M., & Rose, N. (2011). Multiblock redundancy analysis: Interpretation tools and application in epidemiology. *Journal of Chemometrics*, *25*(9), 467–475. https://doi.org/10.1002/cem.1392

Breitenwieser, F., Doll, E. V., Clavel, T., Scherer, S., & Wenning, M. (2020). Complementary use of cultivation and high-throughput amplicon sequencing reveals high biodiversity within raw milk microbiota. *Frontiers in Microbiology*, *11*, 1557. https://doi.org/10.3389/fmicb.2020.01557

Bücher, C., Burtsher, J., Rudavsky, T., Zitz, U., & Domig, K. J. (2024). Sources of propionic acid bacteria contamination in the milking parlor environment on Alpine Dairy Farms. *Journal of Dairy Science*. https://doi.org/10.3168/jds.2024-24752

Burgain, J., Scher, J., Francius, G., Borges, F., Corgneau, M., Revol-Junelles, A. M., Cailliez-Grimal, C., & Gaiani, C. (2014). Lactic acid bacteria in dairy food: Surface characterization and interactions with food matrix components. *Advances in Colloid and Interface Science*, *213*, 21–35. https://doi.org/10.1016/j.cis.2014.09.005

Busetta, G., Garofalo, G., Claps, S., Sardina, M. T., Franciosi, E., Alfonzo, A., Francesca, N., Moschetti, G., Settanni, L., & Gaglio, R. (2024). The wooden shelf surface and cheese rind mutually exchange microbiota during the traditional ripening process. *International Journal of Food Microbiology*, *409*, Article 110478. https://doi.org/10.1016/j.ijfoodmicro.2023.110478

Castellanos-Rozo, J., Pérez Pulido, R., Grande, M. J., Lucas, R., & Gálvez, A. (2020). Analysis of the bacterial diversity of paipa cheese (a Traditional Raw Cow's Milk Cheese from Colombia) by high-throughput sequencing. *Microorganisms*, *8*(2). https://doi.org/10.3390/microorganisms8020218. Article 2.

Castellote, J., Fraud, S., Irlinger, F., Swennen, D., Fer, F., Bonnarme, P., & Monnet, C. (2015). Investigation of *Geotrichum candidum* gene expression during the ripening of Reblochon-type cheese by reverse transcription-quantitative PCR. *International Journal of Food Microbiology*, *194*, 54–61. https://doi.org/10.1016/j.ijfoodmicro.2014.11.009
Castro, H., Rususen, M., & Lindström, M. (2017). Occurrence and growth of Listeria monocytogenes in packaged raw milk. *International Journal of Food Microbiology*, **261**, 1-10. https://doi.org/10.1016/j.ijfoodmicro.2017.08.017

Ceugniez, A., Coucheny, F., Jacques, P., Daube, G., Delcenserie, V., & Drider, D (2017). Anti- Salmonella activity and probiotic trends of Kluyveromyces marxianus S-2-05 and Kluyveromyces lactis S-3-05 isolated from a French cheese, Tomme d'Orchies. *Research in Microbiology*, **168**(6), 575-582. https://doi.org/10.1016/j.resmic.2017.03.004

Chemidlin Prévost-Bouré, N., Karimi, B., Sadet-Bourgeteau, S., Djemiel, C., Brie, M., Dumont, J., Campedelli, M., Nowak, V., Guyot, P., Letourneau, C., Manneville, V., Gillet, F., & Bouton, Y (2021). Microbial transfers from permanent grassland ecosystems to milk in dairy farms in the Comté cheese area. *Scientific Reports*, **11**(1). https://doi.org/10.1038/s41598-021-97373-6. Article 1.

CNIEL. (2022). Chiffres clés des fromages au lait cru 2021. CNIEL- Centre National Interprofessionnel de l'Economie Laitière. https://presse.filiere-laitiere.fr/photos/c/iffres-cles-des-fromages-au-lait-cru-2022-35ed-ef05e.html.

Cogan, T. M., Goerges, S., Gelsomino, R., Larpin, S., Hohenegger, M., Bora, N., Jamet, E., Rea, M. C., Mounier, J., Vancanneyt, M., Guéguen, M., Desmasures, N., Swings, J., Goodfellow, S., Ward, A. C., Sebastiani, H., Irlinger, F., Chamba, J.-F., Beduhn, R., & Scherer, S. (2014). Biodiversity of the surface microbial consortia fromburger, reblochon, livarot, tiltsit, and gubben cheeses. *Microbiology Spectrum*, **2**(1). https://doi.org/10.1128/microbiolspec.CM-0010-2012

Cremones, P., Ceccarani, C., Curone, G., Severgnini, M., Pollera, C., Bronzo, V., Riva, F., Addis, M. F., Filipe, J., Amadori, M., Trevisi, E., Vigo, D., Moroni, P., & Castiglioni, B. (2018). Milk microbiome diversity and bacterial group prevalence in a comparison between healthy Holstein Friesian and Rendena cows. *PLOS ONE*, **13**(10), Article e0205054. https://doi.org/10.1371/journal.pone.0205054

Cruciata, M., Gaglio, R., Scatassa, M. L., Sala, G., Cardamone, C., Palmeri, M., Moschetti, G., La Mantia, T., & Settanni, L. (2018). Formation and characterization of early bacterial biofilms on different wood typologies applied in dairy production. *Applied and Environmental Microbiology*, **84**(4), e02107-e02117. https://doi.org/10.1128/AEM.02107-17

Didienne, R., Defargues, C., Gallon, C., Meylheuc, T., Hulin, S., & Montel, M.-C. (2012). Characteristics of microbial biofilm on wooden vats ('gerles') in PDO Salers cheese. *International Journal of Food Microbiology*, **156**(2), 91-101. https://doi.org/10.1016/j.ijfoodmicro.2012.03.007

Donnelly, C. W. (2013). From pasteur to probiotics: A historical overview of cheese and microbes. *Microbiology Spectrum*, **1**(1). https://doi.org/10.1128/microbiolspec.cm-0001-12. 10.1128/microbiolspec.cm-0001-0012.

Dreier, M., Berthoud, H., Shani, N., Wechsler, D., & Junier, P. (2021). Development of a high-throughput microfluidic qpcr system for the quantitative determination of quality-relevant bacteria in cheese. *Frontiers in Biology*, **11**. https://www.frontiersin.org/articles/10.3389/fmicb.2020.619166.

El Jeni, R., Villot, C., Koyun, O. Y., Osorio-Doblado, A., J. J., Lourenco, J. M., Steele, M., & Callaway, T. R (2024). *Invited review: Probiotic* approaches to improving dairy production: Reassessing "magic foo-foo dust". *Journal of Dairy Science*, **107**(4), 1832-1856. https://doi.org/10.3168/jds.2023-23831

Elmoslemany, A. M., Keefe, G. P., Dohoo, I. R., & Jayarao, B. M. (2009). Risk factors for bacteriological quality of bulk tank milk in Prince Edward Island dairy herds. Part 1: Overall risk factors. *Journal of Dairy Science*, **92**(6), 2634-2643. https://doi.org/10.3168/jds.2008-1812

Elmoslemany, A. M., Keefe, G. P., Dohoo, I. R., Wichtel, J. J., Stryhn, H., & Dingwell, R. T. (2010). The association between bulk tank milk analysis for raw milk quality and on-farm management practices. *Preventive Veterinary Medicine*, **95**(1), 32-40. https://doi.org/10.1016/j.prevetmed.2010.03.007

European Commission. (2005). COMMISSION REGULATION (EC) No 2073/2005 of 15 November 2005 on microbiological criteria for foodstuffs [Regulation]. https://data.europa.eu/eli/reg/2005/2073/oj.

European Commission. (2012). Regulation (EU) No 1151/2012 of the European Parliament and of the Council of 21 November 2012on quality schemes for agricultural products and foodstuffs.

Eurostat. (2022). Milk collection (all milks) and dairy products obtained—Annual data. Eurostat. https://doi.org/10.2908/APRO_MK POBTA

Falardeau, J., Trmcíč, A., & Wang, S. (2021). The occurrence, growth, and biocontrol of Listeria monocytogenes in fresh and surface-ripened soft and semisoft cheeses. *Comprehensive Reviews in Food Science and Food Safety*, **20**(4), 4019-4048. https://doi.org/10.1111/1541-4337.12768

FAO. (2022). Codex General Standard for Use of Dairy Products (CXS 206-1999). https://www.fao.org/fao-who-codexalimentarius/sh/proxy/en/?link=1andurl=https%3A%252F%252Fworkspace.fao.org%252Fsites%252Fcodex%252Fstandards%252F CXS%2B2B06-1999%252FCXS_206e.pdf

Farrokh, C., Jordan, K., Auvray, F., Glass, K., Oppegård, H., Raynaud, S., Thevenot, D., Condron, R., De Reu, K., Govaris, A., Heggum, K., Heyndrickx, M., Hummerjohann, J., Lindsay, D., Miszczycha, S., Moussiegt, S., Verstraete, K., & Cerf, O. (2013). Review of Shiga-toxin-producing *Escherichia coli* (STEC) and their significance in dairy production. *International Journal of Food Microbiology*, **162**(2), 190-212. https://doi.org/10.1016/j.ijfoodmicro.2012.08.008

Flint, S., Bremer, P., Brooks, J., Palmer, J., Burgess, S., & Seale, B. (2022). Filmobiont formation and control in the dairy industry. *Encyclopedia of Dairy Sciences*. https://doi.org/10.1016/B978-0-08-100596-5.00639-9

FSANZ. (2022). Microbiological Risk Assessment of Raw Milk Cheese. https://www.foodstdata. ndards.gov.au/sites/default/files/food-standards-code/proposals/Documents/ P1007%20PPPS%20for%20raw%20milk%201AR%20SD3%20Cheese%20Risk% 20Assessment.pdf.

Gaglio, R., Busetta, G., Gannuscio, R., Settanni, L., Licitra, G., & Todaro, M. (2022). A multivariate approach to study the bacterial diversity associated to the wooden

shelves used for aging traditional sicilian cheeses. *Foods* (Basel, Switzerland), **11**(5). https://doi.org/10.3390/foods1105774. Article 5.

Gleeson, D., O'Brien, B., Flynn, J., O'Callaghan, E., & Galli, F. (2009). Effect of pre-milking teat preparation procedures on the microbial count on teats prior to cluster application. *Irish Veterinary Journal*, **62**(7), 461. https://doi.org/10.1186/2046- 0481-62-7-461

Guzzon, R., Carafa, I., Tuohy, K., Cervantes, G., Vernetti, L., Barmaz, A., Larcher, R., & Franciosi, E. (2017). Exploring the microbiota of the red-brown defect in smear- ripened cheese by 454-pyrosequencing and its prevention using different cleaning systems. *Food Microbiology*, **62**, 160-168. https://doi.org/10.1016/j.fm.2016.10.018

Hocine, M. N., & Chavance, M. (2010). La méthode de la série de cas. *Revue d'Épidémiologie et de Santé Publique*, **58**(6), 435-440. https://doi.org/10.1016/j. respe.2010.06.168

Hohmann, M.-F., Wente, N., Zhang, Y., & Krömker, V. (2020). Bacterial load of the teat apex skin and associated factors at herd level. *Animals: An Open Access Journal from MDPI*, **10**(9), 1647. https://doi.org/10.3390/ani10091647

Iminova, L., Delegan, Y., Frantsuzova, E., Bogun, A., Zvonarev, A., Suzina, N., Anbumani, S., & Solyanikova, I. (2022). Physiological and biochemical characterization and genome analysis of *Rhodococcus qinghengii* strain 7B capable of crude oil degradation and plant stimulation. *Biotechnology Reports*, **35**, e00741. https://doi.org/10.1016/j.btre.2022.e00741

INAQ. (2022). Chiffres clés 2021 des produits laitiers aop et igp. September. IINAO-Institut national de l'origine et de la qualité https://www.inao.gouv.fr/Nos-actualites/ch iffres-cles-produits-laitiers-2021.

Innocente, N., Renoldi, N., Moret, E., Maifreni, M., & Marino, M. (2023). Volatilome of brine-related microorganisms in a curd-based medium. *Journal of Dairy Science*, **106** (12), 8404-8414. https://doi.org/10.3168/jds.2022-23051

Johnson, J. S., Spakowicz, D. J., Hong, B.-Y., Petersen L. M., Demkowicz, P., Chen, L., et al. (2019). Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. *Nature Communications*, **10**(1), 5029. https://doi.org/10.1038/ s41467-019-13036-1

Josse, J., & Husson, F. (2016). missMDA: A package for handling missing values in multivariate data analysis. *Journal of Statistical Software*, **70**, 1-31. https://doi.org/ 10.18637/jss.v07.o1i0

Kable, M. E., Srisengfa, Y., Laird, M., Zaragoza, J., McLeod, J., Heidenreich, J., & Marco, M. L. (2016). The core and seasonal microbiota of raw bovine milk in tanker trucks and the impact of transfer to a milk processing facility. *mBio*, **7**(4). https:// doi.org/10.1128/mbio.00836-16, 10.1128/mbio.00836-16.

Klindworth, A., Pruesse, E., Schweer, T., Peplies, J., Quast, C., Horn, M., & Glöckner, F. O. (2013). Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. *Nucleic Acids Research*, **41**(1), e1. https://doi.org/10.1093/nar/gks808

Li, N., Wang, Y., You, C., Ren, J., Chen, W., Zheng, H., & Liu, Z. (2018). Variation in raw milk microbiota throughout 12 months and the impact of weather conditions. *Scientific Reports*, **8**(1), 2371. https://doi.org/10.1038/s41598-018-20862-8

Licitra, G., Caccamo, M., Valence, F., & Lortal, S. (2017). Traditional wooden equipment used for cheesemaking and their effect on quality. *Global cheesemaking technology* (pp. 157-172). John Wiley and Sons, Ltd. https://doi.org/10.1002/9781119046165. ch06

Lima, S. F., Bicalho, M. L. de S., & Bicalho, R. C. (2018). Evaluation of milk sample fractions for characterization of milk microbiota from healthy and clinical mastitis cows. *PLOS ONE*, **13**(3), Article e0193671. https://doi.org/10.1371/journal.pone.0193671

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biology*, **15**(12), 550. https://doi.org/ 10.1186/s13059-014-0550-8

Lundberg, D. S., Yourstone, S., Mieczkowski, P., Jones, C. D., & Dangl, J. L. (2013). Practical innovations for high-throughput amplicon sequencing. *Nature Methods*, **10**. https://doi.org/10.1038/nmeth.2634. Article 10.

Mariani, C. (2007). *Ecologie microbienne des biofilms présents à la surface des planches d'affinage en bois de l'AOC "Reblochon de Savoie" et effet inhibiteur vis à vis de Listeria monocytogenes* (NNT: 2007AGPT0089 pastel-00005194). ENSIA (AgroParisTech).

Mariani, C., Oulahal, N., Chamba, J.-F., Dubois-Brissonnet, F., Notz, E., & Briandet, R. (2011). Inhibition of Listeria monocytogenes by resident biofilms present on wooden shelves used for cheese ripening. *Food Control*, **22**(8), 1357-1362. https://doi.org/ 10.1016/j foodcont.2011.02.012

McMurdie, P. J., & Holmes, S. (2014). Waste want, not want: Why rarefying microbiome data is inadmissible. *PLOS Computational Biology*, **10**(4), Article e1003531. https:// doi.org/10.1371/journal.pcbi.1003531

Mendonça, R., Furtado, R., Coelho, A., Correia, C. B., Suyarko, E., Borges, V., Gomes, J. P., Pista, A., & Batista, R. (2024). Raw milk cheeses from Beira Baixa, Portugal—A contributive study for the microbiological hygiene and safety assessment. *Brazilian Journal of Microbiology*. https://doi.org/10.1007/s42770-024- 01332-y

Meola, M., Rifa, E., Shani, N., Delbès, C., Berthoud, H., & Chassard, C. (2019). DAIRYdb: a manually curated reference database for improved taxonomy annotation of 16S rRNA gene sequences from dairy products. *BMC Genomics*, **20**(1), 560. https://doi.org/10.1186/s12864-019-5914-8

Monnet, C., Dugat-Bony, E., Swennen, D., Beckerich, J.-M., Irlinger, F., Fraud, S., & Bonnarme, P. (2016). Investigation of the activity of the microorganisms in a reblochon-style cheese by metatranscriptomic analysis. *Frontiers in Microbiology*, **7**. https://doi.org/10.3389/fmicb.2016.00536

Montel, M.-C., Buchin, S., Mallet, A., Delbes-Paus, C., Vuitton, D. A., Desmasures, N., & Berthier, F. (2014). Traditional cheeses: Rich and diverse microbiota with associated benefits. *International Journal of Food Microbiology*, **177**, 136-154. https://doi.org/ 10.1016/j.ijfoodmicro.2014.02.019
Mōkus, K., Viidu, D.-A., Rilanto, T., Niine, T., Orro, T., Viltrop, A., & Bougeard, S. (2024). Application of multiblock analysis to identify key areas and risk factors for dairy cow persistence. *Preventive Veterinary Medicine, **222**, Article 106081. https://doi.org/10.1016/j.prevetmed.2023.106081

Oliver, S. P., Jayarao, B. M., & Almeida, R. A. (2005). *Foodborne Pathogens and Disease, **2*** (2), 115–129. https://doi.org/10.1089/fpd.2005.2.115

Oliver, S. P., & Murinda, S. E. (2011). Milk and raw milk consumption as a vector for human disease. In D. O. Krause, & S. Hendrick (Eds.), *Zoonotic pathogens in the food chain* (pp. 99–118). CABI. https://doi.org/10.1079/9781845936815.0099.

O'Mahony, P. A., & Austin, F. H. (1991). A study of on-farm factors affecting bacteriological quality of bulk farm milk. *Irish Journal of Agricultural Research, **30*** (2), 113–119.

Panthi, R. R., Kelly, A. L., Sheehan, J. J., Bulbul, K., Vollmer, A. H., & McMahon, D. J. (2019). Influence of protein concentration and coagulation temperature on rennet-induced gelation characteristics and curd microstructure. *Journal of Dairy Science, **102*** (1), 177–189. https://doi.org/10.3168/jds.2018-15039

Perrin, F., Tenenhaus-Aziza, F., Michel, V., Miszczycha, S., Bel, N., & Sanaa, M. (2015). Quantitative risk assessment of haemolytic and uremic syndrome linked to O157:H7 and Non-O157:H7 shiga-toxin producing *escherichia coli* strains in raw milk soft cheeses: Quantitative risk assessment of HUS linked to pathogenic STEC in cheese. *Risk Analysis, **35*** (1), 109–128. https://doi.org/10.1111/risa.12267

Ponzoni, E., Mastromauro, F., Giani, S., & Breviario, D (2009). Traceability of plant diet contents in raw cow milk samples. *Nutrients, **1*** (2). https://doi.org/10.3390/nu1020251. Article 2.

Pozzi, A., Nazzicari, N., Capoferri, R., Radovic, S., & Bongioni, G. (2023). Assessment of residual plant DNA in bulk milk for Grana Padano PDO production by a metabarcoding approach. *PLOS ONE, **18*** (7), Article e0289108. https://doi.org/10.1371/journal.pone.0289108

Quigley, L., O'Sullivan, O., Stanton, C., Beresford, T. P., Ross, R. P., Fitzgerald, G. F., & Cotter, P. D. (2013). The complex microbiota of raw milk. *FEMS Microbiology Reviews, **37*** (5), 664–698. https://doi.org/10.1111/1574-6976.12030

Rey-Cadilhac, L., Botreau, R., Ferlay, A., Hulin, S., Hurtaud, C., Lardy, R., Martin, B., & Laurent, C. (2021). Co-construction of a method for evaluating the intrinsic quality of bovine milk in relation to its fate. *Animal: An International Journal of Animal Bioscience, **15*** (7), Article 100264. https://doi.org/10.1016/j.animal.2021.100264

Rey-Cadilhac, L., Ferlay, A., Gelé, M., Léger, S., & Laurent, C. (2023). Regression trees to identify combinations of farming practices that achieve the best overall intrinsic quality of milk. *Journal of Dairy Science, **106*** (2), 1026–1038. https://doi.org/10.3168/jds.2022-22486

Ricci, A., Capello, K., Cibin, V., Pozza, G., Ferrè, N., Barrucci, F., Menin, R., Farina, R., & Marangon, S. (2013). Raw milk-associated foodborne infections: A scoring system for the risk-based categorisation of raw dairy farms. *Research in Veterinary Science, **95*** (1), 69–75. https://doi.org/10.1016/j.rvsc.2013.02.007

Rubiola, S., Chiesa, F., Dalmasso, A., Di Ciccio, P., & Civera, T. (2020). Detection of antimicrobial resistance genes in the milk production environment: Impact of host DNA and sequencing depth. *Frontiers in Microbiology, **11***. https://www.frontiersin.org/article/10.3389/fmicb.2020.01983.

Rubiola, S., Macori, G., Civera, T., Fanning, S., Mitchell, M., & Chiesa, F. (2022). Comparison between full-length 16S rRNA metabarcoding and whole metagenome sequencing suggests the use of either is suitable for large-scale microbiome studies. *Foodborne Pathogens and Disease, **19*** (7), 495–504. https://doi.org/10.1089/fpd.2022.0027

Schvartzman, M. S., Gonzalez-Barron, U., Butler, F., & Jordan, K. (2014). Modeling the growth of *Listeria monocytogenes* on the surface of smear- or mold-ripened cheese. *Frontiers in Cellular and Infection Microbiology, **4***. https://doi.org/10.3389/fcimb.2014.00090

Secchi, G., Amalfitano, N., Carafa, I., Franciosi, E., Gallo, L., Schiavon, S., Sturaro, E., Tagliapietra, F., & Bittante, G. (2023). Milk metagenomics and cheese-making properties as affected by indoor farming and summer highland grazing. *Journal of Dairy Science, **106*** (1), 96–116. https://doi.org/10.3168/jds.2022-22449

Settanni, L., Busetta, G., Puccio, V., Licitra, G., Franciosi, E., Botta, L., Di Gerlando, R., Todaro, M., & Gaglio, R. (2021). In-depth investigation of the safety of wooden

shelves used for traditional cheese ripening. *Applied and Environmental Microbiology, **87*** (23), https://doi.org/10.1128/AEM.01524-21. e01524-21.

Shannon, C. E. (1948). A mathematical theory of communication. *BELL System Technical Journal*, 379-423.

Siebert, A., Hofmann, K., Staib, L., Doll, E. V., Scherer, S., & Wenning, M. (2021). Amplicon-sequencing of raw milk microbiota: Impact of DNA extraction and library-PCR. *Applied Microbiology and Biotechnology, **105*** (11), 4761–4773. https://doi.org/10.1007/s00253-021-11353-4

Skeie, S. B., Håland, M., Thorsen, I. M., Narvhus, J., & Porcellato, D (2019). Bulk tank raw milk microbiota differs within and between farms: A moving goalpost challenging quality control. *Journal of Dairy Science, **102*** (3), 1959–1971. https://doi.org/10.3168/jds.2017-14083

Smith, R. P., Andres, V., Cheney, T. E., Martelli, F., Gosling, R., Marier, E., Rabie, A., Gilson, D., & Davies, R. H. (2018). How do pig farms maintain low *Salmonella* prevalence: A case-control study. *Epidemiology and Infection, **146*** (15), 1909–1915. https://doi.org/10.1017/S0950268818002248

Stinson, L. F., Ma, J., Rea, A., Dymock, M., & Geddes, D. T. (2021). Centrifugation does not remove bacteria from the fat fraction of human milk. *Scientific Reports, **11*** (1). https://doi.org/10.1038/s41598-020-79793-y. Article 1.

Sun, L., & D'Amico, D. J. (2021). Composition, succession, and source tracking of microbial communities throughout the traditional production of a farmstead cheese. *mSystems, **6*** (5), https://doi.org/10.1128/msystems.00830-21, 10.1128/msystems.00830-21.

Sundberg, M., Christiansson, A., Lindahl, C., Wahlund, L., & Birgersson, C. (2011). Cleaning effectiveness of chlorine-free detergents for use on dairy farms. *The Journal of Dairy Research, **78*** (1), 105–110. https://doi.org/10.1017/S0022029910000762

Thomas, M., Webb, M., Ghimire, S., Blair, A., Olson, K., Fenske, G. J., Fonder, A. T., Christopher-Hennings, J., Brake, D., & Scaria, J. (2017). Metagenomic characterization of the effect of feed additives on the gut microbiome and antibiotic resistome of feedlot cattle. *Scientific Reports, **7***, 12257. https://doi.org/10.1038/s41598-017-12481-6

Ueki, A., Shibuya, T., Kaku, N., & Ueki, K. (2015). Aminocella lysinolytica gen. nov., sp. nov., a L-lysine-degrading, strictly anaerobic bacterium in the class *Clostridia* isolated from a methanogenic reactor of cattle farms. *Archives of Microbiology, **197*** (1), 97–104. https://doi.org/10.1007/s00203-014-1066-z

Wang, B., Tan, X., Du, R., Zhao, F., Zhang, L., Han, Y., & Zhou, Z. (2019). Bacterial composition of biofilms formed on dairy-processing equipment. *Preparative Biochemistry and Biotechnology, **49*** (5), 477–484. https://doi.org/10.1080/10826068.2019.1587623

Weber, M., Liedtke, J., Plattes, S., & Lipski, A. (2019). Bacterial community composition of biofilms in milking machines of two dairy farms assessed by a combination of culture-dependent and -independent methods. *PLOS ONE, **14*** (9), Article e0222238. https://doi.org/10.1371/journal.pone.0222238

Williamson, J. H., & Lacy-Hulbert, S. J. (2013). Effect of disinfecting teats post-milking or pre- and post-milking on intramammary infection and somatic cell count. *New Zealand Veterinary Journal, **61*** (5), 262–268. https://doi.org/10.1080/00480169.2012.751576

Wold, S. (1966). *Estimation of principal components and related models by iterative least squares. multivariate analysis* (pp. 391–420). Academic Press.

Wolfe, B. E., Button, J. E., Santarelli, M., & Dutton, R. J. (2014). Cheese rind communities provide tractable systems for in situ and in vitro studies of microbial diversity. *Cell, **158*** (2), 422–433. https://doi.org/10.1016/j.cell.2014.05.041

Yamaguchi, E., Hayama, Y., Murato, Y., Sawai, K., Kondo, S., & Yamamoto, T. (2024). A case-control study of the infection risk of H5N8 highly pathogenic avian influenza in Japan during the winter of 2020-2021. *Research in Veterinary Science, **168*, Article 105149. https://doi.org/10.1016/j.rvsc.2024.105149

Yoon, Y., Lee, S., & Choi, K.-H. (2016). Microbial benefits and risks of raw milk cheese. *Food Control, **63***, 201–215. https://doi.org/10.1016/j.foodcont.2015.11.013

Zago, M., Bardelli, T., Rossetti, L., Nazzicari, N., Carminati, D., Galli, A., & Giraffa, G. (2021). Evaluation of bacterial communities of Grana Padano cheese by DNA metabarcoding and DNA fingerprinting analysis. *Food Microbiology, **93*, Article 103613. https://doi.org/10.1016/j.fm.2020.103613