# Overview of a surface-ripened cheese community functioning by meta-omics analyses

Eric Dugat-Bony, Cecile Straub, Aurélie Teissandier, Djamila Onesime, Valentin Loux, Christophe Monnet, Francoise Irlinger, Sophie Landaud-Liautaud, Marie Noelle Leclercq-Perlat, Pascal Bento, et al.

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

## ► To cite this version:

Eric Dugat-Bony, Cecile Straub, Aurélie Teissandier, Djamila Onesime, Valentin Loux, et al. Overview of a surface-ripened cheese community functioning by meta-omics analyses. PLoS ONE, 2015, 10 (4), pp.e0124360. 10.1371/journal.pone.0124360. hal-01195505

HAL Id: hal-01195505

https://hal.science/hal-01195505v1

Submitted on 24 Apr 2020

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# Overview of a Surface-Ripened Cheese Community Functioning by Meta-Omics Analyses

Eric Dugat-Bony¹,²*, Cécile Straub¹,²*, Aurélie Teissandier³,⁴, Djamila Onésime⁵,⁶, Valentin Loux⁷, Christophe Monnet¹,², Françoise Irlinger¹,², Sophie Landaud²,¹, Marie-Noëlle Leclercq-Perlat¹,², Pascal Bento⁷, Sébastien Fraud⁸, Jean-François Gibrat⁷, Julie Aubert³,⁴, Frédéric Fer¹,²,³,⁴, Eric Guédon⁵,⁶, Nicolas Pons⁹, Sean Kennedy⁹, Jean-Marie Beckerich¹,², Dominique Swennen¹,², Pascal Bonnarme¹,²*

1 INRA, UMR 782 Génie et Microbiologie des Procédés Alimentaires, F-78850, Thiverval-Grignon, France,

2 AgroParisTech, UMR 782 Génie et microbiologie des procédés alimentaires, F-78850, Thiverval-Grignon,

France, 3 AgroParisTech, UMR 518 Mathématiques et Informatiques Appliquées, F-75231, Paris, France,

4 INRA, UMR 518 Mathématiques et Informatiques Appliquées, F-75231, Paris, France, 5 INRA, Institut

Micalis, F-78352, Jouy-en-Josas, France, 6 AgroParisTech, Institut Micalis, F-78352, Jouy-en-Josas,

France, 7 INRA, UR1404 Mathématiques et Informatique Appliquées du Génome à l'Environnement,

F-78352, Jouy-en-Josas, France, 8 Actalia, F-74801, La Roche sur Foron, France, 9 INRA, US 1367

Metagenopolis, F-78352, Jouy-en-Josas, France

* These authors contributed equally to this work.

* Pascal.Bonnarme@grignon.inra.fr

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

OPEN ACCESS

**Citation:** Dugat-Bony E, Straub C, Teissandier A, Onésime D, Loux V, Monnet C, et al. (2015) Overview of a Surface-Ripened Cheese Community Functioning by Meta-Omics Analyses. PLoS ONE 10(4): e0124360. doi:10.1371/journal.pone.0124360

**Academic Editor:** Danilo Ercolini, University of Naples Federico II, ITALY

**Received:** December 20, 2014

**Accepted:** March 2, 2015

**Published:** April 13, 2015

**Copyright:** © 2015 Dugat-Bony et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was supported by the ExEco program (a joint metatranscriptomic and biochemical approach to the cheese ecosystem: for an improved monitoring of the expression of a complex food ecosystem) (ANR-09-ALIA-012-01), funded by the French National Research Agency (ANR). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Abstract

Cheese ripening is a complex biochemical process driven by microbial communities composed of both eukaryotes and prokaryotes. Surface-ripened cheeses are widely consumed all over the world and are appreciated for their characteristic flavor. Microbial community composition has been studied for a long time on surface-ripened cheeses, but only limited knowledge has been acquired about its *in situ* metabolic activities. We applied metagenomic, metatranscriptomic and biochemical analyses to an experimental surface-ripened cheese composed of nine microbial species during four weeks of ripening. By combining all of the data, we were able to obtain an overview of the cheese maturation process and to better understand the metabolic activities of the different community members and their possible interactions. Furthermore, differential expression analysis was used to select a set of biomarker genes, providing a valuable tool that can be used to monitor the cheese-making process.

## Introduction

Microbial communities are of major importance in the fermentation of food products. Fermentation remains a widespread means for food processing and preservation, and fermented foods (including cheese) are widely consumed worldwide. The composition and behavior of the microbial communities in a cheese are important for its characteristic organoleptic properties, shelf life and safety [1]. These communities are involved in the generation of a wide
**Competing Interests:** The authors have declared that no competing interests exist.

range of diverse beneficial functions as a result of individual metabolism and/or complex eco-
logical interactions [2].

To date, a plethora of work related to functions of technological interest have been pub-
lished on cheese-inhabiting microorganism. Such studies mainly concern the ability of those
microorganisms to generate functions such as proteolysis [3], lipolysis [4] and/or catabolic
routes leading to aroma compound production [5-7]. However, although some microorgan-
isms that inhabit cheese are known to be key drivers of the ripening process, our understanding
of how individual microbes and microbial groups change over time within the cheese matrix
and contribute to the structure and function of specific communities remains incomplete.

With the recent advances in high-throughput sequencing technologies (HTS), sensitive profiling of microbial communities from fermented food products can now be performed on an unprecedented scale via the massive sequencing of short DNA fragments [8,9]. Metagenomic studies, including both meta-barcoding (e.g., the deep-sequencing of variable regions of the prokaryotic SSU rRNA gene or of the fungal ITS) and whole metagenome sequencing projects, have made it possible to characterize the microbial community composition of many cheese varieties and to access the diversity of sub-dominant populations [10–13]. Furthermore, genome sequencing of several representative strains isolated from cheese or used as starter culture in the cheese-making process has allowed us to access their metabolic arsenal [14–16]. The next step towards a better understanding of how the cheese ecosystem functions would be to evaluate the expression of these genes *in situ*. As recently demonstrated for a Camembert-type cheese, this is feasible through metatranscriptomic analyses using RNA sequencing [17]. In this example, the authors followed the metabolic activity of a relatively simple community essentially composed of a yeast, *Geotrichum candidum*, and a fungus, *Penicillium camemberti*, and highlighted key functions and metabolic pathways that are expressed during the ripening process. Because the genome of these organisms were not available, *de novo* assembly of the RNA-Seq data (long reads) was required prior to functional assignment of the resulting contigs. This approach could be applied to more complex cheese microbial communities comprising both fungal and bacterial species. In this case, short reads analysis offering a higher sequencing depth is preferable, but it would be highly desirable to have the annotated reference genomes of all of the species.

In the present work, we combined microbiological, biochemical, metagenomic (DNA-Seq) and metatranscriptomic (RNA-Seq) data collected from a simplified microbial community capable of reproducing the complex metabolic pattern of cheese maturation [18,19]. To facilitate these analyses, we established a reference database of all the genomes of the studied community, onto which sequence reads could be mapped. The main objective of the study was to obtain a global view of the dynamics of the microbial community structure as well as the expression profiles of its metabolic potential throughout a ripening cycle at different scales—whole microbial community down to the gene level. Moreover, differential analysis of the ecosystem’s metatranscriptome was performed which should enable us to propose a set of biomarker genes that are representative of the most active species at various stages of ripening. Thus, we expect to reveal the sequential development and/or metabolic features of microbial species, and possibly highlight metabolic complementarities and possible interaction phenomena that sustain the expression of important functionalities of technological interest.

**Materials and Methods**

Cheese production

Full details on microorganisms used for cheese ripening and cheese production are given in the
[S1 File]. Briefly, cheese production was performed with 120 L of pasteurized milk, under aseptic
conditions in a sterilized chamber [20]. A lactic starter culture containing *Lactococcus lactis* subsp. *lactis* S3+ and S3- inoculated at concentrations equivalent to 2 x 10⁶ and 4 x 10⁶ CFU/mL, respectively, was used in combination with a mix of *Kluyveromyces lactis* 3550 (10⁴ CFU/mL), *Debaryomyces hansenii* 304 (10⁴ CFU/mL) and *Geotrichum candidum* ATCC 204307 (10³ CFU/mL). Next, 120 mL of a filter-sterilized CaCl₂ solution (10%) and 40 mL of rennet (20 mg/L of chymosin (Chr. Hansen, Arpajon, France)) were added to allow the milk to coagulate. The curd was cut into small cheeses (diameter: 5 cm; height: 1.5 cm; weight: 26 g) and immersed in sterile brine (270 g/L NaCl, pH 5.5 measured using a contact electrode) to obtain a salt concentration of 1.7%. The five ripening bacteria (*Corynebacterium casei* UCMA 3821, *Brevibacterium aurantiacum* ATCC 9174, *Arthrobacter arilaitensis* CIP 108037, *Staphylococcus equorum* Mu2 and *Hafnia alvei* GB001) were inoculated onto the surface of the cheese at a rate of 2 x 10⁵ CFU/g. The inoculated cheeses were ripened for four weeks at 14°C and 97% relative humidity in sterile crystallizing basins.

## Cheese sampling and microbial analyses

Samples were collected on days 1, 7, 14, 21 and 31. Day 1 corresponds to the cheese curd (before immersion in brine). Three crystallizing basins (corresponding to three replicates) were analyzed at each time-point. As all cheeses were produced from the same batch of milk and inoculum, these replicates should not been considered as true biological replicates. The four small cheeses in each crystallizing basin were crushed and homogenized with sterile forks and knives. Serial dilutions were performed in 9 g/L NaCl from one gram of cheese and plated in triplicate on agar plates. Three selective culture media were used: brain heart infusion (Biokar Diagnostics) with amphotericin (50 mg/L) for cheese-surface bacteria, de Man-Rogosa-Sharpe (pH 6.5, Biokar Diagnostics) with amphotericin (50 mg/L) for lactic acid bacteria, and yeast extract-glucose-chloramphenicol (Biokar Diagnostics) with 2,3,5-triphenyltetrazolium chloride (10 mg/L) for yeasts. The strains could be selectively counted on these media because each had a distinct morphotype.

## Biochemical analyses

**Lactose and lactate concentration.** The lactose and lactate contents of cheeses were determined by high-performance liquid chromatography (HPLC), as previously described [20]. Briefly, a cheese suspension was prepared with finely ground cheese (7.5 g and 10 g for lactose and lactate, respectively) and 10 mL of distilled water. After incubation for 1 h at 50°C, it was homogenized for 2 min at 25000 rpm using a mechanical blender (Ultra-Turrax model T25, Ika Laortechnik).

For lactose measurement, 12.5 mL of Carrez I solution (150 g K₄(Fe(CN)₆), 3H₂O/L water), 12.5 mL of Carrez II solution (240 g Zn(CH₃COO)₂, 2 H₂O/L water) and 2.5 mL of 1 M NaOH were added to the suspension. After homogenization and incubation for 1 h at 25°C, the suspension was filtered through Whatman paper n°42 (GE Healthcare) and through a cellulose filter of 0.22 µm porosity (Minisart SP25, Sartorius). Finally, the filtrate was analyzed by HPLC.

For lactate measurement, 10 mL of 240 g/L trichloroacetic acid (TCA) and 10 mL of water were added to the suspension after cooling at 25°C. The mixture was incubated for 1h at 25°C and filtered through Whatman paper n°42 (GE Healthcare). Finally, the filtrate was analyzed by HPLC.

**Free amino acids concentration.** A cheese suspension was prepared from 5 g of cheese and 45 mL of distilled water. It was homogenized for 2 min at 25000 rpm using a mechanical blender (Ultra-Turrax model T25, Ika Laortechnik) and incubated for 30 min at 50°C. After centrifugation for 30 min at 5000 x g and 4°C, the supernatant was filtered through Whatman
paper n°42 (GE Healthcare). The free amino acids concentration was measured on the filtrate according to the ninhydrin method [21] and quantified using leucine as a standard.

**Proteolysis index.** A cheese suspension was prepared from 10 g of cheese and 90 mL of distilled water. It was homogenized for 2 min at 25000 rpm using a mechanical blender (Ultra-Turrax model T25, Ika Laortechnik). This suspension was used for both total nitrogen and non-casein nitrogen content measure.

Total nitrogen content was measured directly on the cheese suspension using the Kjeldahl method according to standard NF EN ISO 8968-1.

