Correlations between bulk tank milk analysis with weather conditions in dairy farms under tropical environments

Abstract

The consequences of climate change on agriculture have generated concern among researchers and decision-makers, especially regarding the effects it will have on this sector. The phenomenon is expected to affect the productivity of livestock systems, even more so in grazing livestock, where cattle and pastures are directly impacted by climatic variables. This study evaluated the climatic influence on bulk tank milk production and quality in 38 tropical dairy farms, using Pearson’s partial correlation and canonical correlation analyses. Farm production level and milk quality traits were merged with meteorological information obtained from climatological stations distributed around the Valle del Cauca department, Colombia. According to farms’ milk production levels, four evaluation groups (EG) were established for the analysis, and within each group, the available information between 2.0 and 8.5 years was included. Pearson’s partial correlations among productive and climatic variables were scarce and low (r <± 0.10) within each EG. On the other hand, canonical correlation analyses between productive and climatic variable sets presented a linear increase since the establishment of EGs from 1 to 4 (r ranging from 0.39 to 0.59). The magnitude of canonical correlation coefficients depended on bulk tank dairy farm production levels, being the most productive systems the most susceptible, especially for meteorological variables related to temperature and relative humidity. As a consequence, dairy farms with higher performances must define mitigation strategies to reduce the weather effects. Multivariate correlations are recommended to evaluate the relationship between milk production, composition, and environmental variables in grazing dairy systems of tropical areas.

Highlights

• Livestock farming in the tropics are more vulnerable to climate change.

• Heat stress has an impact on the productivity of dairy cows.

• Dairy production systems not intensive under grazing conditions and adapted breeds of ecology zones are not affected by adverse meteorological conditions.

Keywords:

• Dairy systems
• tropical conditions
• heat stress
• grazing livestock

Introduction

Milk production under the current conditions of climate change constitutes a profound challenge for dairy systems, especially for those in livestock systems under grazing conditions (Scharf et al. 2011; Lopes et al. 2016). Climate is one of the main factors that affect the behaviour, physiology, and productive efficiency of animals (Hammami et al. 2013). The climatic variables that most influence the comfort of animals in the tropics are ambient temperature, solar radiation, wind speed, relative humidity, and precipitation (Polsky and von Keyserlingk 2017; Mylostyvyi and Chernenko 2019). The effects of single variables or combinations can lead to unfavourable situations known as heat stress, negatively impacting ingestion, digestion, respiration, and blood circulation, among other physiological processes (Herbut and Angrecka 2018), which, in turn, impact milk production (productivity and quality) and dairy cow reproduction (Pavani et al. 2015; Liu et al. 2019).

It was estimated that climate change has caused an increase of 35% in the costs of livestock production in intertropical climate zones (Thamo et al. 2017). In tropical Latin America, livestock farming represents one of the most important lands uses from a social, economic, and environmental point of view which are usually developed in conditions characterised by high air temperatures and humidity, making livestock production more vulnerable to climate change, hence generating a loss of employment, animals, and profitability (Landers 2007; Oosting et al. 2014). In Colombia, more than 95% of production systems are based on grazing (Carulla and Ortega 2016; Molina et al. 2016).

The exposure of cattle to environmental conditions is influenced by climatic variables, generating direct effects on their performance and behaviour (Scharf et al. 2011), and indirectly impacting the quantity and quality of forage, resulting in consequences such as variability in milk bulk composition (Stürmer et al. 2018).

Although different works have demonstrated the relationship between climate and milk production and quality (Barrera et al. 2015; Mylostyvyi and Chernenko 2019), few studies have jointly addressed the relationships between meteorological variables and milk production and characteristics (Milani et al. 2015b; Gabbi et al. 2017; Stürmer et al. 2018).

Multivariate analysis methodologies such as symmetric canonical correlation are used to analyse functions, relations, or expressions that can allow comprehensive form, to reduce without loss of generality, variables of different types, making your use wide and diverse, thus helping the study of the variables with reliability in the complex production systems (Costa et al. 2018). Canonical correlation and other multivariate methods had been used to analyse the relationship between milk composition and climatic conditions (Stürmer et al. 2018).

The objective of this study was to evaluate the relationship between weather variables and milk tank volume and quality (fat, protein, bacterial count in terms of colony-forming units, and somatic cell count) of grazing dairy farms located in the geographical zone of Valle del Cauca, Colombia, under low tropical conditions, using partial correlation and canonical correlation analysis.

Material and methods

Production systems

