# Analysis of Microbiota Persistence in Quebec's Terroir Cheese Using a Metabarcoding Approach

Annick Raymond-Fleury, Marie-Hélène Lessard, Julien Chamberland, Yves

Pouliot, Eric Dugat-Bony, Sylvie Turgeon, Daniel St-Gelais, Steve Labrie

► To cite this version:

Annick Raymond-Fleury, Marie-Hélène Lessard, Julien Chamberland, Yves Pouliot, Eric Dugat-Bony, et al.. Analysis of Microbiota Persistence in Quebec's Terroir Cheese Using a Metabarcoding Approach. Microorganisms, 2022, 10 (7), pp.1381. 10.3390/microorganisms10071381. hal-03721977

HAL Id: hal-03721977

https://hal.inrae.fr/hal-03721977

Submitted on 13 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Analysis of Microbiota Persistence in Quebec's Terroir Cheese Using a Metabarcoding Approach

Annick Raymond-Fleury <fnref n="1" />,†, Marie-Hélène Lessard <fnref n="1" />,†, Julien Chamberland <fnref n="1" />,†, Yves Pouliot <fnref n="1" />,†, Eric Dugat-Bony <fnref n="2" />,†, Sylvie L. Turgeon <fnref n="1" />,†, Daniel St-Gelais <fnref n="1" />,<fnref n="3" /> and Steve Labrie <fnref n="1" />,* †,†

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

check for
updates

**Citation:** Raymond-Fleury, A.; Lessard, M.-H.; Chamberland, J.; Pouliot, Y.; Dugat-Bony, E.; Turgeon, S.L.; St-Gelais, D.; Labrie, S. Analysis of Microbiota Persistence in Quebec's Terroir Cheese Using a Metabarcoding Approach.
*Microorganisms* **2022**, *10*, 1381.
https://doi.org/10.3390/microorganisms10071381

Academic Editor: Alexander I. Netrusov

Received: 4 June 2022

Accepted: 5 July 2022

Published: 9 July 2022

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

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

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

**Abstract:** Environmental short amplicon sequencing, or metabarcoding, is commonly used to characterize the bacterial and fungal microbiota of cheese. Comparisons between different metabarcoding studies are complicated by the use of different gene markers. Here, we systematically compare different metabarcoding molecular targets using V3–V4 and V6–V8 regions of the bacterial 16S rDNA and fungal ITS1 and ITS2 regions. Taxonomic profiles varied depending on the molecular markers used. Based on data quality and detection capacity of the markers toward microorganisms usually associated with the dairy environment, the ribosomal regions V3–V4 and ITS2 were selected and further used to evaluate variability in the microbial ecosystem of terroir cheeses from the province of Quebec in Canada. Both fungal and bacterial ecosystem profiles were described for 32 different ready-to-eat bloomy-, washed- and natural-rind specialty cheese varieties. Among them, 15 were studied over two different production years. Using the Bray–Curtis dissimilarity index as an indicator of microbial shifts, we found that most variations could be explained by either a voluntary change in starter or ripening culture composition, or by changes in the cheesemaking technology. Overall, our results suggest the persistence of the microbiota between the two years studied—these data aid understanding of cheese microbiota composition and persistence during cheese ripening.

**Keywords:** metabarcoding; amplicon sequencing; cheese; Quebec; terroir; fungi; bacteria; microbiota

## 1. Introduction

Variations in cheese quality may be explained by fluctuations in milk composition [1,2] that affect coagulation properties [3], curd drainage [3] and yield [4]. Fluctuations in the activity of complex microbial cheese ecosystems can also influence product consistency [5–7]. Bacteria and fungi are responsible for the typical organoleptic properties of cheese. Characterization of the bacteria and fungi contained in cheese can therefore help reduce variations in cheese quality and reduce potential associated economic losses [8].

In 2012, over 600 different yeast and mold species were isolated from raw milk produced in Quebec and destined for cheesemaking [9], highlighting that fungi naturally present in milk can develop in cheese, potentially contributing to ripening. However, this study did not analyze bacteria or use the deep characterization that high-throughput sequencing (HTS) methods now allow. The decreasing cost of HTS profiling methods provides an unprecedented opportunity to describe and/or monitor complex microbial

<footnote n="1">¹ Department of Food Sciences and Nutrition, Institute of Nutrition and Functional Foods (INAF), STELA Dairy Research Centre, Université Laval, 2425 rue de l'Agriculture, Quebec City, QC G1V 0A6, Canada; annickRaymond@hotmail.com (A.R.-F.); marie-helene lessard@fsaa.ulaval.ca (M.-H.L.); julien.chamberland@fsaa.ulaval.ca (J.C.); yves.pouliot@fsaa.ulaval.ca (Y.P.); sylvie.turgeon@fsaa.ulaval.ca (S.L.T.); damilk83@hotmail.com (D.S.-G.)</footnote>

<footnote n="2">² UMR SayFood, INRAE, AgroParisTech, Université Paris-Saclay, Avenue Lucien Brétignières, 78850 Thiverval-Grignon, France; eric.dugat-bony@inrae.fr</footnote>

<footnote n="3">³ Agriculture and Agri-Food Canada, Saint-Hyacinthe Research and Development Center, 3600 Casavant Boulevard West, Saint-Hyacinthe, QC J2S 8E3, Canada</footnote>

<footnote>* Correspondence: steve.labrie@fsaa.ulaval.ca</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
communities. Many cheese ecosystems from different countries have been characterized using such HTS technologies [6,10–12]. This has allowed the detection of new microorganisms not previously associated with cheeses and emphasized the importance of the so-called secondary microbiota, including non-inoculated microbiota that develops during the late stages of ripening [5]. HTS methods such as high-throughput amplicon-based metagenomics, or metabarcoding, have made it possible to identify the dominant and subdominant microorganisms involved in an ecosystem, whether they are cultivable in the laboratory or not [13].

Metabarcoding is performed using phylogenetic markers such as the nine highly variable regions of the 16S rDNA (V1–V9) [6,11], the fungal internal transcribed spacers (ITS1 and ITS2) [14–16], and the 18S rDNA or the D1/D2 domain of the 26S rDNA [17]. To date, no consensus has been made about the best metabarcoding targets for the cheese ecosystem, even though they greatly influence the resulting microbial profiles [18,19]. Target region, amplicon length and possible primer mismatches have a direct impact on taxonomic assignment, therefore influencing the measured composition of the same microbial community [20,21]. For example, primer mismatches may induce a distorted taxonomic assessment due to the inability to detect certain groups or genera and the subsequent overestimation of others [22].

Dairy microbiota have previously been studied using the hypervariable regions of the bacterial 16S rDNA V1–V3 [10,23,24], V3–V4 [11,25,26], V4 [5–7], V5–V6 [13] and V6–V8 [27], and of the fungal ITS1 [6,7,12] and ITS2 [11,28]. Regions V3–V4 or V4 are used more frequently to characterize cheese ecosystems compared to V6–V8. Both ITS regions (ITS1 and ITS2) provided similar results for extended studies but yielded different taxonomic resolutions and genera detections according to the environment studied [18,19, 21]. It is, therefore, important to compare markers for specific ecosystems.

To monitor the bacterial and fungal microbiota of terroir cheese from Quebec, we first compared the performance of primers targeting two bacterial 16S rDNA regions (V3–V4 and V6–V8) and two fungal ITS regions (ITS1 and ITS2). We identified V3–V4 and ITS2 as the most appropriate molecular markers to characterize bacterial and fungal ecosystems, respectively, in bloomy-, natural- and washed-ripened specialty cheeses (soft, semi-hard and hard varieties). We evaluate the consistency of the microbial ecosystems in Quebec’s terroir cheese by comparing production years 2015 and 2018. This article reports the first evaluation of the microbial diversity of Quebec’s terroir cheese ecosystems and their variability across different production years.

## 2. Materials and Methods

**Cheese sampling.** A total of 47 ready-to-eat terroir cheeses from Quebec were sampled and analyzed (Table S1). Twenty-seven cheeses made in 2015 were purchased in a single copy from a specialized shop in Quebec City (Québec, QC, Canada). Twenty cheeses produced in 2018 were received in triplicate directly from the cheesemakers (Table S1). These 47 cheeses belong to 32 cheese varieties (15 cheeses from 2015 and 2018, 12 cheeses only from 2015 and 5 varieties only from 2018). Samples included bloomy-, washed- or natural-rind cheeses, the latter being washed with a diluted salt solution but without the addition of ripening starters.

