I Ammoun

C I Kothe

N Mohellibi

C Beal

R Yaacoub

P Renault

Lebanese fermented goat milk products: From tradition to meta-omics

Keywords: Microbial diversity, kefir, metagenomics, genomics, ecology, VBNC, metagenomeassembled genomes, LAB, Lactobacillus kefiranofaciens, Streptococcus parasuis, Mycoplasmopsis agalactiae

Detailed analysis of the production chain of three emblematic Lebanese fermented goat milk products -Ambriss and Serdaleh have a microbiota composition similar to that of kefir -Adaptation of the zoonotic agent Streptococcus parasuis to milk fermentation through the acquisition of lactose utilization genes -Evidence of herd infection by Mycoplasmopsis agalactiae through MAG analysis 1 Abstract

Ambriss, Serdaleh and Labneh El Darff are traditional Lebanese products made from fermented goat's milk from the Bekaa Valley, the Chouf Mountains and Southern Lebanon, respectively. A questionnaire completed by 50 producers of these products showed that they are prepared by periodic percolation either by milk or by Laban in amphorae or goatskins during the lactation season. Production is carried out on a small scale and in a limited number of production units, often by elderly people, resulting in a real risk of disappearance of these products and loss of the corresponding microbial resources will be lost. In this study, 34 samples from 18 producers were characterized by culture-dependent and -independent analyses. The results obtained by these two methods were radically different, especially for Ambriss and Serdaleh, the latter revealing the co-dominance of Lactobacillus kefiranofaciens, a fastidious-growing species, and Lactococcus lactis in a viable but not culturable state. Overall, their composition is reminiscent of kefir grains. However, phylogenomic and functional analyses of the genomes of the key species Lb. kefiranofaciens show that they differed from those of kefir, for example in their polysaccharide genes, explaining the absence of grains. Labneh El Darff displays a dominance of Lactobacillus delbrueckii, probably due to the addition of Laban. In addition, several zoonotic pathogens were found, including Streptococcus parasuis, which was dominant in one sample and which MAG analysis revealed it acquired lactose utilization genes through horizontal gene transfer. The contamination of the herd with Mycoplasmopsis agalactiae in the Chouf region was also revealed by MAG analysis of the Serdaleh samples. Antibiotic resistance genes were detected in most of the samples, particularly in the Serdaleh samples, where the dominant L. lactis strains possessed a plasmid with a multi-resistance island. Finally, this study paves the way for further analyses to shed light on the resilience of these ecosystems established in amphora or in goatskins and to improve hygiene practices for milk production.

Introduction

The Mediterranean region is well known for its goat farming, providing about 21% of the world's goat milk supply (Miller & Lu, 2019;[START_REF] Dubeuf | The Mediterranean sheep and goat sectors between constants and changes over the last decade future challenges and prospects[END_REF]. Due to the nature of its composition, organoleptic characteristics and healthy attributes, goat milk has been the focus of interest in recent years [START_REF] Prosser | Compositional and functional characteristics of goat milk and relevance as a base for infant formula[END_REF][START_REF] Clark | A 100-Year Review: Advances in goat milk research[END_REF]. This sector has received increasing attention, particularly in Lebanon, and continues to grow [START_REF] Kayouli | LEBANON RECOVERY FUND Recovery and Rehabilitation of the Dairy Sector in Bekâa Valley and Hermel-Akkar Uplands1-FINAL PROGRAMME-LRF-21 OSRO/LEB/901/UNJ[END_REF]. Goat milk products are highly appreciated by local consumers [START_REF] Balaa | Sustainability of the Lebanese small ruminant dairy products supply chain[END_REF], especially for their specific "goaty taste" [START_REF] Semaan | Caractérisation chimique et Qualité Bactériologique de Produits Laitiers Caprins Traditionnels Libanais[END_REF]. Compared to cow milk, its casein proteins are more easily digested and it is better tolerated by people with mild lactose intolerance. It has a higher content of short-chain fatty acids in milk fat [START_REF] Shu | Fermentation optimization of goat milk with Lactobacillus acidophilus and Bifidobacterium bifidum by Box-Behnken design[END_REF][START_REF] Zenebe | Review on Medicinal and Nutritional Values of Goat Milk[END_REF] and contains minerals such as magnesium, calcium and phosphorus [START_REF] References Abbas | Physicochemical Characteristics of Goat's Milk[END_REF].

These goat milk products are made either from standardized milk at the industrial scale, such as Laban (Chammas et al., 2006a), Labneh [START_REF] Serhan | Concentrated yogurt (Labneh) made of a mixture of goats' and cows' milk: Physicochemical, microbiological and sensory analysis[END_REF], Shankleesh [START_REF] Serhan | Characterization of four lebanese artisanal goat milk cheeses: Darfiyeh, Aricheh, shankleesh and Serdale by physico-chemical, microbiological and sensory analyses[END_REF][START_REF] Toufeili | Ripening changes and sensory properties of bovine, caprine and ovine shankleesh[END_REF], Kishk [START_REF] Salameh | An overview of recent studies on the main traditional fermented milks and white cheeses in the Mediterranean region[END_REF], or artisanally from raw milk, such as Darfiyeh [START_REF] Hosri | Valoriser le fromage de chèvre traditionnel " Darfiyeh " pour aider au developpement de la région montagnarde nord libanaise[END_REF], Aricheh [START_REF] Serhan | Characterization of four lebanese artisanal goat milk cheeses: Darfiyeh, Aricheh, shankleesh and Serdale by physico-chemical, microbiological and sensory analyses[END_REF], Serdaleh [START_REF] Tabet | Innovation aspects of Serdaleh , a traditional Lebanese cheese produced from raw extensive goat ' s milk. Innovations for Sustainability in Sheep and Goats[END_REF], Ambriss [START_REF] Dimassi | Production and physicochemical properties of labneh anbaris, a traditional fermented cheese like product, in Lebanon[END_REF] and Labneh El Darff [START_REF] Semaan | Caractérisation chimique et Qualité Bactériologique de Produits Laitiers Caprins Traditionnels Libanais[END_REF]. The diversity of these dairy products not only enriches the Lebanese culinary culture and pantry [START_REF] Massaad | Mouneh: Preserving Foods for the Lebanese Pantry[END_REF], but also represents an important source of empowerment for women in poor rural areas. Female-headed households effectively contribute to the development of the goat rearing sector and the manufacturing of goat milk products [START_REF] Kayouli | LEBANON RECOVERY FUND Recovery and Rehabilitation of the Dairy Sector in Bekâa Valley and Hermel-Akkar Uplands1-FINAL PROGRAMME-LRF-21 OSRO/LEB/901/UNJ[END_REF]. The market share of this sector is not well characterized; however, 70% of all small ruminant industries are located in rural areas with limited market access [START_REF] Saadeh | Composition of Dairy Potential Output per Year Lebanese Dairy Sector[END_REF]. Among these fermented goat milk products, Ambriss, Serdaleh and Labneh El Darff, made with milk from the Black Baladi race that grazes in the Bekaa Valley, the Chouf Mountains and Southern Lebanon, respectively (Figure 1A), are real candidates to be awarded the PGI (Protected Geographical Origin) label according to the report of the SRVA (Service Romand de Vulgarisation Agricole; IDEAS [START_REF] Centre | Proposal for the protection of Geographical indication in Lebanon Report on the intention of IDEAS Centre 1204[END_REF]; Ministère de l'Economie et du Commerce de la République Libanaise, 2006).

However, it appears that these products are in danger of disappearing, mainly due to commercial, regulatory and food safety concerns [START_REF] Serhan | The Goat Dairy Sector in Lebanon[END_REF][START_REF] Semaan | Caractérisation chimique et Qualité Bactériologique de Produits Laitiers Caprins Traditionnels Libanais[END_REF]. In fact, these products are made on a seasonal basis, at small or medium production scales and according to traditional methods (Figure 1B), by successive fermentation of raw goat milk in amphorae for Ambriss [START_REF] Dimassi | Production and physicochemical properties of labneh anbaris, a traditional fermented cheese like product, in Lebanon[END_REF] and Serdaleh [START_REF] Serhan | Characterization of four lebanese artisanal goat milk cheeses: Darfiyeh, Aricheh, shankleesh and Serdale by physico-chemical, microbiological and sensory analyses[END_REF] and Laban in goatskins for Labneh El Darff [START_REF] Serhan | Concentrated yogurt (Labneh) made of a mixture of goats' and cows' milk: Physicochemical, microbiological and sensory analysis[END_REF]Nasbimana et al. 2005).

