Udder Ultrasonography of Dairy Cows: Investigating the Relationship between Echotexture, Blood Flow, Somatic Cell Count and Milk Yield during Dry Period and Lactation

Simple Summary Udder health is important for dairy cows’ productivity and welfare. Ultrasonography has been proven a practical method to examine udder health on farms. In this study, we investigated the relationship between ultrasonographic features and indicators of udder health and productivity. Twenty-one Holstein cows were examined repeatedly during a crucial period starting from the end of the late lactation stage, continuing throughout the dry period and ending in the early stage of the consecutive lactation period. Statistical models presented significant associations between blood flow in the milk vein and udder echotexture parameters with daily milk production and milk somatic cell count. Udder ultrasonography is a useful and practical tool for the comprehensive assessment of udder health. Abstract Udder health of dairy cows is related to their productivity and welfare. The period from dry-off to calving and early lactation is crucial. Ultrasonography is a useful and practical tool for the examination of the mammary parenchyma and blood flow. This observational study investigated the relationship between udder echotexture features, blood flow volume (BFVol) in the milk vein, milk somatic cell count (SCC) and daily milk yield (DMY) from late lactation, throughout the dry period and consecutive early lactation. Seventeen repeated measurements were performed on twenty-one Holstein cows. The udder parenchyma was examined with B-mode ultrasonography. Udder echotexture was studied using 15 features: Numerical Pixel Value (NPV), Pixel Standard Deviation (PSD), Skewness, Excess, Contrast, Homogeneity, Correlation, Entropy, Run Percentage, Long-Run Emphasis, Grey Value Distribution, Runlength Distribution, Gradient Mean Value, Gradient Variance and Percentage of Non-zero Gradients. Blood flow in the milk vein was examined with spectral Doppler. Linear mixed-effects models were employed to investigate relationships between BFVol, udder echotexture features, SCC and DMY throughout the study period. Our models showed that a 1 kg increase in DMY was associated with a significant increase of 0.25 L/min in the expected BFVol and that a 1,000,000-cells/mL increase in SCC was associated with a significant BFVol decrease of 0.49 L/min, keeping all other variables constant. Multivariable models showed significant associations between DMY and NPV, between PSD and Long-Run Emphasis, and between SCC and NPV, PSD, Gradient Mean Value, Homogeneity, Gradient Variance and Entropy. In conclusion, udder echotexture and BFVol in the milk vein are related to SCC and milk yield. Ultrasonography can be used for the comprehensive assessment of udder health in support of precision dairy farming.


