Serra da Estrela PDO Cheese Microbiome as Revealed by Next Generation Sequencing

Abstract

Serra da Estrela PDO cheese is the oldest traditional cheese manufactured in Portugal. In this work, its microbiome as well as the main raw materials used in cheese production, raw ewes’ milk and thistle flowers (Cynara cardunculus L.), were characterized using next generation sequencing. Samples were accordingly retrieved from a local producer over two consecutive production campaigns and at different time periods within each campaign. The bacterial and fungi communities associated with each matrix were accessed through sequencing of V3–V4 and Internal Transcribed Spacer 2 regions of rRNA gene amplicons, respectively. A high microbial diversity was found associated to each matrix, differing significantly (p < 0.05) from each other. Over 500 taxa were identified in each analyzed matrix, ranging from dominant (relative abundance > 1%), sub-dominant (0.01–1%) and rare taxa (<0.01%). Specifically, in cheese, 30 taxa were present in all analyzed samples (core taxa), including species of Leuconostoc spp. and Lactococcus spp. for bacteria and Candida spp., Debaryomyces spp. and Yarrowia spp. for fungi, that were cumulatively the most prevalent genera in Serra da Estrela PDO cheese (average relative abundance ≥10%). Ultimately, this characterization study may contribute to a better understanding of the microbial dynamics of this traditional PDO product, namely the influence of raw materials on cheese microbiome, and could assist producers interested in preserving the identity, quality and safety of Serra da Estrela PDO cheese.

1. Introduction

Serra da Estrela cheese (SEC) is the most renowned Portuguese traditional cheese [1]. Produced in the foothills of the Serra da Estrela mountain, to which it owes its name, it follows artisanal protocols passed down between generations of shepherds, tracing back to the Roman occupation of the Iberian Peninsula [1,2,3]. Nowadays, SEC is the primary source of income of local shepherds and farmers, and is one of the most important economic activities in the region [1,2,3].

SEC manufacturing uses ovine raw milk from pasture feed females of “Bordaleira Serra da Estrela” and/or “Churra Mondegueira” breeds, coagulated with the addition of crude extracts of dried Cynara cardunculus L. flower and salt [4]. After a ripening period of 30–45 days, a flat cylinder, highly aromatic with a clean, smooth and slightly acidically flavored cheese is obtained. Featuring a thin, soft, uniform straw–yellowish rind and an ivory–white, buttery, creamy and unctuous paste with none too few little eyes [1,2,3,4]. These intrinsic characteristics earned it, in 1996, the Protected Designation of Origin (PDO) status.

Over the years, extensive fundamental and applied research has been carried out to unveil SEC unique organoleptic characteristics (reviewed in Inácio et al. [4] and Macedo et al. [2]). One critical factor is the use of raw untreated milk over the deliberate addition of starter cultures. This means that the native microbiota of the milk ends up playing a critical role in the manufacture of SEC [2,4,5]. The microbial consortia associated with SEC has been extensively studied, from the raw materials up to the final product [4]. However, these studies were limited to data gathered from cultured dependent methods [3,5,6,7,8], which are unable to truly unveil the diversity associated with these types of product [9,10,11]. Over the last decade, culture-independent methods, such as high-throughput sequencing technologies, have emerged, allowing an in-depth characterization and the identification of both high and low abundance microbiota at a species and, even in some cases, at a subspecies level [12,13,14,15,16]. Hence, Next Generation Sequencing (NGS) is been applied in the characterization of the microbial communities of traditional cheeses [9,17,18,19,20,21]. However, as far as the authors know, this is the first time that is used to characterize SEC microbiome. Therefore, the aim of the present work was to characterize the microbial community, of both bacteria and fungi, associated with the SEC and raw materials using NGS.

2. Materials and Methods

2.1. Cheese Manufacture and Sampling Strategy

For microbial characterization of SEC, raw ewes’ milk (REM), cardoon and cheese samples were analyzed from two consecutive production campaigns (1) 2018/2019; (2) 2019/2020. These were subdivided into three periods: (1) November–January; (2) February–March; and (3) May–June, as depicted in Figure 1.

The sampling strategy allowed a systematic analysis of SEC production, including autumn, winter and spring manufactures, reflecting distinct weather and pasture conditions that, on previous reports, posed a significant influence on microbial abundance [2,6,7]. Furthermore, to account for farmhouse variability REM from the two producers were used and analyzed throughout the study, A in 2018/2019 and B in 2019/2020 [2]. Moreover, four types of cardoon from Cynara cardunculus L. were used and analyzed. One was provided by the cheese producer and is regularly used in SEC manufacture, identified as “commercial”, and three genotypes of Cynara cardunculus L.,1 M, 3 M and 6 M [22].

Finally, SEC was manufactured in a certified PDO cheese producer from Gouveia, Portugal. REM and cardoon samples used for cheese manufacture were collected immediately before use. Serra da Estrela cheeses were produced following standard procedures and a ripening period of 38 days. All samples were transported and kept under refrigeration (≤4 °C) before processing.

2.2. Total DNA Extraction

Upon arrival, milk, cardoon and cheese samples were prepared for DNA extraction. For REM, three replicate tubes containing each 6 mL of milk were pelleted at 13,000 rpm (Hettich, Tuttlingen, Germany) for 10 min [23]. The supernatant was discarded and the pellet was resuspended in DNA/RNA Shield™ solution (ZymoResearch, Irvine, CA, USA). For cardoon, 2 to 10 g of dried thistle flower samples were suspended in 100 mL of a sterile solution containing 10% (v/v) Buffered Peptone Water (Liofilchem srl, Roseto degli Abruzzi, Italy) and 0.01% (v/v) Tween 80^® (PanReac AppliChem). Cardoon suspension was subjected to an ultrasonic bath for 30 s (Jet Program option) (Soltec, Milan, Italy) followed by 150 rpm orbital agitation for 30 min at room temperature (B. Braun Biotech International, Melsungen, Germany). Plant material was grossly filtered out using a stomacher bag filter and the remaining solution was divided into three 20 mL replicate aliquots. Cells were pelleted at 10,000 rpm for 5 min, the supernatant was discarded, and the pellet was resuspended in a DNA/RNA Shield™ solution. For cheese, three replicate portions of 400 mg were aseptically retrieved from the paste side and center of each cheese and were transferred to a tube containing a DNA/RNA Shield™ solution. All sample replicates of milk, cardoon and cheese in DNA/RNA Shield™ solution were homogenized in a bead beater (Benchmark Scientific, Sayreville, USA) for 6 min at maximum speed and stored at −20 °C until DNA extraction.

Total DNA was extracted from each sample replicate using ZymoBIOMICS™ DNA Miniprep Kit (ZymoResearch, Irvine, CA, USA), following the manufacturer’s instructions. The eluted DNAs were cleaned and concentrated using a DNA Clean and Concentrator™ Kit (ZymoResearch, Irvine, CA, USA), following the manufacturer’s instructions. Finally, the purified DNA from three independent extractions of each sample was pooled together and used in downstream procedures.

2.3. Sequencing Preparation, Run and Processing

The pooled DNA from each sample was subjected to two separated PCR amplification runs, one for the V3–V4 hypervariable region of the 16S rRNA gene and the other for the Internal Transcribed Spacer (ITS) region of the 18S rRNA. The PCR amplification runs were performed using KAPA HiFi HotStart PCR Kit (Roche, Indianapolis, IN, USA) according to manufacturer instructions, using 0.3 µM of each primer and 12.5 ng of template DNA, in a total volume of 25 µL. For the amplification of the V3–V4 region, the forward Bakt_341F 5′–CTACGGGNGGCWGCAG–3′ and reverse primers Bakt_805R 5′–GACTACHVGGGTATCTAATCC–3′ were used [24,25]. For the amplification of the ITS region, a pool of forward primers, ITS3NGS1_F 5′–CATCGATGAAGAACGCAG–3′, ITS3NGS2_F 5′–CAACGATGAAGAACGCAG–3′, ITS3NGS3_F 5′–CACCGATGAAGAACGCAG–3′, ITS3NGS4_F 5′–CATCGATGAAGAACGTAG–3′, ITS3NGS5_F 5′–CATCGATGAAGAACGTGG–3′, ITS3NGS10_F 5′–CATCGATGAAGAACGCTG–3′ and reverse primer ITS4NGS001_R 5′–TCCTSCGCTTATTGATATGC–3′ were used [26]. PCR conditions included a 3 min denaturation at 95 °C, followed by 25 cycles of 98 °C for 20 s, 55 °C (V3–V4 region)/60 °C (ITS region) for 30 s and 72 °C for 30 s and a final extension at 72 °C for 5 min. Then, a second PCR reaction added indexes and sequencing adapters to both ends of the amplified target region according to the manufacturer’s recommendations [27]. Negative PCR controls were included for all amplification procedures. PCR products were finally purified and normalized using SequalPrep 96-well plate kit (ThermoFisher Scientific, Waltham, MA, USA) [28], pooled and pair-end sequenced in the Illumina MiSeq^® sequencer with the V3 chemistry, according to the manufacturer’s instructions (Illumina, San Diego, CA, USA) at Genoinseq (Cantanhede, Portugal).