The core (10 g) and rind (25 cm²) cheeses were sampled with sterile scalpels and analyzed separately (98 samples in total). Cryogenic grinding was performed using the CryoMill apparatus (Restch®, Haan, Germany) with liquid nitrogen to obtain a fine frozen cheese powder. The grinding step consisted of an automatic precooling step (>3 min) at a 5 Hz frequency followed by two cycles at 25 Hz for 2 min and a 30-sec cooling step at 5 Hz between the cycles. The samples were named according to the section of the cheese analyzed (c for center, r for rind) and the year of the sampling (e.g., c01y15; center of cheese 01, year 2015).

**DNA extraction, library construction and sequencing.** DNA extraction was performed using the PureLink™ Viral RNA/DNA Mini Kit (Invitrogen, Life Technologies,
Carlsbad, CA, USA), according to the manufacturer's protocol, with minor modifications. A 30–40 mg sample of ground cheese was homogenized in 200 µL of NaCl 0.9% before proteinase K and the lysis buffer were added, following the Purelink protocol. DNA concentration and quality were determined using the Nanodrop (ND-1000) spectrophotometer (ThermoFisher Scientific, Wilmington, NC, USA).

Targeted amplification from total genomic DNA, metabarcoding library construction using specific rDNA primers (Table S2), equimolar pooling, MiSeq high-throughput sequencing (Illumina, San Diego, CA, USA) and paired-end reads (2 × 300 bp) demultiplexing was performed at the Institut de Biologie Intégrative et des Systèmes (IBIS, Université Laval). The bacterial and fungal composition of the core and the rind were determined by targeting the rDNA 16S region V3–V4, respectively [29], or V6–V8 [30] and the Internal Transcribed Spacer 1 or 2 (ITS1, ITS2) [31].

**Bioinformatic analysis, metrics and figures.** Metabarcoding analysis was performed using FROGS with most default parameters [32,33]. Briefly, paired-end sequences were clustered using the swarm clustering method [34], with a maximum of 3 differences (d = 3) allowed between amplicons. All clusters were filtered for chimeras using VSEARCH [35]. Operational Taxonomic Units (OTUs) with an abundance inferior to 0.005% were discarded because they are more likely to be chimeras than rare clusters [36] and may lead to a false overestimation of sample richness. Finally, the affiliation of all bacterial OTUs was validated using the EzBioCloud database [37] and the sequences from type material in 16S rRNA and ITS databases for BlastN in NCBI [38]. The remaining chimeras identified were manually filtered, and these cured OTU tables were used to calculate the biostatistics using PhyloSeq R package v1.28.0, implemented in FROGS [39]. Alpha-diversity metrics (Chao1 and inverse Simpson indexes) were computed after rarefaction. The abundance-based dissimilarity (beta-diversity) index, the Bray–Curtis index [40], was also calculated after rarefaction. Histograms presented in this study were drawn using the R package ggplot2 v.2.2.1 library [41].

## 3. Results

### 3.1. Sequencing Data and Taxonomic Assignation

rDNA libraries were prepared separately for the core and the rind from a total of 32 different ready-to-eat cheeses (Table S1). Two different DNA targets were used for bacteria (V3–V4 and V6–V8 regions from the rDNA 16S) and fungi (ITS1 and ITS2 of the rDNA ITS region) to compare their ability to characterize cheese microbiota using different primer sets (Table S2). A total of 21,616,072 sequence reads were obtained through the Illumina MiSeq sequencing platform.

For the bacterial profiles, 4,505,810 and 6,156,446 paired-end reads were sequenced for the V3–V4 and V6–V8 regions, respectively. Overall, 2,588,028 and 3,451,507 non-chimeric reads were assigned, respectively, to 108 and 154 OTUs, having a relative abundance >0.005% (Table 1). The rDNA target V3–V4 allowed the identification of 57 genera, 22 of which were highly abundant (>10% relative abundance in one or more samples). The V6–V8 region allowed the identification of 82 genera, including 21 highly abundant genera. In total, 38 bacterial genera were detected by both 16S rDNA targets (Table 1).

<table><caption>**Table 1.** Sequencing data of the fungal (ITS1 and ITS2) and the bacterial (V3–V4 and V6–V8) ecosystems.</caption><thead><tr><th rowspan="2">Target Region</th><th colspan="2">Bacteria</th><th colspan="2">Fungi</th></tr><tr><th>V3–V4</th><th>V6–V8</th><th>ITS1</th><th>ITS2</th></tr></thead><tbody><tr><td>Nb. reads sequenced <sup>1</sup></td><td>4,505,810</td><td>6,156,446</td><td>5,207,694</td><td>5,746,122</td></tr><tr><td>Nb. trimmed and assembled sequences</td><td>3,904,825</td><td>5,482,054</td><td>4,931,241</td><td>5,675,978</td></tr><tr><td>Nb. non-chimeric OTU sequences <sup>1</sup></td><td>2,588,028</td><td>3,451,507</td><td>4,683,355</td><td>5,096,678</td></tr><tr><td>Non-chimeric sequences</td><td>66.28%</td><td>62.96%</td><td>94.97%</td><td>89.79%</td></tr><tr><td>Nb. non-chimeric OTUs <sup>1</sup></td><td>108</td><td>154</td><td>67</td><td>78</td></tr><tr><td>Nb. genera assigned <sup>1</sup>/Nb. shared genera assignation</td><td>57/37</td><td>83/37</td><td>34/18</td><td>35/18</td></tr><tr><td>Nb. abundant genera <sup>2</sup>/Nb. shared genera assignation</td><td>22/17</td><td>20/17</td><td>9/6</td><td>9/6</td></tr></tbody></table>

<footnote n="1">¹ with relative abundance >0.005%. ² abundant genera with relative abundance >10%.</footnote>
For the fungal profiles, 5,207,694 and 5,746,122 reads were sequenced for intergenic regions ITS1 and ITS2, respectively. Of this total, 4,683,355 and 5,096,678 non-chimeric sequences were assigned to 67 and 78 OTUs (>0.005%), respectively (Table 1). Among these OTUs, ITS1 and ITS2 allowed the identification of 31 and 34 fungal genera, including 18 genera identified using both regions (Table 1).

## 3.2. Choice of the Best Targets for Metabarcoding

A comparison of different rDNA regions for monitoring the bacterial microbiota in cheese provided similar profiles (Figure 1). At the order taxonomic level, the major taxa identified were *Lactobacillales* (70%), *Micrococcales* (8–10%), *Oceanospirillales*, *Pseudomonadales*, *Enterobacterales*, *Corynebacteriales*, *Bacillales* and *Rhizobiales* (Dataset S1). For bacteria, V3–V4 generated fewer chimeric sequences and fewer low-abundant OTUs than V6–V8 (33.7% compared to 37.0%) (Table 1). In addition, V6–V8 provided ten ambiguous profiles out of the 94 samples sequenced (c06y15, c28y15, r01y15, r03y15, r06y15, r18y15, r21y15, r22y15, r23y15, r15y15; see Section 2 for the sample nomenclature) that could be associated with background or unspecific sequencing. In contrast, only two ambiguous profiles were present using V3–V4 (Figure 1; r01y15, r15y15), therefore suggesting that V3–V4 was a better bacterial target for metabarcoding studies in cheese ecosystems. The two samples suspected to correspond to background or unspecific amplification, namely r01y15 and r15y15, were removed for further analysis.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Comparative distribution of the most abundant bacteria in cheese ecosystems. Each column shows the relative abundance of the bacterial microbiota (16S rDNA), representing over 1% of the average abundance (order <1% are combined and shown in black). Vertical sections show different 16S DNA targets: the top row shows V3–V4 data, and the bottom row shows V6–V8 data. Horizontal sections show different parts of the cheese according to the rind type (bloomy, natural or washed): the left row shows cheese core data, and the right row shows cheese rind data. The last number of the cheese ID refers to the year of production (2015 or 2018) for the 32 cheese varieties.</figcaption></figure>
Fungal metabarcoding performed using ITS1 and ITS2 regions (Figure 2) generated a low proportion of chimeric sequences and rare OTUs, corresponding to 5.0% and 10.2% of the sequences, respectively (Table 1). The fungal profiles were different between ITS1 and ITS2. For ITS1, *Saccharomycetales* (52%), *Hypocreales* (21%), *Eurotiales* (17%) and *Microascales* (9%) were detected as the most abundant orders, while in ITS2 profiles, most sequence reads were affiliated to *Saccharomycetales* (73%), *Eurotiales* (12%), *Hypocreales* (10%) and *Microascales* (4%) (Dataset S1). The large dominance of *Saccharomycetales* using ITS2 is in accordance with specific dairy fungi and supports other evidence indicating that it is the most appropriate target for the characterization of fungal communities in cheeses [31].

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Comparative distribution of the most abundant fungi in cheese ecosystems. Each column shows the relative abundance of the fungal (ITS) microbiota representing over 1% of the average abundance (order <1% are combined and shown in black). Vertical sections show different ITS targets: the top row shows ITS1 data, and the bottom row shows ITS2 data. Horizontal sections show different parts of the cheese according to the rind type (bloomy, natural or washed): the left row shows cheese core data, and the right row shows cheese rind data. The last number of the cheese ID refers to the year of production (2015 or 2018) for the 32 cheese varieties.</figcaption></figure>