Introduction
Ultrasonography is an easy, non-invasive method to visualize the udder parenchyma and blood flow on farms [1]. First, amplitude-mode (A-mode) ultrasonography was used to examine internal structures of the udder. Two decades later, brightness-mode (B-mode) was employed to visualize the bovine udder and teat [2]. Ever since, the physiological appearance of the glandular parenchyma and specific ultrasonographic findings under pathological conditions have been reported [3][4][5]. Ultrasound has also been used to assess parenchymal development in dairy heifers [6,7]. Finally, the applicability of three-dimensional ultrasonography to the examination of the bovine udder has been tested and produced good-quality perspective images of the mammary gland [8].
Echotexture analysis of udder B-mode sonograms has been used to evaluate mammary parenchyma and increase the acquired information. Echotexture refers to the appearance, structure and arrangement of the parts of an object within an ultrasonographic image [9]. Echotexture analysis employs statistical features, such as mean Numerical Pixel Values (NPV) and Pixel Standard Deviation (PSD)-which are first-order echotexture features-to quantify parenchymal changes depicted in grey-scale images (B-mode sonograms) due to the presence of milk or under pathological conditions [10,11]. In buffaloes, first-order echotexture features were significantly correlated with milk somatic cell count (SCC) [12]. Furthermore, second-and higher-order echotexture features have been used in udder echotexture analysis [13]. These features take into account the position of a pixel, as well as its spatial relationship with its neighboring pixels, and reveal the spatial organization of the image echotexture [14][15][16]. Recently, a deep learning algorithm was trained to process udder echotexture features and predict productivity traits of dairy cows [17].
Blood flow of the udder has been studied with invasive techniques [18][19][20] and with non-invasive color Doppler ultrasonography [21][22][23]. The milk vein (subcutaneous abdominal vein) is the main route of blood drainage from the udder. It plays a key role in the mammary function and is commonly used for drug administration [24,25]. The morphology and blood flow of the milk vein have been studied with color Doppler ultrasonography in lactating Brown Swiss cows [25,26]. Another study compared blood flow of cows in the dry period (DP), cows with 10 kg of daily milk yield (DMY) and cows with 20 kg of DMY and found that milk yield has a profound effect on blood flow variables of the milk vein [27].
Research has shown that udder health and productivity traits affect blood flow and echotexture. Despite being the breed with the highest milk production and worldwide spread, only a few studies have focused on Holstein cows' udder echotexture, blood flow and morphological parameters of the milk vein [11,28]. There are also limited longitudinal data on changes in udder echotexture features and blood flow of the same cows throughout the period of late lactation, dry-off, udder involution during the DP, calving and early lactation. Moreover, in most studies, only first-order statistics (NPV and PSD) have been used to quantify udder echotexture. According to recent developments [13,17], secondand higher-order statistics could optimize echotexture analysis of the mammary gland, but further research on this topic is needed.
Finally, there is a lack of a comprehensive approach to the ultrasonographic evaluation of udder health and mammary function in the same study. Such an approach could lead to statistical models taking into account the relationships between udder parenchymal echotexture, blood flow, milk SCC-as an indicator of subclinical mastitis-and milk production. This comprehensive approach could be a valuable contribution to clinical practice, as it could aid the data-driven assessment of mammary parenchyma and blood flow, and the detection of alterations due to pathological conditions. Utilizing this approach, veterinarians and farmers could apply precision management strategies to udder health and productivity, based on the ultrasonographic evaluation of the mammary function.
This longitudinal study of clinically healthy Holstein dairy cows, starting from late lactation, continuing throughout dry period and ending in early lactation, aimed to investigate the following: (i) changes in morphological and blood flow parameters of the right and left milk vein, (ii) udder echotexture changes with the employment of first-, secondand higher-order statistics, and (iii) the relationships between blood flow volume in the milk vein, udder echotexture features, milk somatic cell count and milk production.

Animals
The study was approved by the Research and Ethics Committee of the Faculty of Veterinary Medicine, Aristotle University of Thessaloniki (protocol number: 1399/05-04-2019). All operations were conducted in accordance with the EU Directive 2010/63/EU for animal experiments and the university's guidelines for animal research. Reporting of this observational study was conducted under STROBE-Vet reporting guidelines.
In total, 23 clinically healthy, lactating, pregnant Holstein dairy cows (11 primiparous and 12 multiparous) were initially selected from one commercial herd located in the region of Thessaloniki, Greece. Cows in late and early lactation stages were housed in the same freestall barn, in the low-and high-lactation groups, respectively. They were milked twice daily with a conventional milking system, and DMY was automatically recorded. They were fed a total mixed ration designed specifically for either the low-or high-lactation group, offered once a day, every morning. Cows in DP were housed as far-off and close-up groups and fed total mixed rations designed for these groups, offered every morning. Detailed information on the nutritional content of the rations is provided in Supplementary Materials.
Median DMY and interquartile range (IQR) at the peak of lactation was 47.3 (10.3) kg, and on the last day of lactation, it was 24.3 (8.5) kg. Median duration of the dry period (IQR) was 56 (6) days. Abrupt dry-off with blanket dry cow treatment was applied using an intramammary antibiotic suspension (Nafpenzal ® , Intervet International B.V., Boxmeer, The Netherlands) (300 mg procaine benzylpenicillin, 100 mg nafcillin and 100 mg dihydrostreptomycin). On entrance day in the study, median age (IQR) of the cows was 4.2 (1.5) years old; 11 cows were at the end of 1st lactation, 8 were at the end of 2nd lactation, 3 were at the end of 3rd lactation, and 1 was at the end of 4th lactation.