The bulk tank milk production and composition information were obtained from 38 dairy farms located in the geographical zone of Valle del Cauca, Colombia, in low tropic conditions between 3’30 to 4’78 N and 76’21 to 76’46 W. These production systems usually base their animal feeding directly on rotational grazing of available pastures. The predominant pastures are Cynodon plectostachyus, Megathyrsus maximus, Brachiaria decumbens or Brachiaria Brizantha (to a lesser extent), and native pasture such as Paspalum spp. (Molina et al. 2017; Morales-Vallecilla and Ortiz-Grisales 2018). Pure breeds and multiracial cows were found among dairy farms. Bos taurus × Bos indicus crosses represented 70% of the animal population, whereas Holstein, Jersey, Normando, Lucerna, and Gyrolando breeds represented the remaining 30%. All farms used machine milking twice a day, and the milk collected was cooled. Adequate livestock farming practices were carried out with different levels of technology, according to the production system facilities.

The information used in this work was obtained from the database of one of the milk-purchasing companies of the region. The milk production information was obtained from fortnightly records of total bulk tank milk production (MP; in kg), milk fat (as %), milk protein (as %), milk total solids (TS; as %), colony-forming units (CFU/ml), and somatic cell count (SCC/ml) variables, analysed since 2010 to mid of 2018. Every Once every fifteen days, samples of bulk tank milks were taken from each farm for laboratory analysis, to estimate quality and microbiological variables. The sample was fresh milk free of preservatives. On the same day, the average MP of each farm was estimated by dividing the amount of milk delivered by the number of milked cows in each farm. Fat and protein were determined following the Ekomilk ultrasound technology (Bulteh, Stara Zagora, Bulgary). Total solids were estimated using the Richmond formula $$(T=\frac{L}{4}+1.2F+0.14),$$ where T: total solids; L: lactodensimeter reading, and F: Fat (Patel and Boghra 2018). CFU/ml was analysed through the Standard Plate Count SPC method by surface seeding, with a previous dilution of 10⁻³ (Laird et al. 2004). Somatic cell count was determined using the Delaval optical fluorescence counter (DeLaval International AB, Tumba, Sweden).

Climatic information

Climate information was obtained from the automated meteorological network of sugarcane research centre Cenicaña; it has approximately 30 meteorological stations distributed around the geographical zone. The closest meteorological station was assigned to each farm, ensuring that the distance between them did not exceed 20 km. Daily values of mean temperature (Med_T; °C), minimum temperature (Min_T; °C), maximum temperature (Max_T; °C), mean relative humidity (Med_RH; %), maximum relative humidity (Max_RH; %), minimum relative humidity (Min_RH; %), precipitation (PP; mm), solar radiation (SR; cal/cm²), wind speed (WS; m/s), and temperature and humidity index (THI) were used. The THI was estimated from daily mean temperature and relative humidity as suggested by Gomes da Silva et al. (2007). Historic climatic information from 2010 to mid-2018 was obtained and averaged fortnightly to match them with bulk tank milk production and composition variables of the farms. During this process, some climatic information was excluded due to missing values in productive variables. In the final database, 39% of the farms presented information from January 2010 to June 2018, and the remaining had at least two consecutive years of information.

Statistical analysis

Fortnightly data of bulk tank milk production and its chemical and microbiological composition (MP, fat, protein, TS, CFU, and SCC) were gathered with environmental variables (Max_T, Med_T, Min_T, Max_RH, Med_RH, Min_RH, PP, SR, WS, and THI) in each farm. Data outliers were removed according to box plot representation and biological meaning. Later, descriptive statistics were estimated for each variable within the farms, using the InfoStat software (InfoStat 2019). Fortnightly, bulk tank milk production presented a high variability among dairy systems, with values between 3440 and 21197 kg. To reduce the variation and establish homogeneous groups of farms according to with the average fortnightly bulk tank milk production level, four evaluation groups (EG) were created: EG 1 (1,000–5,000 kg/fortnight), EG 2 (5,001–10,000 kg/fortnight), EG 3 (10,001–15,000 kg/fortnight), and EG 4 (>15,000 kg/fortnight).

To evaluate the association between bulk tank milk production, its composition (chemical and microbiological) and environmental variables of each EG, Pearson’s partial correlation, and canonical correlation was used, as suggested by Stürmer et al. (2018). To avoid confusing the information patterns, the four EGs were analysed independently, and in each case, variables such as MP, protein, TS, CFU, SCC, PP, SR and WS were transformed to the square root or log10 to meet the assumption of multivariate normality (Sherry and Henson 2005). Variable histogram plots were utilised to guarantee the normal distribution. The Max_RH variable was eliminated in all the EGs, as it does not fulfil the assumption.

Pearson’s partial correlation allows estimating the association between two quantitative variables after eliminating the influence of other variables (Vargha et al. 2013). In this case, the procedure allows to control the effect of categorical farm and fortnight variables, and remove their residual correlation (Stürmer et al. 2018). The estimated partial correlations were used to identify bivariate associations among bulk tank milk production, its composition, and meteorological variables. To identify a significant linear trend, only correlation values equal or greater than ± 0.30 were selected for further analysis. Partial correlations were estimated using PROC GLM of SAS 9.4 (SAS Inc SI 2011).