In the Mediterranean region, traditional fermented milk products are closely related to the area in which they are produced and the cultural history and heritage [START_REF] Boyazoglu | Mediterranean dairy sheep and goat products and their quality: A critical review[END_REF]. This fact has drawn the attention of many studies due to the unique microbial diversity of these products [START_REF] Rhaiem | Identification and characterization of lactic acid bacteria isolated from cow milk and olives brine[END_REF][START_REF] Abdalla | Enumeration and Identification of Microflora in Roub, A Sudanese Traditional Fermented Dairy Product[END_REF][START_REF] Ayad | Selection of wild lactic acid bacteria isolated from traditional Egyptian dairy products according to production and technological criteria[END_REF][START_REF] El Soda | Isolation, identification and selection of lactic acid bacteria cultures for cheesemaking[END_REF]. Numerous genera of lactic acid bacteria (LAB) such as Lactobacillus, Lactococcus, Enterococcus, Leuconostoc and Streptococcus are generally involved in these transformations. However, only a few studies have described the native microbiota in the Lebanese fermented milk products, and the characterization of the isolates has been restricted to phenotypic and biochemical analyses. The first attempts were recorded on Laban, where 96 strains of S. thermophilus and L. bulgaricus were characterized (Chammas et al., 2006b). [START_REF] Serhan | Bacterial diversity of Darfiyeh, a Lebanese artisanal raw goat's milk cheese[END_REF] then revealed the presence of the LAB Streptococcus thermophilus, Enterococcus spp., Lactococcus lactis subsp. lactis and cremoris, Lactiplantibacillus plantarum, Lactobacillus curvatus in Darfiyeh cheese using culture-dependent methods. Moreover, [START_REF] Dib | Identification et évaluation de l'effet probiotique des bactéries lactiques isolées dans des fromages caprins traditionnels[END_REF] isolated 25 different strains in Serdaleh, Laban and Darfiyeh cheese belonging to the species Lp. plantarum, Lacticaseibacillus paracasei, Lacticaseibacillus casei and Leuconostoc mesenteroides Until now, these traditional Lebanese products have only been characterized using cultural methods, which provide a patchy view of the microbiota, whereas new methods for characterizing food ecosystems are now available, based on techniques such as 16S and shotgun metagenomic analyses [START_REF] Walsh | Bioinformatic approaches for studying the microbiome of fermented food[END_REF]. Approaches based on these techniques have already been used in a study on traditional Brazilian cheeses and provided a global picture of their microbiota, with taxonomic resolution close to the strain level for dominant populations and functional data at the gene level [START_REF] Kothe | Unraveling the world of halophilic and halotolerant bacteria in cheese by combining cultural, genomic and metagenomic approaches[END_REF]. We therefore decided to develop a better understanding of the microbiota of Ambriss, Serdaleh and Labneh El Darff using these new methods to improve our knowledge of their processes, and to thus develop strategies to preserve them and boost their development. To do this, 34 samples of fermented goat milk products from 18 producers were characterized using a combination of culture-dependent and cultureindependent analyses, including genomic analyses of isolated strains and metagenomic assemblies, as well as functional analyses based on deep metagenomic profiling. The results were interpreted with regard to the manufacturing processes and physico-chemical characteristics of the products.

Materials and Methods

Development of a producer's survey

A producers' survey was conducted over 50 breeders or producers from the three selected production regions -the Bekaa Valley, the Chouf Mountains and Southern Lebanon -for the artisanal products of Ambriss, Serdaleh and Labneh El Darff, respectively. The questionnaire consisted of 34 questions (Supplementary datafile 1), covering three main areas: (i) sociodemographic aspects; (ii) breeding activities; and (iii) production activities, including commercialization. The questions included general and socio-economic information about producers, information regarding livestock size, goat breed, breeding type and practices, milk production volumes, as well as the type of dairy product produced, the production season, the raw materials used, production volumes, prices and marketing practices. In addition, a detailed observation of the processing was carried out in order to build a common production flow chart for each product type, showing the existing differences between the producers as well.

Sample collection, physic-chemical and microbiological analyses

Producers were first selected on the basis of the survey at the different production sites according to two main criteria: their production was small-scale and they had been producing their product for more than 15 years. In total, 34 artisanal samples from 18 producers were collected during two consecutive production seasons in 2018 and 2019 (Table S1).

Samples of fermented goat milk products (10 g) were taken aseptically and homogenized for 60 s in a laboratory blender (Blender LB20E, Waring, Torrington, CT, USA) with 90 mL and 20 mL of distilled water, for pH and salinity tests, respectively. The pH was measured using a digital pH meter (GLP 21,Crison instrument,Barcelona,Spain). The salinity content was determined using chloride test strips (30-600 mg/L; QuanTabTitrator, Hach, Loveland, CO, USA). All measurements were performed in triplicate. Samples were preserved in phosphatebuffered saline (PBS 15X, pH 8.7), and stored at -20°C for further analysis.

For culture-dependent analyses, samples were diluted and homogenized in buffer peptone water (2%). They were spread in triplicate on MRS (pH 6.2, Sigma-Aldrich, St. Louis, MO, USA) and M17 (pH 6.5, Sigma-Aldrich, St. Louis, MO, USA) agar plates. Plates were incubated for 48 h, anaerobically for MRS at 30°C and 42°C to allow mesophilic and thermophilic Lactobacillus counts, respectively, and aerobically for M17 at 30°C for lactococci. Specific procedures were applied to isolate L. lactis and Lactobacillus kefiranofaciens strains. Culture media were supplemented with 2% salt and 10% goat milk for L. lactis, or with 25% goat milk and adjusted at pH 5.3 for Lb. kefiranofaciens, and incubated anaerobically at 30°C for 10 days.

Based on their distinct morphology (shape, color and size), colonies were selected and purified by repetitive streaking on the same medium. Isolates were assessed for their catalase activity and Gram coloration according to the criteria described by Kandler and Weiss (1986), as well as for their ability to produce gas from glucose by inverted tubes (Zúñiga et al., 1993).

