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

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 in-vestigated 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)


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 et al., 2013;Boor et al., 2017;Issa and Tahergorabi, 2019). 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, 2007;Ardö et al., 2017;Nam et al., 2021), and digestion and intestinal functions and integrity in human consumers (prebiotics and probiotics; Aryana and Olson, 2017;Nyanzi et al., 2021). The most interesting unfavorable relationships are the potential effects of some microbial species on the spoilage of milk and dairy products (Quigley et al., 2013;Martin et al., 2021) and on the health of lactating animals (mastitis; Andrews et al., 2019) and of human consumers (pathogenic activities; Verraes et al., 2015).
The recent development of metagenomics is now expanding our knowledge of these aspects of milk microbiota (Addis et al., 2016;Parente et al., 2020). Traditional microbiological studies were based on identifying, isolating, characterizing, and counting individual microbial species or strains (Tilocca et al., 2020). 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 Milk metagenomics and cheese-making properties as affected by indoor farming and summer highland grazing 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 et al., 2017). The difficulties lie mainly in disentangling the confounding effects of environment, management, animal characteristics, season, feedstuffs, and hygiene (Du et al., 2020).
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 et al., 2017). 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 et al., 1999). 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 cheesemaking 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.

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 et al., 2019), and the second reporting some preliminary data on milk microbial counts (Carafa et al., 2020). 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.

98
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 et al., 2020), 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 et al., 2019) 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 , Lactococcus, and Enterococcus (Gagnon et al., 2020). The "other probiotics" category includes all the taxa belonging to the genera Propionibacterium (Rabah et al., 2017) and Bifidobacterium (Prasanna et al., 2014). The "spoilage bacteria" category includes all taxa belonging to the order Clostridiales (Burtscher et al., 2020) and the genera Pseudomonas (Meng et al., 2017), Kocuria (Ribeiro-Júnior et al., 2020), and Alicyclobacillus (Pornpukdeewattana et al., 2020). Finally, the "patho-genic" category includes all taxa belonging to the genus Staphylococcus (Gebremedhin et al., 2022), and the family Enterobacteriaceae (Anand and Griffiths, 2011). 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 et al. (2013) and consisted of the following: rennet coagulation time (RCT, min), curd-firming time (k 20 , min), and curd firmness (a) 30, 45, and 60 min after rennet addition (a 30 , a 45 , and a 60 ; 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 (RCT eq , min), the curd-firming instant rate constant (k CF , %/min), the syneresis instant rate constant (k SR , %/min), maximum curd firmness (CF max , mm), and time to reach CF max (t max , min; Malchiodi et al., 2014).

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 et al., 2018). 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 (%CY CURD ) and cheese solids (%CY SOLIDS ) 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 (REC FAT , %), milk protein (REC PROTEIN , %), total milk solids (REC SOLIDS , %), and milk energy (REC ENERGY , %).

Statistical Analysis
All relative bacterial abundances were log 10 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 log 10 -transformed relative abundances obtained from the milk metagenomic analysis were analyzed with a linear mixed model in the R environment (R Core Team, 2016), 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 et al., 2019) 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) 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 us-ing the same linear mixed model as that used for the metagenomic relative abundances.

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. 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. 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 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.