Associations between variables were also evaluated from a multivariate approach. Bulk tank milk production along with the chemical and microbiological variable sets was related with the respective meteorological variable sets through a canonical correlation analysis applied to each EG, using the PROC CANCORR of SAS 9.4 (SAS Inc SI, 2011). In this case, when high bivariate correlation (r > 0.80) occurred within each variable group, just one was selected for the analysis, to avoid multicollinearity, as a consequence of linear combination between variables within each set (Manly 2004). Canonical correlation analysis creates weighted synthetic variables (also known as canonical functions) within the evaluated groups of characteristics and applies the Pearson correlation method on both synthetic variables. The first canonical function creates the two synthetic variables with the largest possible correlation, thus representing the first canonical correlation. The residual variance left over by the first function is utilised to create a second canonical function, where two more synthetic variables, which are correlated as strongly as possible, are constructed. These new synthetic variables are uncorrelated with the preceding ones (double orthogonality). The process of constructing canonical functions continues until the minimum number of variables in the smaller group are attained or when the total variance of the original variables is explained (Manly 2004; Sherry and Henson 2005).

To estimate the proportion of variance shared between bulk tank milk production—composition and environmental variable sets across all functions—an effect size index was estimated in each EG (Sherry and Henson 2005). Finally, functions that were significant (P < .05) and explained a reasonable amount of the shared variability between the variable sets were interpreted. In these functions, the canonical redundancy statistics available in PROC CANCORR output were analysed to identify how well the environmental canonical variables can predict the original bulk tank milk production and composition variables (SAS Inc SI 2011).

Results and discussion

Descriptive statistics of the bulk tank milk production and composition of the 38 evaluated farms are shown on Table 1. Bulk milk production averaged 3440, 7283, 12022, and 21197 kg/fortnightly for EGs 1, 2, 3, and 4, respectively. EGs 1, 2, and 3 included 42%, 29%, and 16% of the studied farms, respectively, whereas the remaining 13% were in EG 4. This is a common trend since most of the dairy systems in the studied geographical zone belong to small producers, where more productive and technologically advanced systems are scarce (Carulla and Ortega 2016; Molina et al. 2017; Morales-Vallecilla and Ortiz-Grisales 2018). Although the bulk tank milk production was the main difference between EGs, the milk quality parameters were similar. Fat, protein and TS levels ranged between 3.75–3.79%, 3.14–3.22%, and 12.32–12.51%, respectively, whereas CFU and SCC showed a higher variability, with values ranging between 230.18 and 828.11 (×1000 CFU/ml) and between 439.07 and 523.60 (×1000 cells/ml), respectively. In general, the climatic behaviour was characterised by hot days with an average Max_T and Max_RH of 30.6 °C and 97%, respectively, and low precipitation rates (3.4 mm). Productive and climatic conditions found in this analysis were in line with other studies carried out in the same geographical zone (Molina et al. 2016; Morales-Vallecilla and Ortiz–Grisales 2018). The heterogeneity observed among farms related to milk production and composition complies with typically tropical conditions, where different animal management, genetics and physiology states are present (Ramírez et al. 2019).

Table 1. Descriptive statistics of bulk tank milk production, composition and environmental variables of the 38 dairy systems evaluated.