Biochemical identification and carbohydrate fermentation patterns were determined by API 50 CH and API 20 STREP test kits (BioMérieux Inc., Marcy l'Etoile, France) following the manufacturer's recommendations.

In order to identify the different morphotypes, isolates were suspended in 400 µL Tris-EDTA (TE 1X, pH 8.0) buffer in a 2-mL tube. DNAs was lysed using a bead-beater ( 24Microtubes BeadBlaster, Thomas Scientific, NJ., USA), with 100 mg zirconium beads of 0.1 and 0.5-mm diameters and homogenized for 20 s at 4.5 m/s. Once extracted, amplification of the coding gene for 16S rRNA was carried out for 30 cycles on a final mix of 25 µL containing the following universal primers: Bact-27F (5'-AGAGTTTGATCCTGGCTCAG-3') and 1492R (5'-GGTTACCTTGTTACGACTT-3)', using Taqara DNA Taq polymerase [START_REF] Almeida | Construction of a dairy microbial genome catalog opens new perspectives for the metagenomic analysis of dairy fermented products[END_REF].

Products were purified using a purification kit (Gen Elute PCR clean up, Sigma, St. Louis, MO, USA). Genotypic identification was carried out by blasting the different sequences of the isolated strains using the NCBI database.

Genomic DNA extraction and sequencing

DNA from five isolates, chosen as representative of the main species detected in our products, was extracted from the bacterial cells according to the protocol described by [START_REF] Kothe | Unraveling the world of halophilic and halotolerant bacteria in cheese by combining cultural, genomic and metagenomic approaches[END_REF]. DNA sequencing was carried out on an Illumina HiSeq at GATC-Biotech (Konstanz, Germany) to generate about 5-6 million paired-end reads (150 bases in length). For each strain, the paired-end reads were merged and de novo assembly was performed using SPAdes 3.9 [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF]. Only contigs with length > 300 bp and coverage > 100 were considered for further study. Annotations were performed with Rapid Annotation using the Subsystem Technology server [START_REF] Aziz | SEED Servers: High-Performance Access to the SEED Genomes, Annotations, and Metabolic Models[END_REF].

DNA extraction from the products

Samples (0.5 g) were diluted 1:1 (w/v) in guanidinium thiocyanate 4M solution (Sigma-Aldrich, MO, USA) with TrisHCl 0.1M and processed in the same way as for the isolate DNA extraction. The mixture was then supplemented with 10% N-Lauryl sarcosine (Sigma-Aldrich, MO, USA), vortexed and centrifuged at 4°C and 14,000 g for 30 min. Supernatant were discarded and extraction was continued on the emerged pellet, as described by [START_REF] Almeida | Construction of a dairy microbial genome catalog opens new perspectives for the metagenomic analysis of dairy fermented products[END_REF]. DNA was visually evaluated on 0.8% agarose gel and quantified by a fluorometer (Qubit 2.0, ThermoFisher Scientific) supplied with the Qubit dsDNA HS Assay Kit (ThermoFisher Scientific).

Metagenetic analysis

Bacterial diversity was analyzed by sequencing the amplified region V3-V4 of the 16S rRNA gene using primers V3F and V4R (V3F: 5′-CCTACGGGAGGCAGCAG-3′/ V4R: 5′-GGACTACHVGGGTWTCTAAT-3′). The PCR was performed with MTP Taq DNA Polymerase (Sigma-Aldrich, St. Louis, MO, USA), and cycling conditions were: 94°C for 1 min, followed by 30 cycles of amplification at 94°C for 1 min, 65°C for 1 min, and 72°C for 1 min, with a final extension step of 10 min at 72°C. Sequencing was performed with the V3 Illumina MiSeq kit, as described in Poirier et al. (2018), and fasta files were generated at the end of the run (MiSeq Reporter software, Illumina, USA).

The quality of the raw data was evaluated with FastQC (Wingett & Andrews, 2018) and the sequences were imported into the FROGS (Find Rapidly OTUs with Galaxy Solution) pipeline [START_REF] Escudié | FROGS: Find, Rapidly, OTUs with Galaxy Solution[END_REF] to obtain the Operational Taxonomic Units (OTUs). The sequences were filtered by length (480-580 bp) and then pooled into OTUs with SWARM [START_REF] Mahé | Swarm: Robust and fast clustering method for amplicon-based studies[END_REF] with the distance parameter of 1. Chimeras were removed with VSEARCH [START_REF] Rognes | VSEARCH: A versatile open source tool for metagenomics[END_REF] and OTUs with at least 0.005% in the whole dataset were retained [START_REF] Bokulich | Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing[END_REF]. OTUs were affiliated with the 16S SILVA 138.1 database [START_REF] Quast | The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools[END_REF]. Alphadiversity and beta-diversity analyses were performed with R Studio v.3.6.1 using the phyloseq and ggplot2 packages (v1.30.0). Statistics were performed on alpha-diversity using Kruskal-Wallis and linear model regression. DESeq2 was applied to characterize statistically significant differentially abundant OTUs in the different products [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]. The PCoA was plotted with the ward.D2 method and Permanova tests were performed between samples.

Taxonomic composition by shotgun metagenomics study

Sixteen fermented products (six Ambriss, six Serdaleh and four Labneh) were selected on the basis of the metagenetic results for more in-depth analysis. The metagenomic DNA was sequenced using Illumina HiSeq2500 technology at GATC-Biotech (Konstanz, Germany).

Metagenomic reads corresponding to the Capra hircus genome were filtered out by reads mapping against the Capra hircus reference genome (LWLT00000000) using Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF]. The proportion of goat reads in the samples was then assessed using samtools flagstat [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]. On the basis of the fastq file generated after host read filtration, we first estimated microbial composition by mapping the samples. Reads were mapped against the representative clade-specific marker catalogue contained in the MetaPhlAn tool v.

3.0.4 [START_REF] Truong | MetaPhlAn2 for enhanced metagenomic taxonomic profiling[END_REF].

In addition, we performed taxonomic profiling using an assembly-based marker gene analysis, which allows non-supervised binning of metagenomes using metaSPAdes v.3.9 [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF]. Genes were then predicted using Prodigal (v.2.6.3) and marker genes were extracted using fetchMG, v.1.0 [START_REF] Sunagawa | Metagenomic species profiling using universal phylogenetic marker genes[END_REF]. Taxonomic assignation was carried out with the COG0012, ychF gene. Summary species composition plots were created in R v.3.6.1 using the package ggplot2 v.3.3.2.

Metagenome-assembled genome (MAG) analyses

Genome binning was performed using MetaBAT2-2.12.1 [START_REF] Kang | MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities[END_REF], with a minimum contig size of 1,500 nucleotides and the default settings. The quality of the resulting prokaryotic bins was assessed with CheckM [START_REF] Parks | CheckM: Assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes[END_REF], and MAGs < 80% completeness and/or > 10% contamination were excluded. Several MAG assemblies were refined by additional binning, such as blastn analysis with reference genomes of related species. Contigs were then filtered by percentage of the length covered, identity and coverage levels. Finally, manual curing of questionable contigs (new, low coverage of homology on references, etc.) were performed by blastn and blastx analysis against the nr/nt and nr NCBI databases, respectively. Annotations were performed using the Rapid Annotation with Subsystem Technology server.

Phylogenomic and functional analyses

The Automatic Multi-Locus Species Tree web server [START_REF] Alanjary | AutoMLST: An automated web server for generating multi-locus species trees highlighting natural product potential[END_REF] was used to determine closely related genomes based on core gene alignments of the sequenced genomes and recovered MAGs. The closest species were inferred based on the percentage of average nucleotide identity (ANI) calculated using FastANI, v.1.31 (Jain et al., 2018). The kSNP (v.3.0) was used to perform the phylogenomic analysis using distances based on single nucleotide polymorphism variations for the most representative species present in MAGs [START_REF] Gardner | kSNP3.0: SNP detection and phylogenetic analysis of genomes without genome alignment or reference genome[END_REF] with the maximum likelihood option and kmer size of 21. Genomes of the corresponding species were downloaded from Genbank and combined with our MAG to compute this analysis.

Gene content analysis and metabolic reconstruction was performed on SEED using sequence and function-based comparison tools [START_REF] Aziz | SEED Servers: High-Performance Access to the SEED Genomes, Annotations, and Metabolic Models[END_REF]. Functional analysis on metagenomes was performed by read mapping profiling with the sequences of the 16 metagenomes from which the Cellulosimicrobium cellulans and goat reads were removed to avoid bias in the gene quantification. The analysis of microbial pathways was performed by HUMAnN 3.0 [START_REF] Beghini | Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with biobakery 3[END_REF](REF) in order to evidence bacteriocin production genes, potentially health-associated gene clusters or aroma pathways. To screen the potentially health-associated gene clusters, we followed essentially a methodology described previously [START_REF] Leech | Fermented-Food Metagenomics Reveals Substrate-Associated Differences in Taxonomy and Health-Associated and Antibiotic Resistance Determinants[END_REF]. except the functional profiles were established with HumanN3.

These genes were binned into one of the following categories: survival, modulation, and colonization in order to create a Heatmap using Pheatmap. Clusters counts displayed in the heatmap were not absolute counts but scaled so that the values could be compared between the samples. The genes potentially implied in bacteriocin production were searched with the term "bacteriocin" and those of the diacetyl pathways with "Acetoin, Diacetyl, Acetolactate, Citrate synthase, Citrate lyase, Citrate \[pro-3S\]-lyase and Oxaloacetate decarboxylase".

The search for antibiotic resistance (ABR) genes was performed by read mapping against the CARD database [START_REF] Jia | CARD 2017: Expansion and model-centric curation of the comprehensive antibiotic resistance database[END_REF][START_REF] Mcarthur | The comprehensive antibiotic resistance database[END_REF] using the PATRIC web server (Antonopoulos et al., 2019). The presence of particular genes in reference genomes, ABR and genes of technological interest was determined using the PATRIC web server, the FoodMicrobiome tool (https://migale.jouy.inra.fr/foodMicrobiome/, and blast analysis.

Additionally, FoodMicrobiome was used to detect the subdominant populations of given species with high sensitivity and reliability and to discriminate between subspecies. 10 3.9 Data availability Raw sequences were deposited on the European Nucleotide Archive (ENA) under the BioProject ID PRJNA860091 for shotgun metagenomic, amplicon and genomic reads. The MAGs are available in doi: 10.17632/7w4rd76ndj.2.

Results

Goat farmer and dairy producer's survey

According to the results of the questionnaire, 25% of the respondents were identified as both producers and farmers, and 75% as producers only. Production depended on female housekeepers (60%) aged between 40 and 70 years old, and had limited access to the local market place (38%). Most of producers (88%) used only goat milk compared to 4% who used a mixture of goat and cow milk, while the rest used goat Laban in addition to goat milk as the primary dairy material. Producers who reared their own goat herd were able to process larger milk quantities for the three fermented products (135 to 7,500 L/season) than producers who purchased their goat milk (140 to 1,890 L/season). Farming and production parameters for each product are summarized in Table S1.

Ambriss, Serdaleh and Labneh El Darff were made with milk from the Black Baladi race that grazes in the Bekaa Valley, the Chouf Mountains and Southern Lebanon, respectively (Figure 1A). All production phases were completely manual, the production scale was small to medium (Third Quartile = 216.5 kg/season) and both tools and materials were artisanal.

Production takes place during the period when the goats produce milk, i.e., from June to November. The manufacturing process of Serdaleh and Ambriss lasts up to 6 months and consists of fermenting raw goat milk on a daily basis in terracotta jars, except for 20% of Serdaleh producers who used plastic barrels. The jars were first soaked with salt, periodically refilled with fresh raw milk and salt, with periodic whey evacuation from a hole at their bottom (Figure 1B andS1). The manufacturing of Labneh El Darff is conducted in a three-stage process consisting of: (i) Laban production, which is prepared from raw goat milk inoculated with a previously homemade inoculum produced by backslopping; (ii) filling the "Darff" (bag formed by the whole goatskin, but sometimes cloth fabric) with the produced Laban; and (iii) straining of the whey from the suspended Darff during a period of 6 weeks (Figure 1B andS2). The manufacturing processes of these products are carried out at ambient temperature.

Physico-chemical and culture-dependent analyses

Temperature, pH and salinity

Temperature, pH and salinity levels were measured at the time the samples were collected (Table S2). The ambient temperature at which the products were developed varied between 28°C and 33°C for Ambriss and Labneh El Darff, but was lower (p < 0.05) for Serdaleh (between 20°C to 25°C). The pH values showed little variability within the same product type, unlike salinity, which displayed greater intra-group variability (Figure 2A). The pH values were comprised between 3.01 and 4.80. They were higher (p < 0.05) in Ambriss (4.08 to 4.80) when compared to Serdaleh (3.01 to 4.12) and Labneh El Darff (3.62-4.62). The Serdaleh product had the highest level (p < 0.05) of salinity (2.69 to 5.04%) compared to Ambriss (0.87 to 2.24%) and Labneh El Darff (1.10 to 1.75%).

Culture-dependent analyses

The examination of culture results shows a significant variability in the different microbial enumerations (Figure 2B). We observed a dominance of mesophilic lactobacilli (MRS 30°C, 5.4-6.8 Log CFU/g) over lactococci (M17 30°C, 0.5-5.20 Log CFU/g) in the three products.

Thermophilic lactobacilli counts were highly variable in Labneh el Darff (MRS 42°C, 1-6.2 Log CFU/g) compared to its other counterparts.

Overall, 279 strains belonging to 13 species groups were genetically identified by 16S rRNA from 380 isolates. The most commonly cultivated species in the three fermented milks was Lactiplantibacillus plantarum. This species appeared to be dominant in Ambriss (38% of isolates in this product) and Serdaleh (77%), and co-dominant with other species in Labneh El Darff (16%). In Ambriss, we also identified a strong prevalence of the species Enterococcus durans and Lacticaseibacillus rhamnosus (18% of each), in Labneh El Darff, Levilactobacillus brevis was the most frequently obtained cultured species (32% of isolates), followed by Lb. delbrueckii and L. lactis (16% of each) (Figure 2C).

Taxonomic composition using amplicon sequencing

The 16S rRNA amplicon-sequencing analyses were clustered in 209 OTUs, whose abundance was over the threshold of 0.005% in the whole dataset. The rarefaction curves of most of the analyzed samples reached the saturation plateau, indicating that the sequencing depth was sufficient (Figure S3).

The alpha-diversity analysis revealed lower richness of species in Labneh El Darff (p < 0.05) compared to the other two products (Figure 3A). Moreover, the Chao1, Shannon and InvSimpson indices, which measure evenness, confirmed this low diversity. The Principal Coordinates Analysis (PCoA) showed that while Labneh El Darff samples were clearly differentiated, samples of Serdaleh and Ambriss could not be clearly separated and thus likely shared a similar bacterial microbiota (Figure 3B). Indeed, the analysis of their taxonomical composition (Figure 3C and Table S3) showed that Lactobacillus delbrueckii was the main species present in Labneh El Darff samples, representing 55 to 82% of the counts, whereas it was undetected or at levels below 0.4% in the samples of Ambriss and Serdaleh. Moreover, the samples of the latter products contained L. lactis, Lb. kefiranofaciens and Lentilactobacillus kefiri as the major species and, to a lesser extent, Lentilactobacillus diolivorans-group. Differential analysis of OTUs showed that Serdaleh samples had greater L. lactis abundance than Ambriss (padj < 0.01). Remarkably, in the PCoA (Figure 3B), these two products appeared to be spread along Axis 1 as a function of the proportion between Lactococcaceae and Lactobacillaceae species. In particular, a group of 12 Serdaleh samples contained L. lactis at levels of 58-83%, while the other samples were intermixed on the other part of this axis.

Moreover, the 13 samples of Serdaleh collected in 2019 contained low levels (0.67-2.2%) of Mycoplasmopsis agalactiae or M. bovis (Table S3). Finally, three Ambriss samples showed divergent compositions compared to the others: A3-b was dominated by Lentilactobacillus buchneri, A4 by the S. parasuis-sanguinis group (> 53% of the reads), and A6 by high amounts of species poorly present in other samples such as L. gasseri (~36%), Lc rhamnosus (~23%) and Enterobacteriaceae (> 18% of reads). The other species, usually detected at low levels in the various samples, mainly belonged to the Gammaproteobacteria class, such as Aeromonas sp., Citrobacter sp., Klebsiella sp. and Escherichia-Shigella.

Refined microbiota by metagenomics

For a refined taxonomical analysis, 16 samples were selected to perform shotgun metagenomic approaches. Such analyses allow a taxonomic identification up to the species-strain level and the simultaneous assessment of the level of reads corresponding to bacteria and fungi.

Moreover, these analyses make it possible to detect eventual virulence factors in pathogenic species, antibiotic resistance genes and genes of technological interest.

Food metagenomic samples may contain animal reads, which may significantly bias the analysis when they are in high abundance. Since the present products are manufactured from goat milk, we tested for the presence of Capra hircus reads that were detected at levels ranging from 2.3 to 72% (Table S4). They were notably abundant in Labneh El Darff samples where they accounted for more than 50%, as a consequence of the process mode where the product is drained in a goatskin.

Species-level assignments by shotgun metagenomics

In a first approach, the MetaPhlAn tool, which identifies species by mapping the metagenomic reads against a genomic database, was used to determine the taxonomic composition of the samples. It yielded 98 bacterial and one eukaryotic species (Table S5).

Species composition obtained by this analysis was relatively in accordance with those obtained with the 16S rRNA amplicon. Interestingly, it resolved most taxonomical assignation ambiguities for lactic acid bacteria such as Lt. kefiri, Lt. diolivorans, Lactiplantibacillus plantarum, Lactobacillus helveticus and S. thermophilus, and also for contaminants such as S. parauberis and M. agalactiae. It also detected the yeast Saccharomyces cerevisiae, but at low levels and in only three samples.

In order to characterize microbiota independently of fixed references, we applied a marker gene analysis from the assembled metagenomes, which yielded 104 bacterial and 15 fungal species (Table S6). The overall composition is shown in Figure 4. Although the main bacterial species identified here were consistent with the previous analysis made by read profiling, some differences merit highlighting. Whereas MetaPhlAn attributed 35% of the reads of A4 to S. parauberis, the marker gene analysis detected only 6.9% of this species and further detected another streptococcal species, S. parasuis, which retained 55% of the reads. This species was also detected at low levels in four other samples, while S. parauberis was found in eight samples, generally at low levels. This analysis also confirmed the frequent presence at low levels of several species of Enterobacteriaceae such as E. coli, K. pneumoniae and C. freundii, and of Enterococcaceae such as E. faecalis and E. faecium. Lastly, it also established that L. lactis subsp. lactis was detected in all samples, whereas the cremoris subspecies was detected at concentrations ten times lower in S8 and S10-b. A specific analysis by read mapping confirmed the presence of L. lactis subsp. cremoris at 5-15% lower abundance than that of the lactis subspecies in most samples.

Finally, 14 species of yeasts whose genomes are absent from the MetaPhlan database were identified by marker gene analysis. The most frequent and abundant were Pichia kudriavzevii and, to a lesser extent, S. cerevisiae and Kluyveromyces marxianus. Nevertheless, yeast relative abundance was generally below 1%, except in Labneh El Darff where it reached 8%. A more sensitive analysis performed by read mapping (Table S7) confirmed the presence of these yeasts at low levels in more samples, as well as that of Torulaspora delbrueckii.

MAGs and strain-level analyses of LAB

A total of 62 prokaryotic Metagenome Assembled Genomes (MAGs) of high quality were reconstituted from the metagenome assemblies (Table S8). The majority of the MAGs correspond to Lactobacillaceae such as Lt. diolivorans (12), Lt. kefiri (10), Lb. kefiranofaciens (9), Lp. plantarum (5) and L. gasseri [START_REF]) for providing computing and storage resources[END_REF]. The two yoghurt species Lb. delbrueckii subsp.

bulgaricus [START_REF]) for providing computing and storage resources[END_REF] and S. thermophilus [START_REF]) for providing computing and storage resources[END_REF] were also recovered from Labneh El Darff and MAGs of L. lactis subsp. lactis (5) were reconstituted from Ambriss and Serdaleh products. In addition, the results included MAGs of the goat pathogen M. agalactiae (6) from Serdaleh samples, and of milk contaminants such as S. parauberis (3), S. parasuis [START_REF]) for providing computing and storage resources[END_REF], Escherichia coli (1), Klebsiella pneumoniae [START_REF]) for providing computing and storage resources[END_REF], Aeromonas caviae (2) and Acetobacter spp. (2).

