### Experimental design

Because animal interactions occurred only between the farmers and their own animals, ethical approval for this study was not required as per institutional or local legislation, however, informed consent was obtained from the farmers enrolled within the study. Before the start of the trial, the sampling protocol was provided and explained to each farmer, highlighting the study aims and methodology. According to the experimental design, only the first colostrum had to be sampled for the study; therefore, the calf was separated from the dam as soon as possible after birth and was not allowed to suckle. Colostrum from unsupervised calvings was not collected for the study.

The sampling procedure of this study is the same as that described in [12, 13]; briefly, collection of individual colostrum samples (n = 548) took place on nine commercial farms of northern Italy from 2019 to 2020, covering the four calendar seasons. Farms involved in the study had between 60 and 190 lactating Holstein cows and used a total mixed ration method to feed animals, free stall barns, twice-a-day milking, with access to pasture. Moreover, none of the enrolled farms underwent a vaccination protocol against rotavirus, coronavirus, or E. coli. To reduce the impact of time on colostrum composition, only colostrum samples collected within 6 h after calving were included in the study. A representative colostrum sample was taken from a cow based on the pooled colostrum from all functional quarters at first milking.

For each cow, 120 mL were collected in a single sterile plastic tube without preservative (SMIPA srl, Vicenza, Italy). Farmers were in charge of colostrum collection and were instructed to annotate the cow ID on the tube and to freeze (− 20 °C) the samples as soon as possible. Periodically, colostrum samples were retrieved from the farms, transferred to the laboratory of the Department of Agronomy, Food, Natural resources, Animals and Environment of the University of Padova (Legnaro, Italy), and stored at − 20 °C until analysis. Date of birth, date of calving, and parity of sampled cows were retrieved and farm-specific general information about colostrum management was registered.

### Quantification of Ig and BRIX

Samples were thawed overnight at 4 °C in a waterbath and RID analyses were carried out using bovine-specific assays [12]. For this study, 34 ‘Bovine IgG RID Kit’, 34 ‘Bovine IgM RID Kit’, and 30 ‘Bovine IgA RID Kit’ were purchased in advance from Triple J Farms (Bellingham, WA, US). Each kit included a plate and reference sera. Dilution in pure deionized water 1:5 (v/v) for IgG and 1:3 (v/v) for IgA and IgM was performed to reduce the concentration of the target component in such a way that the concentration would fall within the assay detection range [12]. Five µL of diluted colostrum were subsequently pipetted in the wells of the three RID plates (IgG, IgA, and IgM, respectively). After incubation at 20 °C for 24 h, plates were scanned at high resolution to measure the diameter of precipitated rings using the image processing program ImageJ (Laboratory for Optical and Computational Instrumentation, University of Wisconsin-Madison, WI). For each well, the precipitated ring was measured in duplicate to calculate the final diameter (mm) as the average of the two measurements. The final diameter was used to calculate the concentration (g/L) of the target component (IgG, IgA, or IgM) by using the known concentrations and measured diameter of the three reference sera to calculate a standard curve for each individual plate. Concentrations of the reference sera for each Ig were: 1.80, 14.72, and 28.03 g/L for IgG, 0.53, 1.94, and 3.87 g/L for IgA, and 0.62, 2.00, and 3.81 g/L for IgM.

Non-readable or non-circular RID rings due to issues during colostrum dilution or pipetting were treated as missing values, and the same rule was applied to samples with concentration values that were beyond the detection range of the kit. Repeatability of this method was preliminarily tested in cow colostrum as described by [12] and [13]. Based on the indications given by the US Department of Health and Human Services, Food and Drug Administration [16], the repeatability of RID, measured through the intra-assay coefficient of variation, was sufficiently robust (< 10%) to refrain from duplicate analyses [12, 13]. For the duration of the trial, RID analyses were carried out by the same laboratory technician, to reduce interpersonal variation.

Colostral BRIX values were measured using an optical refractometer manufactured by Euro Horse Line (San Marino) with a working range from 0 to 32%. After tube inversion, 0.25 mL of colostrum were placed on the daylight plate of the device and the BRIX values were recorded in duplicate for each sample by the same operator for the duration of the study. The average of the two measurements was used in later analyses.