Combined Effects of Group of Cow and Month of Sampling on Milk Metagenomic Relative Abundances
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).
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 Patterns 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 R 2 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).
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.
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.) 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.
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.
The spoilage bacteria taxa exhibited low correlations with each other and negative correlations with the  . Patterns 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 R 2 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). 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 cat-  . Patterns 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 R 2 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). egories 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.

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).
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).
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 . Patterns 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 R 2 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). 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.

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  Table 2. Significance 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 log 10 relative abundance 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. common factors (environment, diet, hygiene practices, milking routine, etc.) or individual, nonpermanent factors (temporary diseases, cleanness of teat surface, feed selection, etc.; Du et al., 2020;Parente et al., 2020). 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 et al., 2015;Bobbo et al., 2017;Keane, 2019). 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 et al., 2014) and the level of the mastitis indicator represented by the somatic cell content in the milk (Urioste et al., 2010;Pegolo et al., 2021). 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 et al., 2021).

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 et al., 2006;Bobbo et al., 2017).
Multivariate statistical analyses are often used to study associations among different bacterial taxa (Ro- . Heat plot of the correlations between the bacterial and chemical-technological traits included in the factor analysis. LAB = lactic acid bacteria; RCT = rennet coagulation time; k 20 = curd-firming time; a 30, a 45 , a 60 = curd firmness 30, 45, and 60 min after rennet addition, respectively; RCT eq = rennet coagulation time by equation; k CF = curd-firming instant rate constant; k SR = syneresis instant rate constant; CF max = maximum curd firmness; t max = time to reach CF max ; %CY = percentage cheese yields; REC = nutrient recovery traits. drigues et al., 2017) and between these and other milk or cheese characteristics (Nyman et al., 2014). 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  Table 3. Loadings of the latent explanatory factors of the relative abundances of milk bacterial taxa, milk and whey constituents, milk coagulation, curd-firming, and cheese yield 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 RCT eq ), and the time to reach a given (k 20 ) or maximum (t max ) curd firmness. Early gelation is obviously positively correlated with an increase in the traits measuring curd firmness (a 30 , a 45 , a 60 , and CF max ) and allows more time for estimating curd syneresis (k SR ; 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 . 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 et al., 2021). 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  *P < 0.05; **P < 0.01; ***P < 0.001. Figure 8. Patterns 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 R 2 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). 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 et al., 2018b). 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 et al., 2016).
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 . This is why we named this the "udder health factor," as others have done for cattle Macciotta et al., 2012;Dadousis et al., 2018a,b;Cecchinato et al., 2019), as well as for goats and sheep (Manca et al., 2016;Vacca et al., 2016). 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 et al., 2017). 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 et al., 2017). More information on the relationships between milk SCC or SCS and milk properties can be found in Bobbo et al. (2016). 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 (k 20 , a 60 , and CF max ), confirmation that caseins play an important role in curd firming, much more so than in milk coagulation time (Jõudu et al., 2008). Protein is also more correlated than fat with the retention of water in the curd, which explains the positive loading of curd cheese yield . In any case, we should bear in mind that different protein fractions have different effects on milk coagulation and curd-firming traits ), as well as cheese yield and milk nutrient recovery in cheese . It is worth noting that "other LAB taxa" is the only microbiological trait positively included in this factor.

The fifth factor listed in
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 et al., 2019). 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 (k CF , k 20 , and t max ) and syneresis processes (k SR ) is central, and explains the positive loading of a 30 , the small loading of a 45 , and the negative loading of a 60 . 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 et al. (2012). It is worth noting that 2 independent latent factors representing coagulation traits and curd-firming traits were also obtained for goat milk (Todaro et al., 2005).
The eighth factor in Table 3 is characterized by positive loading of the recoveries of milk fat (REC FAT ) and total milk solids (REC ENERGY ) 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 et al., 2020). Other authors have observed large difference in the microbiota of milk produced by cows kept indoors and cows at pasture (Bonizzi et al., 2009;Doyle et al., 2017). 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 . 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 et al., 2020), 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 and Bittante, 2018;Bittante et al., 2021) is due to pasturing. It is well known that different farming systems, as well as individual farms (Bokulich and Mills, 2013;Skeie et al., 2019;Priyashantha and Lundh, 2021) and dairy plants, affect cheese-making efficiency and product quality (Falardeau et al., 2019;Nam et al., 2021). This also supports the claimed specificity of cheeses produced on temporary highland summer farms in Alpine regions (Bittante et al., 2011a,b) and signals the possibility of authenticating the origin of dairy products according to farming system (Bergamaschi et al., 2020).
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 et al., 2018b; 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 et al., 2019).

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.