Variable	n	Mean	Median	S.D	Mínimum	Máximum	n	Mean	Median	S.D	Mínimum	Máximum
EG −1 (No. of farms = 16 – average milking cows 22 ± 7)							EG − 2 (No. of farms = 11 average milking cows 33 ± 8)
MP (Kg)	2048	3440.47	3310.69	1165.62	1190.98	6639.88	968	7283.41	7236.68	1726.15	3571.45	11217.18
Fat (%)	2044	3.79	3.81	0.29	3.06	4.66	962	3.80	3.81	0.23	3.15	4.48
Protein (%)	1971	3.19	3.17	0.12	2.91	3.51	949	3.22	3.23	0.11	2.92	3.47
TS (%)	2044	12.41	12.40	0.38	10.94	13.45	962	12.51	12.47	0.30	11.74	13.35
CFU (×1000 cell/mL)	2036	324.62	102.41	890.48	5.63	6919.50	962	620.11	233.50	1064.38	4.45	5671.36
SCC (×1000 cell/mL)	1507	476.86	413.84	292.50	83.88	1679.63	839	523602.47	471.14	301.05	138.27	1913.55
Max_T (°C)	2045	30.35	30.26	1.05	28.16	33.05	968	30.52	1.02	30.40	28.62	33.02
Med_T (°C)	2045	23.43	23.37	0.77	21.72	25.29	968	23.45	23.38	0.68	22.02	25.11
Min_T (°C)	2045	18.91	18.89	0.62	17.22	20.53	968	18.85	18.84	0.58	17.36	20.27
Max_HR (%)	2045	96.05	96.28	2.50	88.90	99.73	968	97.51	97.85	1.66	92.75	99.81
Med_RH (%)	2044	77.14	77.54	4.57	66.35	86.66	968	78.46	78.96	3.99	69.58	86.85
Min_RH (%)	2045	46.45	47.04	4.84	34.77	56.20	968	46.86	47.60	4.65	36.10	56.10
PP (mm)	2048	3.39	2.80	2.87	0.02	13.28	968	3.21	2.61	2.54	0.05	11.06
SR (cal/cm^2)	2045	403.90	402.27	38.24	307.82	515.91	968	410.10	409.55	37.44	323.46	502.97
WS (m/s)	2045	1.51	1.51	0.17	1.13	2.00	965	1.53	1.52	0.17	1.21	2.03
THI	2048	72.04	72.02	1.07	68.02	74.35	968	72.22	72.16	0.86	70.22	74.31
EG − 3 (No. of farms = 6 – average milking cows 55 ± 12)							EG − 4 (No. of farms = 5 – average milking cows 122±24)
MP (Kg)	670	12022.57	12533.08	4260.77	3204.50	22025.58	635	21197.58	21467	7331.21	2583	46429
Fat (%)	670	3.78	3.79	0.20	3.07	4.34	636	3.75	3.72	0.25	3.09	4.57
Protein (%)	642	3.14	0.11	3.14	2.90	3.44	621	3.17	3.16	0.12	2.6	3.63
TS (%)	670	12.32	12.32	0.26	11.57	12.98	636	12.35	12.3	0.35	11.48	13.8
CFU (×1000 cell/mL)	670	230.18	58.33	743.76	6.17	5325.33	637	828.11	128.82	2072.02	1.00	19498.45
SCC (×1000 cell/mL)	526	520.00	510.83	319.19	77.00	1601.33	515	439.07	400.01	234.73	20.00	1905.46
Max_T (°C)	670	30.81	30.64	1.13	28.68	33.62	639	30.82	30.7	1.32	27.48	34.67
Med_T (°C)	670	23.43	23.36	0.71	22.04	25.19	640	23.58	23.53	0.81	21.51	26.39
Min_T (°C)	670	18.75	18.76	0.63	17.11	20.41	640	18.89	18.88	0.71	16.8	21.31
Max_HR (%)	670	98.23	98.91	1.82	92.95	100.00	640	95.84	96.33	3.68	84.2	100
Med_RH (%)	670	79.37	79.80	4.76	67.53	88.79	640	76.92	77.15	5.05	62.93	89.2
Min_RH (%)	670	45.75	46.14	5.85	31.94	57.19	640	44.75	45.33	5.8	27.4	58.07
PP (mm)	670	3.61	2.99	2.86	0.01	13.04	640	3.22	2.61	2.62	0	15.92
SR (cal/cm^2)	670	401.76	400.66	42.04	267.04	518.00	640	400.77	396.17	54.01	262.99	602.56
WS (m/s)	670	1.33	1.32	0.18	0.91	1.78	640	1.5	1.49	0.22	0.62	2.33
THI	670	72.27	72.19	0.91	70.38	74.61	640	72.31	72.44	1.1	68.94	77.62

EG: evaluated group; n: number of observations; S.D: Standard deviation; MP: milk production; TS: total solids; CFU: colony forming units; SCC: somatic cell count; Max_T: maximum temperature; Med_T: mean temperature; Min_T: minimum temperature; Max_RH: maximum relative humidity; Med_RH: mean relative humidity; Min_RH: minimum relative humidity; PP: precipitation; SR: solar radiation; WS: wind speed; THI: temperature and humidity index.

Scarce and low levels of univariate correlations were found between bulk tank milk production, composition, and climatic variables. Among milk chemical and microbiological traits, a positive relationship between fat and TS was found in all the EGs (0.58–0.72). Higher correlation coefficients of 0.88 and 0.92 were reported by Sourabh et al. (2017) and Bondan et al. (2018), respectively. The same positive tendency was observed for protein and TS in EGs 1 and 3 (0.39–0.41, respectively). These are expected associations since milk fat and protein are the main components of total solids.

Correlations between climatic variables varied within the groups. Broadly, all evaluated groups presented significant and negative relationships between Min_T and Max_T (−0.30 to −0.65), Med_T and Med_RH (−0.66 to −0.97), and Max_T and Min_RH (−0.60 to −0.77); meanwhile, positive associations between Med_T and THI (0.61 to 0.99), and Med_RH and THI (0.35–0.96) were observed. Within the EG 4, direct associations between Min_T with PP and Med_T (0.45 and 0.62, respectively), Max_T with Med_RH and Med_T (0.66 and 0.70, respectively), Med_T and Min_RH (0.45), Med_RH and Min_RH (0.85), and SR and WS (0.33) were also found (Picinin et al. 2019).