### 3.3. Microbial Landscape of Quebec's Terroir Cheeses

Cheese core ecosystems were largely dominated by starter and ripening cultures [11]. Overall, these cheese cores ecosystems were composed of OTUs affiliated with 47 bacterial genera where the most abundant were *Lactococcus* (54%), *Streptococcus* (37%) and the *Lactobacillaceae* family (7%), the latter being present mostly in washed cheese cores (Dataset S1 and Figure 3a). The fungal OTUs were mostly affiliated with the genera *Geotrichum*, *Debaryomyces*, *Penicillium*, *Fusarium*, *Kluyveromyces*, *Cyberlindnera* and *Scopulariopsis* (Dataset S1 and Figure 3b). Fungal core ecosystems showed specific profiles related to rind type. *Geotrichum* (46%) and *Debaryomyces* (20%) were the most abundant yeasts in washed cheese cores, while the surface of the bloomy-rind cheeses was dominated by *Geotrichum* (52%) and *Penicillium* (28%) (Dataset S1 and Figure 3b).
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Comparative distribution of the most abundant genera in cheese ecosystems. Horizontal sections show different parts of the cheese according to the rind type (bloomy, natural or washed): the left row shows cheese core data, and the right row shows cheese rind data. The last number of the cheese ID refers to the year of production (2015 or 2018) for the 32 cheese varieties. (a) Each column shows the relative abundance of the bacterial (V3–V4) microbiota representing over 1% of the average abundance (order <1% are combined and shown in black). (b) Each column shows the relative abundance of the fungal (ITS2) microbiota representing over 1% of the average abundance (order <1% are combined and shown in black).</figcaption></figure>

α-diversity was assessed for all cheese samples studied (Table 2). The Chao index evaluates richness by estimating the number of OTUs in a community [42], and the inverse of the Simpson index is an α-diversity measure (within sample diversity) based on the number of OTU and their abundance [43]. In cheese cores, both the average richness and diversity indexes were similar for fungi and for bacteria, regardless of rind type. The maximum estimated OTU number reached 30 for fungi and 42 for bacteria, while
the maximum inverse of the Simpson α-diversity index reached 3.183 and 2.516 for both kingdoms, respectively (Tables 2 and S1).

<table><caption>**Table 2.** Richness and α-diversity of cheese microbiota.</caption><thead><tr><th rowspan="2">Cheese Section</th><th rowspan="2">Rind Type</th><th colspan="2">Bacteria</th><th colspan="2">Fungi</th></tr><tr><th>Richness <sup>1</sup></th><th>α-Diversity <sup>1</sup></th><th>Richness</th><th>α-Diversity</th></tr></thead><tbody><tr><td rowspan="3">Cheese core</td><td>Bloomy</td><td>16 ± 11</td><td>1.497 ± 0.517</td><td>23 ± 6</td><td>2.020 ± 0.829</td></tr><tr><td>Natural</td><td>16 ± 7</td><td>1.783 ± 0.430</td><td>21 ± 4</td><td>2.251 ± 0.932</td></tr><tr><td>Washed</td><td>24 ± 18</td><td>1.794 ± 0.722</td><td>24 ± 6</td><td>2.317 ± 0.800</td></tr><tr><td rowspan="3">Cheese rind</td><td>Bloomy</td><td>32 ± 8</td><td>3.026 ± 1.683</td><td>19 ± 10</td><td>1.745 ± 0.496</td></tr><tr><td>Natural</td><td>34 ± 16</td><td>2.732 ± 0.372</td><td>18 ± 10</td><td>2.039 ± 1.139</td></tr><tr><td>Washed</td><td>36 ± 10</td><td>5.045 ± 2.223</td><td>17 ± 5</td><td>1.974 ± 0.737</td></tr></tbody></table>

Cheese rind ecosystems were dominated by lactic acid bacteria (LAB), ripening strains and members of the secondary microbiota found in cheeses. On the surface of the cheese sampled, the average fungal richness (18 OTUs) and inverse Simpson index (1.919) were lower than the bacterial richness (34 OTUs) and inverse Simpson index (3.601), no matter the type of rind (Table 2). Overall, the cheese rind's ecosystems were composed of OTUs affiliated to 54 bacterial genera where the most abundant were *Lactococcus*, *Brevibacterium*, *Glutamicibacter*, *Corynebacterium*, *Halomonas*, *Brachybacterium*, *Staphylococcus*, *Psychrobacter*, *Serratia*, *Streptococcus*, *Pseudoalteromonas* and several other less abundant genera (Figure 3a, Dataset S1). The fungal OTUs were mostly affiliated with the genera *Geotrichum*, *Penicillium*, *Debaryomyces*, *Fusarium* and *Scopulariopsis* (Figure 3b, Dataset S1).

These rind ecosystems showed, at the genus taxonomic level, specific populations depending on the type of rind. The highest bacterial richness and α-diversity were observed in washed rinds (maximum richness of 46 OTUs and inverse Simpson index of 7.268), followed by bloomy and natural rinds (Table 2). The richness and inverse Simpson index of fungal cheese rinds were similar in all rinds and to those observed in cheese cores (Table 2).

Washed rinds were dominated by halotolerant psychrophilic and coryneform bacteria belonging to the genera *Halomonas* (20%), *Corynebacterium* (13%), *Glutamicibacter* (13%), *Psychrobacter* (12%), *Staphylococcus* (10%) and *Brevibacterium* (10%) (Figure 3a, Dataset S1). Their fungal ecosystem was mostly composed of *Geotrichum* (35%), *Fusarium* (28%), *Debaryomyces* (22%) and *Scopulariopsis* (14%) (Figure 3b, Dataset S1).

Bloomy rinds showed a different microbiota composed mostly of the fungi *Geotrichum* (58%), *Penicillium* (34%) and *Debaryomyces* (7%), and the bacteria *Lactococcus* (44%), *Serratia* (13%), *Streptococcus* (6%), *Glutamicibacter* (6%) and many less abundant genera (Figure 3a, Dataset S1).

### 3.4. Assessing the Stability of Quebec's Terroir Cheese Ecosystems

The bacterial (V3–V4) and fungal (ITS2) constancy of cheese ecosystems were evaluated by calculating the β-diversity of 15 cheeses produced in 2015 and in 2018. The β-diversity, which is the diversity observed between samples, was calculated separately for the cheese rinds and cores using the Bray–Curtis dissimilarity index [40] to compare these ecosystems between the two years studied. A 0.5 threshold value was set to distinguish similar (0.0–0.5) from dissimilar (0.5–1.0) communities (Figure 4, Table S1). Because it seemed that background or unspecific sequences were obtained for samples r01y15 and r15y15, they were not included in the comparisons. A total of 58 comparisons were considered.

More than half of the Bray–Curtis dissimilarity index calculated between the cheeses sampled in 2015 and 2018 (~65%) could be interpreted as similar ecosystems. Indeed, for the core ecosystems, 9/15 fungi and 13/15 bacteria, while 10/15 fungi and 6/13 bacteria of the rind ecosystems show a Bray–Curtis β-diversity index <0.5. A few ecosystems tend to present a greater dissimilarity. A total of three cheeses had a Bray–Curtis index >0.9 (ITS2c13, V3–V4r12 and V3–V4r13), and two cheeses had a Bray–Curtis index between 0.75 and 0.89 (ITS2c16 and V3–V4r06) (Figure 4, Table S1).

<footnote n="1">¹ Mean values (±Standard Deviation) of richness (estimated number of OTUs) and α-diversity calculated respectively using Chao1 and inverse of Simpson indexes.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 4. Similarity of cheese microbiota from different years according to the Bray–Curtis dis-similarity index.** The Bray–Curtis index shows the similarity (0–0.75) or dissimilarity (0.76–1) of the cheeses from different productions (2015 and 2018). The index has been calculated for fungi (ITS2) and bacteria (V3–V4) from cheese rind (▲) and cheese core (○). Numbers shown are cheese identification numbers for the bloomy- (blue), natural- (green) and washed-rind cheeses (orange).</figcaption></figure>