Further detailed analysis was performed with the lactococci and lactobacilli MAGs and genome sequences of the most frequent and sometimes dominant species in the three products (Table S9). Phylogenomic analysis of the L. lactis FME150 genome and the five L. lactis MAGs was performed with 215 reference genomes of this species, representing strains from its two main subspecies (lactis and cremoris) and isolated from different environments (dairy, vegetable, animal, etc.) (Figure 5A). In this analysis, five groups could be highlighted, one of which included all the representatives of the cremoris subspecies (brown branches), and three were formed by (i) industrial dairy strains (in red), (ii) diacetylactis biovariant strains (in blue), and (iii) dairy strains isolated from traditional products (in violet). The last group (green branches), which was the most diverse, contained environmental strains (plants, animals, etc.). The five L. lactis MAGs from Ambriss sand Serdaleh were included in this last group (highlighted by external red arrows).

A similar study was done with Lb. kefiranofaciens, Lt. kefiri Lt. diolivorans and Lp. plantarum. The Lb. kefiranofaciens FME151 genome and the eight Lb. kefiranofaciens MAGs assembled from Serdaleh (6) and Ambriss (3) metagenomes were compared to seven genomes present in the NCBI database, including those of the two type strains for the subspecies kefiranofaciens and kefirgranum. These strains were isolated from kefirs of different origins.