Bivariate correlations among bulk tank milk production, composition, and climatic variables were low and non-significant (P > .05), with most of the values below ± 0.10. These results contrast with the study of Picinin et al. (2019), where rainfall and high temperatures exert an influence on raw milk volume and quality of dairy farms located under similar environmental conditions of the present study.

According to this univariate approach, environmental variables are not related to the volume and composition of the raw milk collected fortnightly within the evaluated production systems. These results suggest, to some extent, that climatic conditions do not influence the performances of grazing dairy cows under tropical conditions. Low correlations may be linked to the particular management, breeds, and microclimate condition of the farms within each EG (Ruiz et al. 2019). Similar weak correlations between milk production and climate variables were reported in other studies (Milani et al. 2015; Gabbi et al. 2017; Mylostyvyi and Chernenko 2019).

The relationship between bulk tank milk production, composition (chemical and microbiological), and climatic variables was explored deeply in this study with a second multivariate approach. Previous correlation analysis among climatic variables in each EG indicated high associations with THI and Med_T variables (>0.90). In those groups, to avoid multicollinearity, only THI and Med_T were considered for subsequent analysis, since they have a well-known effect on dairy cows’ performance and heat load (Gomes da Silva et al. 2007; Hammami et al. 2013; Gabbi et al. 2017; Mylostyvyi and Chernenko 2019).

Canonical correlation results obtained in the present work are shown in Table 2. All multivariate tests (Wilks’ Lambda, Pillai’s Trace, the Hotelling-Lawley Trace, and Roy’s greatest root) were significant (p < .0001 in all cases) within each of the EGs, suggesting that a latent structure could be identified using the multivariate approach. These results were accompanied by an analysis of effect size indexes, to guarantee the practical significance of the study. Effect size represents the variance proportion shared between the group of studied variables, and large values are always desirable. The effect size was estimated as 1–Wilks’ Lambda coefficient (Sherry and Henson 2005), obtaining values of 0.21, 0.40, 0.43, and 0.47 for EGs 1, 2, 3, and 4, respectively. Even when all the EGs presented significant multivariate test results, only two of them attained a relevant effect size, EGs 3 and 4.

Table 2. Canonical correlation analysis and associated statistics between bulk tank milk production, composition and environmental variables.

EG	Canonical functions	Rc	$R_{\mathrm{c}}^2$	Eigen value proportion	Cumulative	Approximate F value	Pr > F
1	1	0.39	0.15	0.69	0.69	8.64	<.0001
	2	0.22	0.05	0.19	0.88	3.86	<.0001
	3	0.11	0.01	0.05	0.93	2.18	0.0018
	4	0.10	0.01	0.04	0.97	2.20	0.0097
	5	0.08	0.01	0.02	0.99	1.96	0.0685
	6	0.04	0.00	0.01	1.00	1.49	0.2268
2	1	0.44	0.20	0.44	0.44	10.54	<.0001
	2	0.37	0.13	0.28	0.73	8.27	<.0001
	3	0.32	0.10	0.20	0.92	6.14	<.0001
	4	0.14	0.02	0.04	0.96	2.85	0.0007
	5	0.11	0.01	0.02	0.99	2.90	0.0081
	6	0.09	0.01	0.01	1.00	3.39	0.0343
3	1	0.55	0.30	0.66	0.66	7.36	<.0001
	2	0.36	0.13	0.23	0.89	3.68	<.0001
	3	0.24	0.06	0.09	0.98	1.76	0.0204
	4	0.07	0.01	0.01	0.99	0.47	0.9304
	5	0.06	0.00	0.01	1.00	0.48	0.8253
	6	0.04	0.00	0.00	1.00	0.37	0.6885
4	1	0.59	0.35	0.71	0.71	7.19	<.0001
	2	0.28	0.07	0.11	0.82	3.10	<.0001
	3	0.25	0.06	0.09	0.91	2.79	<.0001
	4	0.22	0.04	0.07	0.98	2.22	0.0047
	5	0.11	0.01	0.02	0.99	1.03	0.4143
	6	0.07	0.00	0.00	1.00	0.72	0.5390

EG: evaluated group; R_c: canonical correlation, $R_{\mathrm{c}}^2:$ squared canonical correlation. Canonical correlations were declared significant at 5% level.