A small representative aliquot of colostrum was lyophilized through a freeze-drying process and subsequently analysed with the Kjeldhal method according to AOAC [17] in order to determine the total protein content (PC_{C}, %). For each cow for which the information was available, the milk protein content (PC_{M}, %) and yield (PY_{M}, kg/d) at the first test-day, i.e., between 5 to 75 days in milk, were retrieved from the national database of Italian Holstein, Brown Swiss and Jersey Association (ANAFIBJ, Cremona, Italy). After merging the available test-day milk data and colostrum traits, information on 500 out of the 548 cows was available for further analyses.

### Data analysis

The trial was designed to ensure that only one sample per cow was present in the final dataset. Prior to the multi-trait genetic analysis, significance of the fixed effects of parity (1, 2, 3, 4, and ≥ 5), calving season (January–March, April–June, July–September, and October–December), calving year (2019 and 2020), and herd (n = 9), and least squares means of BRIX for these effects were tested using the software ASReml v4.1 [18]. The model also included random additive genetic animal effects \(\sim N(0, \mathbf{A}{\sigma }_{\text{ag}}^{2}\)) and a residual \(\sim N(0, \mathbf{I}{\sigma }_{\text{e}}^{2})\), where \(\mathbf{A}\) is the additive genetic relationship matrix, \({{\sigma}_{\text{ag}}^{2}}\) is the additive genetic variance, \(\mathbf{I}\) is an identity matrix of appropriate order, and \({{\sigma}_{\text{e}}^{2}}\) is the residual variance. The pedigree file (n = 6714 individuals) used for setting up \(\mathbf{A}\) included all sampled cows and their ancestors up to six generations back.

Subsequently, a four-trait animal model was run with the same software [18] to estimate genetic parameters for BRIX, IgG, IgA, and IgM using the aforementioned fixed and random effects. In matrix notation, the model can be represented as:

$$\left[\begin{array}{c}\begin{array}{c}{\mathbf{y}}_{\mathbf{1}}\\ \dots \end{array}\\ {\mathbf{y}}_{\mathbf{4}}\end{array}\right]=\left[\begin{array}{ccc}{\mathbf{\mathbf{X}}}_{\mathbf{1}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \ddots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{X}}_{\mathbf{4}}\end{array}\right]\left[\begin{array}{c}{\mathbf{b}}_{\mathbf{1}}\\ \dots \\ {\mathbf{b}}_{\mathbf{4}}\end{array}\right]\text{+}\left[\begin{array}{ccc}{\mathbf{Z}}_{\mathbf{1}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \ddots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{Z}}_{\mathbf{4}}\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_{\mathbf{1}}\\ \dots \\ {\mathbf{a}}_{\mathbf{4}}\end{array}\right]\text{+}\left[\begin{array}{c}{\mathbf{e}}_{\mathbf{1}}\\ \dots \\ {\mathbf{e}}_{\mathbf{4}}\end{array}\right],$$

(1)

where \({\mathbf{y}}_{\mathbf{1}}\), \({\mathbf{y}}_{\mathbf{2}}\), \({\mathbf{y}}_{\mathbf{3}}\), and \({\mathbf{y}}_{\mathbf{4}}\) are the vectors of phenotypic observations of the dependent variables BRIX, IgG, IgA, and IgM, respectively; \(\mathbf{b}\) is the vector of fixed effects; \(\mathbf{a}\) is the vector of random additive genetic effects; \(\mathbf{e}\) is the vector of random residuals; and \(\mathbf{X}\) and \(\mathbf{Z}\) are incidence matrices relating the corresponding effects to the dependent variable. In the four-trait animal model, the variance–covariance structure of the random effects was assumed equal to:

$$\mathbf{var}\left[\begin{array}{l}{\mathbf{a}}_{\mathbf{1}}\\ \dots \\ \begin{array}{l}{\mathbf{a}}_{\mathbf{4}}\\ {\mathbf{e}}_{\mathbf{1}}\\ \begin{array}{l}\dots \\ {\mathbf{e}}_{\mathbf{4}}\end{array}\end{array}\end{array}\right]=\left[\begin{array}{l}\begin{array}{llllll}\mathbf{A}{{\upsigma}}_{{{\text{a}}{\text{g}}}_{1}}^{2}& \cdots & \mathbf{A}{{\upsigma}}_{{{\text{ag}}}_{14}} & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \vdots & \ddots & \vdots & \mathbf{0}& \mathbf{0}&\mathbf{0} \\ \mathbf{A}{{\upsigma}}_{{{\text{ag}}}_{14}}& \cdots & \mathbf{A}{{\upsigma}}_{{\text{ag}}_{{\text{4}}}}^{2} & \mathbf{0}& \mathbf{0}& \mathbf{0} \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}{{\upsigma}}_{{{\text{e}}}_{1}}^{2}& \cdots & \mathbf{I}{\upsigma}_{{\text{e}}_{14}} \\ \mathbf{0}& \mathbf{0}&\mathbf{0}& \vdots & \ddots & \vdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{I}{\upsigma}_{{\text{e}}_{14}}& \cdots & \mathbf{I}{\upsigma}_{{\text{e}}_{4}}^{2}\end{array}\end{array}\right].$$

(2)

The phenotypic variance (\({{{\upsigma}}}_{{\text{p}}}^{2}\)), heritability (\({{{h}}}^{2}\)), coefficient of phenotypic (\({\mathrm{CV}}_{\mathrm{p}}\)) and additive genetic variation (CV_{ag}), and the phenotypic (r_{p}) and additive genetic correlations (r_{ag}) were calculated by means of conventional formulas [18]. To test if \({{{h}}}^{2}\) estimates differed significantly from 0, the log-likelihood ratio test was used [19]. Briefly, the log-likelihood of the control model was compared to a model where the additive genetic variance was constrained to be approximately zero (= 1 \(\times\) 10^{–3}). The significance threshold was set at P < 0.05.

The r_{p} and r_{ag} of BRIX and IgG with PC_{C,} PC_{M}, and PY_{M} were also estimated using the approach described above. Two multivariate models, one with colostral BRIX [Eq. (3)] and another one with colostral IgG [Eq. (4)], were run including the same fixed and random effects described for Eq. (1) and using the following phenotypic vectors \(\mathbf{y}\):

$$\left[\begin{array}{c}\begin{array}{c}{\mathbf{y}}_{\mathbf{B}\mathbf{R}\mathbf{I}\mathbf{X}}\\ {\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{C}}}\end{array}\\ \begin{array}{c}{\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{M}}}\\ {\mathbf{y}}_{{\mathbf{P}\mathbf{Y}}_{\mathbf{M}}}\end{array}\end{array}\right],$$

(3)

and

$$\left[\begin{array}{c}\begin{array}{c}{\mathbf{y}}_{\mathbf{I}\mathbf{g}\mathbf{G}}\\ {\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{C}}}\end{array}\\ \begin{array}{c}{\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{M}}}\\ {\mathbf{y}}_{{\mathbf{P}\mathbf{Y}}_{\mathbf{M}}}\end{array}\end{array}\right],$$

(4)

where \({\mathbf{y}}_{\mathbf{B}\mathbf{R}\mathbf{I}\mathbf{X}}\), \({\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{C}}}\), \({\mathbf{y}}_{{\mathbf{P}\mathbf{C}}_{\mathbf{M}}}\), \({\mathbf{y}}_{{\mathbf{P}\mathbf{Y}}_{\mathbf{M}}}\) and \({\mathbf{y}}_{\mathbf{I}\mathbf{g}\mathbf{G}}\) contain phenotypic observations available for BRIX, PC_{C}, PC_{M}, PY_{M}, and IgG, respectively. The use of two separate analyses was needed to avoid convergence issues, which usually arise with the REML algorithm when performing multivariate models to estimate unique solutions and (co)variance components for two (or more) traits that are strongly correlated, such as BRIX and IgG.

### Refractive index cut-off

The R package ‘OptimalCutpoints’ with the Youden method [20, 21] was used to calculate the BRIX value, which is useful to discriminate good- from low-quality colostrum samples, by simultaneously reducing specificity (true negative cases) and improving sensitivity (true positive cases). Following Šimundić [22], the BRIX classification accuracy was evaluated by the area under the curve (AUC) obtained. Receiver operating characteristic (ROC) analyses were carried out considering a RID IgG threshold at 50 and 60 g/L.