### Herd description and milking machine settings

Our study was performed in the research herd at the Norwegian University of Life Sciences. The herd consisted of 91 Norwegian Red cows divided in two groups, each milked by one AMS (DeLaval VMS, Tumba, Sweden) with identical settings. These 91 cows formed the source population. The two groups were situated in immediate proximity and had identical housing conditions. The teatcups were equipped with DeLaval 20 M VMS liner (Product No. 92725901). The system vacuum was set to 45 kPa, the pulsator ratio was 65% and the pulsation rate was 60 cycles/min. Quarter take-off limit (switch-level) was set to 0.24 kg/min. The system was set with a low-vacuum period of 6 s and a delay from initiation of detachment to detachment of 4 s. The herd, including housing conditions, management and equipment, is comparable to commercial AMS farms in Norway. Data on parity and days in milk (DIM) in the study herd were obtained from the Norwegian Dairy Herd Recording System [21]. Except from a period of 3–5 days immediately after calving, the included cows had not been milked by CMS in the current lactation. Cows in second or later lactations had been milked by CMS in previous lactations.

### Teat-end scoring and teat dimensions

Teats on all 91 lactating cows in the herd were evaluated using a scoring system where the thickness and roughness of the callous ring of the teat orifice was categorized into one of eight classes [22]. Teat length (Length) and width 0.5 cm from apex (Apex) and base (Base) was measured by placing the teat between a white background and a transparent 0.5 cm grid (DeLaval, Tumba, Sweden). In addition, the shape of the teat-end was classified in the following groups; pointed, round, flat or inverted, and the teat position was registered. All registrations were performed once per cow during a 2-day period. The udder and teats were cleaned and stripped prior to the evaluation. The same person performed all scoring and registrations. The time from last milking to teat evaluation was not standardized.

Two outcome variables for logistic regression models were established by transforming the results from the teat-end scoring. The variable THICKNESS was given the value 0 if the thickness of the callosity ring was thin or not visible, and the value 1 if it was classified as medium, thick or extreme according to the scoring system by Neijenhuis et al. [22]. The variable ROUGHNESS was given the value 0 if a smooth callosity ring was registered, and the value 1 if a rough callosity ring was registered. This approach was chosen because previous research has shown that severe degrees of teat-end callosity is significantly related to mastitis risk [2, 23].

### Vacuum recordings

VaDia vacuum recorders (BioControl, Rakkestad, Norway) were used to record vacuum at a rate of 200 Hz in the short milk tube, pulsation tube, and MPC in the four individual teatcups during milking in both milking stations. The data were collected during three herd visits, between 1 and 4 weeks after the initial teat-end scoring was performed. We used data from a convenience sample of cows that entered the milking stations voluntarily without interference from the herdsmen. For cases of duplicate vacuum registrations of the same quarter, the recording taken closest to the day of teat-end scoring was used. Cows not entering the milking stations during our visits were excluded from the study. Teats that had been dried-off prior to our observational period and teats where MTT variables could not be calculated due to missing vacuum recordings were also excluded.

### Calculation of MTT variables per quarter

In accordance with common procedures for MTT under field conditions, the different periods of the milking were found by evaluating the vacuum recordings [24]. For each milking two main periods were identified, based on the vacuum registrations from the short milk tube and the MPC: (a) the main milking period and (b) the overmilking period. The main milking period was characterized by high milk flow and stable vacuum conditions in the short milk tube. The start of the main milking period was identified by monitoring the average short milk tube vacuum in 10 s periods, until the decline from one period to the next was less than 0.3 kPa. The end of the main milking period, which coincides with the start of the overmilking period, was defined as the point where a marked change in MPC vacuum became apparent as described by Borkhus and Rønningen [19]. Automatic detection of these MPC vacuum changes was set at the point where the MPC vacuum increased at least 30% plus 2 kPa above a weighted running average. The result of the automated procedure was controlled manually and if necessary adjusted taking changes in short milk tube vacuum into account as described previously [19]. The end of the overmilking period was set to the initiation of detachment, i.e. the point where the short milk tube vacuum started to fall markedly (≥ 5 kPa) shortly before the end of milking.

The average vacuum in the short milk tube during the main milking period (MTVAC) was calculated as the mean of all vacuum recordings from the short milk tube within the main milking period. The average vacuum levels in the MPC during the main (MPCVAC) and the overmilking periods (MPCOM) were calculated accordingly.