For the fungi monitored using the ITS2 region, two washed-rind cheese cores (c13 and c16) were among the most dissimilar ecosystems. The changes observed involved only a few yeast and mold genera. The most dissimilar fungal core ecosystem (Bray–Curtis index of 0.908) was observed in a washed-rind cheese core (c13), where proportions of *Cyberlindnera*, *Kluyveromyces* and *Geotrichum* changed from 64%, 22% and 11% to 84%, 14% and 2% respectively, from 2015 to 2018 (Figure 3b, Dataset S1). For c16 (Bray–Curtis 0.772), the major yeast genera found in 2015 were *Fusarium* (49%), *Kluyveromyces* (19%), *Geotrichum* (18%), *Penicillium* (8%) and *Debaryomyces* (3%). In 2018, these proportions switched to *Geotrichum* (87%), *Penicillium* (12%) and *Kluyveromyces* (1%) (Figure 3b, Dataset S1).

High dissimilarities in the bacterial ecosystems were observed on the surface of two washed-rind cheeses (r12 and r13) and one bloomy cheese (r06). For the bloomy cheese rind (r06), the Bray–Curtis index was 0.769, explained by the high abundance of the starter cultures *Lactococcus* (51%) and *Streptococcus* (38%) in 2015 and by halotolerant psychrophilic bacteria (*Psychrobacter*, *Glutamicibacter* and *Staphylococcus*; totaling 79%) in 2018 (Figure 3a, Dataset S1). The dissimilarities of the washed-rind cheeses involved mostly changes in proportions of the halotolerant psychrophilic bacteria and coryneform bacteria such as *Halomonas*, *Psychrobacter*, *Glutamicibacter*, *Corynebacterium*, *Brevibacterium* and *Staphylococcus* (Figure 3a, Dataset S1).

**4. Discussion**

**Optimal target choice for metabarcoding.** MiSeq sequencing technology was used in this study to analyze Quebec's specialty cheese microbiota composition. We first compared different molecular markers targeting bacteria (regions V3–V4 and V6–V8 of 16S rDNA) and fungi (intergenic regions ITS1 and ITS2). The accuracy of cheese microbial profiling of both regions targeted in rDNA 16S and ITS was evaluated based on several criteria. Overall, the total number of genera detected and the proportions of dominant genera (relative abundance > 10%) were very similar (Table 1). For bacteria, V3–V4 seems to generate 34% of chimeric sequences and low-abundant OTUs, compared to 37% of chimeric sequences for V6–V8 (Table 1). The chimera sequence rate is higher for bacteria than fungi, perhaps because 16S rDNA is highly conserved, promoting recombination during the PCR amplification step prior to library preparation [44,45]. Despite comparable results for
V3–V4 and V6–V8 with other cheese metabarcoding studies [5,6], a lower rate of chimeric
sequences favors V3–V4 as a preferable target for bacterial metabarcoding. Moreover,
V6–V8 gave 10 samples showing only background amplification profiles, compared to
2 samples for V3–V4 (Figure 1) [46,47]. Aside from possible technical problems in the
library preparation, the high number of bloomy-rind cheese rind samples without relevant
metabarcoding profiles could be explained by the generally low abundance of subdominant
bacteria on the surface of bloomy-rind cheeses, which would have to be verified using a
quantification method such as qPCR [7,32,48].

Fungal metabarcoding using ITS1 and ITS2 generated, respectively, 5% and 10% chimeric sequences and low abundant OTUs (Table 1). However, taxonomic assignment favored ITS2 as the best target due to its ability, contrarily to ITS1, to detect genera previously associated with dairy ecosystems such as *Cyberlindnera*, *Cryptococcus*, *Kazachstania*, *Pichia* and *Yarrowia* [9]. Moreover, a high number of polymorphisms have been described previously in *Geotrichum candidum* rDNA, specifically in the ITS1 region [49]. This observation, combined with the fact that the primers used for ITS1 amplification might not allow the detection of the taxonomic class *Saccharomycetes* due to 3' terminal mismatches with several other fungal genera [31,50,51], leads to the conclusion that ITS2 could provide a more accurate representation of the fungal microbiota for cheeses and the dairy environment [52,53].

**The microbial landscape of Quebec's terroir cheeses.** For the first time, the bacterial and fungal microbiota of Quebec's terroir cheeses were studied using metabarcoding. Cheese rinds and cores were analyzed separately because of their different microbial composition [54–56]. Due to their contact with the external environment, cheese rind ecosystems generally have a higher α-diversity than cores [57,58]. Here, when comparing ecosystems, the bacterial α-diversity (calculated using Chao1 and inverse Simpson indexes) was generally higher in rinds than in cores, whereas these indexes were similar for fungal populations [11]. Moreover, cheeses cores and rinds showed a higher α-diversity for washed-ripened cheeses than natural- and bloomy-ripened cheeses (Tables 2 and S1).

The predominant bacteria found in cheese cores belonged to the genera *Lactococcus*, *Streptococcus* and the *Lactobacillaceae* family (Figure 3a, Dataset S1). These lactic acid bacteria are commonly used as a starter or adjunct cultures and were found in similar relative abundances for natural- and washed-rind technologies (Figure 3a). In the cheese samples studied, *Lactobacillaceae* were mostly present in washed-rind cheese cores (19%). They might have been added as adjunct cultures for their sugar metabolism and autolytic and proteolytic capacities that change the texture and provide cheese flavor enhancements or as a starter in some hard cheese varieties [59]. They did not seem to be used either as starter or adjunct cultures in bloomy-rind soft cheeses [60,61]. In these samples, *Lactococcus* and *Streptococcus* were observed in various proportions. Soft bloomy cheeses can be produced using either traditional or stabilized cheesemaking technology, the latter requiring the use of thermophilic starters [62,63].

Cheese surface ecosystems are known to have specific microbial compositions depending on cheesemaking technology [5,6,11]. Bloomy rinds are mostly colonized by yeasts and molds [64], and the bacterial microbiota is generally not considered [6]. The high abundance of *Lactococcus* and other bacteria in the bloomy-cheese rinds analyzed may be explained by an inadvertent sampling of the cheese subsurface [65]. This would need additional and complementary characterization, such as qPCR quantification, to provide a more realistic profile of the surface bacterial ecosystems [32]. Washed rind cheese surfaces have been described abundantly in the literature. As other authors found, the washed-ripened cheese ecosystems contained coryneform bacteria (*Brachybacterium*, *Brevibacterium*, *Corynebacterium* and *Glutamicibacter*) and staphylococci, which were important contributors to the typicity of several washed-rind cheeses [66–68]. They may have been introduced through the personal and the indigenous ecosystem established within the cheesemaking plants [7,69]. As revealed by other studies, despite the common addition of smear bacteria starters in washed-ripened cheese (*Brevibacterium*, *Glutamicibacter*), the cheese rind microbiota was
rather composed of environmental Actinobacteria (Corynebacterium, Brachybacterium) and Gammaproteobacteria (Halomonas, Psychrobacter, Pseudoalteromonas). Their occurrence and relative proportions are known to be influenced by many factors such as salt concentration, pH and low temperature of ripened cheese [70–72]. Halotolerant psychrophilic bacteria are most likely introduced in cheese ecosystems through the brine, tools and surfaces rather than in raw milk [5,73–76].

Fungal communities of cheese cores and rinds showed similar composition for a given variety and were mainly composed of Geotrichum, Penicillium, Debaryomyces, Fusarium, Kluyveromyces, Scopulariopsis and Cyberlindnera (Dataset S1). These dairy yeasts and molds frequently occur in cheeses because they are inoculated as ripening cultures or originate from the indigenous microbiota of the milk or of the dairy plant environment [73,77–79]. They were found in various proportions, without obvious predominance in either rind type, except for Penicillium, which was more abundant in bloomy-rind cheeses and was well known to contribute to the development of their characteristic appearance and flavor (Figure 3b) [78]. In the washed-rind cheese ecosystems, Fusarium (28%) was the second most abundant fungi identified in this study (Dataset S1). Although most commonly found in bloomy-rind cheeses [65], Fusarium can be added to the smearing solutions to reduce rind stickiness [80]. Scopulariopsis (14%, Dataset S1), usually associated with wrapped washed- and natural-rind cheeses [81], can contribute to cheese flavor through ammonia production [82]. Although generally considered spoilage molds [83–85], Fusarium and Scopulariopsis could have been used as ripening agents for their technological properties [78].

The metabarcoding method used in this study allowed the detection of the major contributors to cheese ecosystems and the influence of technological practices, such as the choice of optimal ripening and adjunct cultures, on each cheese variety studied (bloomy-, natural- and washed-rind). These data also provide information on the natural microbiota that could develop during cheese ripening. However, using metabarcoding does not give a quantitative overview of the microbiota nor its impact on organoleptic properties. Further studies using complementary techniques, such as qPCR, texture profile analysis, proteolysis and lipolysis assessment, would be necessary to give a complete overview of these cheeses.