The first canonical correlation coefficients observed between EG presented a linear increase as bulk tank milk production improved within the groups (Figure 1), with values ranging between 0.39 and 0.59 (Table 2). Although the sample sizes were smaller in EGs 3 and 4 compared to EGs 1 and 2, these values still remain significant and higher than those reported in other studies (Aliç Ur Al and Baritçi 2013; Stürmer et al. 2018). Converse to what was found in previous correlations analysis, these results demonstrated how production systems with higher bulk tank milk volumes are more influenced by meteorological conditions; similar findings were reported by Haygert-Velho et al. (2018). This affirmation agrees with other studies where milk production was reduced per unit increase in THI, and the magnitude varies among dairy systems characteristics, environmental conditions, and raised breeds (Promket et al. 2020; Ruiz et al. 2019). Usually, in intensive systems under tropical conditions, high productive dairy cows are more sensitive to hot temperatures, mainly because of the high metabolic rate (Rhoads et al. 2011; Gantner et al. 2017). Cows exposed to high temperatures in tropical ambiances reduce feed intake and increase their metabolic maintenance requirements by 7–25%, causing a decrease in both nutrient supply and nutrient available used for milk synthesis (Ramírez et al. 2019). Maintenance requirements increase due to the rise in basal metabolism caused by activation of the thermoregulatory system (Gaughan and Mader 2014; McCafferty et al. 2017). All these changes negatively impact the individual productive capacity of dairy cows, which can be directly reflected in the amount of milk collected in the production system during a specific time period (Haygert-Velho et al. 2018; Ramírez et al. 2019).

Figure 1. Relationship between a linear combination of bulk tank milk production – composition and environmental variables obtained in the first canonical functions of each EG.

In this study, significant canonical functions explaining a representative amount of variance between bulk tank milk production, composition, and climatic variable sets (the highest squared canonical correlation coefficients ${\mathit{R}}_{\mathit{c}}^2$) were chosen for further interpretation. This criterion was established because some functions may not explain the relationship between the variable sets enough, making their interpretation questionable and increasing the risk of highlighting effects not important or replicable in other studies (Sherry and Henson 2005). This is the case of EGs 1 and 2: even when several canonical functions were significant, most of the variance explained within each group was low, just like their effect size (0.21 and 0.40 for groups 1 and 2, respectively). Canonical correlations coefficients and associated statistics are found in Table 2. According to the squared canonical correlation coefficients, the first canonical functions among all groups explained most of the variance shared between bulk tank milk production, composition, and climatic variables, with values between 0.15 and 0.35. To avoid redundant results interpretation, only these functions were selected for further analysis. Table 3 shows the standardised canonical coefficients of the first function associated to each group. In this case, it is preferable to interpret the standardised coefficients since variables are not measured in the same units (SAS Inc SI 2011). The analysis of these coefficients together with the correlations between the original variables (milk characteristics and climatic conditions), with their respective canonical variable, allows for the identification of the parameters that most contributed to each canonical function; however, their interpretation must be carried out carefully since, as suggested by Moraes et al. (2015), where a correlation is found, its causality cannot always be predicted.

Table 3. Standardised canonical coefficients, correlation of original variables and their canonical function and canonical redundancy statistic obtained within the first canonical function of each Evaluated Group (EG).

	EG 1			EG 2			EG 3			EG 4
	V1	rv,cv	CRS	V1	rv,cv	CRS	V1	rv,cv	CRS	V1	rv,cv	CRS
Production variables
MP	0.285	0.36	0.03	−0.653	−0.71	0.04	0.113	0.24	0.06	−0.122	−0.12	0.09
Fat	0.222	0.15		−0.382	−0.06		0.103	−0.43		0.275	0.61
Protein	−0.390	−0.50		0.316	0.49		0.123	−0.03		0.353	0.71
TS	−0.036	−0.15		0.479	0.35		−0.661	−0.48		0.430	0.67
CFU	−0.585	−0.60		0.396	0.49		0.832	0.84		0.213	0.13
SCC	0.555	0.57		−0.058	0.02		0.043	0.13		0.518	0.48
Environmental variables
Max_T	1.841	0.16	0.005	−1.109	−0.30	0.02	−0.119	0.01	0.04	0.058	−0.01	0.05
Med_T	–			–			–			−0.150	−0.24
Min_T	−0.734	−0.23		0.369	−0.06		0.165	−0.08		−0.418	−0.40
Med_RH	–			–			–			0.863	0.68
Min_RH	1.60	0.35		−0.835	−0.17		−0.308	−0.11		–
PP	−0.065	0.01		0.001	0.03		0.048	−0.21		−0.160	0.23
SR	−0.361	0.02		0.389	0.23		−0.263	−0.21		0.364	0.57
WS	0.035	−0.08		0.650	0.56		0.989	0.93		0.161	0.09
THI	−0.103	0.20		−0.223	−0.38		−0.194	−0.06		0.384	0.07