Study Design and Clinical Examination
This observational study consisted of 17 repeated measurements performed on each of the 23 cows throughout 3 consecutive production stages starting from late lactation (7 days before dry-off and at dry-off day), continuing during the DP (on the 3rd, 7th, 21st and 35th day of the dry period, and 7 and 3 days prior to calving) and ending in early lactation (on calving day and on the 3rd, 7th, 14th, 30th, 45th, 60th, 75th, and 90th days of the new lactation period).
All measurements took place in the middle of the time interval between morning and afternoon milking (10:00-11:00 a.m.). A thorough general physical examination was performed prior to any handling [29]. Udder health was assessed clinically, milk was inspected visually for abnormalities, and California Mastitis Test was performed on farm using a commercially available reagent (BOVIVET CMT Liquid, Jørgen Kruuse A/S, Denmark) [30]. On each measurement day, individual cow DMY was obtained from farm data records and quarter milk SCC (qSCC) and individual cow SCC (iSCC) were recorded on farm with DeLaval Cell Counter (DeLaval International AB, Tumba, Sweden). SCC was logtransformed to somatic cell score (SCS) using the formula SCS = 3 + log 2 (SCC/100.000) [31].
Cows with clinical symptoms, moderate or strong California Mastitis Test reaction and/or major pathogens isolated from their milk samples were excluded after the appearance of these conditions. Among the 23 cows initially enrolled in this study, 2 cows were excluded due to the appearance of clinical conditions (clinical mastitis and left displaced abomasum), and all were diagnosed promptly and treated successfully. The final number of cows included in the study was 21.

Ultrasonography of the Milk Vein (B-Mode and Spectral Doppler)
Cows were examined at standing position, non-sedated, restricted in a cattle crush and handled gently to minimize stress. Ultrasonographic examination of all cows was performed consistently by the first author. Morphology and blood flow of the milk veins were examined bilaterally with B-mode and spectral Doppler (triplex) ultrasonography, respectively. Cows had normal resting heart rate during the vascular flow examination [32]. A portable ultrasound scanner (MyLab™ OneVET, Esaote S.p.A., Genoa, Italy) was used, equipped with a broad-bandwidth multi-frequency linear transducer (SV3513, Esaote S.p.A., Genoa, Italy; 2.5-10 MHz). B-mode frequency was 10 MHz, scanning depth was 4-6 cm, gain 58% and time-gain compensation in neutral position.
A straight part of the milk vein was selected for measurement located at the midpoint between its cranial and caudal part. The area was shaved, washed, degreased with 70% alcohol and covered with coupling gel (Figure 1). To avoid compression of the vein while ensuring sufficient contact, an appropriate amount of coupling gel was applied on the transducer, without the latter directly touching the skin. First, the transducer was positioned in cross-section of the vein. Distance of the milk vein from skin surface (D1) (cm), its vertical diameter from intima to intima (D2) (cm) and vein area in cross-section (A) (cm 2 ) were measured from B-mode images.
appearance of these conditions. Among the 23 cows initially enrolled in this study, 2 cows were excluded due to the appearance of clinical conditions (clinical mastitis and left displaced abomasum), and all were diagnosed promptly and treated successfully. The final number of cows included in the study was 21.

Ultrasonography of the Milk Vein (B-Mode and Spectral Doppler)
Cows were examined at standing position, non-sedated, restricted in a cattle crush and handled gently to minimize stress. Ultrasonographic examination of all cows was performed consistently by the first author. Morphology and blood flow of the milk veins were examined bilaterally with B-mode and spectral Doppler (triplex) ultrasonography, respectively. Cows had normal resting heart rate during the vascular flow examination [32]. A portable ultrasound scanner (MyLab™ OneVET, Esaote S.p.A., Genoa, Italy) was used, equipped with a broad-bandwidth multi-frequency linear transducer (SV3513, Esaote S.p.A., Genoa, Italy; 2.5-10 MHz). B-mode frequency was 10 MHz, scanning depth was 4-6 cm, gain 58% and time-gain compensation in neutral position.
A straight part of the milk vein was selected for measurement located at the midpoint between its cranial and caudal part. The area was shaved, washed, degreased with 70% alcohol and covered with coupling gel (Figure 1). To avoid compression of the vein while ensuring sufficient contact, an appropriate amount of coupling gel was applied on the transducer, without the latter directly touching the skin. First, the transducer was positioned in cross-section of the vein. Distance of the milk vein from skin surface (D1) (cm), its vertical diameter from intima to intima (D2) (cm) and vein area in cross-section (A) (cm 2 ) were measured from B-mode images.  Then, the transducer was positioned longitudinally to the vein, and color Doppler gate was activated and placed at the center of the milk vein. Blood flow sample volume cursor included at least 2/3rds of the vein diameter, 16-20 mm wide, and theta angle was 60 • (Figure 1). Blood flow was examined for 2 min per vein, and 5 optimal spectral Doppler images were obtained from each milk vein measurement for further processing with MyLab™_Desk v. 9.0 software (Esaote S.p.A., Genoa, Italy). Blood flow parameters used were time-averaged mean velocity (TAMV) (cm/s), peak velocity (Vpeak) (cm/s) and milk vein blood flow volume (BFVol) (L/min) [27]. The total number of spectral Doppler sonograms used was 690 [345 measurement days × 2 milk veins (right/left) per cow].