**Assessing the stability of Quebec's terroir cheese ecosystems.** This study analyzed the same cheese varieties (numbered 01 to 16) produced in different years for the first time. The persistence of the microbial communities was thus evaluated using the Bray–Curtis distance [40]. This dissimilarity index, which is a reliable index to point out changed environments, was calculated for 58 comparisons (bacteria and fungi located in the core or on the rind) based on abundance data [86]. Most cheeses showed high similarity of the microbiota from year to year. This may be explained in part by the high abundance of the starter and ripening cultures inoculated during the cheesemaking process. However, the optimized metabarcoding method is sensitive enough to monitor changes in some ecosystems (cheeses 06, 12, 13 and 16). Cheese ecosystems can be modulated by numerous factors, including variation in the indigenous microbiota composition of the milk, the use of heat treatment [5,87], the cheesemaking process (starter, ripening and adjunct cultures used, smearing solutions, brining, utensils and staff) [6,65,76] and the processing/ripening conditions (shelves, temperature, relative humidity, cheese packaging, ripening time and air composition) [1,5–7,65,88,89].

The highest variations between the two years studied were observed in the cheese ecosystems ITS2c13, V3–V4r12 and V3–V4r13, as indicated by a Bray–Curtis index >0.9 (Figure 4). The most dissimilar rind ecosystems were the washed-rind cheeses 12 and 13 (Bray–Curtis indexes 0.958 and 0.952, respectively). The surface of the washed-rind cheeses had a high bacterial α-diversity [6,11,67,77], and the dissimilarity observed between 2015 and 2018 was characterized by shifts in the relative abundance of halotolerant psychrophilic Gammaproteobacteria (Pseudoalteromonas, Halomonas and Psychrobacter), which, to our knowledge, are not added to cheese as ripening starters. The shifts observed between 2015 and 2018 were possibly influenced by a change in the humidity content of the cheese rind [6], affected by possible variations in the environment of the ripening chambers
(temperature, relative humidity and air velocity), the salinity of the washing (smear) solution used [75,76,90], the quality and renewal frequency of brine solutions [76,90,91], the environment in the ripening chamber and/or of the cheesemaking plant [7].

The changes monitored in the dissimilar fungal core ecosystem of cheese 16 (Bray–Curtis index 0.908) were mostly explained by the addition of ripening cultures. The fungal ecosystem described for the 2015 production of cheese 16 was mostly composed of *Cyberlindnera* (64%), *Kluyveromyces* (22%) and *Geotrichum* (11%), which could have been deliberately added as starters culture or could originate from the indigenous milk microflora [92]. *Cyberlindnera* is a facultative fermenting yeast rarely documented in dairy products and mostly found in cheese cores (Figure 3b) [93]. Its occurrence in Quebec’s terroir milk and cheese have been previously described using traditional microbiology [9] and could be explained by milk or brine contamination [94]. *Cyberlindera* produces ethyl acetate and acetaldehyde that could contribute to the cheese flavor [95]. In 2018, the core of the same cheese was largely dominated by the fungal starters *Geotrichum* and *Debaryomyces*, two yeasts frequently inoculated in washed-rind cheeses for their technological properties [77]. *Debaryomyces* can also be found in cheese from involuntary inoculation since it is widespread in the cheesemaking plant environment, especially in brine, due to its tolerance to high salt concentrations and low temperature [7,76,96,97]. These yeasts are therefore highly abundant in cheese cores.

**5. Conclusions**

In conclusion, metabarcoding provided for the first time the microbial profile of many different bloomy-, washed- and natural-rind terroir cheeses from Quebec, evaluating the richness and diversity of both the cores and rinds. This large-scale metabarcoding inventory proposes that V3–V4 and ITS2 are adequate molecular identifiers to explore cheese microbial diversity in depth. This approach leads to an accurate identification (genus-level) within a few hours of analysis using FROGS and allows an understanding of the dissimilarities observed in cheeses through time. This study highlights the persistence of key bacteria and fungi in the Quebec terroir cheeses. More in-depth microbiological analysis could provide biological indicators to help keep high-quality products and reduce specialty cheese spoilage. The level of metabolic activities or quantitative determination of the microbial load was not evaluated and should be pursued in the future.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/microorganisms10071381/s1, Dataset S1: Relative abundance of bacteria and fungi in all cheese core and rind samples. Compared distribution of the bacteria and fungi using all metabarcoding targets (V3–V4, V6–V8, ITS1, ITS2); Table S1: Cheese samples description and biostatistics describing cheese microflora. The richness (estimated OTUs) has been calculated with the Chao1 index, the α-diversity is evaluated with the inverse of Simpson index and the β-diversity is evaluated for the same sample, at different year using Bray–Curtis dissimilarity index, Table S2: Primers sequences for the preparation of the samples. Name of primers or adapters, target, DNA region and sequences.

**Author Contributions:** Conceptualization, S.L.T., D.S.-G. and S.L.; methodology, A.R.-F.; validation, M.-H.L., E.D.-B. and S.L.; formal analysis, A.R.-F. and M.-H.L.; investigation, A.R.-F.; resources, S.L.; writing—original draft preparation, A.R.-F. and M.-H.L.; writing—review and editing, M.-H.L., J.C., E.D.-B., D.S.-G. and S.L.; visualization, A.R.-F. and J.C.; supervision, M.-H.L., E.D.-B., D.S.-G. and S.L.; project administration, M.-H.L. and S.L.; funding acquisition, S.L.T., D.S.-G., Y.P. and S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by Novalait Inc., the Fonds de recherche du Québec—Nature et technologies (FRQ-NT), the Quebec Consortium for Industrial Bioprocess Research and Innovation (CRIBIQ) (Grant ID 2016-LG-195927) and the Natural Sciences and Engineering Research Council of Canada Discovery Grant (NSERC- RGPIN-2017-06388). Eric Dugat-Bony was supported by the AgreenSkills + fellowship program which received funding from the EU’s Seventh Framework Programme (Grant ID FP7-609398, AgreenSkills + contract).
**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Raw demultiplexed sequence data were deposited in the NCBI Sequence Read Archive under the accession number SRP159168.

**Acknowledgments:** The authors are thankful to Laval University's Genomic Analysis Platform (IBIS) and Calcul Quebec for their assistance.

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

## References

1. O'Sullivan, D.J.; Cotter, P.D.; O'Sullivan, O.; Giblin, L.; McSweeney, P.L.; Sheehan, J.J. Temporal and spatial differences in microbial composition during the manufacture of a continental-type cheese. *Appl. Environ. Microbiol.* **2015**, *81*, 2525–2533. [CrossRef] [PubMed]

2. Amenu, B.; Deeth, H.C. The impact of milk composition on Cheddar cheese manufacture. *Aust. J. Dairy Technol.* **2007**, *62*, 171–180.

3. Daviau, C.; Famelart, M.-H.; Pierre, A.; Goudédranche, H.; Maubois, J.-L. Rennet coagulation of skim milk and curd drainage: Effect of pH, casein concentration, ionic strength and heat treatment. *Le Lait* **2000**, *80*, 397–415. [CrossRef]

4. Fenelon, M.A.; Guinee, T.P. The effect of milk fat on Cheddar cheese yield and its prediction, using modifications of the Van Slyke cheese yield formula. *J. Dairy Sci.* **1999**, *82*, 2287–2299. [CrossRef]

5. Quigley, L.; O'Sullivan, O.; Beresford, T.P.; Ross, R.P.; Fitzgerald, G.F.; Cotter, P.D. High-throughput sequencing for detection of subpopulations of bacteria not previously associated with artisanal cheeses. *Appl. Environ. Microbiol.* **2012**, *78*, 5717–5723. [CrossRef] [PubMed]

6. Wolfe, B.E.; Button, J.E.; Santarelli, M.; Dutton, R.J. Cheese rind communities provide tractable systems for in situ and in vitro studies of microbial diversity. *Cell* **2014**, *158*, 422–433. [CrossRef]

7. Bokulich, N.A.; Mills, D.A. Facility-specific “house” microbiome drives microbial landscapes of artisanal cheesemaking plants. *Appl. Environ. Microbiol.* **2013**, *79*, 5214–5223. [CrossRef]

8. Sperber, W.H. Introduction to the microbiological spoilage of foods and beverages. In *Compendium of the Microbiological Spoilage of Foods and Beverages*, 1st ed.; Sperber, W., Doyle, M., Eds.; Springer: New York, NY, USA, 2009; pp. 1–40. [CrossRef]

9. Lavoie, K.; Touchette, M.; St-Gelais, D.; Labrie, S. Characterization of the fungal microflora in raw milk and specialty cheeses of the province of Québec. *Dairy Sci. Technol.* **2012**, *92*, 455–468. [CrossRef]