The non-casein nitrogen content (soluble at pH 4.6) was measured according to the procedure described in standard NF ISO 27871. The cheese suspension was adjusted to pH 4.6 using 5N HCl and filtered through Whatman paper n°42 (GE Healthcare). The nitrogen content of this soluble fraction was measured using the Kjeldahl method according to standard NF EN ISO 8968-1.

The proteolysis index was calculated by dividing the soluble nitrogen content by the total nitrogen content and multiplying by 100, as previously described [22].

**Lipolysis index.** The lipolysis index corresponds to the quantity of KOH (in mg) required to neutralize the free fatty acids contained in 1 g of cheese. It was determined by adapting the trimetric measurement described by Mouillet et al. [23]. Briefly, 2.5 g of cheese was homogenized for 2 min at 25000 rpm using a mechanical blender (Ultra-Turrax model T25, Ika Laortechnik) in 22.5 mL of ethanol/diethyl ether (1:1 v/v) in order to extract fat content. The suspension was filtered through Watman paper No. 42 (GE Healthcare) and the acidity of the filtrate was measured by standard titration with 1M KOH in ethanol using phenolphthalein as a color indicator.

## Metagenomic and metatranscriptomic analyses

**Extraction of DNA from cheese samples.** Total genomic DNA was obtained after casein solubilization and cell recovery from two grams of cheese by adapting the method described by Baruzzi et al. [24]. Cheese samples were mixed with 18 mL of sodium citrate solution (20 g/L trisodium citrate dihydrate) and the mixture was dispersed for 2 min at 24,000 rpm with a mechanical blender (Ultra-Turrax model T25, Ika Laortechnik). A second dispersion was performed after 10 min of incubation at room temperature. The mixture was then centrifuged at 6,400 x g for 10 min at 4°C, and the supernatant was removed. The casein pellet (containing the microbial cells) was resuspended in 5 mL of a Triton X-100 aqueous solution (2.5% v/v), vigorously shaken, heated in a water bath at 70°C for 10 min, centrifuged at 6,400 x g for 10 min at 4°C, and rinsed twice in physiological saline solution. The pellet was dissolved in a mixture of 270 µL of guanidium thiocyanate (4 M) in Tris-HCl (pH 7.5, 0.1 M) and 30 µL of sodium lauroyl sarcosinate (100 g/L), and transferred to a 2 mL tube containing 250 mg of 0.1 mm-diameter zirconium beads and 250 mg of 0.5 mm-diameter zirconium beads (Sigma, St-Quentin-Fallavier, France). Proteinase K treatment, bead-beating, phenol-chloroform extraction, RNase treatment and ethanol precipitation were then performed as previously described by Leclercq-Perlat et al. [25]. The DNA pellet was dissolved in 120 µL of Tris EDTA buffer (10 mM Tris-HCl, pH 8.0, 1 mM EDTA).

**Extraction of RNA from cheese samples.** Total RNA was extracted from 500 mg cheese samples without the prior separation of microbial cells, as previously described [26]. Three separate extractions were performed for each cheese sample and pooled before library construction. The quality of total RNA was analyzed with a 2100 Bioanalyzer and RNA 6000 NANO chips (Agilent, Palo Alto, CA, USA).

**DNA and cDNA library construction and SOLiD sequencing.** Library construction and SOLiD sequencing were performed at the INRA MetaQuant facility (Jouy-en-Josas, France).
The DNA libraries were constructed from 3 µg of total DNA by using the SOLiD Fragment Library Construction Kit (Applied Biosystems, Bedford, MA, USA) and then barcoded with the SOLiD Fragment Library Barcoding Kit. DNA samples collected at day 1 did not allow to construct libraries that meet the quantity and quality criteria enabling to perform sequencing and, thus, were excluded from this study. The cDNA libraries were constructed from 200 to 500 ng total RNA using a SOLiD Whole Transcriptome Analysis Kit and were barcoded with the SOLiD Transcriptome Multiplexing Kit. The SOLiD ePCR kit and SOLiD Bead Enrichment Kit were used to process DNA and cDNA samples for sequencing, and the SOLiD 4 System was used for sequencing.

**Reads processing and mapping.** Short sequence reads of 35 bp were mapped (a maximum of three mismatches were allowed) onto a reference composed of the nine sequenced strains of microorganisms (Table 1) using Bowtie software (version 0.12.7) [27].

After the smart filtering of multiple reads, the numbers of reads mapped onto the reference were counted with METEOR software [28]. Only uniquely mapped reads were further analyzed. General library statistics are available in S1 Table.

**DNA-Seq data analyses.** For each microorganism, sequencing coverage was estimated with the Lander/Waterman equation [29]: $C = LN / G$, where C is the coverage, L is the read length, N is the number of reads and G is the haploid genome length.

**RNA-Seq data analyses.** Sequence reads mapping to the coding DNA sequence (CDS) features were retrieved from the raw dataset.