Ultrasonography and Echotexture Analysis of the Mammary Parenchyma
Udder parenchyma was examined with the same ultrasound scanner, equipped with a broad-bandwidth multi-frequency convex probe (SC3421, Esaote S.p.A., Genoa, Italy; 2.5-6.6 MHz). All four udder quarters were examined at each measurement day, and two optimal B-mode images per quarter were acquired. The probe was placed in coronal plane to examine the front quarters and in sagittal plane for the rear quarters, 15-20 cm above the base of each teat ( Figure 2). Ultrasound settings remained constant throughout the scans: frequency 6.6 MHz, scanning depth of 15 cm, gain 58% and time-gain compensation in neutral position.
Then, the transducer was positioned longitudinally to the vein, and color Doppler gate was activated and placed at the center of the milk vein. Blood flow sample volume cursor included at least 2/3rds of the vein diameter, 16-20 mm wide, and theta angle was 60° (Figure 1). Blood flow was examined for 2 min per vein, and 5 optimal spectral Doppler images were obtained from each milk vein measurement for further processing with MyLab™_Desk v. 9.0 software (Esaote S.p.A., Genoa, Italy). Blood flow parameters used were time-averaged mean velocity (TAMV) (cm/s), peak velocity (Vpeak) (cm/s) and milk vein blood flow volume (BFVol) (L/min) [27]. The total number of spectral Doppler sonograms used was 690 [345 measurement days × 2 milk veins (right/left) per cow].

Ultrasonography and Echotexture Analysis of the Mammary Parenchyma
Udder parenchyma was examined with the same ultrasound scanner, equipped with a broad-bandwidth multi-frequency convex probe (SC3421, Esaote S.p.A., Genoa, Italy; 2.5-6.6 MHz). All four udder quarters were examined at each measurement day, and two optimal B-mode images per quarter were acquired. The probe was placed in coronal plane to examine the front quarters and in sagittal plane for the rear quarters, 15-20 cm above the base of each teat ( Figure 2). Ultrasound settings remained constant throughout the scans: frequency 6.6 MHz, scanning depth of 15 cm, gain 58% and time-gain compensation in neutral position. The total number of sonograms used was 2760 (345 measurements × 4 udder quarters × 2 B-mode images obtained from each udder quarter). The ultrasound scanner produced images with a resolution of 720 by 540 pixels. Each pixel had one of the possible shades of grey (0 to 255). Echotexture analysis procedure presented in Figure 2 was performed with EchoVet v.2.0 software (Aristotle University of Thessaloniki, Greece) [13]. Each image was divided into four imaginary quarters [33]. A circular region of approximately 5000 pixels was manually selected within each of the image quarters. Selection criteria were The total number of sonograms used was 2760 (345 measurements × 4 udder quarters × 2 B-mode images obtained from each udder quarter). The ultrasound scanner produced images with a resolution of 720 by 540 pixels. Each pixel had one of the possible shades of grey (0 to 255). Echotexture analysis procedure presented in Figure 2 was performed with EchoVet v.2.0 software (Aristotle University of Thessaloniki, Greece) [13]. Each image was divided into four imaginary quarters [33]. A circular region of approximately 5000 pixels was manually selected within each of the image quarters. Selection criteria were physiological gland parenchyma, with medium homogeneous echogenicity, avoiding mammary vessels [3]. The overall region of interest selected from each B-mode image was 20,000 pixels [17].