Raw reads were extracted from Illumina MiSeq^® System in FASTQ format and quality-filtered with PRINSEQ version 0.20.4 [29] to remove sequencing adapters, trim bases with an average quality lower than Q25 in a window of 5 bases and reads with less than 100 and 150 bases in length for ITS and V3–V4 files, respectively. Forward and reverse reads were merged by overlapping paired end reads with AdapterRemoval version 2.3.0 [30] using default parameters. ITSx version 1.1.2 [31] was used on ITS files to extract the highly variable fungal ITS2 subregion from merged reads and were then filtered to remove ITS2 reads below 100 bases in length. Sample IDs were assigned to the merged reads and were converted to a FASTA format. Chimeric merged reads were detected and removed using VSearch [32], an implementation of UCHIME [33] against Greengenes database version 13_8 [34] for V3–V4 and UNITE/QIIME ITS database version 8.2 [35] for ITS2 files. Operational Taxonomic Units (OTU) were generated in the Quantitative Insights into Microbial Ecology (QIIME) software version 2020.2 [36]. OTUs were selected at a similarity cut-off of 97% using the open reference strategy and those with less than two reads were discarded.

2.4. Bioinformatics and Statistical Analysis

All files generated in the previous section were analyzed in QIIME software version 2020.8 [36]. Diversity indexes namely, Shannon [37], Simpson [38] and Goods Coverage [39], as well as richness estimators, Chao1 [40] and Abundance Coverage Estimator (ACE) [41], were determined using q2–diversity plugin and alpha pipeline.

A phylogenetic tree was generated using q2–phylogeny plugin and align–to–tree–mafft–fasttree pipeline [42,43]. From there, OTU, Faith Phylogenetic diversity [44] and Shannon [37] rarefaction plots were constructed using q2–diversity plugin and alpha–rarefaction pipeline (Figures S1 and S2 of supplementary material). Furthermore, q2–diversity plugin and core–metrics pipeline was used to: (1) calculate α–diversity metrics, namely Shannon, observed features, Faith phylogenetic diversity and Evenness; (2) calculate β–diversity metrics, namely Jaccard distance [45], Bray–Curtis dissimilarity [46] and unweighted and weighted UniFrac [47,48]; (3) generate principal coordinates analysis (PCoA) plots. The sampling depth parameters used in this step were defined for each group of samples to be the highest possible, while retaining all samples (Tables S1 and S2 of supplementary material). However, due to a very low number of reads (8832), the 6M cardoon sample from the 2019/2020 campaign was intentionally dropped from the ITS diversity analysis between matrices. The generated files were further analyzed within q2–diversity plugin and alpha–group–significance and beta–group–significance pipelines to explore and disclose α and β–diversity differences, using Kruskal–Wallis [49] and Permutational multivariate analysis of variance tests [50], respectively, between matrices, campaigns, periods and cardoon genotypes. Differences with p < 0.05 were considered significant.