ANI analysis showed that all strains shared at least 98.1% identity and could be divided into two main groups, one made by kefir strains and the second by Ambriss and Serdaleh strains, where ANI were ~99.2-99.4% within each group and ~98.2-98.4% between groups. Further phylogenomic analyses (Figure 5B) showed that two groups could be differentiated within kefir strains, a first one with three strains sharing over 99.9% identity and including the Lb. kefiranofaciens subsp. kefiranofaciens type strain, and a second one with four strains that were more divergent between themselves, sharing 99.0 to 99.8% identity and including the Lb. kefiranofaciens subsp. kefirgranum type strain. Lastly, Ambriss and Serdaleh strains also formed two groups sharing 99.4 to 99.8% ANI within and 99.1 to 99.3% outside of their group. Analysis of gene composition confirmed and refined these results. First, it showed a very high relatedness of three Lb. kefiranofaciens subsp. kefiranofaciens strains from kefir that shared more than 95% of their genes (> 2,325 out of ~2,440 genes), while within each of the other groups, the strains were more divergent between themselves (87% of shared genes in the kefirgranum group, 88% in the Serdaleh group and 92% in the Ambriss group). Seventeen EPS-related genes involved in kefiran production present in a 14.4-kb region were present only in the three former strains, whereas only the first part of the operon containing the genes of enzymes involved in EPS regulation and priming was conserved in the other strains. Furthermore, in the latter strains, the genes involved in polymerization, chain length determination and export that could be characterized appeared to be different in each strain, indicating that they produce polysaccharides other than kefiran.

The Lt. kefiri FME71 genome and the ten Lt. kefiri MAGs assembled from Ambriss (3), Serdaleh (6) and Labneh El Darff (1) were compared to the ten reference genomes present in the NCBI database, including seven strains isolated from kefir, two from human gut and one from yoghurt. This analysis showed that all strains shared very high ANI values (> 99.55%), indicating that they were closely related (Table S10). Analysis of their phylogenomic tree did not indicate clustering as a function of the origin of the strains (Figure S4A). Furthermore, the analysis of the core/pan genome showed that all strains shared ~2,148 genes out of ~2400. The analysis of the specific genes in the different strains indicated that they code mostly for "hypothetical proteins" (50%) and mobile element proteins (25%). ANI analysis of Lb. diolivorans FME74 and of our 12 MAG strains with that of the type strain (isolated from maize silage) and one MAG from an ultra-filtered-milk-permeate fermenting bioreactor showed that these strains could be clustered into four groups (Figure S4B), displaying at least ~97.7 to 98.4% ANI . Serdaleh MAGs formed a group of very closely related strains (ANI > 99.9%), Ambriss and Labneh MAGs formed a second homogenous group (ANI > 99.1%), whereas the bioreactor MAG and the type strain were isolated.

The sequences of Lp. plantarum FME91 strain and of the five MAGs obtained from Serdaleh were analyzed with a set of 145 Lp. plantarum strains isolated from different ecosystems, including several types of food, animals, gut, etc. The five MAGs from Serdaleh appeared to be almost identical (ANI > 99.97%), but unrelated to that of FME91 from Ambriss (ANI ~92.7%), and they ranged between 95.7 and 99.7-99.9, with an average of 98.7 with the strains of the set (Table S10). An overview of the phylogenic relationship between these strains is given by the tree obtained by kSNP analysis (Figure S4C). Strains were dispersed irrespective of their environmental origin. The Ambriss strain FME91 was closely related to those isolated from milk or artisanal dairy products, whereas the Serdaleh strains were closely related to those isolated from meat, fruits and vegetables.

Detection of metabolic features by functional profiling

Microbial pathway analysis was carried out by HUMAnN 3.0 in order to highlight interesting features in the microbial communities such as bacteriocin production, clusters of genes potentially associated with health or aromatic pathways. Bacteriocin genes were searched in the 16 metagenomes, and 29 potential bacteriocin-encoding genes were identified (Table S11). The most frequent and abundant were Lactococcin 972 family bacteriocin and Helveticin, which are mainly from L. lactis and Lb. kefiranofaciens, respectively. Other frequent potential bacteriocin-encoding genes are Plantaricin and Thiazolpeptide-type bacteriocin, likely present in Lp plantarum and Lt diolivorans, respectively.

Potentially health-associated gene were divided into three broad categories, "survival", "colonization" and "modulation", as presented previously [START_REF] Leech | Fermented-Food Metagenomics Reveals Substrate-Associated Differences in Taxonomy and Health-Associated and Antibiotic Resistance Determinants[END_REF]. Briefly, these gene clusters were shown to be important for surviving in the gut conditions, colonizing the gut microbiome and stimulating immune system or microbiota products. While most samples displayed relatively similar profiles, several samples appeared to have different profiles (Figure S5). For example, Ambriss A6 presented the highest ratio of survival clusters, while Ambriss A4 presented the lowest ratio for all clusters. A search for overrepresented genes in Ambriss A6 highlighted the presence of genes from Lb. gasseri, a species that was found only in this sample.

The low scores in Ambriss A4 could be due to the dominance of S. parasuis and the low abundance of LAB in this sample.

Finally, genes involved in diacetyl metabolism were searched through citrate utilization, acetolactate and acetoin-butanol production. Since genomic references were known for dominant species of these microbiota, the corresponding genes could be assigned at the species level. It confirmed that the citrate-lyase cluster was absent from L. lactis strains but present in Lb. kefiranofaciens in our samples. Finally, the diacetyl degrading enzymes acetolactate decarboxylase (EC 4.1.1.5) and acetoin (diacetyl) reductase (EC 1.1.1.304) were found in L. lactis, Lt. kefiri and Lt. diolivorans, respectively.

Detection of pathogens by metagenomic analysis

Marker gene analysis showed that undesirable species were present in several Ambriss and Serdaleh samples. All Serdaleh samples collected in 2019 contained low levels (1-4%) of the goat pathogen M. agalactiae. Moreover, zoonotic streptococcal species, namely S. parauberis and S. parasuis, were detected at 0.5-2% in several samples of Ambriss and Serdaleh and, remarkably, at 7 and 55%, respectively, in Ambriss A4. S. aureus was detected at low levels in Serdaleh S10-b. Finally, Enterobacteria such as E. coli and K. pneumoniae were detected in about half of the samples at levels < 1% (except for K. pneumoniae present at 2.7% in Ambriss A1). Deep profiling analysis in the metagenome samples was performed to detect the presence at low abundance of pathogen species potentially found in dairy products (Table S7). No Brucella abortus, Listeria monocytogenes, Mycobacterium vaccae or Staphylococcus agalactiae were detected. However, S. aureus was detected at low levels in five Serdaleh samples.

Further analyses were performed on these undesirable species, which were sufficiently covered to yield MAGs. An ANI analysis of these MAGs with representative genomes of these species confirmed the taxonomic classification previously obtained by the ychF marker. For example, M. agalactiae, S. parauberis and K. pneumoniae MAGs share 99.8, 98.5 and 98.4% ANI with the type strains of their species, E. coli MAG; 99.1% with E. coli K12 laboratory strain, S. parasuis MAG; and 97-98% with four genome references of this species present in the NCBI database.