Statistical Analysis
Data quality control was performed prior to any analysis. Descriptive analyses were performed to summarize the animals' characteristics, compared by production stage (late lactation/dry period/early lactation). Quantitative variables were described with means (standard deviation (SD)) and median (IQR). One-way ANOVA for repeated measurements was used to test for differences between production stages. Linear mixedeffects (LME) models under the Maximum Likelihood method [34] were employed to describe the relationship between and the evolution over time of variables of the parenchymal echotexture, BFVol, DMY and SCC, taking into consideration correlations within repeated measurements.
We included random effects (random intercept and slope) for the time of measurement to improve model fitting by accounting for sources of variability that are not explained by the fixed effects in the models. We accounted for individual differences in baseline levels of the outcome variable (individual variability) and for individual differences in the rate of change over time. This is important because individuals may have different trajectories of change over time, even if they started at the same baseline level. Time was included in the models both as fixed and as random effect. Restricted cubic splines with four knots were used to estimate the time evolution curves of udder BFVol, DMY and SCS. The location of the knots was determined based on Harrell's recommended percentages for a 4-knot RCS [35] with knots at (1) 7 days before dry-off, (2) 36th day of the DP and (3) 30th and (4) 90th day of consecutive lactation. Interactions with time terms were investigated for all models.
All models were minimally adjusted for lactation period and RCS of time to address potential sources of bias due to confounding and capture the outcome's non-linearity. To improve interpretability of our results, we log e transformed Grey Value Distribution and Runlength Distribution. Additionally, we rescaled Homogeneity, Correlation, Gradient Variance and Long-Run Emphasis values. Model selection procedure with backwards evaluation was performed to determine the final multivariable models, keeping variables with p ≤ 0.10. The most significant predictor variables were identified with this method, while controlling for other variables in the models. Model selection was based on the Akaike information criterion and log-likelihood. Multicollinearity was checked for every multivariable model, and collinear variables were excluded based on the highest value of Generalized Variance Inflation Factor [36].
The general form of the LME model with a 4-knot RCS that investigated the relationship of BFVol with DMY and iSCC (Supplementary Materials) was as follows: where Y ij is the BFVol in the milk vein of the ith cow at the j th measurement day, RCS(time) i is the restricted cubic spline with 4 knots for time for the i th cow, LPi is the lactation period of the i th cow, iSCC ij is the individual somatic cell count for the i th cow at the j th measurement day, DMY ij is the daily milk yield of the i th cow at the j th measurement day, MV_Side ij represents the right/left milk vein of the i th cow at the j th measurement day, random(intercept, time) i is the random intercept and random slope for time of the i th cow and e ij is the model error for the ith cow at the j th measurement day. The multivariable LME model formulas regarding the relationships of DMY with each of the 15 echotexture variables and of qSCS with the 15 echotexture variables are presented in Supplementary Materials. The general form of the LME models with a 4-knot RCS that investigated the relationship of DMY with each of the 15 echotexture variables and the relationship of qSCS with the 15 echotexture variables was as follows: where Y ij is the daily milk yield or the quarter somatic cell score of the i th cow at the j th measurement day, RCS(time) i is the restricted cubic spline with 4 knots for time for the i th cow, LPi is the lactation period of the i th cow, Quarter ij is the udder quarter of the i th cow at the j th measurement day, X ij is one of the echotexture variables' measurement of the i th cow at the j th measurement day, random[intercept, RCS(time)] i is the random intercept and random slope for time of the i th cow and e ij is the model error for the i th cow at the j th measurement day. The analytical versions of the models' general equations are presented in the Appendix A. Normality of distribution regarding residuals and random effects of the final multivariable models were graphically checked. All p-values were two-sided, and significance level was set at p ≤ 0.05. All statistical analyses, data manipulation and data visualization were performed using R programming software (version 4.2.2) (R Core Team, 2022). The "lme4" package in R programming language was used to assess convergence of the LME models.