10. Delcenserie, V.; Taminiau, B.; Delhalle, L.; Nezer, C.; Doyen, P.; Crevecoeur, S.; Roussey, D.; Korsak, N.; Daube, G. Microbiota characterization of a Belgian protected designation of origin cheese, Herve cheese, using metagenomic analysis. *J. Dairy Sci.* **2014**, *97*, 6046–6056. [CrossRef]

11. Dugat-Bony, E.; Garnier, L.; Denonfoux, J.; Ferreira, S.; Sarthou, A.S.; Bonnarme, P.; Irlinger, F. Highlighting the microbial diversity of 12 French cheese varieties. *Int. J. Food Microbiol.* **2016**, *238*, 265–273. [CrossRef]

12. Zheng, X.; Liu, F.; Shi, X.; Wang, B.; Li, K.; Li, B.; Zhuge, B. Dynamic correlations between microbiota succession and flavor development involved in the ripening of Kazak artisanal cheese. *Food Res. Int.* **2018**, *105*, 733–742. [CrossRef] [PubMed]

13. Alegria, A.; Szczesny, P.; Mayo, B.; Bardowski, J.; Kowalczyk, M. Biodiversity in Oscypek, a traditional Polish cheese, determined by culture-dependent and -independent approaches. *Appl. Environ. Microbiol.* **2012**, *78*, 1890–1898. [CrossRef] [PubMed]

14. Schoch, C.L.; Seifert, K.A.; Huhndorf, S.; Robert, V.; Spouge, J.L.; Levesque, C.A.; Chen, W.; Fungal Barcoding, C. Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. *Proc. Natl. Acad. Sci. USA* **2012**, *109*, 6241–6246. [CrossRef] [PubMed]

15. Bazzicalupo, A.L.; Bálint, M.; Schmitt, I. Comparison of ITS1 and ITS2 rDNA in 454 sequencing of hyperdiverse fungal communities. *Fungal Ecol.* **2013**, *6*, 102–109. [CrossRef]

16. Blaalid, R.; Kumar, S.; Nilsson, R.H.; Abarenkov, K.; Kirk, P.M.; Kauserud, H. ITS1 versus ITS2 as DNA metabarcodes for fungi. *Mol. Ecol. Resour.* **2013**, *13*, 218–224. [CrossRef]

17. De Filippis, F.; Laiola, M.; Blaiotta, G.; Ercolini, D. Different amplicon targets for sequencing-based studies of fungal diversity. *Appl. Environ. Microbiol.* **2017**, *83*, e00905–e00917. [CrossRef]

18. Mbareche, H.; Veillette, M.; Bilodeau, G.J.; Duchaine, C. An amplicon-based sequencing approach for the study of aeromycology. *J. Xenobiotics* **2018**, *8*, 8–7810. [CrossRef]

19. Mbareche, H.; Veillette, M.; Bilodeau, G.J.; Duchaine, C. Comparison of the performance of ITS1 and ITS2 as barcodes in amplicon-based sequencing of bioaerosols. *PeerJ* **2020**, *8*, e8523. [CrossRef]

20. Mello, A.; Napoli, C.; Murat, C.; Morin, E.; Marceddu, G.; Bonfante, P. ITS-1 versus ITS-2 pyrosequencing: A comparison of fungal populations in truffle grounds. *Mycologia* **2011**, *103*, 1184–1193. [CrossRef]

21. Monard, C.; Gantner, S.; Stenlid, J. Utilizing ITS1 and ITS2 to study environmental fungal diversity using pyrosequencing. *FEMS Microbiol. Ecol.* **2013**, *84*, 165–175. [CrossRef]

22. Poirier, S.; Rué, O.; Coeuret, G.; Champomier-Vergès, M.-C.; Loux, V.; Chaillou, S. Detection of an amplification bias associated to *Leuconostocaceae* family with a universal primer routinely used for monitoring microbial community structures within food products. *BMC Res. Notes* **2018**, *11*, 802. [CrossRef] [PubMed]
23. Ceugniez, A.; Taminiau, B.; Coucheney, F.; Jacques, P.; Delcenserie, V.; Daube, G.; Drider, D. Use of a metagenetic approach to monitor the bacterial microbiota of "Tomme d'Orchies" cheese during the ripening process. *Int. J. Food Microbiol.* **2017**, *247*, 65–69. [CrossRef] [PubMed]

24. Alessandria, V.; Ferrocino, I.; De Filippis, F.; Fontana, M.; Rantsiou, K.; Ercolini, D.; Cocolin, L. Microbiota of an italian Grana-like cheese during manufacture and ripening, unraveled by 16S rRNA-based approaches. *Appl. Environ. Microbiol.* **2016**, *82*, 3988–3995. [CrossRef] [PubMed]

25. Masoud, W.; Wik, M.; Vogensen, F.; Lillevang, S.; Abu Al-Soud, W.; Sørensen, S.; Jakobsen, M. Characterization of bacterial populations in Danish raw milk cheeses made with different starter cultures by denaturing gradient gel electrophoresis and pyrosequencing. *Int. Dairy J.* **2011**, *21*, 142–148. [CrossRef]

26. Riquelme, C.; Camara, S.; Dapkevicius, M.; Vinuesa, P.; da Silva, C.C.; Malcata, F.X.; Rego, O.A. Characterization of the bacterial biodiversity in Pico cheese (an artisanal Azorean food). *Int. J. Food Microbiol.* **2015**, *192*, 86–94. [CrossRef]

27. Chamberland, J.; Lessard, M.-H.; Doyen, A.; Labrie, S.; Pouliot, Y. A sequencing approach targeting the 16S rRNA gene unravels the biofilm composition of spiral-wound membranes used in the dairy industry. *Dairy Sci. Technol.* **2017**, *96*, 827–843. [CrossRef]

28. Ceugniez, A.; Taminiau, B.; Coucheney, F.; Jacques, P.; Delcenserie, V.; Daube, G.; Drider, D. Fungal diversity of "Tomme d'Orchies" cheese during the ripening process as revealed by a metagenomic study. *Int. J. Food Microbiol.* **2017**, *258*, 89–93. [CrossRef]

29. Klindworth, A.; Pruesse, E.; Schweer, T.; Peplies, J.; Quast, C.; Horn, M.; Glockner, F.O. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. *Nucleic Acids Res.* **2013**, *41*, e1. [CrossRef]

30. Comeau, A.M.; Li, W.K.; Tremblay, J.E.; Carmack, E.C.; Lovejoy, C. Arctic Ocean microbial community structure before and after the 2007 record sea ice minimum. *PLoS ONE* **2011**, *6*, e27492. [CrossRef]

31. Tedersoo, L.; Anslan, S.; Bahram, M.; Polme, S.; Riit, T.; Liiv, I.; Koljalg, U.; Kisand, V.; Nilsson, R.H.; Hildebrand, F.; et al. Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi. *MycoKeys* **2015**, *10*, 1–43. [CrossRef]

32. Chamberland, J.; Messier, T.; Dugat-Bony, E.; Lessard, M.-H.; Labrie, S.; Doyen, A.; Pouliot, Y. Influence of feed temperature to biofouling of ultrafiltration membrane during skim milk processing. *Int. Dairy J.* **2019**, *93*, 99–105. [CrossRef]

33. Escudié, F.; Auer, L.; Bernard, M.; Mariadassou, M.; Cauquil, L.; Vidal, K.; Maman, S.; Hernandez-Raquet, G.; Combes, S.; Pascal, G. FROGS: Find, Rapidly, OTUs with Galaxy Solution. *Bioinformatics* **2018**, *34*, 1287–1294. [CrossRef] [PubMed]

34. Mahé, F.; Rognes, T.; Quince, C.; de Vargas, C.; Dunthorn, M. Swarm: Robust and fast clustering method for amplicon-based studies. *PeerJ* **2014**, *2*, e593. [CrossRef]

35. Rognes, T.; Flouri, T.; Nichols, B.; Quince, C.; Mahé, F. VSEARCH: A versatile open source tool for metagenomics. *PeerJ* **2016**, *4*, e2584. [CrossRef]

36. Bokulich, N.A.; Subramanian, S.; Faith, J.J.; Gevers, D.; Gordon, J.I.; Knight, R.; Mills, D.A.; Caporaso, J.G. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. *Nat. Methods* **2013**, *10*, 57–59. [CrossRef] [PubMed]

37. Yoon, S.-H.; Ha, S.-M.; Kwon, S.; Lim, J.; Kim, Y.; Seo, H.; Chun, J. Introducing EzBioCloud: A taxonomically united database of 16S rRNA gene sequences and whole-genome assemblies. *Int. J. Syst. Evol. Microbiol.* **2017**, *67*, 1613–1617. [CrossRef]

38. Altschul, S.F.; Gish, W.; Miller, W.; Myers, E.W.; Lipman, D.J. Basic local alignment search tool. *J. Mol. Biol.* **1990**, *215*, 403–410. [CrossRef]