OTU taxonomy assignment was accomplished using q2–feature–classifier plugin [51]. A series of trial classification runs were made to evaluate the best taxonomy assignment in sequence files obtained from a mock community sample (ZymoResearch, Irvine, CA, USA) sequenced in parallel with the samples. To that end, classify–sklearn and classify–consensus–vsearch pipelines were tested using GreenGenes (version 13_8) [34], NCBI (BioProject 33175) (U.S. National Library of Medicine, Bethesda, USA), EZBioCloud (version of January_2021) (ChunLab, Seoul, Republic of Korea), ARB Silva (version SSURef_NR99, release 138) [52] and UNITE (version 8.2, dynamic fungi release) [35] databases. The best taxonomic outcome for V3–V4 mock file was obtained using the classify–consensus–vsearch pipeline coupled with ARB Silva database and a confidence threshold for taxonomic assignment ≥90%. For mock ITS files, the best performance was obtained using classify–sklearn pipeline coupled with a UNITE database and a confidence threshold for taxonomic assignment ≥70%. These taxonomic annotations were used throughout the work for the construction of stacked bar charts, identification of core taxa by matrix and to identify the taxa responsible for α and β–diversity differences. A Linear Discriminant Analysis (LDA) of Effect Size (LEfSe) was used to estimate the effect size on differentially abundant taxa between and within matrices [53]. LEfSe was performed on a Galaxy computational tool (http://huttenhower.sph.harvard.edu/galaxy/, accessed on 7th of April 2021) using a 0.05 alpha value for the factorial Kruskal–Wallis test among classes, a 0.05 alpha value for the pairwise Wilcoxon test between subclasses, a “one against all” strategy for multi-class analysis, a pairwise comparison only among subclasses with the same name and a logarithmic LDA score threshold of 3.5 and 6.0 for bacteria and fungi, respectively. Observed abundances, at any given taxonomic level, of the total reads were then used to define the dominant (>1%), sub-dominant (between 0.01 and 1%) and rare taxa (<0.01%) [16].

2.5. Nucleotide Sequences Accession Number

Raw reads were deposited in SRA database under BioProject PRJNA723623.

3. Results and Discussion

3.1. Microbial Community of Raw Ewes’ Milk

Milk present in the upper part of a healthy lactation female is generally considered sterile [54,55]. Milk colonization occurs from direct and/or indirect transfer from reservoirs. Direct reservoirs are locations and/or surfaces that at some point contact directly with the previous sterile raw milk, such as the teat canal, teat surface and milking equipment. Indirect reservoirs can be, for example, feed material, grass, soil, air, litter, water (drinking and/or washing), stable, milker and milk parlour [15,54,55].

In this study, over 500 genera (213 of bacteria and 293 of fungi) were identified in the microbial community of REM (Table S3 of supplementary material). This composition is in general, analogous to previous reports on the microbiology of REM from Serra da Estrela region [5], from other locations [56,57,58,59,60,61] and raw milk in general [15,54,55,62].

Despite a wide taxa variability between milk samples, α and β group significance analysis showed that none of the possible factors, such lactation period/season and producer/campaign, posed a significant impact on sample diversity (p > 0.05). Except for the number of fungi OTUs, which were significantly higher in milk retrieved from producer A (2018/2019 campaigns) in comparison to the one from producer B (2019/2020 campaigns) (p < 0.05).

The genus Candida, frequently related to yeast udder inflammation [63,64], with 56% of the relative frequency on average, and C. zeylanoides with 40% were the predominant fungi taxa in REM (Figure 2 and Figure S3 of supplementary material), despite the procedures for the production of certified SEC banning the use of ovine milk obtained from unhealthy females [65]. This high prevalence, without any other information, could be indicative of latent infection, which corresponds to 95–98% of all mastitis cases [63]. Hence, a constant monitorization of the health condition of lactating females is advised to prevent economic losses throughout the production chain.

Psychrotolerant microflora accounted for ≈40% of total bacterial counts (Figure 2 and Figure S3 of supplementary material). This is most likely the result of a combination of milk refrigeration, as recommended in the manufacture manual of SEC, for transportation and storage purposes and milk handling in cold weather, typical from the region of Serra da Estrela [65]. These conditions are known to alter the balance of milk microflora, from a gram-positive to a gram-negative dominance [15,54,55,66]. Finally, Lactic Acid Bacteria (LAB) accounts for a quarter of bacterial counts, with Lactococcus lactis, Enterococcus faecium and Leuconostoc mesenteroides with 12, 4 and 3% of the relative frequency on average, respectively, as the most prevalent species (Figure 2 and Figure S3 of supplementary material).

3.2. Microbial Community of Cardoon

The dried thistle flower of Cynara cardunculus L., cardoon is frequently used as a clotting agent in the production of goat and ovine milk cheeses in Portugal, Spain and Italy [67]. The addition of cardoon in cheese production is not standardized, but is still an empirical exercise of the cheese maker, ranging from 0.15 to 0.6 g of flower in L⁻¹ of milk [4,67]. Furthermore, several methods are reported for the addition of cardoon in cheese making [2]. One of the most popular, and the one used in the selected cheese factory, consists of the maceration of cardoon flowers in water, filtering this extract with the help of a fine clean cloth, which is poured directly into the warm milk. Extensive research was carried out to study the influence of cardoon in cheese production, from a chemical, biochemical, microbiological and sensorial point of view (reviewed in Conceição et al. [67]). However, only a few reports accessed the microbial composition of cardoon extracts [68,69,70].

Over 500 genera (265 of bacteria and 263 of fungi) were identified in cardoon samples, comprised mostly by taxa with a ubiquitous distribution within the environment and frequently isolated from vegetal material such Aspergillus, Aureobasidium, Cladosporium, Lactobacillus, Mucor, Pantoea, Pseudomonas and Rhizopus (Figure 3 and Figure S4 and Table S4 of supplementary material) [55,71,72,73,74,75].

Despite a wide taxa variability between cardoon samples, α and β group significance analysis showed that none of the possible factors, such cardoon genotype and/or production campaign, posed a significant impact on sample diversity (p > 0.05). Except for the number of bacterial OTUs, which were significantly higher in cardoon samples from 2019/2020 campaign in comparison to the ones from 2018/2019 (p < 0.05).

3.3. Microbial Community of Serra da Estrela Cheese

The microbial community of cheese is composed of variable ranges of bacteria, yeasts and molds [54,62,76,77]. Their access to cheese occurs either by deliberate addition, transference from cheese ingredients (e.g., milk, rennet, salt, among others) and/or transference from direct contact with equipment and the cheese maker [76,78].

The microbiological characterization of SEC samples revealed the presence of almost 500 genera (138 of bacteria and 338 of fungi—Table S5 of supplementary material). A composition that is in general, analogous to previous SEC characterizations [3,5,6,7,8,79], other ovine milk cheeses [18,21,79,80,81,82,83,84], and cheese in general [54,62,77].

As expected, the bacterial community is mainly composed of LAB, ranging from 50 up to 90% of all counts (Figure 4 and Figure S5 of supplementary material). Of those, Leuconostoc mesenteroides and Lactococcus lactis stand out as the major constituents. This is in contrast with previous characterizations, that have shown a limited contribution of these species, in particular to the overall microbial composition of SEC [5,6]. However, recent publications concerning other ovine milk cheeses show a high prevalence of L. mesenteroides and L. lactis, namely Azeitão from Portugal [80], Krčki, Istrian and Paski from Croatia [18], Feta from Greece [84] and Fiore Sardo from Italy [82].

The fungi community was dominated by the presence of Candida sake, Candida zeylanoides and Debaryomyces hansenii yeasts (Figure 5 and Figure S6 of supplementary material). D. hansenii is, in fact, one of the most prevalent yeasts found in cheese, previously described as a major constituent of the SEC yeast community [5,6] and other ovine milk cheeses [21,62,77,83,85]. On the contrary, the appearance and role of Candida spp. in cheese is a controversial topic. Generally regarded as a cheese contaminant, some species are known to be opportunistic pathogens (e.g., C. parapsilosis, C. tropicalis, C. albicans, among others), while others are nowadays used as adjunct cultures (e.g., C. krusei and C. colliculosa) [62]. Candida spp. have been previously identified as part of SEC microbiota, namely C. rugosa, C. zeylanoides and C. etchellsii [5]. Moreover, recent microbial characterization studies have highlighted the widespread presence of Candida species in cheese, which has led some authors to raise a yet-undiscovered role during cheese ripening stage [21,62,81,86,87]. In this case, its high prevalence results, most likely, from the high prevalence observed in ovine milk samples, as previously detailed.

Concerning food safety, the analyzed SEC samples were mostly colonized by beneficial and probiotic microbiota. Although, foodborne pathogens, responsible for occasional foodborne disease outbreaks derived from the consumption of contaminated cheese products, have also been detected (reviewed in Fox et al. [88]). Pathogenic bacteria such as Clostridium perfringens, Listeria monocytogenes, Yersinia enterocolitica, Salmonella enterica and Staphylococcus aureus were found, but at rare (<0.01%) or near the rare range (<0.09%) (Table S5 of supplementary material), attesting the overall safety of SEC.

Interestingly, an analysis of the abundance histograms points towards a higher diversity of 2019/2020 vs. 2018/2019 cheeses, while some taxa appear to be season- and campaign-specific (Figure 4 and Figure 5). In fact, α and β group significance analysis showed that the production campaign (2018/2019 vs. 2019/2020), as well as the manufacture period (November–January vs. February–March vs. May–June), poses a significant impact on SEC sample diversity (p < 0.05). Similarly, previous studies have shown the influence of the dairy farm, production campaign, season and axial location factors on overall microbial composition of SEC [2,3,5,6,7]. An LEfSe analysis on the data revealed 31 divergent taxa between campaigns (Figure 6). Of those, most of the divergent taxa were more associated to autumn manufactures, except for Kocuria spp., Lactococcus spp., Carnobacterium spp., C. curvatus and K. servazii, which were more abundant in winter and L. lactis, Lactobacillus spp., Acinetobacter spp., Citrobacter spp., S. stellimalicola and C. metapsilosis in spring cheeses (Figures S8 and S9 of supplementary material). Overall, these differences are likely to reflect a higher/lower presence of each specific taxon as a resident microbiota in the REM used for each cheese manufacture (Table S3 of supplementary material).

3.4. Matrices Taxa Variability and Core Community

The microbial composition and diversity associated with cheese, curd, whey, milk and the other used ingredients in cheese manufacturing, differ significantly between each other, reflecting each matrix type intrinsic characteristics and/or interactions [14]. Unsurprisingly, α and β group significance analysis showed that overall, REM, cardoon and SEC samples differ significantly in their microbial composition (p < 0.05). α group significance analysis showed that cheese samples present more homogenous microbial communities, while cardoon samples are the most diverse matrix. Finally, milk samples presented a high bacterial and low fungi diversity, similar to those observed in cardoon and cheese samples, respectively. Moreover, β group significance analysis showed that is possible to discriminate between each matrix, based on their specific bacterial community (p < 0.001), whereas the fungi community can only be used to distinguish cardoon from the grouped milk and cheese samples (p < 0.001) (Figure S7 of supplementary material). Overall, through an LEfSe analysis, it was found that psychrophiles, LAB and yeasts and plant-associated microorganisms are the most discriminatory biomarkers for REM, SEC and cardoon, respectively (Figure 7 and Figure 8).

The abundance data allowed the determination of the core microbial community associated with each matrix (Figure 9 and Figure 10), that ranged from dominant, sub-dominant to rare taxa (Tables S3–S5 of supplementary material). As demonstrated by α diversity metrics, cardoon is by far the matrix with the most diverse core microbial community with 61 taxa, followed by REM with 47 and SEC with 30. As expected, SEC shared a higher number of core taxa with REM than with cardoon—18 vs. 11—and only 9 were shared between all matrices. L. lactis, one of the most prevalent species in SEC, is one of the most important LAB species in the dairy industry [80]. Its presence in SEC derives predominantly from REM (accounting up to 12% of total counts) and rarely from cardoon (≤0.5% of total counts). It belongs to a restricted group of microorganisms, denominated starter bacteria, whose main function is the production of lactic acid from lactose metabolization [89]. Furthermore, it contributes to the development of the organoleptic characteristics of cheese during ripening, with the production of flavor compounds (e.g., free amino acids and medium- and small-peptides) from the metabolization of milk proteins and to texture through the secretion of exopolysaccharides [90,91].

Likewise, L. mesenteroides composes a significant portion of the SEC microbial community. In contrast to the low abundance levels observed in the raw materials of origin, REM and cardoon, respectively with 4 and 1% of total counts on average, resulting from a poor growth in milk [55]. To stimulate growth, leuconostocs requires the addition and/or synthesis of peptides and amino acids by other microorganisms and, as result of that, are frequently found in close symbiosis with lactococci in the microbiota of fresh and semi-hard cheeses [55,92]. Active during ripening, L. mesenteroides are involved in the development of cheese organoleptic properties through the production of aromatic compounds from co-metabolization of lactose and citrate [55,78].

Enterococci, namely E. faecalis and E. faecium, are another group of core microbiota in SEC. However, they are present at relatively low abundances when compared to other LABs (overall ≤5% of the relative frequency on average). These are a controversial group within LAB, that, depending on the strain, can be considered as a starter culture, probiotic, spoilage, or pathogenic bacteria [93]. They can grow under different conditions and substrates and survive refrigeration, high-temperature and high-salinity environments. However, SEC seems to provide an unfavorable growth environment, either by microbial competition and/or substrate deficiency, since its load on cheese is similar to the ones observed in the raw materials of origin, REM and cardoon (≈4 vs. ≤1%, respectively). Furthermore, the dominance of E. faecium over other enterococci species is maintained across matrices. Nonetheless, enterococci are known to contribute to the development of the organoleptic characteristics of cheese through their proteolytic, lipolytic, and citrate breakdown activities [55,77,78].

L. curvatus and L. plantarum compose the lactobacilli core microbiota of SEC, present at a sub-dominant level (0.01 ≤ relative abundance ≤ 1%), with the exception of L. curvatus on 2019/2020 cheeses that were found with a relative frequency of ≈5% on average. These lactobacilli are only occasionally found in REM and cardoon, which supports an external source of these taxa. They are active during the cheese ripening stage and are involved in a series of events, namely residual sugar metabolization, peptide hydrolysis, amino acid conversion, production of flavor and aroma compounds and the production of antifungal and antibacterial compounds [77].

Salt tolerant bacteria such Hafnia spp., Vibrio spp. and micrococci were found to range between sub-dominant and dominant levels among SEC samples. These taxa are known to survive the brining process and actively contribute to the development of the organoleptic characteristics of cheese [15,78]. Micrococci, and to a lower extent Hafnia spp., were found to inhabit REM and cardoon microflora, while Vibrio spp. was practically absent from these raw materials. An observation that points towards an equipment/facility-related input, for example during cheese washing and salting steps [13].

Psychrobacter spp. was found at a near rare range within SEC, and only occasionally in REM and cardoon samples. Capable of producing aldehydes, ketones and sulphur-volatile compounds, it contributes to the development of cheese organoleptic characteristics [55]. Overall, Enterobacteriaceae composes of a noteworthy group of REM, cardoon and SEC microflora. Despite their high aromatic potential, the presence of Serratia spp., Enterobacter spp. and Citrobacter spp. is often linked to detrimental effects, such hydrogen gas production and off-flavor development [15,54].

Finally, the fungi core taxa of SEC are mostly composed by yeasts (Figure 10). These are important during cheese ripening, namely in matrix deacidification through lactate consumption and in the development of cheese flavor and aroma. Yeast colonization of cheese occurs predominantly in the rind or in its vicinity, apart from some fermenting species that can grow in the cheese interior [13,62,77]. D. hansenii is a natural inhabitant of cheese, owing to its ability to metabolize lactate, lactose and tolerance to acidic and hypersaline habitats [62,77]. As previously stated, it is the third most prevalent yeast in SEC, whose strains arise predominantly from REM (≈6% of total counts) and to a smaller extent from cardoon (≤1%). For Y. porcina, K. lactis and S. cerevisiae, however, REM acts most likely as the primary source of strains for SEC, as it is virtually absent in cardoon samples. K. lactis plays an important role at the beginning of cheese ripening, namely in lactose metabolization, ethanol production and in the secretion of lytic enzymes, hereby contributing to cheese flavor and aroma [15,55,77]. S. cerevisiae, occasionally found in cheese, is active during the ripening stage, when it metabolizes lipids and proteins, contributing to cheese flavor and aroma [77]. On the contrary, Y. porcina is not commonly isolated from cheese, unlike Y. lipolytica [62,77,94]. Y. lipolytica presents a strong lipolytic and proteolytic activity, thus contributing to the organoleptic characteristics of cheese. Accordingly, it is possible that Y. porcina, similarly to the role of Y. lipolytica in other cheeses, is involved in flavor and aroma development in SEC.

4. Conclusions

To the extent of our knowledge, this study represents the first application of next generation sequencing technology to the characterization of the microbial community associated with the Serra da Estrela PDO cheese. This comprehensive work provided an insight into raw materials, raw ewes’ milk and cardoon, used in cheesemaking and cheeses analyzed over two consecutive production campaigns and three periods within each campaign. The here-disclosed microbial profiles, not just of the cheese but also of the milk and thistle flowers used in cheese manufacture, display a high microbial variability and have offered some insights regarding the provenience and abundance patterns of the microbiota found in cheese. Ovine milk was found as the most probable strain source for the core taxa revealed in SEC, with limited contributions of cardoon. Lactobacilli and Vibrio spp. were virtually absent from the raw materials analyzed, suggesting the contribution of other sources to the overall SEC microbiome. Moreover, the presence a newly described yeast, Y. porcina, not commonly associated with cheese, is revealed. This observation requires further investigation, namely in the evaluation of its impact on cheese microbiota dynamics and its contribution to cheese organoleptic characteristics.

Overall, this study contributes to a better understanding of a complex, valuable and appreciated product that must meet rigorous criteria for authenticity and PDO certification. Nonetheless, future research could focus on the expansion of the scope of this characterization, assessing microbial axial location, dairy-farm to dairy-farm variability, and microbial dynamics during cheese production and ripening stages, potentially addressing some of the unanswered questions raised in this work.

Author Contributions

Conceptualization, R.R., M.V.V., J.S. and P.F.; methodology, R.R., J.S. and P.F.; software, R.R.; validation, M.V.V., J.S. and P.F.; formal analysis, R.R.; investigation, R.R.; resources, M.V.V., J.S. and P.F.; data curation, R.R.; writing—original draft preparation, R.R.; writing—review and editing, R.R., M.V.V., J.S. and P.F.; visualization, R.R.; supervision, M.V.V., J.S. and P.F.; project administration, R.R. and M.V.V.; funding acquisition, M.V.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by: “MobFood—Mobilizing scientific and technological knowledge in response to the challenges of the agri–food market” (POCI–01–0247–FEDER–024524), by “MobFood” Consortium, and financed by European Regional Development Fund (ERDF), through the Incentive System to Research and Technological development, within the Portugal2020 Competitiveness and Internationalization Operational Program. UIDB/05937/2020 and UIDP/05937/2020—Centre for Research and Development in Agrifood Systems and Sustainability—funded by national funds, through FCT—Fundação para a Ciência e a Tecnologia.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Raw reads were deposited in SRA database under BioProject PRJNA723623.

Acknowledgments

The authors are thankful to the partner cheese factory for providing cheese, milk and cardoon samples as well as Professor Paulo Barracosa from Instituto Politécnico de Viseu—Portugal for providing thistle flower ecotypes 1 M, 3 M and 6 M.

Conflicts of Interest

The authors declare no conflict of interest.

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

Figures

Figure 1. Schematic representation of the sampling strategy used for the microbial characterization of ewes’ raw milk, cardoon and Serra da Estrela PDO cheese. Co refers to the commercially used cardoon by the cheese producer and 1 M, 3 M and 6 M refers to specific Cynara cardunculus L. genotypes used in this study in the manufacture of Serra da Estrela cheese.

Figure 2. Relative abundances (%) of dominant sequences (>1%) assigned to genus level identified in raw ewes’ milk DNA samples based on partial sequence analysis of the (a) V3–V4 and (b) Internal Transcribed Spacer 2 regions of the rRNA gene. WTI refers to percentage of sequences without taxonomic attribution to the specified taxonomic level.

Figure 3. Relative abundances (%) of dominant sequences (>1%) assigned to genus level identified in dried flowers of Cynara cardunculus L. (cardoon) DNA samples based on partial sequence analysis of the (a) V3–V4 and (b) Internal Transcribed Spacer 2 regions of the rRNA gene. Co refers to the commercially used cardoon by the cheese producer and 3 M and 6 M refers to specific cardoon genotypes. WTI refers to the percentage of sequences without taxonomic attribution to the specified taxonomic level.

Figure 4. Relative abundances (%) of dominant sequences (>1%) assigned to genus level identified in Serra da Estrela cheese DNA samples based on partial sequence analysis of the V3–V4 region in the rRNA gene. Co refers to the commercially used cardoon by the cheese producer and 1 M, 3 M and 6 M refers to specific cardoon genotypes. WTI refers to the percentage of sequences without taxonomic attribution to the specified taxonomic level.

Figure 5. Relative abundances (%) of dominant sequences (>1%) assigned to genus level identified in Serra da Estrela cheese DNA samples based on partial sequence analysis of the Internal Transcribed Spacer 2 region in the rRNA gene. Co refers to the commercially used cardoon by the cheese producer and 1 M, 3 M and 6 M refers to specific cardoon genotypes. WTI refers to the percentage of sequences without taxonomic attribution to the specified taxonomic level.

Figure 6. Linear Discriminant Analysis (LDA) scores computed for: (a) bacterial and (b) fungi features differentially abundant in Serra da Estrela cheese from: (a) 2018/2019 (red); (b) 2019/2020 (green) campaigns.

Figure 7. Linear Discriminant Analysis (LDA) scores computed for bacterial features differentially abundant in: raw ewes’ milk (blue bars); Serra da Estrela cheese (green bars); and dried flowers of Cynara cardunculus L. (cardoon (red bars)).

Figure 8. Linear Discriminant Analysis (LDA) scores computed for fungi features differentially abundant in: raw ewes’ milk (blue bars); Serra da Estrela cheese (green bars); and dried flowers of Cynara cardunculus L. (cardoon (red bars)).

Figure 9. Venn diagram representing shared core taxa of bacteria between raw ewes’ milk (blue area), cardoon (red area) and Serra da Estrela cheese (green area) matrices. f denotes identified bacterial families without taxonomic attribution to genus and/or species.

Figure 10. Venn diagram representing shared core taxa of fungi between raw ewes’ milk (blue area), cardoon (red area) and Serra da Estrela cheese (green area) matrices. f denotes identified bacterial families without taxonomic attribution to genus and/or species.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/microorganisms9102007/s1, Figure S1: Rarefaction curves obtained from V3−V4 region sequencing, Figure S2: Rarefaction curves obtained from Internal Transcribed Spacer 2 region sequencing, Figure S3: Relative abundance histogram assigned to species level identified in raw ewes’ milk, Figure S4: Relative abundance histogram assigned to species level identified in dried flowers of Cynara cardunculus L., Figure S5: Relative abundance histogram of bacteria assigned to species level identified in Serra da Estrela cheese, Figure S6: Relative abundance histogram of fungi assigned to species level identified in Serra da Estrela cheese, Figure S7: Principal coordinates analysis plots based on a Bray–Curtis diversity metric, Figure S8: Relative abundance histograms of biomarker features found for 2018/2019 Serra da Estrela cheeses, Figure S9: Relative abundance histograms of biomarker features found for 2019/2020 Serra da Estrela cheeses, Table S1: Internal Transcribed Spacer 2 sequencing summary data, Table S2: V3−V4 sequencing summary data, Table S3: Raw ewes’ milk frequency tables, Table S4: Cardoon frequency tables, Table S5: Serra da Estrela cheese frequency tables.

References

1. Morales, M.; Bintsis, T.; Alichanidis, E.; Herian, K.; Jelen, P.; Hynes, E.R.; Perotti, M.C.; Bergamini, C.V.; Attard, E.; Grupetta, A. et al. Soft cheeses (with rennet). Global Cheesemaking Technology; Papademas, P.; Bintsis, T. Jonh Wiley & Sons, Ltd: Hoboken, NJ, USA, 2017; pp. 301-325. [DOI: https://dx.doi.org/10.1002/9781119046165.ch4]

2. Macedo, A.C.; Malcata, F.X.; Oliveira, J.C. The technology, Chemistry and Microbiology of Serra Cheese: A review. J. Dairy Sci.; 1993; 76, pp. 1725-1739. [DOI: https://dx.doi.org/10.3168/jds.S0022-0302(93)77505-0]

3. Tavaria, F.K.; Malcata, F.X. On the microbiology of Serra da Estrela cheese: Geographical and chronological considerations. Food Microbiol.; 2000; 17, pp. 293-304. [DOI: https://dx.doi.org/10.1006/fmic.1999.0315]

4. Inácio, R.S.; Gomes, A.; Saraiva, J.A. Serra da Estrela cheese: A review. J. Food Process. Preserv.; 2000; 44, e14412. [DOI: https://dx.doi.org/10.1111/jfpp.14412]

5. Tavaria, F.K.; Malcata, F.X. Microbiological Characterization of Serra da Estrela Cheese throughout Its Appellation d’Origine Protégée Region. J. Food Prot.; 1998; 61, pp. 601-607. [DOI: https://dx.doi.org/10.4315/0362-028X-61.5.601]

6. Macedo, A.C.; Malcata, F.X.; Hogg, T.A. Microbiological profile in Serra ewes’ cheese during ripening. J. Appl. Bacteriol.; 1995; 79, pp. 1-11. [DOI: https://dx.doi.org/10.1111/j.1365-2672.1995.tb03117.x]

7. Macedo, A.C.; Costa, M.L.; Malcata, F.X. Changes in the Microflora of Serra Cheese Evolution Throughout Ripening Time, Lactation Period and Axial Location. Int. Dairy J.; 1996; 6, pp. 79-94. [DOI: https://dx.doi.org/10.1016/0958-6946(94)00051-4]

8. Dahl, S.; Tavaria, F.K.; Malcata, F.X. Relationships between flavor and microbiological profiles in Serra da Estrela cheese throughout ripening. Int. Dairy J.; 2000; 10, pp. 255-262. [DOI: https://dx.doi.org/10.1016/S0958-6946(00)00042-X]

9. Dalmasso, A.; Rio, M.; Civera, T.; Pattono, D.; Cardazzo, B.; Bottero, M.T. Characterization of microbiota in Plaisentif cheese by high–throughput sequencing. LWT Food Sci. Technol.; 2016; 69, pp. 490-496. [DOI: https://dx.doi.org/10.1016/j.lwt.2016.02.004]

10. Park, W.; Yoo, J.; Oh, S.; Ham, J.; Jeong, S.; Kim, Y. Microbiological Characteristics of Gouda Cheese Manufactured with Pasteurized and Raw Milk during Ripening Using Next Generation Sequencing. Food Sci. Anim. Resour.; 2019; 39, pp. 585-600. [DOI: https://dx.doi.org/10.5851/kosfa.2019.e49]

11. Randazzo, C.L.; Pitino, I.; Ribbera, A.; Caggia, C. Pecorino Crotonese cheese: Study of bacterial population and flavour compounds. Food Microbiol.; 2010; 27, pp. 363-374. [DOI: https://dx.doi.org/10.1016/j.fm.2009.11.010] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20227601]

12. Filippis, F.; Parente, E.; Ercolini, D. Metagenomics insights into food fermentations. Microb. Biotechnol.; 2017; 10, pp. 91-102. [DOI: https://dx.doi.org/10.1111/1751-7915.12421]

13. Irlinger, F.; Layec, S.; Hélinck, S.; Dugat-Bony, E. Cheese rind microbial communities: Diversity, composition and origin. FEMS Microbiol. Lett.; 2015; 362, pp. 1-11. [DOI: https://dx.doi.org/10.1093/femsle/fnu015] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25670699]

14. Jonnala, B.R.; McSweeney, P.; Sheehan, J.J.; Cotter, P.D. Sequencing of the Cheese Microbiome and Its Relevance to Industry. Front. Microbiol.; 2018; 9, 1020. [DOI: https://dx.doi.org/10.3389/fmicb.2018.01020] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29875744]

15. O’Sullivan, O.; Cotter, P.D. Microbiota of Raw Milk and Raw Milk Cheeses. Cheese—Chemistry, Physics and Microbiology; 4th ed. McSweeney, P.; Fox, P.F.; Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 301-316. [DOI: https://dx.doi.org/10.1016/B978-0-12-417012-4.00012-0]

16. Quigley, L.; O’Sullivan, O.; Beresford, T.P.; Ross, R.P.; Fitzgerald, G.F.; Cotter, P.D. High–throughput sequencing for detection of subpopulations of bacteria not previously associated with artisanal cheeses. Appl. Environ. Microbiol.; 2012; 78, pp. 5717-5723. [DOI: https://dx.doi.org/10.1128/AEM.00918-12]

17. Escobar-Zepeda, A.; Sanchez-Flores, A.; Quirasco, B.M. Metagenomic analysis of a Mexican ripened cheese reveals a unique complex microbiota. Food Microbiol.; 2016; 57, pp. 116-127. [DOI: https://dx.doi.org/10.1016/j.fm.2016.02.004]

18. Fuka, M.M.; Wallisch, S.; Engel, M.; Welzl, G.; Havranek, J.; Schloter, M. Dynamics of bacterial communities during the ripening process of different Croatian cheese types derived from raw ewe’s milk cheeses. PLoS ONE; 2013; 8, e80734. [DOI: https://dx.doi.org/10.1371/journal.pone.0080734]

19. Jin, H.; Mo, L.; Pan, L.; Hou, Q.; Li, C.; Darima, I.; Yu, J. Using PacBio sequencing to investigate the bacterial microbiota of traditional Buryatian cottage cheese and comparison with Italian and Kazakhstan artisanal cheeses. J. Dairy Sci.; 2018; 101, pp. 6885-6896. [DOI: https://dx.doi.org/10.3168/jds.2018-14403]

20. Riquelme, C.; Câmara, S.; Dapkevicius, M.; Vinuesa, P.; Silva, C.; Malcata, F.X.; Rego, O.A. Characterization of the bacterial biodiversity in Pico cheese (an artisanal Azorean food). Int. J. Food Microbiol.; 2015; 192, pp. 86-94. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2014.09.031]

21. Santos, M.; Benito, M.; Córdoba, M.; Alvarenga, N.; Herrera, S. Yeast community in traditional Portuguese Serpa cheese by culture–dependent and –independent DNA approaches. Int. J. Food Microbiol.; 2017; 262, pp. 63-70. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2017.09.013]

22. Barracosa, P.; Rosa, N.; Barros, M.; Pires, E. Selected Cardoon (Cynara cardunculus L.) Genotypes Suitable for PDO Cheeses in Mediterranean Regions. Chem. Biodivers.; 2018; 15, e1800110. [DOI: https://dx.doi.org/10.1002/cbdv.201800110]

23. Lima, S.F.; Bicalho, M.; Bicalho, R.C. Evaluation of milk sample fractions for characterization of milk microbiota from healthy and clinical mastitis cows. PLoS ONE; 2018; 13, e0193671. [DOI: https://dx.doi.org/10.1371/journal.pone.0193671]

24. Herlemann, D.P.; Labrenz, M.; Jürgens, K.; Bertilsson, S.; Waniek, J.J.; Andersson, A.F. Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J.; 2011; 5, pp. 1571-1579. [DOI: https://dx.doi.org/10.1038/ismej.2011.41]

25. Klindworth, A.; Pruesse, E.; Schweer, T.; Peplies, J.; Quast, C.; Horn, M.; Glöckner, F.O. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next–generation sequencing–based diversity studies. Nucleic Acids Res.; 2013; 41, e1. [DOI: https://dx.doi.org/10.1093/nar/gks808]

26. Tedersoo, L.; Bahram, M.; Põlme, S.; Kõljalg, U.; Yorou, N.S.; Wijesundera, R.; Ruiz, L.; Vasco–Palacios, A.; Thu, P.; Suija, A. et al. Global diversity and geography of soil fungi. Science; 2014; 346, 1256688. [DOI: https://dx.doi.org/10.1126/science.1256688]

27. 16S Metagenomic Sequencing Library Preparation. Preparing 16S Ribosomal RNA Gene Amplicons for the Illumina MiSeq System. 2013; Available online: https://emea.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf (accessed on 14 December 2020).

28. Comeau, A.M.; Douglas, G.M.; Langille, M.G.I. Microbiome Helper: A Custom and Streamlined Workflow for Microbiome Research. mSystems; 2017; 2, e00127-16. [DOI: https://dx.doi.org/10.1128/mSystems.00127-16]

29. Schmieder, R.; Edwards, R. Quality control and preprocessing of metagenomic datasets. Bioinformatics; 2011; 27, pp. 863-864. [DOI: https://dx.doi.org/10.1093/bioinformatics/btr026]

30. Schubert, M.; Lindgreen, S.; Orlando, L. AdapterRemoval v2: Rapid adapter trimming, identification, and read merging Findings Background. BMC Res. Notes; 2016; 9, 88. [DOI: https://dx.doi.org/10.1186/s13104-016-1900-2]

31. Bengtsson-Palme, J.; Ryberg, M.; Hartmann, M.; Branco, S.; Wang, Z.; Godhe, A.; De Wit, P.; Sánchez-García, M.; Ebersberger, I.; Sousa, F. et al. Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data. Methods Ecol. Evol.; 2013; 4, pp. 914-919. [DOI: https://dx.doi.org/10.1111/2041-210X.12073]

32. Rognes, T.; Flouri, T.; Nichols, B.; Quince, C.; Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ; 2016; 4, e2584. [DOI: https://dx.doi.org/10.7717/peerj.2584]

33. Edgar, R.C.; Haas, B.J.; Clemente, J.C.; Quince, C.; Knight, R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics; 2011; 27, pp. 2194-2200. [DOI: https://dx.doi.org/10.1093/bioinformatics/btr381]

34. McDonald, D.; Price, M.N.; Goodrich, J.; Nawrocki, E.P.; DeSantos, T.Z.; Probst, A.; Andersen, G.L.; Knight, R.; Hugenholtz, P. An improved Greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. ISME J.; 2012; 6, pp. 610-618. [DOI: https://dx.doi.org/10.1038/ismej.2011.139]

35. Abarenkov, K.; Zirk, A.; Piirmann, T.; Pöhönen, R.; Ivanov, F.; Nilsson, R.H.; Kõljalg, U. UNITE QIIME release for Fungi. UNITE Community; 2020; [DOI: https://dx.doi.org/10.15156/BIO/786385]

36. Bolyen, E.; Rideout, J.R.; Dillon, M.R.; Bokulich, N.A.; Abnet, C.C.; Al-Ghalith, G.A.; Alexander, H.; Alm, E.J.; Arumugam, M.; Asnicar, F. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol.; 2019; 37, pp. 852-857. [DOI: https://dx.doi.org/10.1038/s41587-019-0209-9]

37. Shannon, C.E. A mathematical theory of communication. Bell Syst. Tech. J.; 1948; 27, pp. 379-423. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x]

38. Simpson, E.H. Measurement of diversity. Nature; 1949; 163, 688. [DOI: https://dx.doi.org/10.1038/163688a0]

39. Good, I.J. The population frequencies of species and the estimation of population parameters. Biometrika; 1953; 40, pp. 237-264. [DOI: https://dx.doi.org/10.1093/biomet/40.3-4.237]

40. Chao, A. Non–parametric estimation of the number of classes in a population. Scand. J. Stat.; 1984; 11, pp. 265-270.

41. Chao, A.; Lee, S.M. Estimating the number of classes via sample coverage. J. Am. Stat. Assoc.; 1992; 87, pp. 210-217. [DOI: https://dx.doi.org/10.1080/01621459.1992.10475194]

42. Katoh, K.; Misawa, K.; Kuma, K.; Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res.; 2002; 30, pp. 3059-3066. [DOI: https://dx.doi.org/10.1093/nar/gkf436]

43. Price, M.N.; Dehal, P.S.; Arkin, A.P. FastTree 2—Approximately maximum–likelihood trees for large alignments. PLoS ONE; 2010; 5, e9490. [DOI: https://dx.doi.org/10.1371/journal.pone.0009490]

44. Faith, D.P. Conservation evaluation and phylogenetic diversity. Biol. Conserv.; 1992; 61, pp. 1-10. [DOI: https://dx.doi.org/10.1016/0006-3207(92)91201-3]

45. Jaccard, P. Nouvelles recherches sur la distribution floral. Bull. Soc. Vard. Sci. Nat.; 1908; 44, pp. 223-270.

46. Sørensen, T. A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on danish commons. Biol. Skr.; 1948; 5, pp. 1-34.

47. Lozupone, C.; Knight, R. UniFrac: A new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol.; 2005; 71, pp. 8228-8235. [DOI: https://dx.doi.org/10.1128/AEM.71.12.8228-8235.2005]

48. Lozupone, C.A.; Hamady, M.; Kelley, S.T.; Knight, R. Quantitative and qualitative beta diversity measures lead to different insights into factors that structure microbial communities. Appl. Environ. Microbiol.; 2007; 73, pp. 1576-1585. [DOI: https://dx.doi.org/10.1128/AEM.01996-06]

49. Kruskal, W.; Wallis, W. Use of Ranks in One–Criterion Variance Analysis. J. Am. Stat. Assoc.; 1952; 47, pp. 583-621. [DOI: https://dx.doi.org/10.1080/01621459.1952.10483441]

50. Anderson, M.J. A new method for non-parametric multivariate analysis of variance. Austral Ecol.; 2001; 26, pp. 32-46. [DOI: https://dx.doi.org/10.1111/j.1442-9993.2001.01070.pp.x]

51. Bokulich, N.A.; Kaehler, B.D.; Rideout, J.R.; Dillon, M.; Bolyen, E.; Knight, R.; Huttley, G.A.; Caporaso, J.G. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome; 2018; 6, 90. [DOI: https://dx.doi.org/10.1186/s40168-018-0470-z]

52. Quast, C.; Pruesse, E.; Yilmaz, P.; Gerken, J.; Schweer, T.; Yarza, P.; Peplies, J.; Glöckner, F.O. The SILVA ribosomal RNA gene database project: Improved data processing and web–based tools. Nucleic Acids Res.; 2013; 41, pp. D590-D596. [DOI: https://dx.doi.org/10.1093/nar/gks1219] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23193283]

53. Segata, N.; Izard, J.; Waldron, L.; Gevers, D.; Miropolsky, L.; Garrett, W.S.; Huttenhower, C. Metagenomic biomarker discovery and explanation. Genome Biol.; 2011; 12, R60. [DOI: https://dx.doi.org/10.1186/gb-2011-12-6-r60]

54. Montel, M.C.; Buchin, S.; Mallet, A.; Delbes-Paus, C.; Vuitton, D.A.; Desmasures, N.; Berthier, F. Traditional cheeses: Rich and diverse microbiota with associated benefits. Int. J. Food Microbiol.; 2014; 177, pp. 136-154. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2014.02.019] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24642348]

55. Quigley, L.; O’Sullivan, O.; Stanton, C.; Beresford, T.P.; Ross, R.P.; Fitzgerald, G.F.; Cotter, P.D. The complex microbiota of raw milk. FEMS Microbiol. Rev.; 2013; 37, pp. 664-698. [DOI: https://dx.doi.org/10.1111/1574-6976.12030] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23808865]

56. Delavenne, E.; Mounier, J.; Asmani, K.; Jany, J.L.; Barbier, G.; Le Blay, G. Fungal diversity in cow, goat and ewe milk. Int. J. Food Microbiol.; 2011; 151, pp. 247-251. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2011.08.029] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21944758]

57. Esteban-Blanco, C.; Gutiérrez-Gil, B.; Puente-Sánchez, F.; Marina, H.; Tamames, J.; Acedo, A.; Arranz, J.J. Microbiota characterization of sheep milk and its association with somatic cell count using 16s rRNA gene sequencing. J. Anim. Breed Genet.; 2020; 137, pp. 73-83. [DOI: https://dx.doi.org/10.1111/jbg.12446] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31602717]

58. Fadda, M.E.; Mossa, V.; Pisano, M.B.; Deplano, M.; Cosentino, S. Occurrence and characterization of yeasts isolated from artisanal Fiore Sardo cheese. Int. J. Food Microbiol.; 2004; 15, pp. 51-59. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2004.02.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15240074]

59. Gaya, P.; Medina, M.; Nuñez, M. Enterobacteriaceae, coliforms, faecal coliforms and salmonellas in raw ewes’ milk. J. Appl. Bacteriol.; 1987; 62, pp. 321-326. [DOI: https://dx.doi.org/10.1111/j.1365-2672.1987.tb04927.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3597207]

60. Salmerón, J.; Vega, C.; Pérez-Elortondo, F.; Albisu, M.; Barrón, L. Effect of pasteurization and seasonal variations in the microflora of ewe’s milk for cheesemaking. Food Microbiol.; 2002; 19, pp. 167-174. [DOI: https://dx.doi.org/10.1006/fmic.2001.0475]

61. Spanamberg, A.; Fraga, C.; Ferreiro, L.; Aguinsky, M.; Sanches, E.; Roehe, C.; Lautert, C.; Santurio, J. Yeasts in the Raw Ewe’s Milk. Acta Sci. Vet.; 2014; 42, 1236.

62. Fröhlich-Wyder, M.T.; Arias-Roth, E.; Jakob, E. Cheese yeasts. Yeast; 2019; 36, pp. 129-141. [DOI: https://dx.doi.org/10.1002/yea.3368]

63. Akdouche, L.; Aissi, M.; Saadi, A. Prevalence and Identification of Yeasts Responsible for Mastitis in Dairy Cattle Farms in the Sidi Lahcene Region in the Wilaya of Sidi Bel abbes–Algeria. J. Adv. Dairy Res.; 2018; 6, 206. [DOI: https://dx.doi.org/10.4172/2329-888X.1000206]

64. Zaragoza, C.S.; Olivares, R.A.; Watty, A.E.; Moctezuma, L.; Tanaca, L.V. Yeasts isolation from bovine mammary glands under different mastitis status in the Mexican High Plateu. Rev. Iberoam. Micol.; 2011; 28, pp. 79-82. [DOI: https://dx.doi.org/10.1016/j.riam.2011.01.002]

65. Specification book—Serra da Estrela cheese protected designation of origin. 2011; Available online: https://tradicional.dgadr.gov.pt/images/prod_imagens/queijos/docs/CE_Queijo_Serra.pdf (accessed on 15 March 2021).

66. Callon, C.; Duthoit, F.; Delbès, C.; Ferrand, M.; Frileux, Y.; Crémoux, R.; Montel, M. Stability of microbial communities in goat milk during a lactation year: Molecular approaches. Syst. Appl. Microbiol.; 2007; 30, pp. 547-560. [DOI: https://dx.doi.org/10.1016/j.syapm.2007.05.004]

67. Conceição, C.; Martins, P.; Alvarenga, N.; Dias, J.; Lamy, E.; Garrido, L.; Gomes, S.; Freitas, S.; Belo, A.; Brás, T. et al. Cynara cardunculus: Use in Cheesemaking and Pharmaceutical Applications. Technological Approaches for Novel Applications in Dairy Processing; Kuka, N. IntechOpen: London, UK, 2018; pp. 73-107. [DOI: https://dx.doi.org/10.5772/intechopen.76530]

68. Fernández-Salguero, J.; Sánchez, E.; Gómez, R.; Mata, C.; Vioque, M.; Tejada, L. Preliminary study of microbiological quality of cardoons of the genus Cynara L. used in manufacture of traditional cheeses. Milchwissenschaft; 1999; 54, pp. 688-689.

69. Gómez, R.; Sánchez, E.; Vioque, M.; Ferreira, J.; Tejada, L.; Fernández-Salguero, J. Microbiological characteristics of ewes’ milk cheese manufactured using aqueous extracts of flowers from various species of cardoon Cynara L. Milchwissenschaft; 2001; 56, pp. 16-19.

70. Ratão, I. Microbiological and Chemical Characterization of Traditional Cheese Made from Milk Produced by the Algarvian Goat Breed. Ph.D. Thesis; Cranfield University: Cranfield, UK, 2011.

71. Bhadra, B.; Rao, R.S.; Singh, P.K.; Sarkar, P.K.; Shivaji, S. Yeasts and yeast–like fungi associated with tree bark: Diversity and identification of yeasts producing extracellular endoxylanases. Curr. Microbiol.; 2008; 56, pp. 489-494. [DOI: https://dx.doi.org/10.1007/s00284-008-9108-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18219522]

72. Egbuta, M.; Mwanza, M.; Babalola, O. A Review of the Ubiquity of Ascomycetes Filamentous Fungi in Relation to Their Economic and Medical Importance. Adv. Microbiol.; 2016; 6, pp. 1140-1158. [DOI: https://dx.doi.org/10.4236/aim.2016.614103]

73. Fox, P.F.; McSweeney, P. Cheese: An Overview. Cheese—Chemistry, Physics and Microbiology; 4th ed. Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 5-21. [DOI: https://dx.doi.org/10.1016/B978-0-12-417012-4.00001-6]

74. Tuszyński, T.; Satora, P. Microbiological characteristics of the Węgierka Zwykła plum orchard in submontane region. Pol. J. Food Nutr. Sci.; 2003; 53, pp. 43-48.

75. Vorholt, J.A. Microbial life in the phyllosphere. Nat. Rev. Microbiol.; 2012; 10, pp. 828-840. [DOI: https://dx.doi.org/10.1038/nrmicro2910] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23154261]

76. Cotter, P.D.; Beresford, T.P. Microbiome Changes During Ripening. Cheese—Chemistry, Physics and Microbiology; 4th ed. Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 389-409. [DOI: https://dx.doi.org/10.1016/B978-0-12-417012-4.00015-6]

77. Irlinger, F.; Helinck, S.; Jany, J.L. Secondary and Adjunct Cultures. Cheese—Chemistry, Physics and Microbiology; 4th ed. Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 273-300. [DOI: https://dx.doi.org/10.1016/B978-0-12-417012-4.00011-9]

78. Picon, A. Cheese Microbial Ecology and Safety. Global Cheesemaking Technology: Cheese Quality and Characteristics; 1st ed. Papademas, P.; Bintsis, T. John Wiley & Sons Ltd.: Hoboken, NJ, USA, 2017; pp. 71-99. [DOI: https://dx.doi.org/10.1002/9781119046165.ch0d]

79. Guerreiro, J. Molecular Methods for Authentication of Protected Denomination of Origin (PDO) Cheeses. Ph.D. Thesis; University of Nottingham: Nottingham, UK, 2006.

80. Cardinali, F.; Ferrocino, I.; Milanović, V.; Belleggia, L.; Corvaglia, M.R.; Garofalo, C.; Foligni, R.; Mannozzi, C.; Mozzon, M.; Cocolin, L. et al. Microbial communities and volatile profile of Queijo de Azeitão PDO cheese, a traditional Mediterranean thistle-curdled cheese from Portugal. Food Res. Int.; 2021; 147, 110537. [DOI: https://dx.doi.org/10.1016/j.foodres.2021.110537]

81. Kačániová, M.; Terentjeva, M.; Kunová, S.; Haščík, P.; Kowalczewski, P.Ł.; Štefániková, J. Diversity of microbiota in Slovak summer ewes’ cheese “Bryndza”. Open Life Sci.; 2021; 16, pp. 277-286. [DOI: https://dx.doi.org/10.1515/biol-2021-0038] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33817319]

82. Pasquale, I.; Cagno, R.; Buchin, S.; Angelis, M.; Gobbetti, M. Spatial Distribution of the Metabolically Active Microbiota within Italian PDO Ewes’ Milk Cheeses. PLoS ONE; 2016; 11, e0153213. [DOI: https://dx.doi.org/10.1371/journal.pone.0153213]

83. Pereira-Dias, S.; Potes, M.E.; Marinho, A.; Malfeito-Ferreira, M.; Loureiro, V. Characterisation of yeast flora isolated from an artisanal Portuguese ewes’ cheese. Int. J. Food Microbiol.; 2000; 60, pp. 55-63. [DOI: https://dx.doi.org/10.1016/S0168-1605(00)00323-8]

84. Spyrelli, E.; Stamatiou, A.; Tassou, C.; Nychas, G.; Doulgeraki, A. Microbiological and Metagenomic Analysis to Assess the Effect of Container Material on the Microbiota of Feta Cheese during Ripening. Fermentation; 2020; 6, 12. [DOI: https://dx.doi.org/10.3390/fermentation6010012]

85. Padilla, B.; Manzanares, P.; Belloch, C. Yeast species and genetic heterogeneity within Debaryomyces hansenii along the ripening process of traditional ewes’ and goats’ cheeses. Food Microbiol.; 2014; 38, pp. 160-166. [DOI: https://dx.doi.org/10.1016/j.fm.2013.09.002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24290639]

86. Lavoie, K.; Touchette, M.; St-Gelais, D.; Labrie, S. Characterization of the fungal microflora in raw milk and specialty cheeses of the province of Quebec. Dairy Sci. Technol.; 2012; 92, pp. 455-468. [DOI: https://dx.doi.org/10.1007/s13594-011-0051-4]

87. Soliman, N.; Aly, S. Occurrence and identification of yeast species isolated from Egyptian Karish cheese. J. Yeast Fungal Res.; 2011; 2, pp. 59-64.

88. Fox, P.F.; Guinee, T.P.; Cogan, T.M.; McSweeney, P.L.H. Pathogens in Cheese and Foodborne Illnesses. Fundamentals of Cheese Science; Fox, P.F.; Guinee, T.P.; Cogan, T.M.; McSweeney, P.L.H. Springer: Berlin/Heidelberg, Germany, 2017; pp. 681-713. [DOI: https://dx.doi.org/10.1007/978-1-4899-7681-9_19]

89. Parente, E.; Cogan, T.M.; Powell, I.B. Starter Cultures: General Aspects. Cheese—Chemistry, Physics and Microbiology; 4th ed. Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 201-226.

90. Cavanagh, D.; Fitzgerald, G.F.; McAuliffe, O. From field to fermentation: The origins of Lactococcus lactis and its domestication to the dairy environment. Food Microbiol.; 2015; 47, pp. 45-61. [DOI: https://dx.doi.org/10.1016/j.fm.2014.11.001]

91. Surber, G.; Schäper, C.; Wefers, D.; Rohm, H.; Jaros, D. Exopolysaccharides from Lactococcus lactis affect manufacture, texture and sensory properties of concentrated acid milk gel suspensions (fresh cheese). Int. Dairy J.; 2021; 112, 104854. [DOI: https://dx.doi.org/10.1016/j.idairyj.2020.104854]

92. McAuliffe, O. Genetics of Lactic Acid Bacteria. Cheese—Chemistry, Physics and Microbiology; 4th ed. Cotter, P.D.; Everett, D.W. Academic Press: Cambridge, MA, USA, 2017; pp. 227-247. [DOI: https://dx.doi.org/10.1016/B978-0-12-417012-4.00009-0]

93. Bhardwaj, A.; Kapila, S.; Mani, J.; Malik, R.K. Comparison of susceptibility to opsonic killing by in vitro human immune response of Enterococcus strains isolated from dairy products, clinical samples and probiotic preparation. Int. J. Food Microbiol.; 2009; 128, pp. 513-515. [DOI: https://dx.doi.org/10.1016/j.ijfoodmicro.2008.10.010]

94. Nagy, E.; Dlauchy, D.; Medeiros, A.O.; Péter, G.; Rosa, C.A. Yarrowia porcina sp. nov. and Yarrowia bubula f.a. sp. nov., two yeast species from meat and river sediment. Antonie Van Leeuwenhoek; 2014; 105, pp. 697-707. [DOI: https://dx.doi.org/10.1007/s10482-014-0125-4]