Advertisement

Milk metagenomics and cheese-making properties as affected by indoor farming and summer highland grazing

Open AccessPublished:November 15, 2022DOI:https://doi.org/10.3168/jds.2022-22449

      ABSTRACT

      The study of the complex relationships between milk metagenomics and milk composition and cheese-making efficiency as affected by indoor farming and summer highland grazing was the aim of the present work. The experimental design considered monthly sampling (over 5 mo) of the milk produced by 12 Brown Swiss cows divided into 2 groups: the first remained on a lowland indoor farm from June to October, and the second was moved to highland pastures in July and then returned to the lowland farm in September. The resulting 60 milk samples (2 kg each) were used to analyze milk composition, milk coagulation, curd firming, and syneresis processes, and to make individual model cheeses to measure cheese yields and nutrient recoveries in the cheese. After DNA extraction and Illumina Miseq sequencing, milk microbiota amplicons were also processed by means of an open-source pipeline called Quantitative Insights Into Microbial Ecology (Qiime2, version 2018.2; https://qiime2.org ). Out of a total of 44 taxa analyzed, 13 bacterial taxa were considered important for the dairy industry (lactic acid bacteria, LAB, 5 taxa; and spoilage bacteria, 4) and for human (other probiotics, 2) and animal health (pathogenic bacteria, 2). The results revealed the transhumant group of cows transferred to summer highland pastures showed an increase in almost all the LAB taxa, bifidobacteria, and propionibacteria, and a reduction in spoilage taxa. All the metagenomic changes disappeared when the transhumant cows were moved back to the permanent indoor farm. The relationships between 17 microbial traits and 30 compositional and technological milk traits were investigated through analysis of correlation and latent explanatory factor analysis. Eight latent factors were identified, explaining 75.3% of the total variance, 2 of which were mainly based on microbial traits: pro-dairy bacteria (14% of total variance, improving during summer pasturing) and pathogenic bacteria (6.0% of total variance). Some bacterial traits contributed to other compositional-technological latent factors (gelation, udder health, and caseins).

      Key words

      INTRODUCTION

      The microbiota of milk has been studied for many decades because of the important relationships between milk microorganisms on one side, and milk characteristics, end product, economic impact, and nutritional and health values on the other (
      • 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.
      ;
      • Boor K.J.
      • Wiedmann M.
      • Murphy S.
      • Alcaine S.
      A 100-Year Review: Microbiology and safety of milk handling.
      ;
      • Issa A.T.
      • Tahergorabi R.
      Chapter 22: Milk bacteria and gastrointestinal tract: Microbial composition of milk.
      ). The most interesting of the favorable relationships between microbes and the various milk characteristics are those that concern the role of microbial species, especially lactic acid bacteria (LAB), in relation to milk end products, particularly cheese (
      • Skeie S.
      Characteristics in milk influencing the cheese yield and cheese quality.
      ;
      • Ardö Y.
      • McSweeney P.L.H.
      • Magboul A.A.A.
      • Upadhyay V.K.
      • Fox . F.
      Chapter 18—Biochemistry of cheese ripening: Proteolysis.
      ;
      • Nam J.H.
      • Cho Y.S.
      • Rackerby B.
      • Goddik L.
      • Park S.H.
      Shifts of microbiota during cheese production: Impact on production and quality.
      ), and digestion and intestinal functions and integrity in human consumers (prebiotics and probiotics;
      • Aryana K.J.
      • Olson D.W.
      A 100-Year Review: Yogurt and other cultured dairy products.
      ;
      • Nyanzi R.
      • Jooste P.J.
      • Buys E.M.
      Invited review: Probiotic yogurt quality criteria, regulatory framework, clinical evidence, and analytical aspects.
      ). The most interesting unfavorable relationships are the potential effects of some microbial species on the spoilage of milk and dairy products (
      • 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.
      ;
      • Martin N.H.
      • Torres-Frenzel P.
      • Wiedmann M.
      Invited review: Controlling dairy product spoilage to reduce food loss and waste.
      ) and on the health of lactating animals (mastitis;
      • Andrews T.
      • Neher D.A.
      • Weicht T.R.
      • Barlow J.W.
      Mammary microbiome of lactating organic dairy cows varies by time, tissue site, and infection status.
      ) and of human consumers (pathogenic activities;
      • Verraes C.
      • Vlaemynck G.
      • Van Weyenberg S.
      • De Zutter L.
      • Daube G.
      • Sindic M.
      • Uyttendaele M.
      • Herman L.
      A review of the microbiological hazards of dairy products made from raw milk.
      ).
      The recent development of metagenomics is now expanding our knowledge of these aspects of milk microbiota (
      • Addis M.F.
      • Tanca A.
      • Uzzau S.
      • Oikonomou G.
      • Bicalho R.C.
      • Moroni P.
      The bovine milk microbiota: Insights and perspectives from -omics studies.
      ;
      • Parente E.
      • Ricciardi A.
      • Zotta T.
      The microbiota of dairy milk: A review.
      ). Traditional microbiological studies were based on identifying, isolating, characterizing, and counting individual microbial species or strains (
      • Tilocca B.
      • Costanzo N.
      • Morittu V.M.
      • Spina A.A.
      • Soggiu A.
      • Britti D.
      • Roncada P.
      • Piras C.
      Milk microbiota: Characterization methods and role in cheese production.
      ). With metagenomics, the entire milk microbiota composition can be identified and characterized. Alongside ecological studies of milk microbial communities, we can now gather new information on many microbial taxa involved in different compositional, technological, and nutritional properties of milk.
      One of the most important, but also difficult to study, issues is the effect of dairy system, particularly pasture grazing, on the milk microbiota, and the effect of the microbiota on the properties associated with milk processing and end products, nutritional value, and consumer health (
      • Doyle C.J.
      • Gleeson D.
      • O’Toole P.W.
      • Cotter P.D.
      Impacts of seasonal housing and teat preparation on raw milk microbiota: A high-throughput sequencing study.
      ). The difficulties lie mainly in disentangling the confounding effects of environment, management, animal characteristics, season, feedstuffs, and hygiene (
      • Du B.
      • Meng L.
      • Liu H.
      • Zheng N.
      • Zhang Y.
      • Guo X.
      • Zhao S.
      • Li F.
      • Wang J.
      Impact of milking and housing environment on milk microbiota.
      ).
      The opposite extremes of dairy farming are represented by the intensive indoor system, with year-round use of TMR, and forage-based systems, where the cows are kept at pasture day and night, and have only limited access to compound feed during milking (
      • O’Callaghan T.F.
      • Mannion D.T.
      • Hennessy D.
      • McAuliffe S.
      • O’Sullivan M.G.
      • Leeuwendaal N.
      • Beresford T.P.
      • Dillon P.
      • Kilcawley K.N.
      • Sheehan J.J.
      • Ross R.P.
      • Stanton C.
      Effect of pasture versus indoor feeding systems on quality characteristics, nutritional composition, and sensory and volatile properties of full-fat Cheddar cheese.
      ). Among the latter, farms that practice transhumance to temporary farms on highland summer pastures are very distinctive for both the extreme environmental conditions the animals face and the renowned quality and nutritional value of their dairy products (
      • Buchin S.
      • Martin B.
      • Dupont D.
      • Bornard A.
      • Achilleos C.
      Influence of the composition of Alpine highland pasture on the chemical, rheological and sensory properties of cheese.
      ). Little is known of the extent to which the specificity of mountain dairy products is due to the milk microbiota. We hypothesized that the summer transhumance to highland summer pastures would alter the microbial population of the milk, and that milk microbiota could affect the cheese-making process.
      The general aim of this research, therefore, was to study the milk microbiota in indoor housing versus summer highland grazing and its relationships with milk quality and technological properties, with particular emphasis on the bacterial taxa related to various specific activities [cheese-making, health maintenance (probiotics), milk spoilage, and pathogeny], and the effects of moving the cows from indoor conditions to the summer highland pasture, and then back to indoor conditions.

      MATERIALS AND METHODS

      Experimental Design and Milk Sampling

      This study is part of a larger project studying the effects of the transhumance of cows to summer highland pastures on their productivity, and on the chemical, technological, and microbiological characteristics of the milk produced. In this project, 2 groups of cows (one kept solely indoors, the other moved to summer pasture) were monitored before, during, and after summer transhumance. Details on the environmental conditions and methodology can be found in 2 previous studies: the first dealing with the cows' body condition and milk yield, milk composition, and cheese-making efficiency (
      • Saha S.
      • Amalfitano N.
      • Sturaro E.
      • Schiavon S.
      • Tagliapietra F.
      • Bittante G.
      • Carafa I.
      • Franciosi E.
      • Gallo L.
      Effects of summer transhumance of dairy cows to alpine pastures on body condition, milk yield and composition, and cheese making efficiency.
      ), and the second reporting some preliminary data on milk microbial counts (
      • Carafa I.
      • Navarro I.C.
      • Bittante G.
      • Tagliapietra F.
      • Gallo L.
      • Tuohy K.
      • Franciosi E.
      Shift in the cow milk microbiota during alpine pasture as analyzed by culture dependent and high-throughput sequencing techniques.
      ). All samples and measurements were obtained during the farms' normal milking procedures; therefore ethics commission approval was not required.
      In line with the aims of this project, the present study was carried out at 2 farms in a mountain area (Trentino Province, northeastern Italian Alps): (1) a modern, permanent farm in a valley (Malè, Trento, Italy; 737 m above sea level), where lactating cows are loose-housed indoors, fed TMR, and milked in a milking parlor; and (2) a temporary summer highland farm (Malga Juribello, within the “Paneveggio – Pale di San Martino” Nature Reserve, Passo Rolle, Trento, Italy; 1,860 m above sea level), where cows are kept at pasture day and night, and are milked and given a supplementary compound feed (3 to 6 kg/d, according to milk production) in a milking parlor in an old barn.
      Briefly, the experimental design consisted of the following steps: the selection at the end of May of 12 healthy, multiparous, early-lactation Brown Swiss cows on the permanent lowland farm, and their random division into 2 groups of 6 cows each. The cows in both groups had similar (P > 0.05, based on t-test) parity numbers (2.5 and 2.8, respectively) and DIM (143 and 120, respectively); all 12 cows were kept together in the same indoor pen before the start of the experiment (beginning of June) and during the first (June) and last month of sampling (October). During summer (July, August, and September) 1 of the 2 groups was moved to the temporary highland farm (high group) at the beginning of July, and returned to the permanent farm at the end of September; the other group remained indoors on the permanent lowland farm (low group). Monthly samples of milk from each cow at the evening milking, from mid-June to mid-October (5 samples per cow) were taken in both farms (60 samples in total).
      The samples taken in June represent the initial condition of the 2 groups: having been reared together, nonsignificant differences for all traits were expected. The samples taken in July, August, and September represent the effects of the 2 farming conditions: indoors in the valley versus on highland pasture. The fifth sample, taken in October after the 2 groups had been together again for a month, represented potential carryover effects of summer pasturing on indoor rearing.
      Each sample was divided into 2 aliquots: the first (50 mL) was immediately frozen in liquid nitrogen, taken to the Research and Innovation Centre, Food Quality and Nutrition Department of the Fondazione Edmund Mach (San Michele all'Adige, Trento, Italy), and stored at −80°C before microbiological analyses within 3 mo. The second aliquot (2 L) was immediately refrigerated at 4°C and transported to the Milk Laboratory of the Department of Agronomy, Food, Natural Resources, Animals and Environment of the University of Padova (Legnaro, Padua, Italy) for evaluation of milk quality, cheese-making aptitude, and cheese yield.

      Metagenomic Analyses

      In a previous study (
      • Carafa I.
      • Navarro I.C.
      • Bittante G.
      • Tagliapietra F.
      • Gallo L.
      • Tuohy K.
      • Franciosi E.
      Shift in the cow milk microbiota during alpine pasture as analyzed by culture dependent and high-throughput sequencing techniques.
      ), we analyzed in detail the bacterial counts of milk regarding a preliminary comparison of the samples taken during summer pasturing (n = 18) with those taken indoors on the permanent farm (n = 42), without taking into account the effects of group and month. In the present study, the principal bacterial taxa identified by Qiime2 (version 2018.2; https://qiime2.org ;
      • 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.
      • Bai Y.
      • Bisanz J.E.
      • Bittinger K.
      • Brejnrod A.
      • Brislawn C.J.
      • Brown C.T.
      • Callahan B.J.
      • Caraballo-Rodríguez A.M.
      • Chase J.
      • Cope E.K.
      • Da Silva R.
      • Diener C.
      • Dorrestein P.C.
      • Douglas G.M.
      • Durall D.M.
      • Duvallet C.
      • Edwardson C.F.
      • Ernst M.
      • Estaki M.
      • Fouquier J.
      • Gauglitz J.M.
      • Gibbons S.M.
      • Gibson D.L.
      • Gonzalez A.
      • Gorlick K.
      • Guo J.
      • Hillmann B.
      • Holmes S.
      • Holste H.
      • Huttenhower C.
      • Huttley G.A.
      • Janssen S.
      • Jarmusch A.K.
      • Jiang L.
      • Kaehler B.D.
      • Bin Kang K.
      • Keefe C.R.
      • Keim P.
      • Kelley S.T.
      • Knights D.
      • Koester I.
      • Kosciolek T.
      • Kreps J.
      • Langille M.G.I.
      • Lee J.
      • Ley R.
      • Liu Y.X.
      • Loftfield E.
      • Lozupone C.
      • Maher M.
      • Marotz C.
      • Martin B.D.
      • McDonald D.
      • McIver L.J.
      • Melnik A.V.
      • Metcalf J.L.
      • Morgan S.C.
      • Morton J.T.
      • Naimey A.T.
      • Navas-Molina J.A.
      • Nothias L.F.
      • Orchanian S.B.
      • Pearson T.
      • Peoples S.L.
      • Petras D.
      • Preuss M.L.
      • Pruesse E.
      • Rasmussen L.B.
      • Rivers A.
      • Robeson M.S.
      • Rosenthal P.
      • Segata N.
      • Shaffer M.
      • Shiffer A.
      • Sinha R.
      • Song S.J.
      • Spear J.R.
      • Swafford A.D.
      • Thompson L.R.
      • Torres P.J.
      • Trinh P.
      • Tripathi A.
      • Turnbaugh P.J.
      • Ul-Hasan S.
      • van der Hooft J.J.J.
      • Vargas F.
      • Vázquez-Baeza Y.
      • Vogtmann E.
      • von Hippel M.
      • Walters W.
      • Wan Y.
      • Wang M.
      • Warren J.
      • Weber K.C.
      • Williamson C.H.D.
      • Willis A.D.
      • Xu Z.Z.
      • Zaneveld J.R.
      • Zhang Y.
      • Zhu Q.
      • Knight R.
      • Caporaso J.G.
      Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2.
      ) were classified into 4 categories according to their potential activity in relation to cheese-making (LAB), other probiotics, spoilage, and pathogenic properties of milk, and were statistically analyzed, disentangling the effects of individual cows, groups of cows, month of sampling, and environmental and feeding conditions. The relationships between the metagenomic information and the qualitative and cheese-making properties of milk were also explored.
      In brief, genomic DNA was extracted using the DNeasyPower Food Microbial Kit (Qiagen) and quantified using a Nanodrop8800 Fluorospectrometer (Thermo Scientific). The Miseq Library (Illumina) was prepared according to the authors' recommendations and followed by Illumina sequencing. All the sequencing data were processed using Qiime2, and the final data were deposited in the NCBI Sequence Read Archive ( https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA528228 ), where they can be accessed under accession number PRJNA528228.

      Bacterial Categories

      The 4 categories of the identified taxa are here briefly reported. The LAB category includes the taxa belonging to the Lactobacillaceae family, including Lactobacillus, Leuconostoc(
      • Zheng J.
      • Wittouck S.
      • Salvetti E.
      • Franz C.M.A.P.
      • Harris H.M.B.
      • Mattarelli P.
      • O’Toole P.W.
      • Pot B.
      • Vandamme P.
      • Walter J.
      • Watanabe K.
      • Wuyts S.
      • Felis G.E.
      • Gänzle M.G.
      • Lebeer S.
      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.
      ), Lactococcus, and Enterococcus(
      • Gagnon M.
      • Ouamba A.J.K.
      • LaPointe G.
      • Chouinard P.Y.
      • Roy D.
      Prevalence and abundance of lactic acid bacteria in raw milk associated with forage types in dairy cow feeding.
      ). The “other probiotics” category includes all the taxa belonging to the genera Propionibacterium (
      • Rabah H.
      • Rosa do Carmo F.
      • Jan G.
      Dairy propionibacteria: Versatile probiotics.
      ) and Bifidobacterium (
      • Prasanna P.H.P.
      • Grandison A.S.
      • Charalampopoulos D.
      Bifidobacteria in milk products: An overview of physiological and biochemical properties, exopolysaccharide production, selection criteria of milk products and health benefits.
      ). The “spoilage bacteria” category includes all taxa belonging to the order Clostridiales (
      • Burtscher J.
      • Hobl L.
      • Kneifel W.
      • Domig K.J.
      Short communication: Clostridial spore counts in vat milk of Alpine dairies.
      ) and the genera Pseudomonas (
      • Meng L.
      • Zhang Y.
      • Liu H.
      • Zhao S.
      • Wang J.
      • Zheng N.
      Characterization of Pseudomonas spp. and associated proteolytic properties in raw milk stored at low temperatures.
      ), Kocuria (
      • Ribeiro-Júnior J.C.
      • Tamanini R.
      • Alfieri A.A.
      • Beloti V.
      Effect of milk bactofugation on the counts and diversity of thermoduric bacteria.
      ), and Alicyclobacillus (
      • Pornpukdeewattana S.
      • Jindaprasert A.
      • Massa S.
      Alicyclobacillus spoilage and control—A review.
      ). Finally, the “pathogenic” category includes all taxa belonging to the genus Staphylococcus (
      • Gebremedhin E.Z.
      • Ararso A.B.
      • Borana B.M.
      • Kelbesa K.A.
      • Tadese N.D.
      • Marami L.M.
      • Sarba E.J.
      Isolation and identification of Staphylococcus aureus from milk and milk products, associated factors for contamination, and their antibiogram in Holeta, Central Ethiopia.
      ), and the family Enterobacteriaceae (
      • Anand S.K.
      • Griffiths M.W.
      Pathogens in milk: Enterobacteriaceae.
      ). The remaining 31 bacterial taxa were grouped as “other milk bacteria.”

      Milk Composition, Cheese-Making Aptitude, and Cheese Yield

      Within 20 h of collection, the second aliquot was used to evaluate milk composition, traditional coagulation properties, whey composition, cheese yields, and milk nutrient recoveries in the curd, and for curd-firming modeling and the manufacture of model cheeses.

      Milk Composition

      In brief, milk composition traits (TS, fat, nonfat solids, protein, casein, lactose, and urea contents) were evaluated with a Milkoscan FT2 infrared analyzer (Foss A/S). Somatic cell counts were obtained with a Fossomatic Minor FC counter (Foss A/S) and then log-transformed to SCS. The fat/protein ratio and casein number (casein as a percentage of protein) were computed from the fat, protein, and casein contents. In the present study, we used the qualitative and technological characteristics of the milk samples to search for potential relationships with the metagenomic information from the same milk samples.

      Milk Coagulation Properties and Curd-Firming Modeling Parameters

      Traditional milk coagulation properties were obtained using a lactodynamograph (Formagraph; Foss A/S) according to
      • Cecchinato A.
      • Cipolat-Gotet C.
      • Casellas J.
      • Penasa M.
      • Rossoni A.
      • Bittante G.
      Genetic analysis of rennet coagulation time, curd-firming rate, and curd firmness assessed over an extended testing period using mechanical and near-infrared instruments.
      and consisted of the following: rennet coagulation time (RCT, min), curd-firming time (k20, min), and curd firmness (a) 30, 45, and 60 min after rennet addition (a30, a45, and a60; mm). Estimates of the curd-firming and syneresis equation parameters of each individual milk sample were obtained by extracting 240 curd firmness values (one every 15 s for 60 min) from the lactodynamograph. The equation parameters were rennet coagulation time by equation (RCTeq, min), the curd-firming instant rate constant (kCF, %/min), the syneresis instant rate constant (kSR, %/min), maximum curd firmness (CFmax, mm), and time to reach CFmax (tmax, min;
      • Malchiodi F.
      • Cecchinato A.
      • Penasa M.
      • Cipolat-Gotet C.
      • Bittante G.
      Milk quality, coagulation properties, and curd firmness modeling of purebred Holsteins and first- and second-generation crossbred cows from Swedish Red, Montbéliarde, and Brown Swiss bulls.
      ).

      Model Cheese-Making, Cheese Yields, Milk Nutrient Recoveries in the Curd, and Whey Composition

      A larger aliquot (1,500 mL) from each milk sample was used to manufacture a model cheese according to a procedure that replicates the process for making full-fat cheese (
      • Stocco G.
      • Cipolat-Gotet C.
      • Gasparotto V.
      • Cecchinato A.
      • Bittante G.
      Breed of cow and herd productivity affect milk nutrient recovery in curd, and cheese yield, efficiency and daily production.
      ). Briefly, the procedure was as follows. (1) Each milk sample was poured into a stainless-steel laboratory vat. (2) The vat was placed in a water bath and heated to 35°C for 30 min. (3) Rennet solution (8 mL; Hansen Standard 215 with 80 ± 5% chymosin and 20 ± 5% pepsin; Pacovis Amrein AG), freshly diluted to 4.29% (wt/vol) in distilled water, was added. (4) After coagulation, the curd was cut. (5) The curd was drained for 30 min. (6) The resulting whey was collected, weighed, and sampled. (7) The chemical composition (TS, fat, protein, and lactose) of the whey samples was analyzed with a Milkoscan FT2 infrared analyzer (Foss A/S). (8) The curd was pressed for 30 min at 250 kPa in a cheese-pressing machine, turning every 10 min. (9) The pressed curd wheel was soaked in a brine solution (20% NaCl) for 30 min. (10) After brining, the cheese wheels were weighed and the pH measured with a Crison Basic 20 electrode (Crison Instruments SA). The percentage yields of fresh cheese (%CYCURD) and cheese solids (%CYSOLIDS) were determined, as well as the following nutrient recovery traits (REC; the quantity of a given nutrient in the cheese as a percentage of the same nutrient in the milk processed): milk fat (RECFAT, %), milk protein (RECPROTEIN, %), total milk solids (RECSOLIDS, %), and milk energy (RECENERGY, %).

      Statistical Analysis

      All relative bacterial abundances were log10 transformed. Two samples were excluded from the statistical analysis because of a lack of microbiological data in one and of qualitative-technological data in the other. All bacterial and qualitative-technological data were checked to identify and exclude outlier values (outside the interval ±3 SD of the mean).

      Mixed-Model ANOVA

      The log10-transformed relative abundances obtained from the milk metagenomic analysis were analyzed with a linear mixed model in the R environment (

      R Core Team. 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.

      ), which included the fixed effects of the month × group interaction (10 levels: 5 mo, June to October; and 2 groups, high and low), and the random effect of cow within group. Polynomial contrasts were estimated between the 5 least squares means of month within the low group to determine the response curve of each trait (linear, quadratic, and cubic components) during the 5 mo the cows were kept indoors on the permanent farm as a measure of the effect of season and advancing lactation in the control group. Contrasts between the high and low groups were estimated within each month to test for the following: homogeneity of groups in the same environment (indoors) at the beginning of the trial (June); the effect of transhumance to highland pasture during the summer months (July, August, and September) compared with the control indoor group; and the carryover effect of summer pasture on the high group after returning to indoor conditions on the permanent farm (October). A similar model with the month × group interaction treated as a random factor was run to quantify the relative importance of this environmental or diet factor, individual animal within group, and residual factors not accounted for by the model. The variances in these 3 sources of variation were expressed as percentages of their sum (total variance).
      The model used here is the same model that we used in the previous study (
      • Saha S.
      • Amalfitano N.
      • Sturaro E.
      • Schiavon S.
      • Tagliapietra F.
      • Bittante G.
      • Carafa I.
      • Franciosi E.
      • Gallo L.
      Effects of summer transhumance of dairy cows to alpine pastures on body condition, milk yield and composition, and cheese making efficiency.
      ) to analyze the qualitative and technological properties of milk. Thus, those results are not reported and discussed here, except where they are useful for interpreting relationships with the metagenomics data.

      Correlation Analysis and Latent Explanatory Factor Analysis

      The 2 data sets of metagenomic relative abundances (only for the bacterial categories of interest; i.e., LAB, other probiotics, spoilage, and pathogenic bacteria), and the qualitative and technological properties of milk were merged for the correlation and multivariate analyses to explore the relationships between bacterial and chemical or technological traits. Correlations were calculated among the metagenomic relative abundances of the selected taxa and groups, and between the metagenomic relative abundances and the qualitative and technological milk traits.
      Due to the high number and complexity of the relationships among all the traits, we used a multivariate factor analysis (FA) to summarize the interrelated measured traits in a small number of unmeasured latent independent explanatory variables (factors). Factor analysis was performed on the selected traits as follows. First, we performed Kaiser–Meyer–Olkin and Bartlett's tests, which showed that the traits were suitable for FA. The FA was carried out with Varimax rotation in the R environment (

      R Core Team. 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.

      ) using the psych package (available at CRAN, the Comprehensive R Archive Network, version 2.2.9; https://cran.r-project.org/web/packages/psych/index.html ) in 3 steps: (1) extraction of factors such that the minimum number of uncorrelated latent factors explained the greatest proportion of common variance; (2) factor rotation until each factor was defined by a few variables with high loadings; and (3) biological interpretation of the factors based on the strength of the loadings of the variables. The eigenvalues of the factors and the communalities of the variables after rotation were also determined.
      Eight latent explanatory factors were extracted from the 47 milk traits selected (17 bacterial and 30 qualitative or technological traits). To better understand their meaning, the signs of all the loadings of the second, sixth, and eighth factors were inverted. The scores of each milk sample for each factor were analyzed using the same linear mixed model as that used for the metagenomic relative abundances.

      RESULTS

      Factors of Variation in Milk Metagenomic Relative Abundances

      As a first step in examining the factors of variation in the relative abundances of the bacterial taxa of the 4 designated groups, we present in Figure 1 a summary of the percentages of total variance represented by the combined effects of group of cows (low vs. high) and of month of sampling (June to October; dark blue), the effects of individual cows within group (red), and the residual sources of variations (light blue). It is clear that the effects of the group × month interaction represent a major source of variation (>45% of total variance) in the relative abundances of about two-thirds of the bacterial traits, with the exception of Leuconostoc, Enterococcus, other LAB, Clostridiales, pathogenic bacteria, and Staphylococcus. The variability due to individual cows was negligible for 10 out of 17 traits, very important for Staphylococcus and Clostridiales taxa, and moderate for Enterococcus, other LAB, Kocuria, pathogenic bacteria and Enterobacteriaceae.
      Figure thumbnail gr1
      Figure 1Sources of the variation (expressed as percentage of total variance) in individual milk bacterial taxa relative abundances and their categories (in bold): effects of the month × group interaction (dark blue), individual cow within group (red), and residual variability (light blue). LAB = lactic acid bacteria.

      Combined Effects of Group of Cow and Month of Sampling on Milk Metagenomic Relative Abundances

      Table 1 shows the levels of statistical significance of the combined effects of group and month of sampling for the relative abundances of the 44 bacterial taxa identified, and their sums in categories defined by their prevalent activity. The month × group interaction exerted a significant effect on the relative abundances of LAB, other probiotics, and spoilage bacteria, and of all the taxa within these categories except for the Leuconostoc taxa. The pathogenic bacteria group was not significantly affected, although the Staphylococcus and Enterobacteriaceae taxa within this group were; these 2 taxa went in opposite directions. The relative abundances of the 31 bacterial taxa belonging to neither the desired nor the undesired groups were significantly affected in fewer than half of cases.
      Table 1Descriptive statistics, significance levels of the month × group interaction, and root mean square error (RMSE) of the log10 relative abundances of milk bacterial taxa with known dairy (LAB), other probiotic, spoilage, and pathogenic activities, and of other bacteria found in the milk
      TraitSamples, NDescriptive statisticsMonth × group, F-valueRMSE
      Mean±SD
      Lactic acid bacteria (LAB)570.7800.52810.0
      P < 0.001.
      0.335
      Lactobacillus570.2770.45225.0
      P < 0.001.
      0.200
      Leuconostoc570.0820.1651.50.158
      Lactococcus570.2940.48720.6
      P < 0.001.
      0.237
      Enterococcus580.3840.3892.8
      P < 0.05
      0.335
        Other LAB580.1910.2883.9
      P < 0.01
      0.214
      Other probiotics570.2870.42615.8
      P < 0.001.
      0.230
      Propionibacterium550.0590.14010.6
      P < 0.001.
      0.086
      Bifidobacterium570.2450.39914.2
      P < 0.001.
      0.224
      Spoilage bacteria581.1820.58512.9
      P < 0.001.
      0.342
      Clostridiales580.3080.3462.4
      P < 0.05
      0.258
      Pseudomonas580.7450.6315.4
      P < 0.001.
      0.482
      Kocuria580.6590.6547.7
      P < 0.001.
      0.425
      Alicyclobacillus580.1110.1786.6
      P < 0.001.
      0.129
      Pathogenic bacteria580.6560.5221.80.447
      Staphylococcus570.4040.4942.4
      P < 0.05
      0.341
      Enterobacteriaceae570.2290.39511.9
      P < 0.001.
      0.220
      Other bacteria
      Actinomyces560.0190.0963.0
      P < 0.01
      0.083
      Corynebacterium550.0160.0710.90.072
      Rhodococcus570.2150.3142.9
      P < 0.05
      0.250
        Other Actinomycetales560.1530.2434.1
      P < 0.001.
      0.197
      Bacteroidales570.2830.2795.0
      P < 0.001.
      0.207
      Flavobacterium570.0940.1892.2
      P < 0.05
      0.173
      Chryseobacterium580.5360.5153.8
      P < 0.01
      0.378
        Other Flavobacteriales570.0940.1441.40.139
      Sphingobacterium570.1520.1953.6
      P < 0.01
      0.154
        Other Sphingobacteriales570.0280.0771.70.073
      Chitinophagaceae570.0150.0611.00.056
      Solibacillus570.0550.1631.80.153
      Jeotgalicoccus560.0410.1221.90.103
      Exiguobacterium570.0130.0710.90.071
        Other Bacillales570.1120.2833.2
      P < 0.01
      0.220
      Aerococcaceae580.1480.2863.1
      P < 0.01
      0.229
      Carnobacteriaceae560.0450.1511.80.128
        Other Firmicutes560.0040.0201.00.020
      Ochrobactrum570.0680.1593.2
      P < 0.01
      0.118
      Paracoccus570.0800.1947.9
      P < 0.001.
      0.132
      Sphingomonas570.0500.0932.6
      P < 0.05
      0.080
        Other Alphaproteobacteria580.2760.3131.60.278
      Delftia580.1480.2301.70.213
        Other Betaproteobacteria570.1770.2251.00.225
      Deltaproteobacteria550.0010.0070.90.007
      Epsilonproteobacteria560.0150.0420.90.042
      Ruminobacter550.0220.0660.80.066
      Acinetobacter580.7510.5076.7
      P < 0.001.
      0.326
      Enhydrobacter560.0730.1690.70.174
      Xanthomonadaceae580.5990.4654.8
      P < 0.001.
      0.363
        Other Gammaproteobacteria570.1870.2801.90.260
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      Figure 2 shows the plots of the relative abundances of the sum of all the milk LAB taxa having a desired effect on cheese-making, and of the individual taxa (Lactobacillus, Leuconostoc, Lactococcus, Enterococcus, and other LAB). In the first plot, we can see that the relative abundance of LAB exhibits a significant quadratic pattern in the low group during the 5 mo of the experiment, with the lowest value occurring in August. In June, the difference between the high and the low groups of cows was not significant, which is expected because all the cows were housed and fed together indoors on the permanent lowland farm. In contrast, during the 3 summer months (July, August, and September), when the high group was on the summer highland pastures, their values were always significantly higher than the low group. At the end of summer, when the high group returned to join the low group on the permanent farm in the valley, the 2 groups showed no significant differences, indicating the absence of carryover effects of summer transhumance.
      Figure thumbnail gr2
      Figure 2Patterns of the relative abundances of milk lactic acid bacteria (LAB) taxa having desired dairy characteristics (the “LAB” category and individual bacterial taxa) during the experimental period. Blue circles represent LSM of the cows kept solely indoors, green triangles represent the cows moved to summer highland pastures, and blue triangles represent the latter cows when indoors before and after the summer transhumance. Bars represent SE of estimates. Lines and curves represent significant linear, quadratic, or cubic patterns, with their R2 values, for cows kept solely indoors. Asterisks indicate the significance levels of the differences between the 2 groups in each month (*P < 0.05; **P < 0.01; ***P < 0.001).
      Lactobacillus, Lactococcus, and Enterococcus showed a trend very similar to the LAB category, probably because of their high relative abundances during the summer highland grazing period. The situation is very different for Leuconostoc and the other LAB taxa: over the 5 mo of the experiment, the low group exhibited a linear decreasing pattern in the case of Leuconostoc and a linear increasing pattern for other LAB.
      Other probiotics followed a cubic pattern (Figure 3) for the low group, with the highest relative abundances in August and September, mainly due to the pattern of Bifidobacterium. Throughout the study period, the high group exhibited higher abundances than the low group, with the difference increasing month by month. Again, as for LAB, no carryover effect was observed after the high cows returned to the indoor permanent farm (October).
      Figure thumbnail gr3
      Figure 3Patterns of the relative abundances of the other milk bacterial taxa having desired probiotic characteristics (the “other probiotics” category and individual bacteria taxa) during the experimental period. Blue circles represent LSM of the cows kept solely indoors, green triangles represent the cows moved to summer highland pastures, and blue triangles represent the latter cows when indoors before and after summer transhumance. Bars represent SE of the estimates. Lines and curves represent significant linear, quadratic, or cubic patterns, with their R2 values, for the cows kept solely indoors. Asterisks indicate the significance levels of the differences between the 2 groups in each month (*P < 0.05; **P < 0.01; ***P < 0.001).
      Spoilage bacteria showed a complex pattern (Figure 4) and some seasonal variation in the low group according to bacterial taxa. With the exception of Kocuria, the relative abundances of all the spoilage bacteria taxa were lower in the high than in the low group during summer pasturing. This difference was significant in July and August for Alicyclobacillus, in August for Clostridiales, in September for Pseudomonas, and in August and September for the whole group. No carryover effect was observed in October.
      Figure thumbnail gr4
      Figure 4Patterns of the relative abundances of the milk bacterial taxa having undesired spoilage characteristics (the “spoilage bacteria” category and individual bacteria taxa) during the experimental period. Blue circles represent LSM of the cows kept solely indoors, green triangles represent the cows moved to summer highland pastures, and blue triangles represent the latter cows when indoors before and after the summer transhumance. Bars represent SE of estimates. Lines and curves represent significant linear, quadratic, or cubic patterns, with their R2 values, for cows kept solely indoors. Asterisks indicate the significance levels of the differences between the 2 groups in each month (*P < 0.05; **P < 0.01; ***P < 0.001).
      Finally, the relative abundances of the pathogenic bacteria (Figure 5) followed a cubic pattern for the low group of cows (with the highest value in July), due to the pattern of Staphylococcus. The differences between the high and low groups were not significant for these bacteria because of the opposite patterns in Staphylococcus and Enterobacteriaceae: the high group had a lower relative abundance of Staphylococcus than the low group and a higher relative abundance of Enterobacteriaceae during summer pasturing. (No differences were observed before and after transhumance.)
      Figure thumbnail gr5
      Figure 5Patterns of the relative abundances of the milk bacterial taxa having undesired pathogenic characteristics (the “pathogenic bacteria” category and individual bacteria taxa) during the experimental period. Blue circles represent LSM of the cows kept solely indoors, green triangles represent the cows moved to summer highland pastures, and blue triangles represent the latter cows when kept indoors before and after the summer transhumance. Bars represent SE of estimates. Lines and curves represent significant linear, quadratic, or cubic patterns, with their R2 values, for cows kept solely indoors. Asterisks indicate significance levels of the differences between the 2 groups in each month (*P < 0.05; **P < 0.01; ***P < 0.001).
      The results of the mixed-model ANOVA of the other 31 bacterial taxa detected in milk are summarized in Table 2. Regarding the seasonal pattern of the low cows, only 12 of the 31 taxa exhibited a significant trend: Jeotgalicoccus, Aerococcaceae, and Acinetobacter showed a linear increase over time; Rhodococcus and Delftia showed a quadratic pattern with a maximum during summer; Bacteroidales showed a quadratic pattern with a minimum during summer; Chitiniphagaceae, Solibacillus, other Bacillales, Ochrobactrum, Sphingomonas, and other Alphaproteobacteria followed a cubic pattern with the maximum value in July and the minimum in September. The other 19 taxa showed no significant variation over time.
      Table 2Significance levels and order and shape of the patterns observed for microbial abundances of bacterial taxa with no specific dairy, probiotic, spoilage, or pathogenic activity in milk samples collected from cows kept permanently indoors, and differences between these and cows moved to summer highland pastures in the months before (June), during (July, August, and September), and after (October) transhumance, in terms of log10 relative abundance
      ItemPattern for indoor cows
      L: up = linear pattern increasing from June to October; Q: up-down = zenithal quadratic pattern rising to a maximum during summer and then decreasing; Q: down-up = nadir quadratic pattern decreasing to a minimum during summer and then increasing; C: up-down-up = cubic pattern rising to a maximum in July, decreasing to a minimum in September, and then increasing again. Dashes indicate absence of a significant pattern for indoor cows.
      Difference between transhumant vs. indoor cows
      P-valueOrder: shapeJuneJulyAugustSeptemberOctober
      ActinomycesNS0.000.010.000.20
      P < 0.001.
      0.00
      CorynebacteriumNS−0.010.010.000.090.04
      Rhodococcus0.005Q: up-down−0.01−0.34
      P < 0.05
      −0.36
      P < 0.05
      0.12−0.10
      Other ActinomycetalesNS0.110.20−0.070.51
      P < 0.001.
      −0.08
      Bacteroidales0.017Q: down-up−0.36
      P < 0.05
      −0.19−0.260.19−0.04
      FlavobacteriumNS0.000.070.03−0.080.03
      ChryseobacteriumNS0.020.71
      P < 0.01
      0.520.98
      P < 0.001.
      −0.23
      Other FlavobacterialesNS−0.130.08−0.130.01−0.17
      SphingobacteriumNS−0.07−0.15−0.05−0.33
      P < 0.01
      0.11
      Other SphingobacterialesNS−0.020.10
      P < 0.05
      0.000.07−0.02
      Chitinophagaceae0.046C: up-down-up−0.01−0.020.000.00−0.03
      Solibacillus0.015C: up-down-up−0.11−0.18
      P < 0.05
      0.000.00−0.09
      Jeotgalicoccus0.002L: up0.05−0.010.00−0.15
      P < 0.05
      −0.11
      ExiguobacteriumNS0.000.000.000.00−0.02
      Other Bacillales<0.001C: up-down-up0.02−0.47
      P < 0.01
      −0.02−0.06−0.27
      Aerococcaceae0.001L: up0.10−0.13−0.09−0.32
      P < 0.05
      −0.04
      CarnobacteriaceaeNS−0.060.000.01−0.100.13
      Other FirmicutesNS0.000.00−0.010.00−0.03
      P < 0.05
      Ochrobactrum0.04C: up-down-up−0.03−0.29
      P < 0.01
      −0.100.010.02
      ParacoccusNS0.000.52
      P < 0.001.
      0.140.17
      P < 0.05
      0.01
      Sphingomonas<0.001C: up-down-up−0.01−0.15
      P < 0.01
      −0.05−0.03−0.06
      Other Alphaproteobacteria0.031C: up-down-up0.08−0.02−0.43
      P < 0.05
      0.23−0.20
      Delftia0.033Q: up-down−0.09−0.22−0.23−0.27
      P < 0.05
      0.14
      Other BetaproteobacteriaNS0.160.11−0.180.070.03
      DeltaproteobacteriaNS0.00−0.01
      P < 0.05
      0.000.000.00
      EpsilonproteobacteriaNS0.000.000.010.02−0.02
      RuminobacterNS−0.04−0.05−0.030.000.00
      Acinetobacter0.003L: up−0.100.210.81
      P < 0.01
      −0.280.72
      EnhydrobacterNS0.07−0.090.00−0.07−0.11
      XanthomonadaceaeNS0.08−0.57
      P < 0.05
      −0.28−0.74
      P < 0.001.
      0.45
      P < 0.05
      Other GammaproteobacteriaNS0.040.30−0.050.43
      P < 0.01
      0.03
      1 L: up = linear pattern increasing from June to October; Q: up-down = zenithal quadratic pattern rising to a maximum during summer and then decreasing; Q: down-up = nadir quadratic pattern decreasing to a minimum during summer and then increasing; C: up-down-up = cubic pattern rising to a maximum in July, decreasing to a minimum in September, and then increasing again. Dashes indicate absence of a significant pattern for indoor cows.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      Comparing the high and low cows within each month, we found an unexpected significant difference for Bacteroidales in June. We detected some difference in 19 of the 31 taxa during the 3 summer months (i.e., when the high group was on highland pastures while the low group remained indoors): in 10 in July, 3 in August, and 10 in September. Higher values of 9 taxa were observed in the high group, and 14 taxa in the low group. Only other Firmicutes and Xanthomonadaceae showed significant carryover effects in October.

      Correlations

      Correlation analyses were carried out on the relative abundances of the bacterial taxa known to have some specific activity in the 4 designated categories, and their Pearson correlations are summarized in a heat plot in Figure 6. The relative abundances of the individual taxa, and the LAB and other probiotic categories, were generally positively correlated with each other, with the exception of the Leuconostoc taxa, which is almost entirely independent of all the other taxa and categories.
      Figure thumbnail gr6
      Figure 6Heat plot of the correlations among the bacterial traits included in the factor analysis. LAB = lactic acid bacteria.
      The spoilage bacteria taxa exhibited low correlations with each other and negative correlations with the relative abundances of the taxa of the LAB and other probiotics categories (Figure 6). The 2 taxa included in the pathogenic bacteria category were negatively correlated with each other, and their correlations with other bacterial categories were in opposite directions: Enterobacteriaceae were positively correlated with the LAB and other probiotics categories and taxa, and Staphylococcus had low correlations with spoilage bacteria taxa.
      The correlations among the constituents and technological properties of milk are not among the objectives of this study, as they are already well known; therefore, they are not illustrated in detail and discussed here.
      The correlations between the bacterial taxa and the constituents and technological traits of milk are summarized in a heat plot in Figure 7. These correlations vary greatly according to taxa and milk trait. It is worth noting that the LAB and other probiotics categories and their individual taxa exhibited correlations with many of the 30 milk composition and technological traits that were often in opposite directions to those shown by the spoilage bacteria category and individual taxa.
      Figure thumbnail gr7
      Figure 7Heat plot of the correlations between the bacterial and chemical-technological traits included in the factor analysis. LAB = lactic acid bacteria; RCT = rennet coagulation time; k20 = curd-firming time; a30, a45, a60 = curd firmness 30, 45, and 60 min after rennet addition, respectively; RCTeq = rennet coagulation time by equation; kCF = curd-firming instant rate constant; kSR = syneresis instant rate constant; CFmax = maximum curd firmness; tmax = time to reach CFmax; %CY = percentage cheese yields; REC = nutrient recovery traits.

      Latent Explanatory Factors of the Bacterial, Compositional, and Technological Traits of Milk

      The multivariate FA carried out on the 47 selected bacterial, compositional, and technological milk traits identified 8 latent explanatory factors. The loadings of each factor, excluding those that were nonrelevant (<0.30), are reported in Table 3, with high loadings (>0.50) indicated by asterisks. All 17 bacterial taxa were included in one (n = 6) or more (n = 11) factors. Leuconostoc, other LAB, Clostridiales, Alicyclobacillus, and Enterobacteriaceae had the lowest loading values (never higher than 0.44) and therefore are not well represented by any latent explanatory factor. In contrast, all 4 bacterial categories presented high communality loading values (0.68 to 0.90).
      Table 3Loadings of the latent explanatory factors of the relative abundances of milk bacterial taxa, milk and whey constituents, milk coagulation, curd-firming, and cheese yield
      Item
      LAB = lactic acid bacteria; RCT = rennet coagulation time; k20 = curd-firming time; a30, a45, a60 = curd firmness 30, 45, and 60 min after rennet addition, respectively; RCTeq = rennet coagulation time by equation; kCF = curd-firming instant rate constant; kSR = syneresis instant rate constant; CFmax = maximum curd firmness; tmax = time to reach CFmax; %CY = percentage cheese yields; REC = nutrient recovery traits.
      Factor 1: GelationFactor 2: Pro-dairyFactor 3: Cheese yieldFactor 4: Udder healthFactor 5: CaseinsFactor 6: PathogensFactor 7: CurdlingFactor 8: Fat recoveryCommunality
      Milk bacterial group
       LAB0.63
      High loading, >0.50.
      −0.390.68
      Lactobacillus0.83
      High loading, >0.50.
      0.82
      Leuconostoc−0.40−0.330.37
      Lactococcus0.65
      High loading, >0.50.
      0.61
      Enterococcus−0.53
      High loading, >0.50.
      −0.370.48
      Other LAB0.330.310.330.44
       Other probiotics0.91
      High loading, >0.50.
      >0.90
      Propionibacterium0.52
      High loading, >0.50.
      −0.460.55
      Bifidobacterium0.92
      High loading, >0.50.
      >0.90
       Spoilage bacteria−0.79
      High loading, >0.50.
      −0.420.85
      Clostridiales−0.310.24
      Pseudomonas−0.57
      High loading, >0.50.
      −0.310.43
      Kocuria−0.48−0.52
      High loading, >0.50.
      0.57
      Alicyclobacillus0.380.410.36
       Pathogenic bacteria0.80
      High loading, >0.50.
      0.74
      Staphylococcus−0.350.78
      High loading, >0.50.
      0.77
      Enterobacteriaceae0.48−0.360.43
      Milk technological trait
       Milk composition
      TS0.89
      High loading, >0.50.
      >0.90
      Milk fat0.95
      High loading, >0.50.
      >0.90
      Nonfat solids−0.360.52
      High loading, >0.50.
      0.70
      High loading, >0.50.
      >0.90
      Milk protein0.35−0.320.81
      High loading, >0.50.
      >0.90
      Fat/protein ratio0.84
      High loading, >0.50.
      −0.44>0.90
      Milk casein0.300.86
      High loading, >0.50.
      >0.90
      Casein number−0.380.330.53
      High loading, >0.50.
      0.61
      Milk urea−0.430.43
       Udder health trait
      SCS−0.400.25
      Milk lactose−0.310.75
      High loading, >0.50.
      >0.90
       Coagulation property
      RCT−0.96
      High loading, >0.50.
      >0.90
      k20−0.47−0.31−0.56
      High loading, >0.50.
      0.69
      a300.87
      High loading, >0.50.
      0.31>0.90
      a450.79
      High loading, >0.50.
      >0.90
      a600.44−0.360.41−0.59
      High loading, >0.50.
      >0.90
       Curd-firming modeling
      RCTeq−0.96
      High loading, >0.50.
      >0.90
      kCF0.90
      High loading, >0.50.
      >0.90
      kSR0.67
      High loading, >0.50.
      0.480.84
      CFmax0.81
      High loading, >0.50.
      0.38>0.90
      tmax−0.75
      High loading, >0.50.
      −0.54
      High loading, >0.50.
      >0.90
       Cheese yield
      %CYCURD0.58
      High loading, >0.50.
      0.50
      High loading, >0.50.
      0.64
      %CYSOLIDS0.85
      High loading, >0.50.
      0.31>0.90
      RECFAT0.87
      High loading, >0.50.
      >0.90
      RECPROTEIN0.62
      High loading, >0.50.
      0.45
      RECSOLIDS0.65
      High loading, >0.50.
      −0.56
      High loading, >0.50.
      >0.90
      RECENERGY0.68
      High loading, >0.50.
      −0.380.47>0.90
       Whey composition
      Whey TS−0.350.75
      High loading, >0.50.
      −0.36>0.90
      Whey fat0.47−0.83
      High loading, >0.50.
      >0.90
      Whey protein0.39−0.410.340.57
      Whey lactose−0.360.88
      High loading, >0.50.
      >0.90
      1 LAB = lactic acid bacteria; RCT = rennet coagulation time; k20 = curd-firming time; a30, a45, a60 = curd firmness 30, 45, and 60 min after rennet addition, respectively; RCTeq = rennet coagulation time by equation; kCF = curd-firming instant rate constant; kSR = syneresis instant rate constant; CFmax = maximum curd firmness; tmax = time to reach CFmax; %CY = percentage cheese yields; REC = nutrient recovery traits.
      * High loading, >0.50.
      The 30 compositional and technological traits of milk, with the exception of milk urea, SCS, and whey protein, often contributed to characterizing the factors and presented high communality values.
      The 8 latent explanatory factors represented 75.3% of the total covariance of the whole matrix, with individual values ranging from 14.2% for the first factor to 5.3% for the eighth (Table 4). It is worth noting that the mixed-model ANOVA of the scores of each factor revealed that all the latent explanatory factors, except factor 7, were significantly affected by the combined effect of group of cows and month of sampling (Table 4).
      Table 4Descriptive statistics, significance levels of the group × month interaction, and root mean square error (RMSE) of the latent explanatory factors of milk bacterial taxa, milk and whey constituents, milk coagulation and curd-firming properties, and cheese yield
      Latent explanatory factorExplained variance (%)Group × month F-valueRMSE
      By each factorCumulative
      Factor 1: Gelation14.214.27.3
      P < 0.001.
      0.49
      Factor 2: Pro-dairy bacteria14.028.223.1
      P < 0.001.
      0.46
      Factor 3: Cheese yield11.739.92.4
      P < 0.05
      0.84
      Factor 4: Udder health9.449.313.0
      P < 0.001.
      0.44
      Factor 5: Casein9.258.56.5
      P < 0.001.
      0.55
      Factor 6: Pathogenic bacteria6.064.62.8
      P < 0.05
      0.69
      Factor 7: Curdling5.570.11.80.58
      Factor 8: Fat recovery5.375.33.5
      P < 0.01
      0.76
      * P < 0.05
      ** P < 0.01
      *** P < 0.001.
      The 10 least squares means values, the standard errors, the seasonal patterns of low cows, and the significance levels of the differences between the 2 groups of cows within each month are illustrated in Figure 8 (one plot for each of the 8 latent explanatory factors). Factor 1 was characterized by a quadratic trend with the maximum value in October for the cows kept solely indoors (low group), and presented no significant differences between the 2 experimental groups of cows in any month of sampling. Factor 2 presented a quadratic seasonal trend with the maximum value in August and the high group having significantly greater values during the 3 summer months. Factor 3 presented a quadratic seasonal pattern for low cows, but with the minimum in August and the high group having significantly higher values in July and August. Factor 4 presented a pattern that was almost the opposite of that of factor 3: low cows followed a cubic pattern with the maximum in August, and high cows had significantly lower values in September. Factor 5 presented a linearly increasing pattern for low cows and significantly lower values in the high group only in August. Factor 6 presented a cubic seasonal pattern with the maximum in July and no differences between the 2 groups of cows. Factor 7 presented a linearly increasing pattern for the low group and no differences between the 2 groups. Factor 8, like factor 7, presented a linearly increasing pattern for the low group and no differences between the 2 groups.
      Figure thumbnail gr8
      Figure 8Patterns of the scores of the latent explanatory factors during the experimental period. Blue circles represent LSM of the cows kept solely indoors, green triangles represent the cows moved to summer highland pastures, and blue triangles represent the latter cows when indoors before and after the summer transhumance. Bars represent SE of estimates. Lines and curves represent significant linear, quadratic, or cubic patterns, with their R2 values, for the cows kept solely indoors. Asterisks indicate the significance levels of the differences between the 2 groups in each month (*P < 0.05; **P < 0.01; ***P < 0.001).

      DISCUSSION

      Effects of Group of Cows and Individual Cows

      In this study, we compared 2 groups of cows kept in the same physiological, environmental, nutritional, and management conditions only in June, when they were all kept indoors on the lowland permanent farm. At P < 0.05 there is a 5% probability that the differences between the high and low are due to chance, and at P < 0.01 the probability is 1%. Of the 56 contrasts tested (the relative abundances of 44 taxa and 4 categories of bacteria shown in Figures 2, 3, 4, and 5 and Table 2, and the scores of 8 explanatory latent factors shown in Figure 8), 3 were significant at P < 0.05 (Bacteroidales taxa, factor 3, factor 5) and 1 at P < 0.01 (Pseudomonas taxa), so we can assume the 2 groups were homogeneous. No information is available in the literature on the variability in bacterial traits in different, randomly composed groups of cows. The increasing variability among groups of cows may be due to permanent differences among different cows observed in subsequent samplings (animal effect). As seen in Figure 1, the permanent animal effect is generally small or not observable for the majority of bacterial taxa. This means that the (high) variability observed among different milk samples is due to the effects of other common factors (environment, diet, hygiene practices, milking routine, etc.) or individual, nonpermanent factors (temporary diseases, cleanness of teat surface, feed selection, etc.;
      • Du B.
      • Meng L.
      • Liu H.
      • Zheng N.
      • Zhang Y.
      • Guo X.
      • Zhao S.
      • Li F.
      • Wang J.
      Impact of milking and housing environment on milk microbiota.
      ;
      • Parente E.
      • Ricciardi A.
      • Zotta T.
      The microbiota of dairy milk: A review.
      ). It is not surprising that the bacterial taxa with the highest animal effect was Staphylococcus—that is, a potential pathogen generally considered to be a cause of clinical mastitis (
      • Verraes C.
      • Vlaemynck G.
      • Van Weyenberg S.
      • De Zutter L.
      • Daube G.
      • Sindic M.
      • Uyttendaele M.
      • Herman L.
      A review of the microbiological hazards of dairy products made from raw milk.
      ;
      • Bobbo T.
      • Ruegg P.L.
      • Stocco G.
      • Fiore E.
      • Gianesella M.
      • Morgante M.
      • Pasotto D.
      • Bittante G.
      • Cecchinato A.
      Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
      ;
      • Keane O.M.
      Symposium review: Intramammary infections—Major pathogens and strain-associated complexity.
      ). It is worth noting that the cows selected for this study were all healthy, and none of them developed clinical mastitis during the study. On the other side, single episodes of clinical mastitis do not influence the (permanent) animal effect, only the residual variance, and repeated cases of clinical mastitis normally result in the cows being culled. The high animal effect observed in Figure 1 for Staphylococcus could be interpreted as the predisposition of some healthy cows to carry greater or lesser quantities of these bacteria. Whether this indicates a predisposition for subclinical mastitis is not known, and is an issue worth investigating with a much larger number of animals. Cows have a (modest) heritability for both the incidence of clinical mastitis (
      • Koeck A.
      • Loker S.
      • Miglior F.
      • Kelton D.F.
      • Jamrozik J.
      • Schenkel F.S.
      Genetic relationships of clinical mastitis, cystic ovaries, and lameness with milk yield and somatic cell score in first-lactation Canadian Holsteins.
      ) and the level of the mastitis indicator represented by the somatic cell content in the milk (
      • Urioste J.I.
      • Franzén J.
      • Strandberg E.
      Phenotypic and genetic characterization of novel somatic cell count traits from weekly or monthly observations.
      ;
      • Pegolo S.
      • Giannuzzi D.
      • Bisutti V.
      • Tessari R.
      • Gelain M.E.
      • Gallo L.
      • Schiavon S.
      • Tagliapietra F.
      • Trevisi E.
      • Ajmone Marsan P.
      • Bittante G.
      • Cecchinato A.
      Associations between differential somatic cell count and milk yield, quality, and technological characteristics in Holstein cows.
      ). This indicates an interaction between the cow's genome and the infectiousness of the pathogens causing mastitis, so the variability in the relative abundances of Staphylococcus taxa among different animals observed here could also depend on their genome.
      Other bacterial taxa shown in Figure 1 exhibiting a nontrivial animal effect are other LAB, Clostridiales, and Kocuria. We found no information in the scientific literature on animal repeatability of the relative abundances of these taxa, so this, too, could be an interesting line of research with respect to cheese-making efficiency or cheese defects (
      • de Paiva Anciens Ramos G.L.
      • Vigoder H.C.
      • dos Santos Nascimento J.
      Kocuria spp. in foods: Biotechnological uses and risks for food safety.
      ).

      Associations Between Milk Microbiota, Milk Composition, and Cheese-Making Properties

      Association studies between milk microbiota and composition and cheese-making properties are infrequent and deal mainly with the potential effects of pathogenic bacteria on udder health, and mastitis in particular (
      • Leitner G.
      • Krifucks O.
      • Merin U.
      • Lavi Y.
      • Silanikove N.
      Interactions between bacteria type, proteolysis of casein and physico-chemical properties of bovine milk.
      ;
      • Bobbo T.
      • Ruegg P.L.
      • Stocco G.
      • Fiore E.
      • Gianesella M.
      • Morgante M.
      • Pasotto D.
      • Bittante G.
      • Cecchinato A.
      Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
      ).
      Multivariate statistical analyses are often used to study associations among different bacterial taxa (
      • Rodrigues M.X.
      • Lima S.F.
      • Canniatti-Brazaca S.G.
      • Bicalho R.C.
      The microbiome of bulk tank milk: Characterization and associations with somatic cell count and bacterial count.
      ) and between these and other milk or cheese characteristics (
      • Nyman A.K.
      • Persson Waller K.
      • Bennedsgaard T.W.
      • Larsen T.
      • Emanuelson U.
      Associations of udder-health indicators with cow factors and with intramammary infection in dairy cows.
      ). The most commonly used method is principal component analysis, as it is very efficient, although the results are not always easy to interpret. Factor analysis has the advantage of better clustering the observed traits so that the latent explanatory factors can be related to a small number of important traits. Unfortunately, this method is seldom used with metagenomics data sets, and we are not aware of any FA combining the microbiological, compositional, and technological properties of milk, which means that it is not possible to compare our results with those of other authors.
      It is worth noting that in this analysis we did not obtain any latent factors based on the simultaneous strong influence of bacterial and other milk characteristics. This means that milk characteristics do not seem to be highly dependent on microbiological populations, nor vice versa. Two of the 8 factors were based mainly on bacterial taxa, and the other 6 on milk composition and cheese-making properties.
      As seen in Table 4, the first latent explanatory factor (14.2% of total variation) is based on traits obtained mainly from lactodynamographic tests and modeling, and particularly, with negative loadings, on the time from rennet addition to coagulation (RCT and RCTeq), and the time to reach a given (k20) or maximum (tmax) curd firmness. Early gelation is obviously positively correlated with an increase in the traits measuring curd firmness (a30, a45, a60, and CFmax) and allows more time for estimating curd syneresis (kSR; Table 3). This is why we named this factor the “milk gelation factor.” It also includes the protein and casein contents of milk and whey, with positive loadings; proteins, especially caseins, are known to have a favorable effect on the rapidity and intensity of coagulation (
      • Amalfitano N.
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Malacarne M.
      • Summer A.
      • Bittante G.
      Milk protein fractions strongly affect the patterns of coagulation, curd firming, and syneresis.
      ). Among the bacterial traits, only 2 LAB taxa are involved in the gelation factor, but with opposite signs: negative in the case of Leuconostoc, positive in the case of other LAB. Finally, milk urea is also included with a negative loading.
      The second most important latent factor obtained from the analysis (14.0% of total variation), shown in Table 4, is substantially based on milk metagenomics, as it includes, with positive loadings, the 2 categories (and their major taxa) that have putative positive effects on the commercial and nutritional value of milk: LAB and other probiotics. However, it also includes the spoilage bacteria taxa, with a negative loading, and spoilage microorganisms are well known to negatively influence the value and quality of dairy products (
      • Martin N.H.
      • Torres-Frenzel P.
      • Wiedmann M.
      Invited review: Controlling dairy product spoilage to reduce food loss and waste.
      ). The pathogenic bacteria category is not included in this factor because of the opposite sign of the loadings of the 2 taxa it includes. It is evident that this important latent factor can be considered an index of the favorableness of the microbiological profile of milk, and for this reason we named it “pro-dairy bacteria.” This factor also includes, with a negative sign, some traits related to milk and whey composition, but none of these characterize the latent factor (Table 3).
      The third latent explanatory factor in Table 4 (11.7% of total variation) is based on milk fat content and cheese yield. It is worth noting that fat is a major component of curd solids but is, in particular, the component with the highest variability, much greater than that of protein. This explains why milk fat is so closely associated with milk solids (and the fat/protein ratio) and with cheese yield, expressed as cheese solids as a percentage of the TS of the processed milk. This last trait is obviously associated with the recovery of milk solids and energy in the curd, and also with water retained in the curd and fresh cheese yield. This is why in Table 3 we named this latent explanatory factor the “cheese yield factor.” It is also associated with, but not characterized by, whey fat, because with the increasing fat content of milk we expect an increase in both fat retained in the curd and fat lost in the whey, although proportionally in favor of the former due to the positive relationships between the fat content of milk and fat recovery in the curd. The other traits (negatively) associated with the cheese yield factor were casein number and milk lactose content. Because lactose is the major component of milk solids and ends up being mainly lost with the whey, it is evident that, as it increases in milk, it causes a reduction in the yield of cheese solids. The meaning of the negative loading of casein number in this factor is less evident, but we should bear in mind that this trait is at the same time also included in 2 other factors, discussed later. A cheese yield factor was also identified in a previous large data set that did not include metagenomic information (
      • Dadousis C.
      • Cipolat-Gotet C.
      • Schiavon S.
      • Bittante G.
      • Cecchinato A.
      Inferring individual cow effects, dairy system effects and feeding effects on latent variables underlying milk protein composition and cheese-making traits in dairy cattle.
      ). It was the most important factor in that study, representing 14% of all variation and, as in this study, was characterized by yield of cheese solids, fat content, and recovery of milk energy in the curd, as well as by protein content. A cheese yield latent factor was also identified from analysis of a large data set obtained from observations on dairy ewes (
      • Manca M.G.
      • Serdino J.
      • Gaspa G.
      • Urgeghe P.
      • Ibba I.
      • Contu M.
      • Fresi P.
      • Macciotta N.P.P.
      Derivation of multivariate indices of milk composition, coagulation properties, and individual cheese yield in dairy sheep.
      ).
      The fourth factor presented in Table 3 is characterized by the lactose content of milk and whey, milk nonfat solids, and whey solids (lactose being the major constituent of both of the latter). As expected, because lactose is retained in very small proportions in the curd, this factor included recovery of milk solids and energy in the curd, with negative loadings. An increase in lactose content is frequently associated with a lower somatic cell content (see the negative loading of SCS) and whey proteins (see the positive loading of casein number), and, taken together, these are interpreted as indicators of good udder health (
      • Macciotta N.P.P.
      • Cecchinato A.
      • Mele M.
      • Bittante G.
      Use of multivariate factor analysis to define new indicator variables for milk composition and coagulation properties in Brown Swiss cows.
      ). This is why we named this the “udder health factor,” as others have done for cattle (
      • Cecchinato A.
      • Ribeca C.
      • Maurmayr A.
      • Penasa M.
      • De Marchi M.
      • Macciotta N.P.P.
      • Mele M.
      • Secchiari P.
      • Pagnacco G.
      • Bittante G.
      Short communication: Effects of β-lactoglobulin, stearoyl-coenzyme A desaturase 1, and sterol regulatory element binding protein gene allelic variants on milk production, composition, acidity, and coagulation properties of Brown Swiss cows.
      ;
      • Macciotta N.P.P.
      • Cecchinato A.
      • Mele M.
      • Bittante G.
      Use of multivariate factor analysis to define new indicator variables for milk composition and coagulation properties in Brown Swiss cows.
      ;
      • Dadousis C.
      • Cipolat-Gotet C.
      • Bittante G.
      • Cecchinato A.
      Inferring genetic parameters on latent variables underlying milk yield and quality, protein composition, curd firmness and cheese-making traits in dairy cattle.
      ,
      • Dadousis C.
      • Cipolat-Gotet C.
      • Schiavon S.
      • Bittante G.
      • Cecchinato A.
      Inferring individual cow effects, dairy system effects and feeding effects on latent variables underlying milk protein composition and cheese-making traits in dairy cattle.
      ;
      • Cecchinato A.
      • Macciotta N.P.P.
      • Mele M.
      • Tagliapietra F.
      • Schiavon S.
      • Bittante G.
      • Pegolo S.
      Genetic and genomic analyses of latent variables related to the milk fatty acid profile, milk composition, and udder health in dairy cattle.
      ), as well as for goats and sheep (
      • Manca M.G.
      • Serdino J.
      • Gaspa G.
      • Urgeghe P.
      • Ibba I.
      • Contu M.
      • Fresi P.
      • Macciotta N.P.P.
      Derivation of multivariate indices of milk composition, coagulation properties, and individual cheese yield in dairy sheep.
      ;
      • Vacca G.M.
      • Paschino P.
      • Dettori M.L.
      • Bergamaschi M.
      • Cipolat-Gotet C.
      • Bittante G.
      • Pazzola M.
      Environmental, morphological, and productive characterization of Sardinian goats and use of latent explanatory factors for population analysis.
      ). Note that the udder health factor is associated with variations in 5 bacterial taxa: the decreases in Enterococcus, Propionibacterium, and Enterobacteriaceae, and the increases in other LAB and Alicyclobacillus. The meaning of these associations is clear for Enterobacteriaceae but not for the other taxa.
      In a previous study we found relationships between the presence of a few bacterial groups causing mastitis and the qualitative and technological properties of milk (
      • Bobbo T.
      • Ruegg P.L.
      • Stocco G.
      • Fiore E.
      • Gianesella M.
      • Morgante M.
      • Pasotto D.
      • Bittante G.
      • Cecchinato A.
      Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
      ). Staphylococcus aureus was the contagious bacteria most frequently isolated in milk samples from individual cows. The presence of this pathogen was associated with decreases in daily milk yield, casein number, and milk lactose (i.e., the udder health factor), but not in other milk constituents. The contagious milk bacteria were also associated with a worsening of milk coagulation and curd firmness properties, and a decrease in cheese-making efficiency and, in particular, recovery of milk fat and protein in the curd (
      • Bobbo T.
      • Ruegg P.L.
      • Stocco G.
      • Fiore E.
      • Gianesella M.
      • Morgante M.
      • Pasotto D.
      • Bittante G.
      • Cecchinato A.
      Associations between pathogen-specific cases of subclinical mastitis and milk yield, quality, protein composition, and cheese-making traits in dairy cows.
      ). More information on the relationships between milk SCC or SCS and milk properties can be found in
      • Bobbo T.
      • Cipolat-Gotet C.
      • Bittante G.
      • Cecchinato A.
      The nonlinear effect of somatic cell count on milk composition, coagulation properties, curd firmness modeling, cheese yield, and curd nutrient recovery.
      .
      The fifth factor listed in Table 3 is characterized by milk casein, milk protein (whey protein), nonfat solids, casein number, and protein recovery in the curd. It is evident that casein content has a central role in this factor, so we named it the “casein factor.” It also includes some traits related to the rate and degree of curd firming (k20, a60, and CFmax), confirmation that caseins play an important role in curd firming, much more so than in milk coagulation time (
      • Jõudu I.
      • Henno M.
      • Kaart T.
      • Püssa T.
      • Kärt O.
      The effect of milk protein contents on the rennet coagulation properties of milk from individual dairy cows.
      ). Protein is also more correlated than fat with the retention of water in the curd, which explains the positive loading of curd cheese yield (
      • Cipolat-Gotet C.
      • Malacarne M.
      • Summer A.
      • Cecchinato A.
      • Bittante G.
      Modeling weight loss of cheese during ripening and the influence of dairy system, parity, stage of lactation, and composition of processed milk.
      ). In any case, we should bear in mind that different protein fractions have different effects on milk coagulation and curd-firming traits (
      • Amalfitano N.
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Malacarne M.
      • Summer A.
      • Bittante G.
      Milk protein fractions strongly affect the patterns of coagulation, curd firming, and syneresis.
      ), as well as cheese yield and milk nutrient recovery in cheese (
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Malacarne M.
      • Bittante G.
      • Summer A.
      Variations in milk protein fractions affect the efficiency of the cheese-making process.
      ). It is worth noting that “other LAB taxa” is the only microbiological trait positively included in this factor.
      As seen in Table 3, the sixth factor, like the pro-dairy bacteria factor, is based only on bacterial traits. It is characterized mainly by the pathogenic bacteria category, specifically by the Staphylococcus taxa, but not by Enterobacteriaceae (included negatively in the udder health factor). This is why we named it “pathogenic bacteria.” It is worth noting that mastitis has been described as a dysbiosis, an imbalance in the healthy mammary gland microbiome (
      • Andrews T.
      • Neher D.A.
      • Weicht T.R.
      • Barlow J.W.
      Mammary microbiome of lactating organic dairy cows varies by time, tissue site, and infection status.
      ). This latent factor is also associated with a negative loading of the LAB (Leuconostoc and Enterococcus taxa) and spoilage bacteria categories (particularly the Kocuria and Pseudomonas taxa), and with a positive loading of Alicyclobacillus.
      The seventh factor (Table 3) is also mainly characterized by traits obtained during the lactodynamographic test, but here the time intervals from rennet addition to coagulation are not included (earliness of coagulation, as in the case of the gelation factor), whereas the rapidity of the curd-firming (kCF, k20, and tmax) and syneresis processes (kSR) is central, and explains the positive loading of a30, the small loading of a45, and the negative loading of a60. This is why we named this the “curdling factor.” Small correlations between gelation time, rapidity of curd firming, and rapidity of syneresis were previously reported by
      • Macciotta N.P.P.
      • Cecchinato A.
      • Mele M.
      • Bittante G.
      Use of multivariate factor analysis to define new indicator variables for milk composition and coagulation properties in Brown Swiss cows.
      . It is worth noting that 2 independent latent factors representing coagulation traits and curd-firming traits were also obtained for goat milk (
      • Todaro M.
      • Scatassa M.
      • Giaccone P.
      Multivariate factor analysis of Girgentana goat milk composition.
      ).
      The eighth factor in Table 3 is characterized by positive loading of the recoveries of milk fat (RECFAT) and total milk solids (RECENERGY) in the curd, and a negative loading of the fat and TS contents in the whey, so we named this the “fat recovery” factor.

      Effects of Season or Lactation Stage, Farming System, and Pasture Carryover

      We were unable in this study to disentangle the effects of season and lactation stage, because, as is often the case in pasture-based farming systems, the cows at the beginning of pasturing are generally in the first half of lactation, and the lactation progresses together with the advancing season. Having initially created 2 homogeneous groups of cows for lactation stage, and having kept their composition constant during the experiment, we found a similar overlap in advancing season and lactation stage in the high (moved to summer pastures) and the low group (kept solely indoors). In the case of the high group, the overlap also went hand in hand with the gradual maturation of the forage on the summer pasture. This should be borne in mind when interpreting the results, as the often significant seasonal pattern of the low cows (Figures 2, 3, 4, and 5, and Table 2) represents the simultaneous change in season and lactation stage, but not feeding regimen (constant TMR).
      We were also unable to model the same pattern in the case of the high group, because we cannot assume continuous evolution along the 5 experimental months. This group, in fact, underwent 2 abrupt changes: the move from the permanent farm in the valley to the highland pastures, and then the return to indoor conditions. In this case, we considered the low group as the control, and compared the high group against them month by month. As seen before, the substantially low incidence of significant contrasts between the 2 groups in June supports the assumption of initial homogeneity of the 2 groups.
      The large number of significant differences observed in July, August, and September for the bacterial traits and the latent factors confirms the hypothesis of very large effects of farming system (permanent indoor housing vs. summer highland pasture), and is consistent with results previously obtained by the same project using bacterial culture-dependent and -independent approaches (
      • Carafa I.
      • Navarro I.C.
      • Bittante G.
      • Tagliapietra F.
      • Gallo L.
      • Tuohy K.
      • Franciosi E.
      Shift in the cow milk microbiota during alpine pasture as analyzed by culture dependent and high-throughput sequencing techniques.
      ). Other authors have observed large difference in the microbiota of milk produced by cows kept indoors and cows at pasture (
      • Bonizzi I.
      • Buffoni J.N.
      • Feligini M.
      • Enne G.
      Investigating the relationship between raw milk bacterial composition, as described by intergenic transcribed spacer-PCR fingerprinting, and pasture altitude.
      ;
      • Doyle C.J.
      • Gleeson D.
      • O’Toole P.W.
      • Cotter P.D.
      Impacts of seasonal housing and teat preparation on raw milk microbiota: A high-throughput sequencing study.
      ). The autochthonous LAB of milk produced on highland pastures are known to be important for the cheese-making process and the final quality of traditional mountain cheeses (
      • Carafa I.
      • Stocco G.
      • Franceschi P.
      • Summer A.
      • Tuohy K.M.
      • Bittante G.
      • Franciosi E.
      Evaluation of autochthonous lactic acid bacteria as starter and non-starter cultures for the production of Traditional Mountain cheese.
      ). As the metagenomics results are expressed as relative abundances (i.e., proportions among different taxa and not their population sizes), it would be useful to draw comparisons with bacterial plate counts on different selective media to obtain a clearer picture of the actual amounts of viable bacterial populations in the milk samples analyzed. In our previous study (
      • Carafa I.
      • Navarro I.C.
      • Bittante G.
      • Tagliapietra F.
      • Gallo L.
      • Tuohy K.
      • Franciosi E.
      Shift in the cow milk microbiota during alpine pasture as analyzed by culture dependent and high-throughput sequencing techniques.
      ), we found that summer highland grazing increased the counts of all bacterial categories. The increases in aerobic (+41%), anaerobic (+54%), and mesophilic lactococci (+45%) and bifidobacteria counts (+47%) from the low to the high group were similar, whereas the increases in mesophilic lactobacilli (+411%), Propionibacteria (+125%), and coliform counts (+631%) were several times greater, largely congruent with the increases in their relative abundances found here using metagenomics. The substantial increases in the relative abundances of the LAB category and its main taxa (Figure 2), and of the pro-dairy factor (Figure 8) during summer highland pasturing, confirm that the improvement in milk chemical composition observed in this project and in several other studies (
      • Bergamaschi M.
      • Cipolat-Gotet C.
      • Stocco G.
      • Valorz C.
      • Bazzoli I.
      • Sturaro E.
      • Ramanzin M.
      • Bittante G.
      Cheesemaking in highland pastures: Milk technological properties, cream, cheese and ricotta yields, milk nutrients recovery, and products composition.
      ;
      • Bergamaschi M.
      • Bittante G.
      From milk to cheese: Evolution of flavor fingerprint of milk, cream, curd, whey, ricotta, scotta, and ripened cheese obtained during summer Alpine pasture.
      ;
      • Bittante G.
      • Cecchinato A.
      • Tagliapietra F.
      • Schiavon S.
      • Toledo-Alvarado H.
      Effects of breed, farm intensiveness and cow productivity level on cheese-making ability predicted using infrared spectral data at the population level.
      ) is due to pasturing. It is well known that different farming systems, as well as individual farms (
      • Bokulich N.A.
      • Mills D.A.
      Facility-specific “house” microbiome drives microbial landscapes of artisan cheesemaking plants.
      ;
      • Skeie S.B.
      • Håland M.
      • Thorsen I.M.
      • Narvhus J.
      • Porcellato D.
      Bulk tank raw milk microbiota differs within and between farms: A moving goalpost challenging quality control.
      ;
      • Priyashantha H.
      • Lundh Å.
      Graduate Student Literature Review: Current understanding of the influence of on-farm factors on bovine raw milk and its suitability for cheesemaking.
      ) and dairy plants, affect cheese-making efficiency and product quality (
      • Falardeau J.
      • Keeney K.
      • Trmčić A.
      • Kitts D.
      • Wang S.
      Farm-to-fork profiling of bacterial communities associated with an artisan cheese production facility.
      ;
      • Nam J.H.
      • Cho Y.S.
      • Rackerby B.
      • Goddik L.
      • Park S.H.
      Shifts of microbiota during cheese production: Impact on production and quality.
      ). This also supports the claimed specificity of cheeses produced on temporary highland summer farms in Alpine regions (
      • Bittante G.
      • Cecchinato A.
      • Cologna N.
      • Penasa M.
      • Tiezzi F.
      • De Marchi M.
      Factors affecting the incidence of first-quality wheels of Trentingrana cheese.
      ,
      • Bittante G.
      • Cologna N.
      • Cecchinato A.
      • de Marchi M.
      • Penasa M.
      • Tiezzi F.
      • Endrizzi I.
      • Gasperi F.
      Monitoring of sensory attributes used in the quality payment system of Trentingrana cheese.
      ) and signals the possibility of authenticating the origin of dairy products according to farming system (
      • Bergamaschi M.
      • Cipolat-Gotet C.
      • Cecchinato A.
      • Schiavon S.
      • Bittante G.
      Chemometric authentication of farming systems of origin of food (milk and ripened cheese) using infrared spectra, fatty acid profiles, flavor fingerprints, and sensory descriptions.
      ).
      It is worth noting that, in a previous large survey, a significant effect of dairy system was found for latent factors named “cheese yield,” “udder health,” and “yield” (
      • Dadousis C.
      • Cipolat-Gotet C.
      • Schiavon S.
      • Bittante G.
      • Cecchinato A.
      Inferring individual cow effects, dairy system effects and feeding effects on latent variables underlying milk protein composition and cheese-making traits in dairy cattle.
      ; production traits were not included in our study). The effect of farming system favored modern indoor farming for cheese yield and yield, but favored traditional farming for udder health.
      The potential carryover effects of summer pasture after cows return to indoor rearing have not been extensively studied at the microbiological level. It is worth noting that, in this study, the carryover effect on milk bacteria was negligible, whereas some effects on milk quality, composition, and cheese-making aptitude have been observed (
      • Saha S.
      • Amalfitano N.
      • Sturaro E.
      • Schiavon S.
      • Tagliapietra F.
      • Bittante G.
      • Carafa I.
      • Franciosi E.
      • Gallo L.
      Effects of summer transhumance of dairy cows to alpine pastures on body condition, milk yield and composition, and cheese making efficiency.
      ).

      CONCLUSIONS

      We can conclude that the milk microbiota is very complex, and the majority of bacterial taxa are strongly influenced by farming system as well as by the advancement of season and lactation stage. Transhumance of dairy cows from indoor conditions to summer highland pastures may increase the relative abundances of LAB and other probiotic bacteria (bifidobacteria and propionibacteria) and decrease the abundances of spoilage bacteria, thereby improving the milk in terms of cheese-making aptitude and benefits to human health. This effect disappears after the cows return indoors in the autumn. Systematic differences in milk microbiota among different cows concern some bacterial taxa, particularly the pathogenic bacteria and Clostridiales, signaling the need for new studies on the relationships between the cow genome and milk microbiota. Metagenomic analysis of the milk microbiota appears to be a powerful tool for studying the complex relationships between farming system, individual cow characteristics, and milk value for cheese-making and for human and animal health.

      ACKNOWLEDGMENTS

      The data presented in this study were compiled as part of the TrentinCLA project, co-funded by the Fondazione CARITRO, Cassa di Risparmio di Trento e Rovereto (Palazzo Calepini Via Calepina 1, 38122, Trento, TN, Italy). The authors have not stated any conflicts of interest.