39. McMurdie, P.J.; Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. *PLoS ONE* **2013**, *8*, e61217. [CrossRef]

40. Bray, J.R.; Curtis, J.T. An ordination of the upland forest communities of southern Wisconsin. *Ecol. Monogr.* **1957**, *27*, 325–349. [CrossRef]

41. Wickham, H. *ggplot2: Elegant Graphics for Data Analysis*; Springer: New York, NY, USA, 2016. [CrossRef]

42. Chao, A. Nonparametric estimation of the number of classes in a population. *Scand. J. Stat.* **1984**, *11*, 265–270.

43. Simpson, E.H. Measurement of diversity. *Nature* **1949**, *163*, 688. [CrossRef]

44. Haas, B.J.; Gevers, D.; Earl, A.M.; Feldgarden, M.; Ward, D.V.; Giannoukos, G.; Ciulla, D.; Tabbaa, D.; Highlander, S.K.; Sodergren, E.; et al. Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons. *Genome Res.* **2011**, *21*, 494–504. [CrossRef] [PubMed]

45. Bjørnsgaard Aas, A.; Davey, M.L.; Kauserud, H. ITS all right mama: Investigating the formation of chimeric sequences in the ITS 2 region by DNA metabarcoding analyses of fungal mock communities of different complexities. *Mol. Ecol. Resour.* **2017**, *17*, 730–741. [CrossRef] [PubMed]

46. Wang, X.; Ren, H.; Zhan, Y. Characterization of microbial community composition and pathogens risk assessment in typical Italian-style salami by high-throughput sequencing technology. *Food Sci. Biotechnol.* **2018**, *27*, 241–249. [CrossRef] [PubMed]

47. Quince, C.; Lanzén, A.; Curtis, T.P.; Davenport, R.J.; Hall, N.; Head, I.M.; Read, L.F.; Sloan, W.T. Accurate determination of microbial diversity from 454 pyrosequencing data. *Nat. Methods* **2009**, *6*, 639–641. [CrossRef]

48. Chamberland, J.; Beaulieu-Carbonneau, G.; Lessard, M.H.; Labrie, S.; Bazinet, L.; Doyen, A.; Pouliot, Y. Effect of membrane material chemistry and properties on biofouling susceptibility during milk and cheese whey ultrafiltration. *J. Membr. Sci.* **2017**, *542*, 208–216. [CrossRef]

49. Alper, I.; Frenette, M.; Labrie, S. Ribosomal DNA polymorphisms in the yeast *Geotrichum candidum*. *Fungal Biol.* **2011**, *115*, 1259–1269. [CrossRef]
50. Tedersoo, L.; Ramirez, K.S.; Nilsson, R.H.; Kaljuvee, A.; Koljalg, U.; Abarenkov, K. Standardizing metadata and taxonomic identification in metabarcoding studies. *Gigascience* **2015**, *4*, 34. [CrossRef]

51. Nilsson, R.H.; Anslan, S.; Bahram, M.; Wurzbacher, C.; Baldrian, P.; Tedersoo, L. Mycobiome diversity: High-throughput sequencing and identification of fungi. *Nat. Rev. Microbiol.* **2019**, *17*, 95–109. [CrossRef]

52. Taylor, D.L.; Walters, W.A.; Lennon, N.J.; Bochicchio, J.; Krohn, A.; Caporaso, J.G.; Pennanen, T. Accurate estimation of fungal diversity and abundance through improved lineage-specific primers optimized for Illumina amplicon sequencing. *Appl. Environ. Microbiol.* **2016**, *82*, 7217–7226. [CrossRef]

53. Bellemain, E.; Carlsen, T.; Brochmann, C.; Coissac, E.; Taberlet, P.; Kauserud, H. ITS as an environmental DNA barcode for fungi: An in silico approach reveals potential PCR biases. *BMC Microbiol.* **2010**, *10*, 189. [CrossRef] [PubMed]

54. Ogier, J.-C.; Son, O.; Gruss, A.; Tailliez, P.; Delacroix-Buchet, A. Identification of the bacterial microflora in dairy products by temporal temperature gradient gel electrophoresis. *Appl. Environ. Microbiol.* **2002**, *68*, 3691–3701. [CrossRef]

55. Leclercq-Perlat, M.N.; Buono, F.; Lambert, D.; Latrille, E.; Spinnler, H.E.; Corrieu, G. Controlled production of Camembert-type cheeses. Part I: Microbiological and physicochemical evolutions. *J. Dairy Res.* **2004**, *71*, 346–354. [CrossRef] [PubMed]

56. Gobbetti, M.; Burzigotti, R.; Smacchi, E.; Corsetti, A.; De Angelis, M. Microbiology and biochemistry of gorgonzola cheese during ripening. *Int. Dairy J.* **1997**, *7*, 519–529. [CrossRef]

57. Giraffa, G.; Rossetti, L.; Mucchetti, G.; Addeo, F.; Neviani, E. Influence of the temperature gradient on the growth of thermophilic Lactobacilli used as natural starters in Grana cheese. *J. Dairy Sci.* **1998**, *81*, 31–36. [CrossRef]

58. Almena-Aliste, M.; Mietton, B. Cheese classification, characterization, and categorization: A global perspective. *Microbiol. Spectr.* **2014**, *2*, CM-0003-2012. [CrossRef]

59. Broadbent, J.R.; Cai, H.; Larsen, R.L.; Hughes, J.E.; Welker, D.L.; De Carvalho, V.G.; Tompkins, T.A.; Ardo, Y.; Vogensen, F.; De Lorentiis, A.; et al. Genetic diversity in proteolytic enzymes and amino acid metabolism among *Lactobacillus helveticus* strains. *J. Dairy Sci.* **2011**, *94*, 4313–4328. [CrossRef]

60. Parente, E.; Cogan, T.M.; Powell, I.B. Chapter 8—Starter Cultures: General Aspects. In *Cheese*, 4th ed.; McSweeney, P.L.H., Fox, P.F., Cotter, P.D., Everett, D.W., Eds.; Academic Press: San Diego, CA, USA, 2017; pp. 201–226.

61. Irlinger, F.; Helinck, S.; Jany, J.L. Chapter 11—Secondary and Adjunct Cultures. In *Cheese*, 4th ed.; McSweeney, P.L.H., Fox, P.F., Cotter, P.D., Everett, D.W., Eds.; Academic Press: San Diego, CA, USA, 2017; pp. 273–300.

62. Auclair, J.; Accolas, J.P. Use of thermophilic lactic starters in the dairy industry. *Antonie Leeuwenhoek* **1983**, *49*, 313–326. [CrossRef]

63. Lawrence, R.C.; Creamer, L.K.; Gilles, J. Texture development during cheese ripening. *J. Dairy Sci.* **1987**, *70*, 1748–1760. [CrossRef]

64. Lessard, M.H.; Belanger, G.; St-Gelais, D.; Labrie, S. The composition of Camembert cheese ripening cultures modulates both mycelial growth and appearance. *Appl. Environ. Microbiol.* **2012**, *78*, 1813–1819. [CrossRef]

65. Irlinger, F.; Layec, S.; Helinck, S.; Dugat-Bony, E. Cheese rind microbial communities: Diversity, composition and origin. *FEMS Microbiol. Lett.* **2015**, *362*, 1–11. [CrossRef] [PubMed]

66. Mounier, J.; Rea, M.C.; O’Connor, P.M.; Fitzgerald, G.F.; Cogan, T.M. Growth characteristics of *Brevibacterium*, *Corynebacterium*, *Microbacterium*, and *Staphylococcus* spp. isolated from surface-ripened cheese. *Appl. Environ. Microbiol.* **2007**, *73*, 7732–7739. [CrossRef] [PubMed]

67. Brennan, N.M.; Brown, R.; Goodfellow, M.; Ward, A.C.; Beresford, T.P.; Simpson, P.J.; Fox, P.F.; Cogan, T.M. *Corynebacterium mooreparkense* sp. nov. and *Corynebacterium casei* sp. nov., isolated from the surface of a smear-ripened cheese. *Int. J. Syst. Evol. Microbiol.* **2001**, *51*, 843–852. [CrossRef]

68. Rea, M.C.; Gorges, S.; Gelsomino, R.; Brennan, N.M.; Mounier, J.; Vancanneyt, M.; Scherer, S.; Swings, J.; Cogan, T.M. Stability of the biodiversity of the surface consortia of Gubbeen, a red-smear cheese. *J. Dairy Sci.* **2007**, *90*, 2200–2210. [CrossRef]

69. Grice, E.A.; Segre, J.A. The skin microbiome. *Nat. Rev. Microbiol.* **2011**, *9*, 244–253. [CrossRef]