To estimate the forces applied on the teat-end by the liner during the closed-phase of the pulsation cycle, the variable teat compression intensity (COMPR) was calculated for each teat. This variable is comparable with “drucksumme”, as described by Spohr and Uhlenbruck [12]. Vacuum records from the pulsation tube and short milk tube of 10 consecutive pulsation cycles from 60 to 70 s into the main milking period were used for the calculation of COMPR. Differential pressure across the liner wall, i.e. the difference between the short milk tube vacuum and the pulsation tube vacuum, was calculated throughout the 10 pulsation cycles. Touchpoint pressure difference (TPPD) is the pressure difference across the liner wall when two sides of the liner achieves or loses contact during closing and opening respectively [9]. TPPD is traditionally measured without a teat in the teatcup [25]. In this study, however, approximated values for TPPD was derived from the pulsation tube vacuum curves at the points of fastest liner wall movement during opening and closing, respectively [26]. The average of the approximated TPPD at opening and closing was used in further calculations. The closed-phase of the pulsation cycles were defined as the period when the differential pressure across the liner wall was higher than the approximated TPPD. For each of the 10 pulsation cycles, an integral of the differential pressure across the liner wall minus approximated TPPD as a function of time over the liner closed period was calculated. The average of the integrals found in the 10 pulsation cycles represented teat compression intensity, COMPR (kPa s).

Quarter average milk flow rate (AVGFLOW) and quarter peak milk flow rate (PEAKFLOW), for the same milkings in which the vacuum measurements were performed, were obtained from the AMS software (DeLaval, Tumba, Sweden). The milking stations were equipped with ICAR-approved milk meters using near-infrared technology providing in-line data on milk flow and milk yield used in the calculation of these variables.

### Statistical analyses

Principal component analysis (PCA) is a multivariate technique that can be used to explore multi-dimensionality of data and to reduce a large set of variables to a small number of latent variables, principal components, which nevertheless retain most of the information in the dataset. PCA is also useful for the analysis of intercorrelation of variables. We applied PCA to study the relationships between the variables obtained from the different registrations, and as suggested by Dohoo et al. [27], we used PCA as a complementary technique in the subsequent model-building process. We used the 2-dimensional scatter plot of loadings for two specified components from PCA, which is most useful for interpreting principal component 1 versus principal component 2, as they contain the most important information in the data. In our application of PCA, we focused on the geometric interpretation of the relationships between variables, plotted as points in the component space using their loadings as coordinates on the “circle of correlations”. In addition to PCA, we also used a linear regression model to describe the mathematical relationship between MTVAC and AVGFLOW.

The following variables were evaluated as potential explanatory variables in logistic regression models describing the likelihood of a teat-end having a rough or thickened callosity ring, respectively: DIM, Parity, MTVAC, COMPR, AVGFLOW, PEAKFLOW, MPCVAC, MPCOM, Length, Apex, Base, teat position and overmilking time. Parity was categorized as; first lactation, second lactation, and greater than second lactation.

Linearity between the outcome variables and the explanatory variables was assessed by inspecting lowess smoothing curves obtained with a logit transformation of the outcome variables; THICKNESS and ROUGHNESS (Stata SE/14, Stata Corp., College Station, TX, USA).

To establish logistic regression models for teat-end callosity roughness and thickness, we initially tested each of the explanatory variables in univariable logistic regression models to detect associations with the two outcome variables, THICKNESS and ROUGHNESS. Variables with a P-value less than 0.2 were further evaluated for inclusion in multivariable models. In order to avoid including highly correlated variables in the same model, variables identified by PCA as belonging either to the same cluster or in clusters aligned on the opposite side of the circle of correlation, were evaluated in separate models. We repeated a backwards selection procedure to build multiple multivariable models for each outcome, and variables with P-value ≥ 0.05 were excluded from the final multivariable models. Because teat-end callosity thickness and roughness has been shown to vary between parities and lactation stages, parity and DIM were forced into all multivariable models [1].

We expected registrations between teats within the same cow to be correlated, and to account for this, we included a random intercept at cow level in both the univariable and multivariable models. Because the multivariable models were non-nested, i.e. the predictors in one model could not be considered as subsets of the predictors in other models, Bayesian Information Criterion (BIC) was calculated to compare model fit [27].

Data from different sources were assembled in SAS 9.4 (SAS Institute Inc., Cary, NC, USA) to form a final dataset. We used the STATA meqrlogit procedure for the logistic regression analysis and the regress procedure for the linear regression analysis (Stata SE/14, Stata Corp., College Station, TX, USA). The PCA was conducted using the statistical software JMP Pro version 12 (SAS Institute Inc., Cary, NC, USA).