First, functional classification of the metatranscriptomic dataset was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations [30]. To do this, data were filtered to remove genes displaying an average of less than five reads per sample across the entire dataset (15 samples) and normalized according to the library size using custom scripts built under the statistical environment R (http://www.r-project.org/).

Second, differential expression analysis over time was done by comparing the number of mapped reads for each gene at two different time-points (n = 3 per ripening time). For this analysis, data (six samples for each comparison) were filtered: (i) genes with less than five mapped reads were eliminated; and (ii) genes whose replicates were heterogeneous were discarded (variation coefficient > minimum mean observed for the two groups of three replicates). Data normalization and determination of differentially expressed genes were then conducted using the Bioconductor DESeq2 package in the statistical environment R [31,32]. Raw p-values were adjusted for multiple testing using the Benjamini-Hochberg procedure [33], which assesses the False Discovery Rate. Gene transcripts with an adjusted p-value < 0.05 were considered to be differentially abundant between two ripening times.

## Sequence accession numbers

The raw SOLiD read data for all samples was deposited in the European Bioinformatics Institute's European Nucleotide Archive under the accession number PRJEB6315. Accession numbers of reference genomes are given in Table 1.

## Results and Discussion

### Overview of the cheese microbial community composition and activity during ripening

In order to estimate the overall contribution of each microorganism throughout the cheese-ripening kinetics, we used three complementary approaches: microbiological counting providing a measurement of the viable cells, metagenomic sequencing (DNA-Seq) indicating the
<table><caption>Table 1. Reference genomes used for mapping of the sequence reads.</caption><thead><tr><th>Species</th><th>Nb of genes (CDS)</th><th>Genome size (Mb)</th><th>Status</th><th>BioProject Accession number</th><th>Reference</th></tr></thead><tbody><tr><td colspan="6"><strong>Six bacteria</strong></td></tr><tr><td>Arthrobacter arilaitensis CIP 108037</td><td>3,689</td><td>3.92</td><td>Complete</td><td>PRJEA50353</td><td>[40]</td></tr><tr><td>Brevibacterium aurantiacum ATCC9174</td><td>4,105</td><td>4.39</td><td>Draft (76 contigs)</td><td>PRJNA405</td><td>no reference</td></tr><tr><td>Corynebacterium casei UCMA3821</td><td>3,015</td><td>3.12</td><td>Draft (106 contigs)</td><td>PRJEA76363</td><td>[41]</td></tr><tr><td>Hafnia alvei GB001</td><td>4,692</td><td>4.89</td><td>Draft (137 contigs)</td><td>PRJEB6257</td><td>to be published</td></tr><tr><td>Lactococcus lactis S3</td><td>2,615</td><td>2.49</td><td>Draft (163 contigs)</td><td>PRJEB6259</td><td>to be published</td></tr><tr><td>Staphylococcus equorum Mu2</td><td>2,932</td><td>2.96</td><td>Draft (30 contigs)</td><td>PRJEA88899</td><td>[76]</td></tr><tr><td colspan="6"><strong>Three yeasts</strong></td></tr><tr><td>Debaryomyces hansenii CBS767</td><td>6,421</td><td>12.2</td><td>Complete</td><td>PRJNA13832</td><td>[14]</td></tr><tr><td>Geotrichum candidum ATCC204307</td><td>6,958</td><td>24.9</td><td>Complete</td><td>PRJEB5752</td><td>to be published</td></tr><tr><td>Kluyveromyces lactis NRRL Y1140</td><td>5,120</td><td>10.8</td><td>Complete</td><td>PRJNA13835</td><td>[14]</td></tr></tbody></table>

doi:10.1371/journal.pone.0124360.t001

proportion of DNA molecules from both viable and dead cells in each sample, and metatranscriptomics (RNA-Seq) reflecting the active populations (Fig 1).

Growth measurements revealed that the early development of yeasts, especially *D. hansenii* and *G. candidum*, contributed to the increase in the pH of the cheese curd (Fig 1A), which, in turn, prompted the development of acid-sensitive bacteria (*A. arilaitensis*, *H. alvei*, *B. aurantiacum*, *S. equorum*, *C. casei*). The lactic acid bacterium *L. lactis* was the most prominent bacterial

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

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

<figure><img src="image_7.png" /><figcaption>**Fig 1.** Changes in the microbial community structure during surface-ripened cheese maturation. (A) Microbiological counts and pH measurements. (B) Distribution of metagenomic data by species. (C) Distribution of metatranscriptomic data by species (only reads mapping CDS features were taken into account). SE: *Staphylococcus equorum*. BA: *Brevibacterium aurantiacum*. AA: *Arthrobacter arilaitensis*. HA: *Hafnia alvei*. CC: *Corynebacterium casei*. DH: *Debaryomyces hansenii*. GC: *Geotrichum candidum*. KL: *Kluyveromyces lactis*. LL: *Lactococcus lactis*. NA: data not available.</figcaption></figure>

doi:10.1371/journal.pone.0124360.g001
<table><caption>Table 2. Sequencing coverage (C) and percentage of genes (P) with at least an average of five uniquely mapped reads in the DNA-Seq dataset across the three replicates for each microbial genome during ripening.</caption><thead><tr><th rowspan="2">Species<sup>a</sup></th><th colspan="2">Day 7</th><th colspan="2">Day 14</th><th colspan="2">Day 21</th><th colspan="2">Day 31</th></tr><tr><th>C</th><th>P</th><th>C</th><th>P</th><th>C</th><th>P</th><th>C</th><th>P</th></tr></thead><tbody><tr><td>AA</td><td>0.52</td><td>81.4%</td><td>11.53</td><td>98%</td><td>8.39</td><td>97.8%</td><td>4.13</td><td>97%</td></tr><tr><td>BA</td><td>0.03</td><td>0.3%</td><td>1.02</td><td>93.1%</td><td>1.05</td><td>93.1%</td><td>0.99</td><td>92.8%</td></tr><tr><td>CC</td><td>1.35</td><td>93.6%</td><td>11.54</td><td>96.3%</td><td>110.59</td><td>96.8%</td><td>175.54</td><td>96.8%</td></tr><tr><td>DH</td><td>10.90</td><td>98.6%</td><td>9.04</td><td>98.5%</td><td>3.63</td><td>98.3%</td><td>2.78</td><td>98%</td></tr><tr><td>GC</td><td>14.39</td><td>100%</td><td>12.16</td><td>100%</td><td>4.99</td><td>100%</td><td>3.78</td><td>100%</td></tr><tr><td>HA</td><td>0.21</td><td>42.1%</td><td>12.25</td><td>99.9%</td><td>8.47</td><td>99.8%</td><td>4.40</td><td>99.5%</td></tr><tr><td>KL</td><td>6.29</td><td>97.1%</td><td>3.54</td><td>97.1%</td><td>0.86</td><td>96.3%</td><td>0.26</td><td>75.5%</td></tr><tr><td>LL</td><td>218.17</td><td>99.5%</td><td>134.74</td><td>99.5%</td><td>43.40</td><td>99.5%</td><td>21.31</td><td>99.3%</td></tr><tr><td>SE</td><td>0.04</td><td>0.2%</td><td>0.12</td><td>13.7%</td><td>0.11</td><td>12.4%</td><td>0.09</td><td>7.8%</td></tr></tbody></table>

doi:10.1371/journal.pone.0124360.t002

population until day 14, but did not grow during the ripening process. *Kluyveromyces lactis* was detected at day 1 but was not detectable thereafter. *C. casei* was the dominant bacterium at the end of the ripening process (~10⁹ CFU/g of cheese). *Geotrichum candidum* was the major yeast at day 7. At the end of the ripening process (day 31), the populations of *G. candidum* and *D. hansenii* were 7 x 10⁷ CFU and 1 x 10⁷ CFU/g of cheese, respectively. The pH of the cheese curd increased from 4.6 at the beginning of ripening to 7.5 at day 31.

Similar trends were observed between microbiological counts and DNA-Seq data (Fig 1B). The proportions of reads attributed to each microorganism revealed a co-dominance of *L. lactis* and yeasts during the first two weeks of the ripening process, which were progressively replaced by surface bacteria, principally *C. casei* (71% of the total reads at day 31). *S. equorum* exhibited the lowest genome coverage in our dataset with a maximum of 0.12 X at day 14 (Table 2). In contrast, abundant species such as *L. lactis* and *C. casei* showed the highest genome coverages with 218 X for the former (day 7) and 176 X for the latter (day 31). The most abundant yeast, *G. candidum*, displayed a maximum coverage of 14 X at day 7. The differences observed between CFU counting (Fig 1A) and the proportions of DNA-Seq reads per species (Fig 1B) might be due to several reasons. First, CFU counting is supposed to measure living cells whereas DNA sequencing target both living and dead cells. Second, CFU counting may provide an altered view of the living community structure because of viable but nonculturable (VBNC) issues, the number of individual cells generating every colonies or the difficulty to count some species (such as *G. candidum* which is a filamentous yeast). Third, DNA extraction bias (not assessed in this study), especially the cell separation step, as well as genome size may influence the proportion of reads by species observed with DNA-Seq.

RNA-Seq data indicated that *L. lactis* and *K. lactis* were the most active species at day 1 (77% and 13% of total CDS reads, respectively) (Fig 1C). The proportion of reads from *G. candidum* accounted for 84% of total CDS reads at day 7 and remained highly dominant over time (with 65% of total CDS reads at day 31). Read counts corresponding to *D. hansenii* transcripts increased from day 1 (1%) to day 14 (9%) and remained stable thereafter. Ripening bacteria, especially *C. casei* and *H. alvei*, were mainly detected at the end of ripening with RNA-Seq (9% and 7% of the CDS reads at day 31, respectively). Again, RNA extraction efficiency might be variable between microorganisms (especially between eukaryotes and prokaryotes) and maybe also between ripening times. This may have influenced the distribution of reads observed with this dataset.

<footnote><sup>a</sup>AA = *Arthrobacter arilaitensis*; BA = *Brevibacterium aurantiacum*; CC = *Corynebacterium casei*; DH = *Debaryomyces hansenii*; GC = *Geotrichum candidum*; HA = *Hafnia alvei*; KL = *Kluyveromyces lactis*; LL = *Lactococcus lactis*; SE = *Staphylococcus equorum*</footnote>
Overall, the results highlighted the successive development and metabolic activity of different microbial groups, *L. lactis* and *K. lactis* at the beginning of cheese maturation, followed by *G. candidum* and *D. hansenii* and, finally, acid-sensitive bacteria. This dynamics is in accordance with other studies conducted on surface-ripened cheese, which depicted a similar development sequence of the microbial species [19,34,35]. However, we provide new information here with respect to the metabolic activity of the different cheese-ripening populations.

## Functional expression of the cheese ecosystem over time

The global expression pattern of the cheese ecosystem throughout the ripening process was evaluated by normalizing all data by library size (S2 Table) and classifying RNA-Seq reads according to KEGG annotations (Fig 2). Data were also mapped onto the KEGG general map to obtain a dynamic view at the ecosystem level and per microbial species (S1 Fig). The KEGG categories—Amino acid metabolism, Carbohydrate metabolism, Energy metabolism, Transport and catabolism, Folding, sorting and degradation, and Translation and Signal transduction—accounted for the most abundant transcripts throughout the ripening period. These data are consistent with those observed on a Camembert-type cheese using a similar approach [17]. Differences in the expression dynamics between the functional classes were observed. However, this possibly reflected changes of several distinct metabolic pathways and required a more detailed analysis. Furthermore, it should be noticed that a great proportion of reads (between 69 and 81% depending on the sample) mapped genes without ortholog in the KEGG database. This included both annotated genes such as those encoding non enzymatic proteins or enzymes not yet referenced in KEGG pathways, as well as unknown genes.

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

**Functional classes**

**Fig 2. Functional classification of the metatranscriptome during surface-ripened cheese maturation.** Functional classes were determined according to KEGG annotations. Read counts corresponding to all species were cumulated. Read numbers were normalized (according to the library size) to 50,000 reads per sampling day.

doi:10.1371/journal.pone.0124360.g002
**Energy metabolism and iron homeostasis.** The great majority of genes involved in energy production detected in our metatranscriptomic dataset were involved in the oxidative phosphorylation pathway (PATH:ko00190) (S2 Fig). As expected, the dynamics in the detection of transcripts involved in this pathway followed the changes in the community structure observed using microbial counts (Fig 1A). *G. candidum* genes encoding ATPases (EC:3.6.3.14 and 3.6.3.6), NADH dehydrogenase (EC:1.6.5.3) and cytochrome C oxidase (EC:1.9.3.1) accounted for the most frequently detected genes (S2 Table). Their maximum detection levels were observed between days 7 and 21. At day 1, it was possible to observe transcripts of genes encoding *K. lactis* ATPase (EC:3.6.3.6) and *L. lactis* ATPase (EC:3.6.3.14), NADH dehydrogenase (EC:1.6.99.3) and cytochrome D ubiquinol oxidase (EC:1.10.3.-), which reflected the early growth of these two microbial species. At the late stage of ripening (days 21 to 31), the detection of gene transcripts involved in the oxidative phosphorylation pathway from *H. alvei*, *C. casei* and, to a lesser extent, *A. arilaitensis* increased.

Metals ions, especially iron, are cofactors and/or components of various enzymatic systems involved in key microbial metabolisms (e.g., respiration). In cheese, iron often forms complexes with various curd components such as proteins (e.g., lactoferrin, ferritin and casein) and peptides, leading to low bioavailability [36]. Eukaryotic and prokaryotic microorganisms have different specific systems to transport metals. Numerous bacteria synthesize and excrete siderophores to trap iron [37]. However, some yeasts, including *K. lactis*, *D. hansenii* and *G. candidum*, express different siderophore transport systems to take advantage of siderophores present in their environment [38,39]. Another interesting feature is that fungal genes encoding iron high-affinity transport systems are generally induced under iron starvation conditions [38]. In our metatranscriptome, gene transcripts related to iron capture and/or transport were frequently detected (S2 Table), which may indicate a need for iron mobilization within the microbial community. For instance, *feoB* transcripts encoding a ferrous iron transport protein in *L. lactis* were mainly detected at day 1, transcripts of SIT1 encoding a ferrioxamine B transporter in *G. candidum* were essentially detected from day 7 to day 21, and transcripts of several iron-siderophore ABC transporters from *C. casei* (locus tag CCAS_01175, CCAS_03565, CCAS_00130, CCAS_05220, CCAS_05225, CCAS_12545, CCAS_12550, CCAS_12555, CCAS_10415) were mainly detected at day 31. Furthermore, sequence reads mapping to genes encoding high-affinity iron transport systems from both *G. candidum* (e.g., FTR1, FET3, FTH1), *D. hansenii* (e.g., FET3) and *K. lactis* (e.g., FTH1) were also detected in our conditions. Surprisingly, we did not observe bacterial transcripts for genes possibly involved in siderophore biosynthetic pathways, although *A. arilaitensis*, *C. casei* and *B. aurantiacum* are known to carry such pathways in their genome [26,40,41]. This could be due to an insufficient sequencing depth for cheese-surface bacterial species. A second hypothesis, which requires complementary analysis to be confirmed, could be that peptides released during the ripening process through casein hydrolysis would act as iron trappers, as suggested by [36], and be used as siderophore-like complexes by members of the microbial community.

**Lactose consumption.** According to Fig 2, carbohydrate metabolism was mainly detected in the early stage of ripening. Lactose is the major carbohydrate compound present in cheese curd and is rapidly consumed during the process (Fig 3A). We identified numerous expressed genes encoding enzymes involved in lactose uptake and degradation (Fig 3B and S3A Fig). For example, *L. lactis* genes encoding components of the lactose-specific phosphoenolpyruvate-dependent transport system (PTS), including *lacEF* (EC:2.7.1.69) and *ptsI* (EC:2.7.3.9), were highly expressed at day 1. The PTS system allows the concomitant translocation and phosphorylation of carbohydrates in *L. lactis*. The detection of components from other carbohydrate-specific PTS systems such as *celB* (EC:2.7.1.69) and *manXYZ* (EC:2.7.1.69) might indicate a versatile use of these systems for lactose uptake in this bacterium, as previously suggested by Aleksandrzak-Piekerczyk
<figure><img src="image_9.png" /><figcaption>**Fig 3.** Lactose metabolism during surface-ripened cheese maturation. (A) Lactose and lactate concentrations. (B) Expression dynamics of lactose degradation pathways in *Lactococcus lactis* and *Kluyveromyces lactis*. Read numbers were normalized (according to the library size) to 50,000 reads per sampling day. For each degradation step, a histogram represents cumulative read numbers when several genes were involved.</figcaption></figure>

doi:10.1371/journal.pone.0124360.g003

[42]. Transcripts of the *lacG* gene encoding the 6P-β-Galactosidase (EC:3.2.1.85), as well as *lacACD* (EC:5.3.1.26; EC:2.7.1.144; EC:4.1.2.40) involved in Lactose-6-P degradation via the D-Tagatose-6-P pathway in *L. lactis*, were also detected, mainly during the early stages of ripening (days 1 to 7). The products from this pathway, i.e., β-D-Glucose and Glyceraldehyde-3-P, then undergo glycolysis and/or the pentose phosphate pathway. Transcripts involved in these pathways were also highly detected in *L. lactis* at the same time (S2B and S2C Fig).

*K. lactis* is also able to metabolize lactose, but this yeast uses a slightly different pathway [43]. Lactose is imported through a specific lactose permease (encoded by the gene LAC12) and is then metabolized by a β-Galactosidase (EC:3.2.1.23) encoded by LAC4 and enzymes involved in the Leloir pathway (EC:2.7.1.6; EC:2.7.7.12) and encoded by genes GAL1 and GAL7. The expression of all genes from *K. lactis* involved in this pathway was also observed in our dataset, their detection being the most important during the early stage of ripening (days 1 to 14) (Fig 3B).

*G. candidum* does not consume lactose [44], whereas *D. hansenii* can efficiently catabolize this disaccharide [22]. Although our analysis indicated the detection of transcripts corresponding to the KEGG galactose pathway (PATH:ko00052) in these species (S3A Fig), this mainly corresponded to the detection of genes encoding enzymes not necessarily specific to this pathway, such as phosphoglucomutase (EC:5.4.2.2), 6-phosphofructokinase (EC:2.7.1.11) and hexokinase (EC:2.7.1.1). Indeed, these enzymes might instead reflect the detection of glycolysis
(PATH:ko00010, S3B Fig) and the pentose phosphate pathway (PATH:ko00030, S3C Fig). Consequently, these species are not likely to significantly contribute to the lactose degradation in our model cheese.

Together, these results highlighted that in our experimental conditions, one of the key functions sustaining the cheese-ripening process, namely lactose biodegradation, involved the active participation of both *L. lactis* and *K. lactis* and revealed a functional redundancy existing for this metabolism within the studied ecosystem, as suggested from previous work [45].

**Lactate metabolism.** Lactate is of major importance in cheese making. It is produced from lactose present in milk by lactic acid bacteria (Fig 3A). Two lactate dehydrogenases transcripts of *L. lactis*, *ldhA* (EC:1.1.1.27) and *dld* (EC:1.1.2.4), were generally detected at day 1 (S2 Table), suggesting that *L. lactis* actively produces lactate at the early stage of ripening. Lactate exporter-(s), responsible for lactate extrusion from the intracellular environment in *L. lactis*, have not yet been characterized. Carvalo et al. [46] proposed llmg_2513 CDS from *L. lactis* MG1363 as a good candidate to encode a lactate transporter based on its predicted protein sequence topology that is similar to known lactate transporters from other species (such as LldP from *E. coli*, Yqkl from *B. subtilis* and JEN1 from *S. cerevisiae*), but failed to confirm this activity using the mutagenesis approach. LLACS3_00055 from *L. lactis* S3, used as starter culture in our experiment, the only CDS presenting a strong homology with llmg_2513 (89% sequence identity at the nucleic level), was not detected in our metatranscriptome and, thus, would probably not be responsible for lactate export. Several gene transcripts encoding transporters in *L. lactis* were although detected, mostly at day 1, including MFS and ABC family transporters that could be good candidates to facilitate the movement of small solutes such as lactate across cell membranes. However, none of them could be identified as a potential lactate transporter on the basis of sequence similarity with known sequences from other organisms.

Lactate degradation by the cheese microflora, principally yeasts, is a key driver for curd deacidification and occurred essentially during the two first weeks of ripening in our conditions (Fig 3A). Concomitantly, we detected reads for genes JEN1, encoding a lactate transporter, together with DLD2 (EC:1.1.2.4), encoding a lactate dehydrogenase, in both *D. hansenii* and *G. candidum* (S2 Table). It should also be mentioned that *G. candidum* DLD1 transcripts encoding a D-lactate dehydrogenase (EC:1.1.2.4) and CYB2 encoding a L-lactate dehydrogenase (EC:1.2.2.3) were also detected, generally between day 7 and day 21. In *Saccharomyces cerevisiae* [47,48] and *D. hansenii*—cultivated in a cheese-like medium [49]—the authors reported that DLD1 and CYB2 were induced by lactate, making them good candidates for lactate degradation.

**Protein and amino acid degradation.** Proteins and, more precisely, caseins, are major carbon and energy sources for microbial species living on cheese, together with lactose, lactate and lipids [6]. Proteolysis, which refers to the cleavage of caseins into small peptides and, ultimately, free amino acids by microbial proteases and peptidases, took place progressively and steadily in cheese, as attested to by the evolution of the proteolysis index and free amino acid concentrations (Fig 4A). Because a specific category for proteolysis does not exist in the KEGG database, we manually built two enzyme classes (Proteases and Peptidases) based on gene product annotations present in the reference genomes, and analyzed their expression profile during cheese ripening (Fig 4B and 4C). Normalized expression data for all these genes along with their annotation is given in S3 Table. These results suggested that *G. candidum* was the major contributor to proteolysis in this simplified surface-ripened cheese and supported other studies that indicated that proteolysis mainly occurs during the first three weeks of soft cheese ripening [17,50].

In the genome of *G. candidum*, only one predicted CDS encoding a putative protease contains a signal sequence indicating a possible extracellular localization according to SignalP [51] and the PSORT server [52]. However, we didn't detect any read mapping to this CDS in our
<figure><img src="image_10.png" /><figcaption>**Fig 4. Protein degradation during surface-ripened cheese maturation.** (A) Proteolysis and free amino acid concentration. Expression data observed for genes encoding proteases (B) and peptidases (C). Read numbers were normalized (according to the library size) to 50,000 reads per sampling day. SE: *Staphylococcus equorum*. BA: *Brevibacterium aurantiacum*. AA: *Arthrobacter arilaitensis*. HA: *Hafnia alvei*. CC: *Corynebacterium casei*. LL: *Lactococcus lactis*. KL: *Kluyveromyces lactis*. DH: *Debaryomyces hansenii*. GC: *Geotrichum candidum*.</figcaption></figure>

doi:10.1371/journal.pone.0124360.g004

metatranscriptome. Similar genome analysis revealed that there is no extracellular protease en-coding genes predicted in the *D. hansenii* and *K. lactis* available genomes, but contrary to *Y. lipolytica*, another yeast commonly found at the surface of cheese. Furthermore, extracellular protease activity assay performed using the strains used for our experimental cheese produc-tion supported those predictions (S4 Fig). Thus, we hypothesized that there were three possible non-exclusive ways for *G. candidum* to utilize caseins. First, it may directly uptake casein-de-rived peptides present in the extracellular environment, freed from the hydrolytic activities of both rennet used for milk coagulation [53] and *L. lactis* used as starter culture [54]. In our metatranscriptomic dataset, this was supported by the high number of reads observed for OPT2 and GAP1 genes of *G. candidum* (S2 Table), mainly between days 7 and 21, encoding an oligopeptide transporter and a general amino acid permease, respectively, which are known to contribute to the amino acid and protein uptake in yeasts [55,56]. Furthermore, several tran-scripts reflecting the proteolytic activity of *L. lactis* were detected earlier in the ripening kinet-ics, between days 1 to 7, namely *pepA*, *pepDB*, *pepC*, *pepN*, *pepM*, *pepT*, *pepV* and *pepX* encoding various peptidases (EC:3.4.11.7; EC:3.4.11.23; EC:3.4.22.40; EC:3.4.13.-; EC:3.4.11.2; EC:3.4.11.18; EC:3.4.11.4; EC:3.4.14.11) involved in the intracellular cleavage of small peptides into free amino acids [54] as well as genes *optA*, *optC* and *optF* encoding an oligopeptide ABC transporter system that has been shown to be responsible for casein-derived peptide transport for amino acid supply in *L. lactis* [57]. Second, *G. candidum* may internalize caseins to the vac-uole by endocytosis and degrade them using vacuolar proteases and peptidases. In agreement
with this hypothesis, expression data revealed the important detection as of day 7, of several genes involved in endocytosis in *G. candidum*, among which the most frequently detected ones were EDE1 encoding an endocytic protein [58], SAC6 encoding fimbrin [59], ALY2 encoding an alpha arrestin [60], and PIL1 and LSP1 encoding primary components of eisosomes [61]. Among the most frequently transcribed genes encoding proteases and peptidases, CPS1 (EC:3.4.17.4), PRC1 (EC:3.4.16.5) and PEP4 (EC:3.4.23.25) encode enzymes known to be active in the vacuole [62]. The third possible way that casein can be used by *G. candidum* is the libera- tion of cellular proteases and peptidases in the extracellular environment during cell lysis, e.g., the highly expressed metalloendopeptidase encoded by the PRD1 gene. Although no data in our metatranscriptomic dataset enable us to support this hypothesis, this mechanism has al- ready been suggested for a cheese isolate of *D. hansenii* [63].

The next step following initial proteolysis is amino acid biodegradation, which is generally linked to the cheese matrix alkalinization and volatile compound production known to play an important role in aroma perception [5, 6, 64–66]. The most dominant amino acids composing caseins are glutamate, proline, leucine, lysine, aspartate, valine, serine, tyrosine and isoleucine. As shown in Fig 5, among the complete list of amino acid metabolic pathways, those responsi- ble for the metabolism of these dominant amino acids are also the most detected in our meta- transcriptome. *G. candidum* accounted for the majority of the expression data observed regarding most amino acid metabolism. For example, genes involved in glutamate catabolism were highly detected in this organism. This included genes encoding the NAD-dependent glu- tamate dehydrogenase (GDH2, EC:1.4.1.2) and the NADP-dependent glutamate dehydroge- nase (GDH3, EC:1.4.1.4) responsible for the deamination of this amino acid to generate 2-oxoglutarate, which is then supplied to the TCA cycle, as well as genes encoding the gluta- mate decarboxylase (GAD1, EC:4.1.1.15), degrading glutamate into 4-aminobutyrate (GABA), which is then converted into succinate by the enzymes 4-aminobutyrate aminotransferase (EC:2.6.1.19) and NAD(P)-dependent succinate semialdehyde dehydrogenase (EC:1.2.1.79), encoded by genes UGA1 and UGA2, respectively, which also feed the TCA cycle. GDH2 induc- tion has already been observed in *G. candidum* by RT-qPCR in Reblochon-type cheese (a French surface-ripened cheese) at the end of ripening, and this gene was proposed as a bio- marker for amino acid catabolism [67]. However, *G. candidum* was not the only microorgan- ism responsible for amino acid catabolism in our experimental surface-ripened cheeses. Transcripts involved in glycine, serine and threonine metabolism (PATH:ko00260) as well as the valine, leucine and isoleuline degradation pathway (PATH:ko00280) were also highly de- tected in *L. lactis* at day 1.

**Lipid metabolism.** Lipolysis refers to the hydrolysis of triglycerides and results in the lib- eration of free fatty acids (FFAs) that are important precursors of catabolic reactions that pro- duce volatile compounds that contribute to cheese quality and flavor [4,66,68]. In our metatranscriptomic data, we observed that lipid metabolism pathway expression approximate- ly followed the same dynamics as amino acid metabolism, i.e., a strong increase from day 1 to day 7 followed by a global stability or a slow decrease during the next three weeks (Fig 2). This is in agreement with the dynamics observed by Lessard et al. [17] on a Camembert-type cheese. Furthermore, in our study, the lipolysis index steadily increased along the ripening process, which indicates a regular lipolysis of the food matrix (S5 Fig). *G. candidum* once again ac- counted for the most highly expressed genes detected for this metabolism and is likely to be the major contributor to cheese lipolysis in our conditions. However, we also detected two genes encoding putative esterases from *L. lactis* (locus tag: LLACS3_11125 and LLACS3_11440), mostly at day 1, which may contribute to the initial lipolysis. Indeed, it has been demonstrated that esterases produced by lactic acid bacteria could degrade milk fat into FFAs in hard cheese [69]. Subsequently, few lipase (EC:3.1.1.3) or esterase (EC:3.1.1.13) encoding genes from G.
Categories: amino acids and other amino acids

Alanine, aspartate and glutamate metabolism [PATH:ko00250]

Glycine, serine and threonine metabolism [PATH:ko00260]

Cysteine and methionine metabolism [PATH:ko00270]

Valine, leucine and isoleucine degradation [PATH:ko00280]

Valine, leucine and isoleucine biosynthesis [PATH:ko00290]

Lysine biosynthesis [PATH:ko00300]

Lysine degradation [PATH:ko00310]

Arginine and proline metabolism [PATH:ko00330]

Histidine metabolism [PATH:ko00340]

Tyrosine metabolism [PATH:ko00350]

Phenylalanine metabolism [PATH:ko00360]

Tryptophan metabolism [PATH:ko00380]

Phenylalanine, tyrosine and tryptophan biosynthesis [PATH:ko00400]

Beta-Alanine metabolism [PATH:ko00410]

Taurine and hypotaurine metabolism [PATH:ko00430]

Phosphonate and phosphinate metabolism [PATH:ko00440]

Selenocompound metabolism [PATH:ko00450]

Cyanoamino acid metabolism [PATH:ko00460]

D-Glutamine and D-glutamate metabolism [PATH:ko00471]

D-Arginine and D-ornithine metabolism [PATH:ko00472]

D-Alanine metabolism [PATH:ko00473]

Glutathione metabolism [PATH:ko00480]

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

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

<figure><img src="image_13.png" /><figcaption>Cysteine and methionine
metabolism [PATH:ko00270]</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Valine, leucine and isoleucine
degradation [PATH:ko00280]</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Tyrosine metabolism
[PATH:ko00350]</figcaption></figure>

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

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

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

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

**Fig 5. Gene expression related to amino acid metabolism.** For each pathway, the heatmap represents the expression dynamics over time (cumulative number of normalized reads per pathway) using a gray scale bar from 0 read in white to 500 reads in black. For seven pathways, histogram charts detail this dynamic per microbial species. CC: Corynebacterium casei, HA: Hafnia alvei, AA: Arthrobacter arilaitensis, BA: Brevibacterium aurantiacum, SE: *Staphylococcus equorum*, LL: *Lactococcus lactis*, KL: *Kluyveromyces lactis*, DH: *Debaryomyces hansenii*, GC: *Geotrichum candidum*.

doi:10.1371/journal.pone.0124360.g005

candidum (i.e., ATG15, TGL1, ROG1) were detected, and their expression level increased along with the ripening time. We also detected several transcripts involved in free fatty acid import including FAA1, FAA2 and FAA4 encoding long-chain fatty acyl-CoA synthetases (EC:6.2.1.3) and PXA1 (synonym: PAT2) and PXA2 (synonym PAT1) encoding peroxisomal fatty acyl-CoA ABC transporters. Furthermore, genes involved in the peroxisomal version of the β-oxidation pathway, namely POX1 (EC:1.3.3.6), CTA1 (EC:1.11.1.6), FOX2 (EC:4.2.1.119), ECI1 (EC:5.3.3.8), SPS19 (EC:1.3.1.34) and POT1 (EC:2.3.1.16), were highly detected from day 7 to day 31.

Identification of possible biomarkers of the cheese-ripening process by differential expression analysis

*In situ* gene expression measurement is considered as a promising tool for improving our understanding of the microflora activity in cheese and for monitoring the ripening process. Indeed, considerable efforts have been made over the last few years to develop consistent and
<table><caption>Table 3. Number of differentially expressed genes according to ripening time.</caption><thead><tr><th>Day 1</th><th>Day 7</th><th>Day 14</th><th>Day 21</th><th>Day 31</th><th></th></tr></thead><tbody><tr><td>0</td><td>302</td><td>314</td><td>257</td><td>203</td><td>Day 1</td></tr><tr><td></td><td>0</td><td>56</td><td>177</td><td>438</td><td>Day 7</td></tr><tr><td></td><td></td><td>0</td><td>79</td><td>482</td><td>Day 14</td></tr><tr><td></td><td></td><td></td><td>0</td><td>330</td><td>Day 21</td></tr><tr><td></td><td></td><td></td><td></td><td>0</td><td>Day 31</td></tr></tbody></table>

doi:10.1371/journal.pone.0124360.t003

repeatable methods for *in situ* quantification of mRNA transcripts based on RT-qPCR techniques. Several examples in the context of cheese ripening are now available in the literature [67,70–73]. However, such studies enable the transcription monitoring of a limited number of genes (generally up to one hundred) and, thus, rely on the previous selection of appropriate biomarkers relevant to the targeted functions.

In the following section, we used our metatranscriptomic dataset to select and propose a set of biomarker genes that cover diverse metabolisms relevant in the cheese-ripening process. To do this, we performed differential expression analysis using the DESeq 2 package [31] and first compared the number of differentially expressed genes between each ripening time (Table 3). Based on this result, we decided to separate the ripening kinetics into two phases, the early and late stage of ripening corresponding to day 1 to day 14 (D1 vs. D14) and day 14 to day 31 (D14 vs. D31), respectively, for which differentially expressed genes were observed. The complete lists of gene transcripts showing a differential abundance in these two comparisons, along with their actual adjusted p-values, are available in S4–S5 Tables, respectively. We then manually selected 70 genes from among this list with a high number of reads and related to a biological and/or technological function according to the genome annotations. For the most prominent microorganisms, selected genes are reported in Table 4 and classified into categories that might be relevant for the monitoring of cheese ripening.

**Brief description of differentially expressed genes between D1 and D14.** On the basis of the comparative analysis of D1 vs. D14, two genes related to glycolysis, encoding glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12) and enolase (EC:4.2.1.11), could be proposed as common biomarkers of microbial species activity at the early stage of the ripening process. Interestingly, they were more abundant at D1 in *K. lactis* and *L. lactis*, in agreement with an intense development and metabolic activity for these facultative anaerobic microorganisms during the first days of ripening. As mentioned above, both species are involved in lactose degradation within the cheese matrix. This was exactly the opposite for the cheese-surface aerobic yeasts *G. candidum* and *D. hansenii* for which both genes were more expressed at D14. Mitochondrial genes of yeasts related to the respiration chain, e.g., COB and COX, followed the same variation pattern.

Regarding *G. candidum*, genes involved in protein degradation (e.g., PRD1, PEP4) and in glutamate—the most abundant amino acid—degradation (e.g., GDH2) were more frequently detected at D14 when compared to D1 and could thus be proposed as biomarkers for proteolytic activity and amino acid catabolism. GUT1 encoding glycerol kinase (EC:2.7.1.30) involved in lipid metabolism followed the same trend and could be used as a biomarker for lipolysis.

A striking feature of *L. lactis* is the occurrence of several stress-related genes (e.g., sodA, dnaK, clpB) with high expression levels at D1 (S4 Table). This might be related to several possible stresses such as the shift from anaerobic (lactic acid production during the milk acidification and coagulation phase) to aerobic conditions during the ripening process, the osmotic pressure induced by salting or variations in the pH. Indeed, transcriptomic analysis of *L. lactis*
<table><caption>Table 4. Selection of genes showing a differential abundance between Day 1 and Day 14 and/or Day 14 and Day 31.</caption><thead><tr><th>Locus tag (gene)</th><th>Product (EC number)</th><th>Sepcies</th><th>Log2 fold change Day 14 vs. Day 1 (mean)</th><th>Log2 fold change Day 31 vs. Day 14 (mean)</th></tr></thead><tbody><tr><td colspan="5"><strong>Central Metabolism</strong></td></tr><tr><td>DEHA2G14058g (ENO1)</td><td>Enolase I (EC:4.2.1.11)</td><td>D. hansenii</td><td>3.20 (22.81)</td><td>-2.36 (21.01)</td></tr><tr><td>DEHA2G18348g (PDC1)</td><td>Pyruvate decarboxylase (EC:4.1.1.1)</td><td>D. hansenii</td><td>2.43 (46.35)</td><td>-</td></tr><tr><td>DEHA2F04796g (TDH3)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>D. hansenii</td><td>2.03 (88.25)</td><td>-</td></tr><tr><td>DEHA2E13530g (MLS1)</td><td>Malate synthase (EC:2.3.3.9)</td><td>D. hansenii</td><td>-</td><td>1.27 (69.92)</td></tr><tr><td>GECA22s00351g (TDH3)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>G. candidum</td><td>3.51 (126.71)</td><td>-</td></tr><tr><td>GECA17s02408g (ENO2)</td><td>Enolase II (EC:4.2.1.11)</td><td>G. candidum</td><td>3.72 (167.68)</td><td>-</td></tr><tr><td>GECA06s01693g (MDH1)</td><td>Malate dehydrogenase (EC:1.1.1.37)</td><td>G. candidum</td><td>4.03 (71.02)</td><td>-2.25 (69.13)</td></tr><tr><td>GECA06s04762g (PYC2)</td><td>Pyruvate carboxylase 2 (EC:6.4.1.1)</td><td>G. candidum</td><td>3.24 (200.63)</td><td>-1.97 (193.23)</td></tr><tr><td>GECA20s01396g (ADH3)</td><td>Alcohol dehydrogenase 3 (EC:1.1.1.1)</td><td>G. candidum</td><td>-</td><td>-1.66 (176.65)</td></tr><tr><td>KLLA0A09185g (ENO1)</td><td>ENO1 Enolase I (EC:4.2.1.11)</td><td>K. lactis</td><td>-6.62 (96.53)</td><td>-</td></tr><tr><td>KLLA0A11011g (PGK1)</td><td>Phosphoglycerate kinase (EC:2.7.2.3)</td><td>K. lactis</td><td>-5.58 (119.43)</td><td>-</td></tr><tr><td>KLLA0A11858g (TDH2)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>K. lactis</td><td>-4.68 (184.66)</td><td>-</td></tr><tr><td>KLLA0E07569g (FBA1)</td><td>Fructose 1 6-bisphosphate aldolase (EC:4.1.2.13)</td><td>K. lactis</td><td>-2.79 (97.16)</td><td>-</td></tr><tr><td>KLLA0E16303g (PDC1)</td><td>Pyruvate decarboxylase (EC:4.1.1.1)</td><td>K. lactis</td><td>-3.38 (163.45)</td><td>-</td></tr><tr><td>KLLA0F20988g (GAP1)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>K. lactis</td><td>-2.47 (80.16)</td><td>-3.41 (11.87)</td></tr><tr><td>KLLA0F23397g (PYK1)</td><td>Pyruvate kinase (EC:2.7.1.40)</td><td>K. lactis</td><td>-3.85 (50.88)</td><td>-</td></tr><tr><td>LLACS3_03095 (gap)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>L. lactis</td><td>-5.00 (5185.82)</td><td>-</td></tr><tr><td>LLACS3_01725 (gap)</td><td>Glyceraldehyde-3-phosphate dehydrogenase (EC:1.2.1.12)</td><td>L. lactis</td><td>-5.80 (5302.58)</td><td>-</td></tr><tr><td>LLACS3_09760 (eno)</td><td>Enolase (EC:4.2.1.11)</td><td>L. lactis</td><td>-6.13 (449.67)</td><td>-</td></tr><tr><td>LLACS3_02295 (adhE)</td><td>Aldehyde-alcohol dehydrogenase (EC:1.1.1.1; EC:1.2.1.10)</td><td>L. lactis</td><td>-3.44 (79.61)</td><td>-</td></tr><tr><td>LLACS3_04185 (gnd)</td><td>Phosphogluconate dehydrogenase (EC:1.1.1.44)</td><td>L. lactis</td><td>-4.87 (333.29)</td><td>-</td></tr><tr><td>LLACS3_06550 (deoB)</td><td>Phosphopentomutase (EC:5.4.2.7)</td><td>L. lactis</td><td>-5.76 (491.42)</td><td>-</td></tr><tr><td>CCAS_01315 (pyc)</td><td>Pyruvate carboxylase (EC:6.4.1.1)</td><td>C. casei</td><td>-</td><td>4.41 (24.97)</td></tr><tr><td>CCAS_05260</td><td>2-methylcitrate dehydratase (EC:4.2.1.79)</td><td>C. casei</td><td>-</td><td>4.96 (118.78)</td></tr><tr><td>HALV_03660 (adhE)</td><td>Aldehyde-alcohol dehydrogenase (EC:1.1.1.1; EC:1.2.1.10)</td><td>H. alvei</td><td>-</td><td>3.35 (36.21)</td></tr><tr><td>HALV_13325 (fba)</td><td>Fructose-bisphosphate aldolase class II (EC:4.1.2.13)</td><td>H. alvei</td><td>-</td><td>3.56 (73.94)</td></tr><tr><td colspan="5"><strong>Carbohydrate uptake</strong></td></tr><tr><td>LLACS3_01460 (manX)</td><td>PTS system mannose-specific EIIAB component (EC:2.7.1.69)</td><td>L. lactis</td><td>-4.54 (72.28)</td><td>-</td></tr><tr><td>LLACS3_01465 (manY)</td><td>PTS system mannose-specific EIIIC component</td><td>L. lactis</td><td>-4.51 (74.66)</td><td>-</td></tr><tr><td>LLACS3_11905 (lacE)</td><td>PTS system lactose-specific IIBC component (EC:2.7.1.69)</td><td>L. lactis</td><td>-3.43 (169.59)</td><td>-</td></tr></tbody></table>

(Continued)
<table><caption>Table 4. (Continued)</caption><thead><tr><th>Locus tag (gene)</th><th>Product (EC number)</th><th>Species</th><th>Log2 fold change Day 14 vs. Day 1 (mean)</th><th>Log2 fold change Day 31 vs. Day 14 (mean)</th></tr></thead><tbody><tr><td>LLACS3_11930 (lacA)</td><td>Galactose-6-phosphate isomerase subunit lacA (EC:5.3.1.26)</td><td>L. lactis</td><td>-3.64 (189.96)</td><td>-</td></tr><tr><td>LLACS3_07715 (ptsl)</td><td>Phosphoenolpyruvate protein phosphotransferase (EC:2.7.3.9)</td><td>L. lactis</td><td>-4.95 (2656.44)</td><td>-</td></tr><tr><td>LLACS3_07720 (ptsh)</td><td>Phosphocarrier protein HP<sup>+</sup> (EC:2.7.11.-)</td><td>L. lactis</td><td>-6.80 (2126.77)</td><td>-</td></tr><tr><td colspan="5"><strong>Protein degradation</strong></td></tr><tr><td>GECA06s04927g (PRD1)</td><td>Zinc metalloendopeptidase</td><td>G. candidum</td><td>6.59 (191.17)</td><td>-2.65 (186.53)</td></tr><tr><td>GECA01s02760g (PEP4)</td><td>Vacuolar aspartyl protease (proteinase A) (EC:3.4.23.25)</td><td>G. candidum</td><td>4.26 (118.86)</td><td>1.02 (293.28)</td></tr><tr><td>GECA10s04102g (CPS1)</td><td>Vacuolar carboxypeptidase yscS (EC:3.4.17.4)</td><td>G. candidum</td><td>-</td><td>-1.58 (62.04)</td></tr><tr><td>GECA01s04322g (OPT2)</td><td>Oligopeptide transporter</td><td>G. candidum</td><td>-</td><td>-4.46 (129.07)</td></tr><tr><td>LLACS3_02215 (ftsh)</td><td>ATP-dependent zinc metalloprotease (EC:3.4.24.-)</td><td>L. lactis</td><td>-3.55 (91.99)</td><td>-</td></tr><tr><td>CCAS_13715</td><td>Xaa-Pro dipeptidase (EC:3.4.13.9)</td><td>C. casei</td><td>-</td><td>2.96 (18.35)</td></tr><tr><td>HALV_05100 (hflB)</td><td>ATP-dependent metallopeptidase (EC:3.4.24.-)</td><td>H. alvei</td><td>-</td><td>2.97 (16.07)</td></tr><tr><td>HALV_13275 (pepP)</td><td>Proline-specific aminopeptidase (EC:3.4.11.5)</td><td>H. alvei</td><td>-</td><td>4.97 (11.34)</td></tr><tr><td colspan="5"><strong>Amino acid degradation</strong></td></tr><tr><td>GECA04s04949g (AAT2)</td><td>Cytosolic aspartate aminotransferase (EC:2.6.1.1)</td><td>G. candidum</td><td>6.10 (76.20)</td><td>2.18 (349.16)</td></tr><tr><td>GECA22s00494g (CAR2)</td><td>L-ornithine transaminase (OTAse) (EC:2.6.1.13)</td><td>G. candidum</td><td>-</td><td>-2.37 (99.20)</td></tr><tr><td>GECA03s01374g (GAP1)</td><td>General amino acid permease</td><td>G. candidum</td><td>-</td><td>-3.74 (113.88)</td></tr><tr><td>GECA04s00373g (PP2)</td><td>Serine/threonine protein phosphatase Z (EC:3.1.3.16)</td><td>G. candidum</td><td>-</td><td>3.61 (481.29)</td></tr><tr><td>GECA05s04498g (DUR1)</td><td>Urea amidolyase (EC:6.3.4.6; EC:3.5.1.54)</td><td>G. candidum</td><td>-</td><td>2.20 (116.57)</td></tr><tr><td>GECA20s00296g (LAP3)</td><td>Cysteine aminopeptidase (EC:3.4.22.40)</td><td>G. candidum</td><td>-</td><td>-4.12 (62.59)</td></tr><tr><td>GECA06s00802g (GDH2)</td><td>NAD(+)-dependent glutamate dehydrogenase (EC:1.4.1.2)</td><td>G. candidum</td><td>5.17 (350.98)</td><td>-</td></tr><tr><td>GECA09s02639g (CYS3)</td><td>Cystathionine gamma-lyase (EC:4.4.1.1)</td><td>G. candidum</td><td>3.11 (47.57)</td><td>-</td></tr><tr><td>HALV_00935 (tdcC)</td><td>Threonine/serine transporter</td><td>H. alvei</td><td>-</td><td>4.45 (52.07)</td></tr><tr><td>HALV_13395 (speA)</td><td>Arginine decarboxylase (EC:4.1.1.19)</td><td>H. alvei</td><td>-</td><td>3.35 (14.06)</td></tr><tr><td colspan="5"><strong>Sulfur metabolism</strong></td></tr><tr><td>GECA07s04267g</td><td>Putative sulfate permease</td><td>G. candidum</td><td>-</td><td>1.64 (253.70)</td></tr><tr><td>GECA18s00835g (MET6)</td><td>Cobalamin-independent methionine synthase (EC:2.1.1.14)</td><td>G. candidum</td><td>-</td><td>2.44 (255.35)</td></tr><tr><td colspan="5"><strong>Respiration and iron transport</strong></td></tr><tr><td>DEHA_mCDS7140 (COB)</td><td>Apocytochrome b</td><td>D. hansenii</td><td>2.96 (99.64)</td><td>-</td></tr><tr><td>DEHA_mCDS20585 (COX1)</td><td>Cytochrome c oxidase subunit 1 (EC:1.9.3.1)</td><td>D. hansenii</td><td>2.52 (359.26)</td><td>1.11 (818.11)</td></tr><tr><td>GECA_mCDS15388 (COB)</td><td>Apocytochrome b</td><td>G. candidum</td><td>2.39 (598.75)</td><td>-0.90 (659.71)</td></tr><tr><td>GECA_mCDS24166 (COX1)</td><td>Cytochrome c oxidase subunit 1 (EC:1.9.3.1)</td><td>G. candidum</td><td>2.77 (841.18)</td><td>-</td></tr><tr><td>GECA05s03805g (FET3)</td><td>Ferro-O2-oxidoreductase</td><td>G. candidum</td><td>7.16 (67.98)</td><td>-6.71 (58.21)</td></tr><tr><td>GECA12s01770g (SIT1)</td><td>Ferrioxamine B transporter</td><td>G. candidum</td><td>6.14 (112.96)</td><td>-4.92 (97.92)</td></tr></tbody></table>

(Continued)
<table><caption>Table 4. (Continued)</caption><thead><tr><th>Locus tag (gene)</th><th>Product (EC number)</th><th>Species</th><th>Log2 fold change Day 14 vs. Day 1 (mean)</th><th>Log2 fold change Day 31 vs. Day 14 (mean)</th></tr></thead><tbody><tr><td>GECA26s00230g (FTR1)</td><td>High affinity iron permease</td><td>G. candidum</td><td>-</td><td>-1.21 (89.06)</td></tr><tr><td>KLLA_mCDS20667 (COB)</td><td>Apocytochrome b</td><td>K. lactis</td><td>-4.74 (431.07)</td><td>-</td></tr><tr><td>KLLA_mCDS27490 (COX1)</td><td>Cytochrome c oxidase subunit 1 (EC:1.9.3.1)</td><td>K. lactis</td><td>-4.72 (813.44)</td><td>-</td></tr><tr><td>CCAS_12555</td><td>Iron-siderophore ABC transporter</td><td>C. casei</td><td>-</td><td>5.99 (53.33)</td></tr><tr><td colspan="5"><strong>Stress response</strong></td></tr><tr><td>LLACS3_06675 (dnaK)</td><td>Chaperone protein DnaK</td><td>L. lactis</td><td>-5.94 (10536.90)</td><td>1.22 (475.75)</td></tr><tr><td>LLACS3_08635 (clpB)</td><td>Chaperone protein ClpB</td><td>L. lactis</td><td>-5.89 (15396.23)</td><td>-</td></tr><tr><td>LLACS3_00650 (sodA)</td><td>Superoxide dismutase (EC:1.15.1.1)</td><td>L. lactis</td><td>-6.40 (6002.79)</td><td>-</td></tr><tr><td>CCAS_01895</td><td>Catalase (EC:1.11.1.6)</td><td>C. casei</td><td>-</td><td>4.91 (47.50)</td></tr><tr><td>GECA07s02177g (RIM101)</td><td>Transcriptional repressor involved in response to pH</td><td>G. candidum</td><td>-</td><td>1.41 (221.55)</td></tr><tr><td colspan="5"><strong>Lipid metabolism</strong></td></tr><tr><td>GECA13s03354g (GUT1)</td><td>Glycerol kinase (EC:2.7.1.30)</td><td>G. candidum</td><td>2.49 (147.08)</td><td>-2.09 (131.76)</td></tr><tr><td colspan="5"><strong>Endocytosis</strong></td></tr><tr><td>GECA06s02496g (ALY2)</td><td>Alpha arrestin, endocytosis</td><td>G. candidum</td><td>-</td><td>1.30 (82.97)</td></tr><tr><td colspan="5"><strong>Hyphal growth</strong></td></tr><tr><td>GECA07s00142g (MIT1)</td><td>Transcriptional regulator of pseudohyphal growth</td><td>G. candidum</td><td>-</td><td>2.17 (82.88)</td></tr></tbody></table>

The log2 fold change and the average number of reads (mean) are indicated only if the gene revealed a significant difference between the two ripening times (adjusted p-value < 0.05).

doi:10.1371/journal.pone.0124360.t004
has already revealed that *sodA* gene overexpression is associated with oxidative stress response during milk fermentation [74], and that genes encoding chaperones induction reflects the growth arrest of *L. lactis* in cheese observed almost 24h post-inoculation [75]. In our study, model cheeses were analyzed as a whole, thus including both surface and core parts. However, local conditions might have a great influence on physiochemical stresses encountered by *L. lactis* and thus should be taken into account for data interpretation.

**Brief description of differentially expressed genes between D14 and D31.** The comparative analysis of D14 vs. D31 revealed a shift of the metabolism of *G. candidum* towards peptide (OPT2, PRD1, LAP3, CPS1) and amino acid (GAP1, AAT2, CAR2) catabolism/transport compared to D1/D14. Another interesting feature is the detection of genes related to iron capture and/or transport, higher in D14 compared to D31 for *G. candidum* (SIT1, FET3, FTR1), but lower in *C. casei* (iron-siderophore ABC transporter). This indicated that iron capture is crucial for microbial species living on the cheese surface, as demonstrated by Monnet et al. [26]. In *C. casei*, we observed the highest abundance at D31 of CCAS\_05260 encoding a 2-methylcitrate dehydratase (EC:4.2.1.79), an important enzyme involved in the methylcitrate pathway. In *H. alvei*, gene transcripts involved in peptide/amino acid transport (*meoA*, *tdcC*) and degradation (*speA*, *pepP*, *hlfB*) were also more abundant at D31, revealing the catabolic activity of this species at the end of ripening. Finally, two transcripts involved in sulfur metabolism, encoding a putative sulfate permease (GECA07s04267g) and a methionine synthase (MET6), were more highly detected at D31 in *G. candidum*. They could be involved in sulfur metabolism through sulfur recycling from various sulfur compounds, including volatile ones, which are released in great quantity at the surface of ripened cheeses [7].

**Conclusion**

In this study, we used a combination of microbiological, biochemical, metagenomic and metatranscriptomic methods to obtain a detailed picture of an experimental surface-ripened cheese ecosystem that functions during the ripening process. Overall, we were able to reveal the major contribution of the most dominant microbial species (e.g. *L. lactis*, *K. lactis*, *G. candidum*, *D. hansenii* and *C. casei*) and possible interactions regarding key functions involved in the dairy matrix degradation. *L. lactis* and *K. lactis* activities during the early stage of ripening enabled the rapid consumption of lactose. Lactate, produced from lactose by *L. lactis*, was then rapidly consumed by *D. hansenii* and *G. candidum* for which we detected high levels of lactate dehydrogenase transcripts. Regarding protein and lipid metabolism, the great majority of RNA-Seq reads mapped *G. candidum* genes, which suggested a strong influence of this species on casein and fat degradation. At the end of ripening, our dataset indicated the expression of amino acid degradation-related genes by *G. candidum* and acid-sensitive bacteria such as *C. casei* and *H. alvei*, which were linked to their late development at the cheese surface. We demonstrated that global gene expression data collected at the ecosystem scale were in good accordance with the observed phenomena (e.g., biochemical and microbiological data) and provided the unique opportunity to simultaneously address questions related to different metabolisms and involving several individual species.

Furthermore, statistical methods based on differential expression analysis made it possible to select a short list of potential biomarkers. This valuable tool might be particularly useful for more precise and in-depth studies aiming at understanding and/or simply following the contribution of different strains or species in the ripening process, sustaining the production of different surface-ripened cheeses.

Altogether, metatranscriptomic analysis revealed the proportion of the genes that are actually expressed within a food microbial community composed of both eukaryotes and prokaryotes. When combined with biochemical data, it may also indicate the microbial populations
that are metabolically active and how they respond to a perturbation. Thus, it may be applied
to detect near-instantaneous responses to environmental perturbations (e.g. biotic and abiotic
constraints) which could occur during the ripening process and other related food fermenta-
tion processes. Currently, the main limitation of this approach remains the difficulty to detect
low-abundant species which could also contribute to the matrix transformation. The combina-
tion of both using efficient rRNA depletion methods and increasing the sequencing depth
might enable to partly overcome this limitation.

**Supporting Information**

**S1 Fig. Metabolic pathways detected in the metatranscriptomic dataset.** Genes exhibiting an average of > 5 normalized reads were mapped in black onto KEGG general metabolic pathways (ko01100). E: complete ecosystem. LL: *Lactococcus lactis*. KL: *Kluyveromyces lactis*. GC: *Geotrichum candidum*. DH: *Debaryomyces hansenii*. CC: *Corynebacterium casei*. HA: *Hafnia alvei*. (TIF)

**S2 Fig. Expression profile of the oxidative phosphorylation pathway during surface-ripened cheese maturation.** Histogram charts represent the cumulative number of normalized reads per sampling day and per microbial species. CC: *Corynebacterium casei*, HA: *Hafnia alvei*, AA: *Arthrobacter arilaitensis*, BA: *Brevibacterium aurantiacum*, SE: *Staphylococcus equorum*, LL: *Lactococcus lactis*, KL: *Kluyveromyces lactis*, DH: *Debaryomyces hansenii*, GC: *Geotrichum candidum*.
(TIF)

**S3 Fig. Expression profile of carbohydrate metabolic pathways.** Histogram charts represent the expression dynamics (cumulative number of normalized reads per sampling day and per microbial species) of the galactose metabolism (A), glycolysis-gluconeogenesis pathway (B) and pentose phosphate pathway (C). CC: *Corynebacterium casei*, HA: *Hafnia alvei*, AA: *Arthrobacter arilaitensis*, BA: *Brevibacterium aurantiacum*, SE: *Staphylococcus equorum*, LL: *Lactococcus lactis*, KL: *Kluyveromyces lactis*, DH: *Debaryomyces hansenii*, GC: *Geotrichum candidum*.
(TIF)

**S4 Fig. Extracellular protease assay.** *Yarrowia lipoltica* 1E07 (Yl), *Kluyveromyces lactis* 3550 (Kl), *Debaryomyces hansenii* 304 (Dh) and *Geotrichum candidum* ATCC 204307 (Gc) were spotted on protease assay medium (0.67% yeast nitrogen base without ammonium sulfate and amino acids (Difco Laboratories), 0.1% glucose, 50 mM phosphate buffer, pH 6.8, 2% skim milk (Difco Laboratories)), and incubated 6 days at 18°C. Clarification zone around the colony indicated extracellular protease activity.
(TIF)

**S5 Fig. Lipolysis index measurement during surface-ripened cheese maturation.**
(TIF)

**S1 Table. Metagenome and metatranscriptome statistics.**
(XLSX)

**S2 Table. Expression data for all CDSs detected in the metatranscriptome (read numbers were normalized to 50,000 reads per sampling day).**
(XLSX)

**S3 Table. Expression data for CDSs encoding peptidases and proteases (read numbers were normalized to 50,000 reads per sampling day).**
(XLSX)
**S4 Table.** Gene transcripts showing a differential abundance between Day 1 and Day 14.
(XLSX)

**S5 Table.** Gene transcripts showing a differential abundance between Day 14 and Day 31.
(XLSX)

**S1 File.** Supplementary methods: experimental surface ripened-cheese production.
(DOC)

**Acknowledgments**

We thank Serge Casaregola and Guillaume Morel (INRA, UMR1319 Micalis, Thiverval-Grignon, France) for providing the genome of *G. candidum* ATCC 204307, and Anne Sophie Sarthou and Jessie Castellote (INRA, UMR782 GMPA, Thiverval-Grignon, France) for technical assistance. We are grateful to our undergraduate bioinformatics students Maxime Ferrier, Roxane Boyer and Jeffrey Boggio (IUT Génie Biologique option bioinfomatique, Aurillac, France) for the development of scripts for functional classification, and Nicolas Parisot (Université d'Auvergne, Clermont-Ferrand, France) for supervising their work.

**Author Contributions**

Conceived and designed the experiments: P. Bonnarme JMB CM FI DS. Performed the experiments: CS CM FI SL SK DO MNLP SF. Analyzed the data: EDB CS P. Bento VL NP AT JA FF EG. Contributed reagents/materials/analysis tools: P. Bonnarme JMB JFG. Wrote the paper: EDB CS P. Bonnarme JMB DS.

**References**

1. Montel MC, Buchin S, Mallet A, Delbes-Paus C, Vuitton DA, Desmesures N, et al. Traditional cheeses: Rich and diverse microbiota with associated benefits. Int J Food Microbiol. 2014; 177: 136–154. doi: 10.1016/j.ijfoodmicro.2014.02.019 PMID: 24642348

2. Irlinger F, Mounier J. Microbial interactions in cheese: implications for cheese quality and safety. Curr Opin Biotechnol. 2009; 20: 142–148. doi: 10.1016/j.copbio.2009.02.016 PMID: 19342218

3. Sousa MJ, Ardö Y, McSweeney PLH. Advances in the study of proteolysis during cheese ripening. Int Dairy J. 2001; 11: 327–345.

4. Collins YF, McSweeney PLH, Wilkinson MG. Lipolysis and free fatty acid catabolism in cheese: a review of current knowledge. Int Dairy J. 2003; 13: 841–866. doi: 10.1016/S0958-6946(03)00109-2

5. Yvon M, Rijnen L. Cheese flavour formation by amino acid catabolism. Int Dairy J. 2001; 11: 185–201.

6. Marilley L. Flavours of cheese products: metabolic pathways, analytical tools and identification of producing strains. Int J Food Microbiol. 2004; 90: 139–159. doi: 10.1016/S0168-1605(03)00304-0 PMID: 14698096

7. Landaud S, Helinck S, Bonnarme P. Formation of volatile sulfur compounds and metabolism of methionine and other sulfur compounds in fermented food. Appl Microbiol Biotechnol. 2008; 77: 1191–1205. doi: 10.1007/s00253-007-1288-y PMID: 18064452

8. Van Hijum SA, Vaughan EE, Vogel RF. Application of state-of-art sequencing technologies to indigeneous food fermentations. Curr Opin Biotechnol. 2013; 24: 178–186. doi: 10.1016/j.copbio.2012.08.004 PMID: 22960050

9. Ercolini D. High-Throughput Sequencing and Metagenomics: Moving Forward in the Culture-Independent Analysis of Food Microbial Ecology. Appl Environ Microbiol. 2013; 79: 3148–3155. doi: 10.1128/AEM.00256-13 PMID: 23475615

10. Quigley L, O'Sullivan O, Beresford TP, Ross RP, Fitzgerald GF, Cotter PD. High-Throughput Sequencing for Detection of Subpopulations of Bacteria Not Previously Associated with Artisanal Cheeses. Appl Environ Microbiol. 2012; 78: 5717–5723. doi: 10.1128/AEM.00918-12 PMID: 22685131

11. Bokulich NA, Mills DA. Facility-Specific “House” Microbiome Drives Microbial Landscapes of Artisan Cheesemaking Plants. Appl Environ Microbiol. 2013; 79: 5214–5223. doi: 10.1128/AEM.00934-13 PMID: 23793641
12. Wolfe BE, Button JE, Santarelli M, Dutton RJ. Cheese Rind Communities Provide Tractable Systems for In Situ and In Vitro Studies of Microbial Diversity. Cell. 2014; 158: 422–433. doi: 10.1016/j.cell.2014.05.041 PMID: 25036636

13. Delcenserie V, Taminiau B, Delhalle L, Nezer C, Doyen P, Crevecoeur S, et al. Microbiota characterization of a Belgian protected designation of origin cheese, Herve cheese, using metagenomic analysis. J Dairy Sci. 2014; 97: 6046–6056. doi: 10.3168/jds.2014-8225 PMID: 25064656

14. Dujon B, Sherman D, Fischer G, Durrens P, Casaregola S, Lafontaine I, et al. Genome evolution in yeasts. Nature. 2004; 430: 35–44. PMID: 15229592

15. Cheeseman K, Ropars J, Renault P, Dupont J, Gouzy J, Branca A, et al. Multiple recent horizontal transfers of a large genomic region in cheese making fungi. Nat Commun. 2014; 5. doi: 10.1038/ncomms3876

16. Almeida M, Hébert A, Abraham AL, Rasmussen S, Monnet C, Pons N, et al. Construction of a dairy mi-crobial genome catalog opens new perspectives for the metagenomic analysis of dairy fermented prod-ucts. BMC Genomics. 2014; 15: 1101. doi: 10.1186/1471-2164-15-1101 PMID: 25496341

17. Lessard MH, Viel C, Boyle B, St-Gelais D, Labrie S. Metatranscriptome analysis of fungal strains Penicillium camemberti and Geotrichum candidum reveal cheese matrix breakdown and potential development of sensory properties of ripened Camembert-type cheese. BMC Genomics. 2014; 15: 235. doi: 10.1186/1471-2164-15-235 PMID: 24670012

18. Bonaïti C, Irlinger F, Spinnler HE, Engel E. An Iterative Sensory Procedure to Select Odor-Active Associations in Complex Consortia of Microorganisms: Application to the Construction of a Cheese Model. J Dairy Sci. 2005; 88: 1671–1684. doi: 10.3168/jds.S0022-0302(05)72839-3 PMID: 15829658

19. Mounier J, Monnet C, Vallaeys T, Ardit i R, Sarthou AS, Helias A, et al. Microbial Interactions within a Cheese Microbial Community. Appl Environ Microbiol. 2008; 74: 172–181. doi: 10.1128/AEM.01338-07 PMID: 17981942

20. Leclercq-Perlat MN, Oumer A, Bergere JL, Spinnler HE, Corrieu G. Growth of Debaryomyces hansenii on a bacterial surface-ripened soft cheese. J Dairy Res. 1999; 66: 271–281. null

21. Doi E, Shibata D, Matoba T. Modified colorimetric ninhydrin methods for peptidase assay. Anal Bio-chem. 1981; 118: 173–184. doi: 10.1016/0003-2697(81)90175-5 PMID: 7039409

22. Leclercq-Perlat MN, Oumer A, Buono F, Bergere JL, Spinnler HE, Corrieu G. Behavior of Brevibacterium linens and Debaryomyces hansenii as ripening flora in controlled production of soft smear cheese from reconstituted milk: protein degradation. J Dairy Sci. 2000; 83: 1674–1683. doi: 10.3168/jds.S0022-0302(00)75036-3 PMID: 10984142

23. Mouillet L, Luquet FM, Nicod H, Boudier JF, Mahieu H. La lipolyse des laits. Etude d'une méthode rapide de mesure. Le Lait. 1981; 61: 171–186.

24. Baruzzi F, Matarante A, Caputo L, Morea M. Development of a Culture-Independent Polymerase Chain Reaction-Based Assay for the Detection of Lactobacilli in Stretched Cheese. J Rapid Methods Autom Microbiol. 2005; 13: 177–192. doi: 10.1111/j.1745-4581.2005.00018.x

25. Leclercq-Perlat MN, Picque D, Martin del Campo Barba ST, Monnet C. Dynamics of Penicillium camemberti growth quantified by real-time PCR on Camembert-type cheeses under different conditions of temperature and relative humidity. J Dairy Sci. 2013; 96: 4031–4040. doi: 10.3168/jds.2012-6372 PMID: 23587382

26. Monnet C, Back A, Irlinger F. Growth of Aerobic Ripening Bacteria at the Cheese Surface Is Limited by the Availability of Iron. Appl Environ Microbiol. 2012; 78: 3185–3192. doi: 10.1128/AEM.00085-12 PMID: 22367081

27. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10: R25. doi: 10.1186/gb-2009-10-3-r25 PMID: 19261174

28. Cotillard A, Kennedy SP, Kong LC, Prifti E, Pons N, Le Chatelier E, et al. Dietary intervention impact on gut microbial gene richness. Nature. 2013; 500: 585–588. doi: 10.1038/nature12480 PMID: 23985875

29. Lander ES, Waterman MS. Genomic mapping by fingerprinting random clones: A mathematical analysis. Genomics. 1988; 2: 231–239. doi: 10.1016/0888-7543(88)90007-9 PMID: 3294162

30. Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012; 40: D109–D114. doi: 10.1093/nar/gkr988 PMID: 22080510

31. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. bioRxiv. 2014; 002832. doi: 10.1101/002832

32. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004; 5: R80. doi: 10.1186/gb-2004-5-10-r80 PMID: 15461798
33. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Ser B Methodol. 1995; 57: 289–300. doi: 10.2307/2346101

34. Larpin S, Mondoloni C, Goerges S, Vernoux JP, Guéguen M, Desmasures N. *Geotrichum candidum* dominates in yeast population dynamics in Livarot, a French red-smear cheese. FEMS Yeast Res. 2006; 6: 1243–1253. doi: 10.1111/j.1567-1364.2006.00127.x PMID: 17156021

35. Larpin-Laborde S, Imran M, Bonaïti C, Bora N, Gelsomino R, Goerges S, et al. Surface microbial con-
sortia from Livarot, a French smear-ripened cheese. Can J Microbiol. 2011; 57: 651–660. doi: 10.1139/
w11-050 PMID: 21815832

36. Vegarud GE, Langsrud T, Svenning C. Mineral-binding milk proteins and peptides; occurrence, bio-
chemical and technological characteristics. Br J Nutr. 2000; 84: 91–98.

37. Miethke M. Molecular strategies of microbial iron assimilation: from high-affinity complexes to cofactor assembly systems. Metalomics. 2013; 5: 15. doi: 10.1039/c2mt20193c PMID: 23192658

38. Blaiseau PL, Seguin A, Camadro JM, Lesuisse E. Iron Uptake in Yeasts. In: Cornelis P, Andrews SC, editors. Iron Uptake and Homeostasis in Microorganisms. Caister Academic Press; 2010. pp. 265–284.

39. Mor H, Barash I. Characterization of siderophore-mediated iron transport in Geotrichum candidum, a non-siderophore producer. Biol Met. 1990; 2: 209–213. PMID: 2143917

40. Monnet C, Loux V, Gibrat JF, Spinnler E, Barbe V, Vacherie B, et al. The Arthrobacter arilaitensis Re117 Genome Sequence Reveals Its Genetic Adaptation to the Surface of Cheese. Ahmed N, editor. PLoS ONE. 2010; 5: e15489. doi: 10.1371/journal.pone.0015489 PMID: 21124797

41. Monnet C, Loux V, Bento P, Gibrat JF, Straub C, Bonnarme P, et al. Genome Sequence of Corynebacterium casei UCMA 3821, Isolated from a Smear-Ripened Cheese. J Bacteriol. 2012; 194: 738–739. doi: 10.1128/JB.06496-11 PMID: 22247534

42. Aleksandrzak-Piekarczyk T. Lactose and β-Glucosides Metabolism and Its Regulation in Lactococcus lactis: A Review. In: Kongo JM, editor. Lactic Acid Bacteria—R & D for Food, Health and Livestock Purposes. InTech; 2013. Available: http://www.intechopen.com/books/lactic-acid-bacteria-r-d-for-food-health-and-livestock-purposes/lactose-and-glucosides-metabolism-and-its-regulation-in-lactococcus-lactis-a-review. Accessed 2015 Mar 19.

43. Rodicio R, Heinisch JJ. Yeast on the milky way: genetics, physiology and biotechnology of *Kluyveromyces lactis*: Kluyveromyces lactis. Yeast. 2013; 30: 165–177. doi: 10.1002/yea.2954 PMID: 23576126

44. Boutrou R, Guéguen M. Interests in Geotrichum candidum for cheese technology. Int J Food Microbiol. 2005; 102: 1–20. doi: 10.1016/j.ijfoodmicro.2004.12.028 PMID: 15924999

45. Arfi K, Leclercq-Perlat MN, Baucher A, Tâche R, Delettre J, Bonnarme P. Contribution of several cheese-ripening microbial associations to aroma compound production. Le Lait. 2004; 84: 435–447. doi: 10.1051/lait:2004016

46. Carvalho AL, Turner DL, Fonseca LL, Solopova A, Catarino T, Kuipers OP, et al. Metabolic and Trans-
scriptional Analysis of Acid Stress in Lactococcus lactis, with a Focus on the Kinetics of Lactic Acid
Pools. Zilberstein D, editor. PLoS ONE. 2013; 8: e68470. doi: 10.1371/journal.pone.0068470 PMID:
23844205

47. Lodi T, Alberti A, Guiard B, Ferrero I. Regulation of the Saccharomyces cerevisiae DLD1 gene encod-
ing the mitochondrial protein D-lactate ferricytochrome c oxidoreductase by HAP1 and HAP2/3/4/5. Mol
Gen Genet MGG. 1999; 262: 623–632.

48. Ramil E, Agrimonti C, Shechter E, Gervais M, Guiard B. Regulation of the CYB2 gene expression: tran-
scriptional co-ordination by the Hap1p, Hap2/3/4/5p and Adr1p transcription factors. Mol Microbiol.
2000; 37: 1116–1132. PMID: 10972830

49. Cholet O, Henaut A, Casaregola S, Bonnarme P. Gene Expression and Biochemical Analysis of Cheese-Ripening Yeasts: Focus on Catabolism of L-Methionine, Lactate, and Lactose. Appl Environ Microbiol. 2007; 73: 2561–2570. doi: 10.1128/AEM.02720-06 PMID: 17308183

50. Boutrou R, Kerriou L, Gassi JY. Contribution of Geotrichum candidum to the proteolysis of soft cheese. Int Dairy J. 2006; 16: 775–783. doi: 10.1016/j.idairyj.2005.07.007

51. Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat Methods. 2011; 8: 785–786. doi: 10.1038/nmeth.1701 PMID: 21959131

52. Nakai K, Horton P, Nakai K, Horton P. PSORT: a program for detecting sorting signals in proteins and predicting their subcellular localization. Trends Biochem Sci. 1999; 24: 34–35. doi: 10.1016/S0968-0004(98)01336-X PMID: 10087920

53. Milesi MM, McSweeney PLH, Hynes ER. Impact of Chymosin- and Plasmin-Mediated Primary Proteol-
ysis on the Growth and Biochemical Activities of Lactobacilli in Miniature Cheddar-Type Cheeses. J
Dairy Sci. 2008; 91: 3277–3290. doi: 10.3168/jds.2008-1197 PMID: 18765587
54. Savijoki K, Ingmer H, Varmanen P. Proteolytic systems of lactic acid bacteria. Appl Microbiol Biotechnol. 2006; 71: 394–406. doi: 10.1007/s00253-006-0427-1 PMID: 16628446

55. Risinger AL, Cain NE, Chen EJ, Kaiser CA. Activity-dependent reversible inactivation of the general amino acid permease. Mol Biol Cell. 2006; 17: 4411–4419. PMID: 16885415

56. Lubkowitz MA, Barnes D, Breslav M, Burchfield A, Naider F, Becker JM. Schizosaccharomyces pombe isp4 encodes a transporter representing a novel family of oligopeptide transporters. Mol Microbiol. 1998; 28: 729–741. PMID: 9643541

57. Lamarque M, Charbonnel P, Aubel D, Piard JC, Atlan D, Juillard V. A Multifunction ABC Transporter (Opt) Contributes to Diversity of Peptide Uptake Specificity within the Genus Lactococcus. J Bacteriol. 2004; 186: 6492–6500. doi: 10.1128/JB.186.19.6492-6500.2004 PMID: 15375130

58. Gagny B, Wiederkehr A, Dumoulin P, Winsor B, Riezman H, Haguenauer-Tsapis R. A novel EH domain protein of Saccharomyces cerevisiae, Ede1p, involved in endocytosis. J Cell Sci. 2000; 113: 3309–3319. PMID: 10954428

59. Kaksonen M, Toret CP, Drubin DG. A Modular Design for the Clathrin- and Actin-Mediated Endocytosis Machinery. Cell. 2005; 123: 305–320. doi: 10.1016/j.cell.2005.09.024 PMID: 16239147

60. Lin CH, MacGurn JA, Chu T, Stefan CJ, Emr SD. Arrestin-Related Ubiquitin-Ligase Adaptors Regulate Endocytosis and Protein Turnover at the Cell Surface. Cell. 2008; 135: 714–725. doi: 10.1016/j.cell.2008.09.025 PMID: 18976803

61. Walther TC, Brickner JH, Aguilar PS, Bernales S, Pantoja C, Walter P. Eisosomes mark static sites of endocytosis. Nature. 2006; 439: 998–1003. doi: 10.1038/nature04472 PMID: 16496001

62. Wilkinson D, Ramsdale M. Proteases and caspase-like activity in the yeast *Saccharomyces cerevisiae*. Biochem Soc Trans. 2011; 39: 1502–1508. doi: 10.1042/BST0391502 PMID: 21936842

63. Kumura H, Takagaki K, Sone T, Tsukahara M, Tanaka T, Shimazaki K. Casein digestion by Debaryomyces hansenii isolated from cheese. Biosci Biotechnol Biochem. 2002; 66: 1370–1373. doi: 10.1271/bbb.66.1370 PMID: 12162560

64. Mansour S, Bailly J, Delettre J, Bonnarme P. A proteomic and transcriptomic view of amino acids catabolism in the yeast *Yarrowia lipolytica*. PROTEOMICS. 2009; 9: 4714–4725. doi: 10.1002/pmic.200900161 PMID: 19764064

65. Gori K, Mortensen HD, Arneborg N, Jespersen L. Ammonia Production and Its Possible Role as a Mediator of Communication for Debaryomyces hansenii and Other Cheese-Relevant Yeast Species. J Dairy Sci. 2007; 90: 5032–5041. doi: 10.3168/jds.2006-750 PMID: 17954742

66. Molimard P, Spinnler HE. Review: Compounds Involved in the Flavor of Surface Mold-Ripened Cheeses: Origins and Properties. J Dairy Sci. 1996; 79: 169–184. doi: 10.3168/jds.S0022-0302(96)76348-8

67. Castellote J, Fraud S, Irlinger F, Swennen D, Fer F, Bonnarme P, et al. Investigation of Geotrichum canidum gene expression during the ripening of Reblochon-type cheese by reverse transcription-quintitative PCR. Int J Food Microbiol. 2015; 194: 54–61. doi: 10.1016/j.ijfoodmicro.2014.11.009 PMID: 25461609

68. Sourabié AM, Spinnler HE, Bourdat-Deschamps M, Tallon R, Landaud S, Bonnarme P. S-methyl thioesters are produced from fatty acids and branched-chain amino acids by brevibacteria: focus on l-leucine catabolic pathway and identification of acyl-CoA intermediates. Appl Microbiol Biotechnol. 2012; 93: 1673–1683. doi: 10.1007/s00253-011-3500-3 PMID: 21858675

69. Hickey DK, Kilcawley KN, Beresford TP, Wilkinson MG. Starter Bacteria Are the Prime Agents of Lipolysis in Cheddar Cheese. J Agric Food Chem. 2006; 54: 8229–8235. doi: 10.1021/jf060819h PMID: 17032033

70. Falentin H, Postollec F, Parayre S, Henaff N, Le Bivic P, Richoux R, et al. Specific metabolic activity of ripening bacteria quantified by real-time reverse transcription PCR throughout Emmental cheese manufacture. Int J Food Microbiol. 2010; 144: 10–19. doi: 10.1016/j.ijfoodmicro.2010.06.003 PMID: 20630608

71. Falentin H, Henaff N, Le Bivic P, Deutsch SM, Parayre S, Richoux R, et al. Reverse transcription quantitative PCR revealed persistency of thermophilic lactic acid bacteria metabolic activity until the end of the ripening of Emmental cheese. Food Microbiol. 2012; 29: 132–140. doi: 10.1016/j.fm.2011.09.009 PMID: 22029927

72. Rossi F, Gardini F, Rizzotti L, Gioia FL, Tabanelli G, Torriani S. Quantitative Analysis of Histidine Deboxylase Gene (hdcA) Transcription and Histamine Production by Streptococcus thermophilus PR160 under Conditions Relevant to Cheese Making. Appl Environ Microbiol. 2011; 77: 2817–2822. doi: 10.1128/AEM.02531-10 PMID: 21378060
73. Monnet C, Straub C, Castellote J, Onesime D, Bonnarme P, Irlinger F. Quantification of Yeast and Bacterial Gene Transcripts in Retail Cheeses by Reverse Transcription-Quantitative PCR. Appl Environ Microbiol. 2013; 79: 469–477. doi: 10.1128/AEM.02360-12 PMID: 23124230

74. Cretenet M, Gall GL, Wegmann U, Even S, Shearman C, Stentz R, et al. Early adaptation to oxygen is key to the industrially important traits of Lactococcus lactis ssp. cremoris during milk fermentation. BMC Genomics. 2014; 15: 1054. doi: 10.1186/1471-2164-15-1054 PMID: 25467604

75. Cretenet M, Laroute V, Ulvé V, Jeanson S, Nouaille S, Even S, et al. Dynamic Analysis of the Lactococcus lactis Transcriptome in Cheeses Made from Milk Concentrated by Ultrafiltration Reveals Multiple Strategies of Adaptation to Stresses. Appl Environ Microbiol. 2011; 77: 247–257. doi: 10.1128/AEM.01174-10 PMID: 21075879

76. Irlinger F, Loux V, Bento P, Gibrat JF, Straub C, Bonnarme P, et al. Genome Sequence of Staphylococcus equorum subsp. equorum Mu2, Isolated from a French Smear-Ripened Cheese. J Bacteriol. 2012; 194: 5141–5142. doi: 10.1128/JB.01038-12 PMID: 22933766