70. Goerges, S.; Mounier, J.; Rea, M.C.; Gelsomino, R.; Heise, V.; Beduhn, R.; Cogan, T.M.; Vancanneyt, M.; Scherer, S. Commercial ripening starter microorganisms inoculated into cheese milk do not successfully establish themselves in the resident microbial ripening consortia of a South german red smear cheese. *Appl. Environ. Microbiol.* **2008**, *74*, 2210–2217. [CrossRef] [PubMed]

71. Feurer, C.; Vallaeys, T.; Corrieu, G.; Irlinger, F. Does smearing inoculum reflect the bacterial composition of the smear at the end of the ripening of a French soft, red-smear cheese? *J. Dairy Sci.* **2004**, *87*, 3189–3197. [CrossRef]

72. Irlinger, F.; Mounier, J. Microbial interactions in cheese: Implications for cheese quality and safety. *Curr. Opin. Biotechnol.* **2009**, *20*, 142–148. [CrossRef]

73. Quijada, N.M.; Mann, E.; Wagner, M.; Rodriguez-Lazaro, D.; Hernandez, M.; Schmitz-Esser, S. Autochthonous facility-specific microbiota dominates washed-rind Austrian hard cheese surfaces and its production environment. *Int. J. Food Microbiol.* **2018**, *267*, 54–61. [CrossRef]

74. Ishikawa, M.; Kodama, K.; Yasuda, H.; Okamoto-Kainuma, A.; Koizumi, K.; Yamasato, K. Presence of halophilic and alkaliphilic lactic acid bacteria in various cheeses. *Lett. Appl. Microbiol.* **2007**, *44*, 308–313. [CrossRef]

75. Monnet, C.; Landaud, S.; Bonnarme, P.; Swennen, D. Growth and adaptation of microorganisms on the cheese surface. *FEMS Microbiol. Lett.* **2015**, *362*, 1–9. [CrossRef] [PubMed]

76. Mounier, J.; Goerges, S.; Gelsomino, R.; Vancanneyt, M.; Vandemeulebroecke, K.; Hoste, B.; Brennan, N.M.; Scherer, S.; Swings, J.; Fitzgerald, G.F.; et al. Sources of the adventitious microflora of a smear-ripened cheese. *J. Appl. Microbiol.* **2006**, *101*, 668–681. [CrossRef] [PubMed]
77. Mounier, J.; Coton, M.; Irlinger, F.; Landaud, S.; Bonnarme, P. Chapter 38—Smear-ripened cheeses. In *Cheese*, 4th ed.; McSweeney, P.L.H., Fox, P.F., Cotter, P.D., Everett, D.W., Eds.; Academic Press: San Diego, CA, USA, 2017; pp. 955–996.

78. Spinnler, H.E. Chapter 36—Surface mold-ripened cheeses. In *Cheese*, 4th ed.; McSweeney, P.L.H., Fox, P.F., Cotter, P.D., Everett, D.W., Eds.; Academic Press: San Diego, CA, USA, 2017; pp. 911–928.

79. Cotter, P.D.; Beresford, T.P. Chapter 15—Microbiome changes during ripening. In *Cheese*, 4th ed.; McSweeney, P.L.H., Fox, P.F., Cotter, P.D., Everett, D.W., Eds.; Academic Press: San Diego, CA, USA, 2017; pp. 389–409.

80. Bachmann, H.P.; Bobst, C.; Bütikofer, U.; Casey, M.G.; Dalla Torre, M.; Fröhlich-Wyder, M.T.; Fürst, M. Occurrence and significance of *Fusarium domesticum* alias *Anticollanti* on smear-ripened cheeses. *LWT-Food Sci. Technol.* **2005**, *38*, 399–407. [CrossRef]

81. Marcellino OSB, S.N.; Benson, D.R. The good, the bad, and the ugly: Tales of mold-ripened cheese. *Microbiol. Spectr.* **2013**, *1*, CM-0005-0012. [CrossRef]

82. Bothast, R.J.; Lancaster, E.B.; Hesseltine, C.W. *Scopulariopsis brevicaulis*: Effect of pH and substrate on growth. *Eur. J. Appl. Microbiol. Biotechnol.* **1975**, *1*, 55–66. [CrossRef]

83. Houbraken, J.; Samson, R.A. Current taxonomy and identification of foodborne fungi. *Curr. Opin. Food Sci.* **2017**, *17*, 84–88. [CrossRef]

84. Woudenberg, J.H.C.; Meijer, M.; Houbraken, J.; Samson, R.A. *Scopulariopsis* and scopulariopsis-like species from indoor environments. *Stud. Mycol.* **2017**, *88*, 1–35. [CrossRef]

85. Geiser, D.M.; Aoki, T.; Bacon, C.W.; Baker, S.E.; Bhattacharyya, M.K.; Brandt, M.E.; Brown, D.W.; Burgess, L.W.; Chulze, S.; Coleman, J.J. One fungus, one name: Defining the genus *Fusarium* in a scientifically robust way that preserves longstanding use. *Phytopathology* **2013**, *103*, 400–408. [CrossRef]

86. Clarke, K.R.; Somerfield, P.J.; Chapman, M.G. On resemblance measures for ecological studies, including taxonomic dissimilarities and a zero-adjusted Bray–Curtis coefficient for denuded assemblages. *J. Exp. Mar. Biol. Ecol.* **2006**, *330*, 55–80. [CrossRef]

87. Quigley, L.; McCarthy, R.; O'Sullivan, O.; Beresford, T.P.; Fitzgerald, G.F.; Ross, R.P.; Stanton, C.; Cotter, P.D. The microbial content of raw and pasteurized cow milk as determined by molecular approaches. *J. Dairy Sci.* **2013**, *96*, 4928–4937. [CrossRef]

88. Dolci, P.; Alessandria, V.; Rantsiou, K.; Bertolino, M.; Cocolin, L. Microbial diversity, dynamics and activity throughout manufacturing and ripening of Castelmagno PDO cheese. *Int. J. Food Microbiol.* **2010**, *143*, 71–75. [CrossRef] [PubMed]

89. Petersen, K.M.; Westall, S.; Jespersen, L. Microbial succession of *Debaryomyces hansenii* strains during the production of Danish surfaced-ripened cheeses. *J. Dairy Sci.* **2002**, *85*, 478–486. [CrossRef]

90. Bockelmann, W.; Willems, K.P.; Neve, H.; Heller, K.H. Cultures for the ripening of smear cheeses. *Int. Dairy J.* **2005**, *15*, 719–732. [CrossRef]

91. Mounier, J.; Gelsomino, R.; Goerges, S.; Vancanneyt, M.; Vandemeulebroecke, K.; Hoste, B.; Scherer, S.; Swings, J.; Fitzgerald, G.F.; Cogan, T.M. Surface microflora of four smear-ripened cheeses. *Appl. Environ. Microbiol.* **2005**, *71*, 6489–6500. [CrossRef]

92. Bertuzzi, A.S.; Walsh, A.M.; Sheehan, J.J.; Cotter, P.D.; Crispie, F.; McSweeney, P.L.H.; Kilcawley, K.N.; Rea, M.C. Omics-based insights into flavor development and microbial succession within surface-ripened cheese. *MSystems* **2018**, *3*, e00211–e00217. [CrossRef]

93. Cogan, T.M.; Goerges, S.; Gelsomino, R.; Larpin, S.; Hohenegger, M.; Bora, N.; Jamet, E.; Rea, M.C.; Mounier, J.; Vancanneyt, M.; et al. Biodiversity of the Surface Microbial Consortia from Limburger, Reblochon, Livarot, Tilsit, and Gubbeen Cheeses. *Microbiol. Spectr.* **2014**, *2*, CM-0010-2012. [CrossRef]

94. Büchl, N.R.; Seiler, H. Yeasts and molds | Yeasts in milk and dairy products. In *Encyclopedia of Dairy Sciences*, 2nd ed.; Fuquay, J.W., Ed.; Academic Press: San Diego, CA, USA, 2011; pp. 744–753.

95. Buerth, C.; Tielker, D.; Ernst, J.F. *Candida utilis* and *Cyberlindnera (Pichia) jadinii*: Yeast relatives with expanding applications. *Appl. Microbiol. Biotechnol.* **2016**, *100*, 6981–6990. [CrossRef]

96. Butinar, L.; Santos, S.; Spencer-Martins, I.; Oren, A.; Gunde-Cimerman, N. Yeast diversity in hypersaline habitats. *FEMS Microbiol. Lett.* **2005**, *244*, 229–234. [CrossRef]

97. Breuer, U.; Harms, H. *Debaryomyces hansenii*—An extremophilic yeast with biotechnological potential. *Yeast* **2006**, *23*, 415–437. [CrossRef]