EG: evaluated group; V1: canonical function 1; r_v,cv: correlation between original variables and their canonical function; CRS: canonical redundancy statistic; MP: milk production; TS: total solids; CFU: colony formig units; SCC: somatic cell count; Max_T: maximum temperature; Med-T: mean temperature; Min_T: minimum temperature; Max_RH: maximum relative humidity; Med_RH: mean relative humidity; Min_RH: minimum relative humidity; PP: precipitation; SR: solar radiation; WS: wind speed; THI: temperature and humidity index.

The first canonical correlations among all groups were positive with increasing values since EGs 1–4 (Table 2), suggesting that the relationship between the linear combination of bulk tank milk production, composition, and environmental variables was strengthened as bulk tank milk volume increased among groups. Stürmer et al. (2018) also found a strong canonical correlation (0.95) between milk production, composition, and price with climatic variables in systems with bulk tank milk production close to 14500/fortnightly, similar to the volumes observed in the farms belonging to EGs 3 and 4, where the highest canonical correlations coefficients were obtained (0.55 and 0.59, respectively).

The EG 1 presented a canonical correlation of 0.39, with an explained variance of 15%. Standardised canonical coefficients showed that the first canonical variable for bulk tank milk characteristics is a weighted difference between CFU/ml (−0.585), protein (−0.390), and SCC (0.555). Environmental standardised coefficients indicated that Max_T (1.841), Min_RH (1.60), and Min_T (−0.734) were the variables that most contributed to this component. Correlations between original variables and their canonical function also presented the highest associations (−0.60–0.57 and −0.23–0.35, for bulk tank milk characteristics and climatic variables, respectively). These results can be interpreted to mean that high Max_T and Min_RH, and lower Min_T are associated with bulk milk tanks with lower protein and CFU, and higher SCC. Previous climatic conditions are comparable to those found in tropical dry seasons, with high ambient temperature and relative humidity. Under summer conditions, Bertocchi et al. (2014) reported lower protein content and high SCC in the bulk milk data of 3,328 dairy farms. Bouraoui et al. (2002) also observed the same tendency in lactating Holstein cows, while Renna et al. (2010) reported a reduction in milk protein content during the hottest months of the year in grazing cows.

The EG 2 presented a 0.44 canonical correlation and 20% explained variance. Variables like MP, TS, and CFU presented the most contributing standardised coefficients associated with the milk characteristics component (−0.653, 0.479, and 0.396, respectively); meanwhile, the Max_T, Min_RH, and WS variables were for the climatic canonical function (−1.109, −0.835, and 0.650, respectively). Correlations between original variables and their canonical function also presented the highest associations (−0.71–0.49 and −0.30–0.56, for bulk tank milk characteristics and climatic variables, respectively). Within the EG 2, results suggest that low Max_T and Min_RH, and high WS conditions are associated with a reduction in bulk tank milk production and an increase in TS and CFU levels. Some studies report the negative effect of temperature, relative humidity, and wind speed on milk production (Zuluaga and Restrepo 2009; Lambertz et al. 2014; Promket et al. 2020), while Ramírez et al. (2019) described an inverse relation between TS and the milk produced, supporting the observed fact that lower milk volumes are associated with higher TS due to the concentration effect of milk solids.

In EG 3, a 0.55 canonical correlation and 30% explained variance were found. The highest standardised coefficients in canonical bulk tank milk characteristics were observed in TS (−0.661) and CFU (0.832), while WS (0.989) and Min_RH (−0.308) were for the environmental canonical function. Additionally, TS, CFU, WS, and Min_RH showed a strong association with their respective canonical functions (−0.48, 0.84, 0.93, and −0.11 respectively). A conclusion that may arise from these results is that a high WS and lower Min _RH have corresponded with lower bulk tank TS and an increase in CFU. A similar inverse association between relative humidity and CFU had been reported in other studies ( Zuluaga and Restrepo 2009; Stürmer et al. 2018). A direct relationship between relative humidity and TS was found by Stürmer et al. (2018), contrary to what was obtained within this group. No associations between wind speed and TS or CFU were identified. The influencing variables within this group are similar to those found in EG 2, reinforcing the effect of relative humidity and wind speed, mainly on milk TS and CFU variables.

Finally, the EG 4 showed a canonical correlation of 0.59 and a 35% explained variance. Bulk tank milk characteristics within this group were represented mainly by SCC (0.518), TS (0.430), and protein (0.353); correlations between these variables and their canonical functions were significant (0.48, 0.67, and 0.71, respectively). The highest standardised coefficients related to the environmental canonical function were observed in Med_RH (0.863), Min_T (−0.418), and SR (0.364), with correlation coefficients between these variables and their canonical function of 0.68, −0.40, and 0.57, respectively. These results suggest that environmental conditions with high Med_RH and SR, and low Min_T are associated with increases in bulk tank milk protein, TS, and SCC levels. Zuluaga and Restrepo (2009) showed positive correlations between relative humidity with milk protein and SCC. Stürmer et al. (2018) found negative associations among minimum temperature and milk protein, and TS and SCC. Relative humidity and temperature variables used to construct temperature and humidity indexes have shown a positive correlation with SCC (Bertocchi et al. 2014). Solar radiation also positively influences milk composition, especially the protein and TS percentages (Hill and Wall 2015;; Stürmer et al. 2018).