Descriptive Statistics
Repeated measurements from 21 clinically healthy Holstein dairy cows corresponding to 345 cow measurement days were included in our analysis. Changes during the study period in each of the 15 echotexture features used for the investigation of udder parenchymal changes are presented in Supplementary Materials. Morphological and blood flow parameters of the milk vein throughout the study period are presented in Table 1 and Figure 3. Distance of the milk vein from skin surface, its diameter and its cross-sectional area remained constant across production stages.
A significant change was observed in the values of blood flow variables (BFVol, TAMV and Vmax) of the milk vein between late lactation, dry period and early lactation (Table 1, Figure 3). The highest value of mean BFVol was recorded at calving day, at 7.57 ± 2.18 L/min, and the lowest at mid-DP (21st day in the DP). The most abrupt changes in BFVol were observed between the measurements "dry-off day" and "3rd day in DP" (decrease) and between the measurements "3 days before calving" and "calving" (increase). TAMV and Vmax followed the same trend as BFVol throughout the study period. . Figure 3. Changes in blood flow (A-C) and morphological (D-F) parameters of the milk vein during the study period: starting from 7 days before the dry period ("DP"), throughout the DP, 7 and 3 days prior to calving, at calving ("C") and from the 3rd to the 90th day in the new lactation period ("DIM": days in milk).
A significant change was observed in the values of blood flow variables (BFVol, TAMV and Vmax) of the milk vein between late lactation, dry period and early lactation (Table 1, Figure 3). The highest value of mean BFVol was recorded at calving day, at 7.57 Figure 3. Changes in blood flow (A-C) and morphological (D-F) parameters of the milk vein during the study period: starting from 7 days before the dry period ("DP"), throughout the DP, 7 and 3 days prior to calving, at calving ("C") and from the 3rd to the 90th day in the new lactation period ("DIM": days in milk).

Blood Flow Volume, Daily Milk Yield and Somatic Cell Count
The results of the univariable and multivariable LME models investigating the relationship between BFVol in the milk vein, DMY and iSCC are presented in Table 2

Udder Echotexture and Somatic Cell Score
The results of the univariable and multivariable LME models investigating the relationship between udder parenchymal echotexture parameters and udder quarter SCS are presented in Table 4. Table 4. Association between udder quarter somatic cell score and udder echotexture parameters using linear mixed models.