Gene content analysis was performed with S. parasuis MAG_A4, which appears to be dominant in Ambriss A4. In particular, genes involved in lactose metabolism were searched to explain the fitness of this strain in milk fermentation. Among the two systems used for lactose utilization, MAG-A4 contains the Lac-PTS system, including ten genes encoding the Lac-PTS transporter, the phospho-beta-galactosidase, the tagatose pathway, glucokinase, a repressor and a gene of unknown function (LacX). A blast search with these genes provided a hit at 97-99% nucleotide identity on its entire length with an ICESsuJH1308-1 mobile element present in S. suis JH1308 and a lower identity with several other S. suis strains. However, none of the nine S. parasuis genomes available contained these genes. These data show that S. parasuis MAG_A4 acquired lactose utilization genes by horizontal transfer, possibly from S. suis.

Interestingly, MAGs from M. agalactiae were of very good quality, probably due to particular features of their sequences, such as genomic guanine-cytosine content (~29.7%) and di-tri-nucleotide composition, allowing a better binning. In particular, MAGs from Serdaleh S8 and S9 consisted of only four and ten contigs, respectively. ANI analysis indicated that the strains from the Serdaleh samples were very closely related to each other (ANI > 99.9%), as well as to the type strain (ANI ~99.85%; Table S10). Phylogenomic analysis confirmed the close relationship of M. agalactiae Serdaleh MAGs together and with PG2 T (Figure 5C). Furthermore, gene content analysis indicates that these two MAGs differed from PG2 T by less than 20 genes encoding hypothetical proteins.

Detection of antibiotic resistance genes

A mapping-based approach against a comprehensive collection of antibiotic resistance genes was performed. A total of 18 ABR genes were identified, belonging to seven different classes of antibiotics (Figure 6). They were detected in almost all samples, and notably those from Serdaleh had higher levels and more diverse genes than those from Ambriss. Concerning the Labneh samples, only sample L18 presented detectable amounts of ABR. However, it cannot be concluded that L15, L16 and L17 samples did not contain ABR because there were not enough useful reads in these samples due to the presence of goat DNA.

Overall, the most abundantly detected antibiotic classes were aminoglycosides, MLS (Macrolide, Lincosamide, Streptogramin), nucleosides and tetracycline. In particular, Serdaleh samples contained high levels of aminoglycoside (aad(6) and APH(3'')-IIIa), MLS (ermB), nucleosides (sat-4) and tetracycline (tetS). The coverage of these genes is similar to that of L. lactis, which was the highly dominant species in these samples. Since this level of coverage allows a reliable level of assembly, we searched for contigs carrying these determinants by BLAST analysis in Serdaleh samples. This search pointed out two contigs of 12,360 and 3,370 nucleotides in sample S10-b, with coverages ~2 times that of L. lactis, and carrying, respectively, the four first determinants and tetS. The first contig is ~98% identical on 76-85% of its length to large plasmids from E. faecium, E. faecalis and Staphylococcus pseudointermedius. The second contig is identical on its full length to regions of plasmids and genome drafts of various food bacteria such as L. raffinolactis, L. lactis, S. thermophilus, Carnobacterium divergence and Listeria innocua. The other Serdaleh metagenomes contained contigs of comparable sizes indicating similar events.

The metagenome of Serdaleh S10-b also contained a ~6.6 kb contig carrying the tetL, tetM and cat genes that was identical on its entire length to regions of plasmids from E. faecalis.

The coverage of this contig is 40 times lower compared to the two former ones, which is a level corresponding to that of Lt. kefiri, Lb. kefiranofaciens and Lt. diolivorans. None of the other metagenomes contained these three genes together on one contig. Nevertheless, samples S8, S7b, S8-b, S9-b and S11-b contained contigs of ~5 and 0.9 kb carrying tetLM and cat genes, respectively. The coverages of these contigs were significantly different in the other assemblies where tetLM coverage was 2 to 10 times higher that of cat genes, indicating that their organization might be different from that found in S10-b.

Lastly, in Ambriss A4, tetM was found on a ~32 kb contig of S. parasuis, ermB on a 7 kb contig sharing 90% identity on 3 kb with plasmids of various Lactobacillaceae, cat on a ~3 kb contig sharing 93% identity to L. garviae plasmids, aad and sat on a ~23 kb contig sharing 93% identity on ~12 kb (including aad and sat) of the S. agalactiae B111 chromosome and 99% identity on 4-6 kb of E. faecalis plasmids and the S. pseudointermedius chromosome.

Discussion

Traditional Lebanese fermented goat milk products are perceived in Lebanon as having an important historical and cultural heritage, but their microbiota are still poorly characterized. A detailed characterization of 50 producers allowed us to select those who used genuine traditional processes for the production of Ambriss, Serdaleh and Labneh El Darff. The manufacturing processes were established and the physico-chemical characteristics of the products were determined. Our survey also made it possible to define the socio-economic profiles of the producers, which confirmed the danger of extinction of these products, thus making it possible to require a Protected Designation of Origin label. This study explored the microbial biodiversity of 20 34 samples of three fermented goat milk products from three regions by combining microbiological, genomic and metagenomic approaches.

The samples analyzed here shared similar technical properties for pH, salinity, temperature and microbial counts as those of the same products previously studied [START_REF] Dimassi | Production and physicochemical properties of labneh anbaris, a traditional fermented cheese like product, in Lebanon[END_REF][START_REF] Serhan | Concentrated yogurt (Labneh) made of a mixture of goats' and cows' milk: Physicochemical, microbiological and sensory analysis[END_REF]Nasbimana et al., 2005;[START_REF] Hajj | Characterisation Chimique Et Quaalité Bactériologique De Produits Laitiers[END_REF], and were therefore probably [START_REF] Hajj | Characterisation Chimique Et Quaalité Bactériologique De Produits Laitiers[END_REF] [START_REF] Hajj | Characterisation Chimique Et Quaalité Bactériologique De Produits Laitiers[END_REF] [START_REF] Hajj | Characterisation Chimique Et Quaalité Bactériologique De Produits Laitiers[END_REF] [START_REF] Hajj | Characterisation Chimique Et Quaalité Bactériologique De Produits Laitiers[END_REF] representative of these products in Lebanon.

Microbial originality of traditional Lebanese fermented goat milk products

An overall view of the microbiota is provided by metagenetic analysis, which deciphers species or lower taxonomic levels. This study clearly separated Labneh samples from those of Ambriss and Serdaleh. Labneh contained a relatively simple microbiota, strongly dominated by Lb. delbrueckii, which was almost absent from the other products. Furthermore, it contained Lt. diolivorans, Lt. kefiri and S. thermophilus to a lesser extent. The dominance of Lb. delbrueckii is due to the way Labneh is inoculated by Laban, whose composition is that of yoghurt made of Lb. delbrueckii and S. thermophilus (Chammas et al., 2006). The dominance of L. bulgaricus is probably due to its better resistance to acidity than S. thermophilus [START_REF] Béal | Production of Laban[END_REF].