According to previous results, climatic conditions from the studied geographical zone differentially influence the bulk tank milk production and composition among the EGs, which can be attributed to differences in dairy systems management, proper milking practices implementation, animal genetics, and microenvironmental conditions of the farms associated to each group. Climatic variables like Max_T, Min_T, Med_RH, Min_RH, WS, and SR, together with bulk tank milk volume, protein, TS, CFU, and SCC variables constituted the main factors that explained the canonical correlations found. It is important to highlight that temperature and relative humidity conditions were the most frequented climatic variables influencing bulk tank milk production and composition between the EGs. This observation agrees with other research reports demonstrating that these climatic variables are the most frequent factors impacting the production level and comfort of grazing animals (Gomes da Silva et al. 2007; Polsky and von Keyserlingk 2017; Mylostyvyi and Chernenko 2019); however, other climatic variables such as wind speed and solar radiation can also have a significant influence, so they should be considered in future studies (Gauly et al. 2013). These results demonstrate how the climatic variables interact to differentially influence the productivity of the farms and the milk physicochemical and microbiological characteristics, even with production systems located under the same study area.

The relevant result observed with the multivariate approach is that the meteorological variables are correlated with bulk tank milk production, composition, and hygienic conditions in tropical dairy systems based on grazing (R_c = 0.39–0.59). The magnitude of the association varied, depending on the bulk tank milk volume of the dairy farm, being the most productive systems, the most susceptible to climatic conditions. In tropical areas, high ambient temperatures and relative humidity are common, increasing the occurrence risk of heat stress cases in dairy cows, affecting the dairy farm productivity directly, especially those with bulk tank milk production levels above 15,000 kg/fortnight. To reduce the heat stress impact and face new climate change scenarios, these systems require establishing heat stress mitigation practices like the use of trees in grazing areas, implementation of cooling systems (e.g. fans, sprinkles, self-controlled showers) during milking, improving water availability and quality, diet manipulation, and animal selection for heat stress resistances (Lopes et al. 2016; Gantner et al. 2017; Polsky and von Keyserlingk 2017).

The canonical redundancy statistics presented in Table 3 examined how well the canonical variables can predict the original ones (SAS Inc SI 2011). In this study, it is useful to identify the amount of variation in bulk tank milk production and composition variables that is accounted for by the linear combination of meteorological variables. Results showed that canonical redundancy statistics increased between 3 and 9% from EG 1 to 4. Only the 9% value found in EG 4 was similar to the 10.2% reported by Stürmer et al. (2018), which can be attributed to the fact that in both situations, the dairy farms included in the analysis presented a bulk tank milk production equal or above 15,000 kg or l/fortnightly. The canonical redundancy coefficients obtained indicate that a linear combination of climatic variables explained at least 3% and a maximum of 9% of the variance of bulk tank milk production and quality variables, depending on the farm production level. These findings are in line with the previous canonical correlation analysis results, where the increasing influence of climatic variables on bulk tank milk characteristics was emphasised.

The positive relationship between bulk tank milk production and composition with environmental variables observed in the present study with canonical correlation analysis is similar to that found with other multivariate methods like cluster, factor, and discriminant analysis (Shiker 2012; Haygert-Velho et al. 2018; Stürmer et al. 2018). These findings corroborate the importance of multivariate approaches to fully understanding the complex influence of several climatic variables on a particular dairy production scenario.

Conclusions

Compared with univariate analysis, multivariate approaches allow for a deeper exploration of factors, contributing to relationship analyses among the climatic variables and bulk tank milk production. Canonical correlation analysis found a significant association between climatic variable sets and bulk tank milk production, and quality in tropical grazing-based dairy systems. The effect and magnitude of climatic influence depended on the bulk tank milk volume of dairy farms, with the most productive systems being the most influenced. Climatic variables like Max_T, Min_T, Med_RH, Min_RH, WS, and SR, along with bulk tank milk volume, protein, TS, CFU/ml, and SCC characteristics were the main factors contributing to the observed canonical correlations. Temperature and relative humidity remained the most important variables affecting bulk tank milk production and quality, even though other variables such as wind speed and solar radiation also deserve attention for further studies. As a result, systems with high production levels must establish practices to mitigate the impact of climatic variables that affect animals’ performance and the herd’s milk yield in tropical environments.

Acknowledgments

The work was carried out as part of the project “Quantification of climate impacts on dairy cattle systems in the low tropic through data mining and modeling”, Announcement 776-2017, MinCiencias.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.