Discussion
Ultrasonography was employed to examine dairy cows' udder parenchyma and blood flow repeatedly from late lactation through to the DP and early lactation. This period is very demanding for the udder, as it transitions from a state of milk synthesis to involution and, through intense growth, back to high milk production [37]. This interval is also crucial for udder health because it is related to high infection risk [38]. In this study, we recorded echotextural changes in the mammary parenchyma, as well as changes in morphology and blood flow features of the milk vein throughout the above period. Furthermore, we investigated the associations between udder parenchymal echotexture, blood flow in the milk vein, somatic cell count and milk yield of Holstein dairy cows.
Venous morphological parameters remained constant throughout late lactation, dry period and early lactation stage, despite significant changes in BFVol (Table 1). Change in BFVol is related to milk synthesis taking place in the mammary gland during late and early lactation. Indeed, in our linear mixed model, BFVol appeared significantly associated both with DMY and with production stage ( Table 2). This is also apparent in the change in BFVol and DMY over time, as BFVol follows the curve of milk production: a rapid decrease within 3 days after the abrupt dry-off, followed by a stationary (non-lactating) period and then a rapid increase in BFVol 1 week before calving (Figure 3), when parenchymal re-growth and milk synthesis take place. Given the calculation formula BFVol = TAMV × 60 × A [22], the increased blood supply during lactation was achieved by increased flow velocities (TAMV, Vmax), while venous diameter and cross-sectional area remained constant. This explains why BFVol, TAMV and Vmax curves throughout the study period presented almost identical trends.
Our results agree with [22], who observed a significant correlation between DMY and BFVol, and with [27], who found that lactating cows had significantly higher BFVol than dry cows. Brown Swiss cows' milk vein diameter recorded in [27] was smaller than Holstein cows' in our study, with no significant difference regarding production status in either study. Even though BFVol of both right and left milk vein did not differ significantly in [27], we included the "side" parameter in our linear mixed model, as it appeared significant. The presence of this parameter in the model formula improves its performance. This is explained by the variability in udder quarter milk yield contribution to the total DMY throughout a lactation period [39], affecting BFVol of each side, too. Finally, mean BFVol values of the milk vein reported by [28] are higher than our findings and those of [26,27]. This could be due to the different technique applied in that study, measuring blood flow in cross-section instead of longitudinal section of the vein.
A significant association was observed in the BFVol-iSCC univariable model, where measurement day (corresponding to production stage) was significant, too. Cows with high iSCC present lower DMY [40][41][42][43], which in turn is significantly associated with BFVol, as discussed above. The lower DMY and its impact on BFVol could possibly explain the reduced BFVol in cows with high iSCC. Further research could use these clinical observations to distinguish cause and consequence in the BFVol-DMY-SCC relationship and to fully understand how subclinical mastitis affects udder blood flow.
Echotexture of the mammary parenchyma is subject to changes from dry-off, throughout the DP and into regrowth and high lactation [44]. Milk quantity, histomorphology of the udder and ultrasound transducer frequency have an impact on mean NPV and PSD of dairy sheep's udder sonograms [45]. Milk quantity affects the echotexture features of dairy cows' udder sonograms, and removal of the hypoechoic milk alters parenchymal echotexture of the mammary gland [11]. In our multivariable model (Table 3), milk yield and production stage were significantly associated with changes in echotexture features (NPV, PSD and Long-Run Emphasis). The variability in milk quantity contained in the mammary parenchyma across different production stages could explain the significant changes in echotexture features throughout the study period. A limitation of our study was the lack of quarter-level DMY data that could aid optimal fitting of the model. Future studies could focus on udder quarter echotextural differences, utilizing quarter milk yield data produced from automatic milking systems to optimize the models' predictive accuracy.
Normally, udder parenchyma is uniformly echogenic, whereas milk appears anechoic. Mastitis milk has high cellular content, leading to increased milk echogenicity. Subclinical mastitis also leads to increased SCC, and, in chronic cases, it may induce changes in the mammary parenchyma [42], resulting in increased presence of echogenic connective tissue [3,4]. In our multivariable model (Table 4), quarter SCS was significantly associated with NPV, PSD, Gradient Mean Value, Homogeneity, Gradient Variance and Entropy, whilst udder quarter and time were significant, too. The authors of [12] studied buffaloes and also found that echotexture features (NPV and PSD) were significantly affected by SCC, lactation stage, scanning plane and udder quarter. These results indicate that echotexture analysis could aid the detection of subclinical mastitis-induced alterations of the mammary parenchyma. To optimize models' accuracy, future research should focus on the impact of udder quarter, examination before/after milking, age/lactation period, ultrasound scanning plane and scanner device technical specifications on udder echotexture models.
Our study attempted to connect the dots between udder echotexture and blood flow in the milk vein with milk production and somatic cell count. Based on these results, echotexture analysis could potentially be applied for early diagnosis of subclinical mastitis in lactation, aid the assessment of parenchymal alterations due to chronic subclinical mastitis and suggest a novel method to detect parenchymal alterations due to subclinical mastitis in the DP. Finally, this comprehensive approach to udder ultrasonography, combined with deep learning algorithms [17], could lead to sophisticated models evaluating udder health and productivity in support of precision dairy farming.

Conclusions
Ultrasonographic examination of the mammary parenchyma and milk veins is a practical, non-invasive method on farm sites. Blood flow parameters in the milk vein change significantly throughout late lactation, dry period and early lactation, corresponding to lac-tation needs, while morphological parameters remain constant. To serve these needs, BFVol in the milk vein increases by 0.25 L/min for every 1 kg increase in DMY, keeping all other variables constant. The presence of high SCC in milk is associated with decreased BFVol in the milk vein. Furthermore, second-and higher-order statistical parameters can provide useful information on udder echotexture. Echotexture parameters significantly associated with DMY are NPV, PSD and Long-Run Emphasis. Echotexture parameters significantly associated with SCC are NPV, PSD, Gradient Mean Value, Homogeneity, Gradient Variance and Entropy. These parameters can potentially be used to detect mastitis-related parenchymal alterations and evaluate udder health and productivity, both at individual cow level and at herd level. This comprehensive approach to ultrasonographic examination of dairy cows' udder echotexture and blood flow can contribute to the data-driven assessment of udder health in support of precision dairy farming.  Data Availability Statement: The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.