Similar observations were made with yoghurt and Doogh (also known as Ayran), a traditional Iranian fermented ewe milk product [START_REF] Davati | 16S rRNA metagenomic analysis reveals significant changes of microbial compositions during fermentation from ewe milk to doogh with antimicrobial activity[END_REF]. The presence of Lt. diolivorans and Lt. kefiri, which are also found in Ambriss and Serdaleh, might be due to the specificity of the local goat milk. Moreover, the conditioning of Labneh in a goatskin may also influence the microbial composition by releasing nutrients or other metabolic products, as indicated by the very high level of goat DNA obtained from the metagenomic data for Labneh. Interestingly, Lt. diolivorans is usually found in vegetable and fruit fermentation, whereas it is occasionally reported in dairy products (Tamang et al. 2005;[START_REF] Zhang | Propionic acid production by cofermentation of Lactobacillus buchneri and Lactobacillus diolivorans in sourdough[END_REF][START_REF] Zheng | A taxonomic note on the genus Lactobacillus: Description of 23 novel genera, emended description of the genus Lactobacillus beijerinck 1901, and union of Lactobacillaceae and Leuconostocaceae[END_REF]. To the best of our knowledge, it was isolated from one raw milk in Brazil [START_REF] Agostini | Characterization of technological and probiotic properties of indigenous Lactobacillus spp. from south Brazil[END_REF], traditional yoghurts in Iran [START_REF] Roushanzadeh | Phenotypic and genotypic diversity of dominant lactic acid bacteria isolated from traditional yoghurts produced by tribes of Iran[END_REF], traditional fermented dairy products in Mongolia [START_REF] Yu | Diversity of lactic acid bacteria associated with traditional fermented dairy products in Mongolia[END_REF], and detected in metagenomic samples of Koumiss [START_REF] You | Metagenomic features of traditional fermented milk products[END_REF][START_REF] Kozhakhmetov | Metagenomic Analysis of Koumiss in Kazakhstan[END_REF]. The present study is thus the first one reporting the genome of Lt. diolivorans isolated from a dairy product.

Ambriss and Serdaleh samples contained, as dominant species, L. lactis and Lb. kefiranofaciens, which were almost absent from Labneh, followed by Lt. kefiri and, to a lesser extent, by Lt. diolivorans and Lp. plantarum. Interestingly, this composition was reminiscent of that of dairy kefir, a fermented milk with health properties, which is produced by the inoculation of a stable and complex microbial community made of lactic acid bacteria and yeasts entrapped in gelatinous grains [START_REF] Blasche | Metabolic cooperation and spatiotemporal niche partitioning in a kefir microbial community[END_REF][START_REF] Prado | Milk kefir: Composition, microbial cultures, biological activities, and related products[END_REF][START_REF] Nalbantoglu | Metagenomic analysis of the microbial community in kefir grains[END_REF]. In this product, the dominance of Lb. kefiranofaciens, which otherwise grows poorly alone in milk, is attributed to the metabolic cooperation it develops with the community, and the synthesis of an exopolysaccharide.

In this respect, growth conditions in Ambriss and Serdaleh, where the microbiota remained essentially in a coagulated milk matrix, displayed strong similarities with the one encountered in the kefir grain, which may thus explain the dominance of this particular Lactobacillus and of several sub-dominant lactic acid bacteria. MAGs and genome comparisons of Lb. kefiranofaciens and Lt. kefiri confirm the very close phylogenomic relationship between kefir, Ambriss and Serdaleh strains. Their gene content mainly differed by unknown functions and mobile elements. In Lb. kefiranofaciens strains isolated from Ambriss and Serdaleh and in Lb. kefirgranum subspecies strains, only the first six genes encoding the general functions of initiating and controlling EPS synthesis were conserved, whereas glycosyl transferase conferring the specific structure of EPS units differed. This fact explains why these strains did not produce kefiran, the specific polysaccharide for kefir grain formation. Lastly, metabolic cooperation in Ambriss and Serdaleh certainly differed from that of kefir, as suggested by the low level of Leuconostoc and Acetobacter species, which play an important role in milk fermentation and citrate metabolism [START_REF] Laëtitia | The Citrate Metabolism in Homo-and Heterofermentative LAB: A Selective Means of Becoming Dominant over Other Microorganisms in Complex Ecosystems[END_REF]. Nevertheless, the low input of milk on a the coagulated milk matrix gradually offers an environment for this "slow growing species" that is Lb. kefiranofaciens [START_REF] Blasche | Metabolic cooperation and spatiotemporal niche partitioning in a kefir microbial community[END_REF][START_REF] Georgalaki | Lactobacillus kefiranofaciens: From isolation and taxonomy to probiotic properties and applications[END_REF]. Since the microbiota of Ambriss and Serdaleh were not inoculated, the environmental origin of Lb. kefiranofaciens in these products still remains unknown. Interestingly, this species was also found to be dominant in several samples of Airag, a traditional Mongolian fermented mare milk made in horse-hide containers [START_REF] Watanabe | Diversity of lactic acid bacteria and yeasts in Airag and Tarag, traditional fermented milk products of Mongolia[END_REF] and in Koumiss (M. [START_REF] Zhang | Comparison of Bacterial Microbiota in Raw Mare's Milk and Koumiss Using PacBio Single Molecule Real-Time Sequencing Technology[END_REF], neither of which contain grain.

In the present set of samples, the balance between L. lactis and Lb. kefiranofaciens appeared to be the main factor differentiating Ambriss and Serdaleh, the former generally containing lower levels of L. lactis and higher levels of Lb. kefiranofaciens. Nevertheless, the samples of these two similarly processed products, originating from two distinct regions, could not be clearly separated since some samples of Serdaleh were also dominated by Lb. kefiranofaciens. Although Serdaleh samples were obtained at lower temperatures, displayed higher salt content, and reached lower pH in comparison to Ambriss samples, these factors do not fully explain this population shift. Alternatively, there may be subtle differences in production methods, such as the quantity of milk added at once, or a potential mixing that modifies the speed of milk infiltration through the product, which may affect the metabolism of the microbial communities in different ways, as discussed earlier. For example, if a significant quantity of the milk poured into the jar is not mixed with the curd, it may favor a rapid planktonic development of L. lactis, compared to the slow solid-state development where Lb. kefiranofaciens and Lt. kefiri compete more effectively.

A second important outcome of this study is the detection of uncultivated species as dominant microbiota in Ambriss and Serdaleh. While the initial failure to isolate Lb. kefiranofaciens from the samples can be explained by the fact that this species is tedious to grow and requires specific media and growth conditions [START_REF] Georgalaki | Lactobacillus kefiranofaciens: From isolation and taxonomy to probiotic properties and applications[END_REF], the lack of recovery of L. lactis is less understandable. Indeed, this bacterium is routinely cultured from different environments such as plant or dairy products [START_REF] Soundharrajan | Isolation of lactococcus lactis from whole crop rice and determining its probiotic and antimicrobial properties towards gastrointestinal associated bacteria[END_REF][START_REF] Li | Fermentation Characteristics of Lactococcus lactis subsp. lactis Isolated From Naturally Fermented Dairy Products and Screening of Potential Starter Isolates[END_REF]. In this work, we thus tested the hypothesis that L. lactis cells were viable but not-cultivable (VBNC) [START_REF] Oliver | Recent findings on the viable but nonculturable state in pathogenic bacteria[END_REF] and would thus require specific conditions to recover growth throughout their transfer from the product to M17 medium. Actually, a transfer of several hours in media containing goat milk allowed a straightforward isolation of strains of this species. This result indicates that the production process of Ambriss and Serdaleh induced a particular physiological state in L. lactis cells that hampers its growth recovery, thus corresponding to the definition of VBNC. Indeed, during the Ambriss and Serdaleh process, various environmental factors (such as pH, salinity and temperature) and the physiological age of the culture could induce the entry of bacterial cells into a VBNC state. Interestingly, the present process is still reversible, which is rarely the case [START_REF] Dong | Induction, detection, formation, and resuscitation of viable but non-culturable state microorganisms[END_REF], and hence offers a reproducible model for the study of the entry into the VBNC state. Further work would be required to understand why L. lactis growth recovery was deprived when directly plated on M17.

Potential safety concerns

A wide metagenomic search for bacteria potentially involved in public health issues indicated the presence of DNA of several kinds of undesirable species. Among hygiene indicators, Enterobacteria were detected at low levels in about half of the samples, as well as S. aureus, which was detected at low levels, mostly in Serdaleh samples. These results are comparable to those of a recent study on raw cow milk where, from a total of 195 samples collected from different Lebanese farms, 65% and 17% of the samples were non-compliant with national standards for Escherichia coli and S. aureus, respectively [START_REF] Joubrane | Microbiological contamination , antimicrobial residues , and antimicrobial resistance in raw bovine milk in Lebanon[END_REF].

Enterobacteriaceae are often reported as a typical bacterial population in raw milk from many mammalian origins [START_REF] Ramos | Antibiotic resistance profile and detection of degradative enzymes by enterobacteriaceae isolated from raw goat milk[END_REF]. They can also be transmitted to milk and dairy products by fecal contamination during the milking process or because of poor hygiene conditions [START_REF] Giraffa | Functionality of enterococci in dairy products[END_REF]. S. aureus was highly prevalent in some products and was the most frequent contaminant of milk and goat milk products [START_REF] Oliveira | Characterisation of Staphylococcus aureus strains from milk and goat cheese and evaluation of their inhibition by gallic acid, nisin and velame of the Brazilian caatinga[END_REF][START_REF] Xing | Prevalence and characterization of Staphylococcus aureus isolated from goat milk powder processing plants[END_REF].

The transmission of S. aureus to milk is due to its ability to form biofilms and to remain on surfaces and utensils, its presence being mainly attributed to poor handling during processing [START_REF] Kümmel | Staphylococcus aureus entrance into the Dairy Chain: Tracking S. aureus from dairy cow to cheese[END_REF].

A striking finding was the presence of several zoonotic streptococcal species such as S. parauberis and S. parasuis, which were detected at 0.5-2% in several samples of Ambriss and Serdaleh, respectively, and even at 7 and 55% in Ambriss A4. S. parasuis is an important pathogen for humans and pigs [START_REF] Nomoto | Reappraisal of the taxonomy of streptococcus suis serotypes 20, 22 and 26: Streptococcus parasuis sp. nov[END_REF] and was already isolated as a contaminant in raw camel milk. Intriguingly, in Ambriss A4, S. parasuis was strongly dominant over the traditional lactic acid bacteria, which accounted for only 10% of the population. The fact that the pH of this sample has a value similar to that of the other samples indicates that the strain of S. parasuis present in this sample is an acidifying strain. Indeed, our analysis revealed the presence of ten genes encoding a Lac-PTS system in its genome, similar to the one that allows efficient lactose utilization in several lactic acid bacteria such as L. lactis (Aleksandrzak-Piekarczyk, 2013). These genes were probably acquired by a horizontal gene transfer from S. suis.

Adaptation to lactose utilization in milk has already been observed for other non-starter streptococci, in particular, S. infantarius in fermented dairy products from East and West Africa [START_REF] Jans | Prevalence and comparison of Streptococcus infantarius subsp. infantarius and Streptococcus gallolyticus subsp. macedonicus in raw and fermented dairy products from East and West Africa[END_REF] and in artisanal cheese from Brazil [START_REF] Kothe | Unraveling the world of halophilic and halotolerant bacteria in cheese by combining cultural, genomic and metagenomic approaches[END_REF].

Serdaleh samples contained low levels (1-4%) of the goat pathogen M. agalactiae. This bacterium is one of the causal agents of contagious agalactia, a serious enzootic disease of small ruminants, which notably occurs in the Mediterranean Basin [START_REF] Kumar | Mycoplasma agalactiae, an etiological agent of contagious agalactia in Small ruminants: A review[END_REF]. Once it has infected the animal, it disseminates though the blood to different mucosa, including the mammary gland, and can consequently be found in milk. M. agalactiae has already been detected in artisanal cheeses via non-cultural techniques [START_REF] Martín-Platero | Polyphasic study of microbial communities of two Spanish farmhouse goats' milk cheeses from Sierra de Aracena[END_REF]. Its presence here is, therefore, a sign that the herd in the Chouf region was infected with this pathogen.

Finally, transmissible antibiotic resistance genes, which are considered to be a growing threat to public health [START_REF] Ghimpețeanu | Antibiotic Use in Livestock and Residues in Food-A Public Health Threat: A Review[END_REF][START_REF] Yadav | Antibiotic resistance: Global health crisis and metagenomics[END_REF], were detected. In particular, Serdaleh samples contained a multi-resistance cassette carrying aad(6) and APH(3'')-IIIa, ermB and sat-4, as well as tetracycline resistance genes. Since these antibiotics are used in prophylaxis for M. agalactiae, it is likely that such treatment contaminated the milk and favored the emergence of antibiotic resistant lactic acid bacteria strains.

Conclusion

In this study, we analyzed three artisanal Lebanese fermented products (Ambriss, Serdaleh and Labneh El Darff) that are processed from the raw milk of the Black Baladi goat, a local breed. The process consists of the daily introduction of salted milk or Laban into a terracotta jar or a goatskin over several months, together with the draining of the whey, until the container is full. This production remains highly artisanal and offers opportunities for the empowerment of women, but it is in danger of extinction due to the low volumes that are produced, the age of the producers and the absence of marketing. In industrial products that are inoculated with starters or by backslopping, the microbiota rapidly develops in large batches of milk, whereas in the current products, it develops slowly during percolation of the milk that is added periodically. During batch growth, fast-growing microorganisms are highly favored, whereas in a dense community, slow-growing bacteria can compete more effectively, as is the case with the well-described kefir grain where Lb. kefiranofaciens clearly predominates despite its poor growth performance in milk. Interestingly, kefir is promoted as a probiotic beverage and Lb. kefiranofaciens is presented as a key factor in this attribute [START_REF] Georgalaki | Lactobacillus kefiranofaciens: From isolation and taxonomy to probiotic properties and applications[END_REF].

Although no studies have yet proven this, traditional Lebanese goat milk products are believed to deliver health benefits, and it would thus be interesting to further test this possibility. These products are also interesting as sources of biodiversity because they contain strains of Lb. kefiranofaciens and Lt. kefiri, which are genetically quite different from those present in kefir that are highly conserved.

Furthermore, progress has to be made to ensure the safety of these products from herd prophylaxis to milk chain management via a global action in order to avoid contamination by potential pathogens and the spread of antibiotic resistance genes. Finally, this study also provides key elements related to the microbiology and microbiota functioning of these products that will facilitate the development of strategies to improve these processes and preserve their identity. Zúñiga, M., Pardo, I., & Ferrer, S. (1993). An improved medium for distinguishing between homofermentative and heterofermentative lactic acid bacteria. International Journal of Food Microbiology,18(1), 37-42. https://doi.org/10.1016/0168-1605(93)90005-2 9 The different groups of L. lactis are highlighted. The colors of the branches make it possible to visualize cremoris subspecies (brown), environmental sp. lactis strains (green), industrial dairy sp. lactis (red), diacetylactis biovariant (blue) and lactis sp from traditional dairy products (violet). The internal circle presents the environment from which the strains were isolated (blue: milk; green: vegetable; purple: animal). Ambriss and Serdaleh MAGs and the isolated strain genome are highlighted by external red arrows and names are shaded in gray Table S2: Physico-chemical parameters and culture counts of Ambriss, Serdaleh and Labneh El Darff samples as a function of their producers and production dates. Table S3. Percentage (%) of bacterial operational taxonomic units (OTUs) in the samples. Table S4. Percentage (%) of Capra hircus reads aligned in each fermented food metagenome. Table S5. Relative abundance of microbial composition of 16 traditional fermented goat milk metagenomes from Lebanon using the MetaPhlAn taxonomic assignment method. Table S6. Relative abundance of microbial composition of 16 traditional fermented goat milk metagenomes from Lebanon using ychF marker gene assignation. Table S7. Relative abundance of microbial composition of 16 traditional fermented goat milk metagenomes from Lebanon using the FoodMicrobiome tool.

Figure 1 .

 1 Figure 1. Geographical origin (A) and schematic processing diagram (B) of Ambriss, Serdaleh and Labneh El Darff. Ambriss and Serdaleh are made in terracotta jars periodically filled with raw goat's milk, and Labneh El Darf in goatskin filled with Laban, which is either a homemade yoghurt of raw goat milk. The detailed flow charts are presented in Figure S1 and S2.

Figure 2 .

 2 Figure 2. Physico-chemical and microbiological characteristics of Ambriss, Serdaleh and Labneh El Darff: Temperature, pH and salinity (A); boxplot representation of microbial counts on M17 and MRS (B); and histogram of isolated species (C).

Figure 3 .

 3 Figure 3. Metagenetic analysis based on OTU data of Ambriss, Serdaleh and Labneh El Darff: Boxplots of alpha-diversity indices (A); Principal coordinate analysis among bacterial communities identified by OTUs (B); bacterial composition plot depicting the relative abundances of the 20 main OTUs (C). Presence of different letters (a, b, c) above the box-plots indicate their statistical significance compared to the others Colors in (A) and (B) refer to the different fermented milk products, i.e., yellow for Ambriss, blue for Serdaleh and red for Labneh El Darff.

Figure 4 .

 4 Figure 4. Relative abundance of microbial composition of Ambriss, Serdaleh and Labneh El Darff. Values were calculated from the coverage of the ychF marker gene assembled in the metagenomes. The asterisk (*) indicates a low depth of the sequences for L15, L16 and L17 samples.

Figure 5 .

 5 Figure 5. Phylogenomic tree for L. lactis (A), Lb. kefiranofaciens (B) and M. agalactiae (C).

Figure 6 .Figure S1 :

 6S1 Figure 6. Heatmap showing distribution of the 18 antibiotic resistance genes detected in reads/kb within the 13 metagenomic samples. MLS: Macrolide, Lincosamide and Streptogramin.

Figure S2 :

 S2 Figure S2: Flow chart of the Labneh El Darff manufacturing process Figure S3. Rarefaction curves depicting the depth of 16S sequencing and species richness for the data obtained from Ambriss, Serdaleh and Labneh El Darff. The x-axis represents the sequencing depth (reads) and the y-axis represents the estimated OTU richness detected at the species level.

Figure S4 :

 S4 Figure S4: Phylogenomic tree forLb. kefirano, Lt. kefiri (A), Lt. diolivorans (B) and Lp. plantarum (C).

Figure S5 :

 S5 Figure S5: Heatmap showing the presence of potentially health-associated gene clusters across the Lebanese fermented goat milk product samples. Gene clusters are binned according to the potential ability of the microbiota to colonize the gastrointestinal tract, survive transit to the gut, and modulate the host phenotype. Each line is standardized across all samples, allowing only products to be compared with each other.Table S1: Producer survey results of production parameters, rearing practices and price indicators related to Ambriss, Serdaleh and Labneh El Darff.

Table S8 .

 S8 Quality of prokaryotic MAGs

Table S9 .

 S9 Genomes sequenced and assembly statistics

Table S10 .

 S10 FastANI values for MAGs, isolated strain genomes and NCBI reference genomes for L. lactis,Lb. kefiranofaciens, Lt. kefiri, Lt. diolivorans, L. plantarum, M. agalactiae and S.

	parasuis.

Table S11 .

 S11 Bacteriocin profile (in gene counts) of 16 traditional fermented goat milk metagenomes from Lebanon.

7 Acknowledgements This work was supported by the PHC CEDRE program 2018 (French-Lebanese campus), project number 40316UJ), the Lebanese University Grant Research program, and the French Dairy Interbranch Organization (CNIEL, Centre National Interprofessionnel de l'Economie Laitière), Paris, France. The authors are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility, doi: