Thermotolerance capabilities, blood metabolomics and mammary gland hemodynamics and transcriptomic profiles of slick-haired Holstein cattle during mid-lactation in Puerto Rico

Holstein cattle carrying a prolactin receptor gene mutation (SLICK) exhibit short and sleek hair coats enhancing thermotolerance and productivity compared with wild type-haired Holstein (WT) under tropical conditions. The objectives were to unravel the physiological and molecular mechanisms that confer an advantage to this SLICK genotype in Puerto Rico and determine potential correlations between metabolites and physiological variables. At 160 ± 3 d in milk (DIM) we compared vaginal temperatures (VT) and voluntary solar radiation exposure (VSRE) during 48 h between 9 SLICK and 9 WT Holsteins, whereas a subsample of 7 SLICK and 7 WT were used to assess udder skin temperature, mammary gland hemodynamics and transcriptomics, and blood plasma untargeted metabolomics at a single time point. The SLICK cattle showed lower vaginal temperatures throughout the day and greater VSRE at 1000 h and 1100 h compared with their WT counterparts. Total mammary blood flow was greater in SLICK Holsteins compared with WT. The metabolite 9-nitrooctadecenoic acid was identified as a potential biomarker for MBF; moreover, SLICK cattle had greater amounts of this metabolite in their plasma. Prostaglandin D2 synthase ( PTGS ) was up-regulated in the SLICK mammary gland, while plasma prostaglandin D2 was positively correlated with milk yield and increased in SLICK Holsteins compared with WT. Interestingly, the arachidonic acid metabolism pathway was enriched in the mammary gland transcriptome and perturbed in the blood metabolome in the SLICK Holsteins. In conclusion, SLICK Holsteins exhibited lower body temperatures, greater voluntary solar radiation exposure, enhanced blood supply to the mammary gland, and alterations in genes and metabolites involved in arachidonic acid metabolism at the mammary gland and blood plasma. The usage of the slick-haired Holstein cattle genetics in dairy operations could be a feasible alternative to mitigate the adverse consequences of heat stress.


INTRODUCTION
The dairy industry is the most affected livestock production system by heat stress in the US, with an estimated annual loss around 1 billion dollars (St. Pierre et al., 2003).The thermoneutral zone of dairy cattle has been estimated to be between 5 to 25°C (Young et al., 1981;Berman et al., 1985).With climate change, ambient temperatures can easily surpass 25°C in tropical, subtropical, and even temperate regions (De Rensis et al., 2015), threatening livestock production and food security.When the ambient temperature and relative humidity increase, cows' capacity to efficiently dissipate body heat to the environment is impaired, resulting in increased core body temperature and a reduction in feed intake and milk production (Kadzere et al., 2002;West, 2003).Laporta et al. (2020) estimated that dairy cattle in southern US are exposed to heat stress during 33 to 60% of the year, with a temperature-humidity index (THI) ≥ 68.However, the aforementioned study did not include environmental data from Puerto Rico, a tropical US territory, where such threshold in THI is chronically exceeded, even during the winter season (Contreras-Correa et al., 2018).
Management practices such as environmental modifications, nutritional schemes, and selection for adapted Thermotolerance capabilities, blood metabolomics and mammary gland hemodynamics and transcriptomic profiles of slickhaired Holstein cattle during mid-lactation in Puerto Rico breeds have been suggested to alleviate the adverse effects of heat stress on lactating dairy cattle (Beede and Collier, 1986).Selection for adapted breeds is an attractive option where the introgression of heat tolerant genes from zebu into susceptible Bos taurus dairy cattle populations has been proposed; but the negative traits of indicine cattle inherited by the crossbreds may limit its use in the dairy industry (Madalena et al., 1990).Nevertheless, Littlejohn et al. (2014) reported the SLICK1 as a prolactin receptor (PRLR) mutation in Senepol cattle which consists of a single base pair deletion in exon 10 introducing a premature stop codon (p.Leu462*), truncating the PRLR protein.The Puerto Rican Holstein cows inheriting this gene mutation are known as "Pelonas" exhibiting a sleek and very short hair coat phenotype resulting from the SLICK1 mutation.Such mutation was presumably inherited from the Criollo cattle that was adapted to the island environmental conditions after the Spanish colonization (Sánchez-Rodríguez, 2019a); while the SLICK1 allele has been deliberately introduced into Holsteins by crossbreeding these animals with Senepol bulls in Florida and California (Olson et al., 2003;Dikmen et al., 2014;Carmickle et al., 2022).An additional 5 SLICK mutations have been described in cattle from Central America and South America (Porto-Neto et al., 2018;Florez-Murillo et al., 2020).
Regardless of its origin, the slick-haired Holstein cows (SLICK) better tolerate the adverse effects of heat stress, which is reflected in lower vaginal temperatures and respiration rates (Dikmen et al., 2008;Castro et al., 2015;Sánchez et al., 2015) when compared with wild-type haired (WT) Holsteins.Also, Puerto Rican SLICK Holsteins have larger sweat glands than their WT counterparts (Contreras-Correa et al., 2017;Sánchez-Rodríguez, 2019b), which may allow for greater heat dissipation through evaporation.Recently, Sosa et al. (2022a) did not observed differences in sweat gland size and reported limited differential gene expression in skin and liver samples from SLICK and WT heifers (Sosa et al., 2022b); however, examination of lactating animals has been unexplored.Heifers and lactating cattle have indisputable differences in energy requirements and metabolism, thereby the mammary gland gene expression and blood metabolite of lactating SLICK Holsteins has yet to be evaluated to elucidate differences in lactogenesis and metabolism between genotypes.Lastly, SLICK cattle have shown superior lactation performance (Dikmen et al., 2014;Contreras-Correa et al., 2016), thereby these animals are expected to have higher demand of nutrients and oxygen by the mammary gland tissue (Berger et al., 2016) which can be determined by examining mammary hemodynamics.
Anecdotally, the SLICK cattle endure more time grazing under direct solar radiation than their WT counterparts; but empirical data to support this claim is limited (Sánchez-Rodríguez and Domenech-Pérez, 2021).In addition, besides a readily apparent heat dissipation benefit through a shorter hair coat to mitigate heat stress, the physiological mechanisms that allow SLICK cattle to have superior thermotolerance and lactation performance under tropical conditions remain unclear.We hypothesized that SLICK Holstein cattle have greater mammary gland blood perfusion, increased transcripts related to angiogenesis, and increased metabolites that act as vasodilators as part of the physiological and molecular mechanisms that confer superior performance.Thus, the objectives of this study were to compare thermotolerance capabilities, as well as the mammary gland hemodynamics, mammary gland transcriptome, and blood plasma metabolome of mid-lactation SLICK and WT Holstein cattle during the summer in Puerto Rico.In addition, we aimed to identify potential correlations between blood metabolites and physiological variables of interest such as milk yield, body temperature, and mammary gland hemodynamic variables in Holstein cattle.

Animals, Housing and Management
This experiment was conducted at the University of Puerto Rico Agricultural Experiment Station (Lajas, P.R.) from July 18th to August 5th, 2022.All experimental procedures were approved by the University of Puerto Rico, Mayagüez campus (UPRM) Institutional Animal Care and Use Committee.All lactating cows at the herd were genotyped for the SLICK1 mutation using the TaqMan based assay as previously described by Sosa et al. (2021).Nine WT (with no SLICK1 allele present) and 9 SLICK (heterozygous for the SLICK1 allele mutation) Holstein cows were selected from the herd based on parity, DIM, and BW and enrolled in the study.The WT and SLICK groups were balanced by parity (1.89 ± 0.31 and 2.00 ± 0.33; P = 0.81), DIM (159.44 ± 14.34 and 160.33 ± 11.37;P = 0.96),and BW (516.76 ± 14.18 and 530.01 ± 14.18 kg; P = 0.52).Cows were milked twice daily (0400 and 1400 h) in a double-6 herringbone parlor with artificial shade.Approximately 1 h before each milking, cows were moved to a barn with artificial shade and no walls, where concentrate feed was provided, and udders were washed using a water hose.In between milkings, all animals had access to a grazing paddock with available natural shade and ad libitum water.The grazing pas-Contreras-Correa et al.: Slick-haired Holstein physiology ture consisted of pangola grass (Digitaria eriantha) and all animals were grouped together in the same pasture.

Environmental conditions
The air temperature and relative humidity were recorded utilizing the HOBO Pro v2 temp/RH Data Loggers (Onset Computer Corporation; Bourne, MA, USA) placed in the barn (n = 2) and the grazing paddock (n = 2).Data loggers were installed approximately 2 m above the ground.Air temperature and relative humidity data were used to calculate the thermal humidity index (THI) using the formula previously published by Vitali et al. (2009): Where AT is air temperature in °C and RH is relative humidity expressed as a decimal.The THI was determined every 15 s and averaged by hour for descriptive statistics of the study period (Figure 1A).Similarly, light intensity data loggers (HOBO Pendant MX Temperature/Light Data Logger, Onset Computer Corporation, Bourne, MA) were used to record environmental light intensity values every 15 s in the barn (n = 3) and in the grazing paddock (n = 3).All data loggers were placed with the recording surface pointing upward at approximately 2 m high from the floor level.As previously validated by Sánchez-Rodríguez and Domenech-Pérez (2021), the obtained light intensity data was used as a solar radiation index.The light intensity values recorded throughout the study were further averaged by hour for descriptive statistics (Figure 1B).

Vaginal Temperatures and Voluntary Solar Radiation Exposure
Vaginal temperatures were recorded similar to Sánchez et al. (2016) and Contreras-Correa et al. (2021).Briefly, TidbiT v2 Water Temperature Data Loggers (Onset Computer Corporation; Bourne, MA) were attached to a non-progesterone CIDR (Zoetis, Exton, PA, USA) via a plastic tie wrap.The CIDR and data logger devices were washed in water and soap, rinsed in a chlorhexidine solution (Durvet Inc., Blue Springs, MO, USA), and lubricated (OB Lube; AgriLabs, St. Joseph, MO, USA) immediately before vaginal insertion with a CIDR applicator.For the insertion, the cows' tails were restricted, and the vulvas were washed with soap and water, rinsed with the chlorhexidine solution, and dried with a clean cloth towel.Immediately after, each cow received a nylon collar (Nasco, Fort Atkinson, WI, USA) with an attached HOBO Pendant MX Temperature/Light data logger (Onset Computer Corporation, Bourne, MA, USA) in the superior part of the collar and a weight (132 g) in the inferior part to keep the logger in place on top of the cows' necks.These light sensors data loggers were previously validated by Sánchez-Rodríguez and Domenech-Pérez (2021) as viable indicators of solar radiation exposure in dairy cattle.Data loggers were inserted and secured in each cow 24 h before data collection began and programed to record vaginal temperatures and light intensity every 15 s for 2 consecutive days.The obtained data were then averaged by hour for further statistical analysis.

Color Doppler Ultrasonography and Udder Skin Temperatures
Data were collected from all 18 cows, but 4 cows (2 WT and 2 SLICK) presented anormal mammary hemodynamics, thereby a subsample of 14 animals (7 WT and 7 SLICK) were used for the mammary gland hemodynamics, udder skin temperatures, mammary gland biopsies, and blood sample collection.These animals had no differences in parity (1.57± 0.29 vs. 2.28 ± 0.29; P = 0.11), DIM (157.7 ± 14.36 vs. 163.4 ± 14.36; P = 0.78) or BW (503.9 ± 15.6 vs. 529.7 ± 15.6 kg; P = 0.29) between WT and SLICK cows, respectively.However, differences in milk yield were observed after analyzing the daily milk productions collected by the Dairy Herd Improvement 5 d before enrollment in the study (19.74 ± 1.24 kg/d vs. 23.77± 1.24, in the WT versus SLICK cows, respectively; P = 0.04).The mammary gland hemodynamics were determined via the color function on the Doppler ultrasound (M-Turbo, Sonosite Inc., Bothell, VA, USA) using a trans-rectal probe (Linear Endorectal L52x probe, Sonosite Inc., Bothell, VA, USA).The left and right external pudendoepigastric arteries blood flow was assessed similar to Gotze et al. (2010) and Sánchez-Rodríguez et al. (2012).Briefly, the abdominal aorta was identified with the transducer probe facing dorsally and moving caudally toward the external iliac arteries.In each side of the cow (left and right), following each external iliac artery, ventrally the femoral artery was found and the external pudendoepigastric artery was located immediately after.The ultrasound transducer probe was aligned to the pudendoepigastric arteries, and 2 to 3 independent cardiac cycle waveforms were used to calculate systolic velocity (Vmax; cm/s), diastolic velocity (MDV; cm/s), s:d ratio (Vmax/MDV), pulsatility index (PI), and resistance index (RI) using preset functions on the Doppler ultrasound.Additionally, the ellipse preset function in the ultrasound was used to determine the pudendoepigastric arteries cross sectional diameter, area, and circumference.The mean velocity was calculated using the equation: MnV = (Vmax -MDV)/PI.Blood flow was calculated using the equation: BF = MnV*vessel area*60 s.These formulas were previously published by Lemley et al. (2018).Total mammary blood flow is reported as the summation of the left and right pudendoepigastric arteries blood flow.Animals were divided into 2 groups, evenly represented by both genotypes and ultrasound examinations were carried out in 2 consecutive days between 0700 to 0900 h lasting approximately 20 min per animal.Researchers were aware of treatment allocation during ultrasonography because of the outward phenotypic appearance; however, researchers were blinded to treatment allocation during ultrasound data calculations and analysis.Cows were evaluated at a random order.Immediately after performing each transrectal ultrasound, udder images were captured by a single operator using a FLIR T540 professional thermal camera (Teledyne FLIR LLC, Oregon, USA) to assess udder skin temperature.The operator stood approximately 1.0 m behind the cow, the udder was centered, and the thermograph of both hindquarters were taken with the camera autofocus.The camera resolution was 464 × 348 pixels, and the accuracy was ± 2°C.Images were evaluated and the udder skin temperature were obtained by utilizing the FLIR Tools+ software (Teledyne FLIR LLC, Oregon, United States).The udder skin temperature reported is the average of the left and right hindquarters.

Physiological Variables Statistical Analysis
The UNIVARIATE procedure of SAS software version 9.4 (SAS Institute, Cary, NC) was used to determined data normality using the Shapiro-Wilks statistics.Each animal was considered as the experimental unit.The GLIMMIX procedure of SAS was used with vaginal temperature and voluntary solar radiation exposure as dependent variables.In these models, time of the day and hair coat type were included as fixed effects, while the cows' identification numbers were used as a random effect.The interaction between time of the day and hair coat type was also tested.The VT and VSRE were considered repeated measurements.The mammary gland hemodynamics, udder skin temperature, and milk yield data were analyzed by one-way ANOVA using the MIXED procedure of SAS.The dependent variables were resistance index, pulsatility index, diameter, total mammary blood flow, udder skin temperature, and milk yield.The hair coat type was the independent variable and the day and order the animals were sampled were added in the model as covariates.Lastly, the CORR procedure of SAS was used to determine the Pearson correlation coefficients between milk yield and mammary gland blood flow.Statistical significance was declared at P < 0.05 and tendencies at 0.05 < P < 0.10.

Mammary Gland Biopsies and RNA Extraction
Mammary biopsies were performed following the protocol published by Daley et al. (2018) with some minor modifications.Briefly, the biopsies were guided using a color Doppler ultrasound with a linear probe to avoid major blood vessels and to prevent hemorrhaging or development of mastitis.The cows underwent standing sedation by injection of 0.01 mg/kg of xylazine hydrochloride (20 mg/mL) into the coccygeal vein.Surgical clippers were used to remove the hair from the left hindquarter and the area was sterilized with povidoneiodine 7.5% surgical scrub and 70% ethanol.The biopsy site received local anesthesia by administrating 5 mL of 2% lidocaine hydrochloride subcutaneously.A 1-to 2-cm vertical incision was performed using a no.22 scalpel blade through the mammary skin and subcutaneous tissue avoiding any major blood vessels.A reusable core biopsy instrument (MG1522; BD, Franklin Lakes, NJ, USA) was used with a sterile disposable 12 g x 16 mm biopsy needle (MN1216; BD, Franklin Lakes, NJ, USA) to collect the mammary tissue sample.The tissue weight was between 15 to 20 mg and 2 samples were collected from the same incision.The mammary tissue was washed in 1x phosphate-buffered saline, snap frozen in liquid nitrogen, and stored at −80°C until further analysis.The incision was closed using a medical grade skin stapler with 35 wide stainless-steel staples (AHS35W; American Health Service, Lake Forest, IL, USA) and aerosol bandage (4% aluminum powder) was sprayed over the closed incision.The described protocol took approximately 10 to 15 min per animal.After the biopsy collection, all animals received 1.1 mg/kg of flunixin meglumine (50 mg/mL) intravenously and 6.6 mg/kg of ceftiofur crystalline free acid (200 mg/mL) subcutaneously into the posterior aspect of the ear as preventive care.Animals were monitored for signs of pain and inflammation twice daily and continued on their regular milking schedule after a 48-h withdrawal period.The staples were removed 12 d post-procedure.No adverse effects were observed.
Total RNA was extracted from the frozen mammary tissue by homogenizing approximately 15 mg of tissue in 350 mL of RLT Lysis buffer (QIAGEN, Hilden, Germany) and RNA was purified using the RNeasy Mini Kit (QIAGEN, Hilden, Germany) following the manufacturer's recommendations.The quality of the total RNA extraction was assessed and quantified using a NanoDrop One spectrophotometer (Thermo Scientific, Waltham, MA).All RNA samples were diluted to 100 ng/µL and 20 µL were shipped on dry ice for unstranded RNA-seq library creation and Illumina Novaseq 2x150bp sequencing (Novogene Corporation Inc., Sacramento, CA).The RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).Samples with an RNA Integrity number (RIN) of ≥7.0 were kept for library construction.

Differential Expression and Gene Set Enrichment Analysis
To quantify the transcript level expression, the raw reads were aligned using salmon [v1.9.0] (Patro et al., 2017) to the Bos taurus (ARS-UCD1.3)RefSeq transcripts with the genome used as decoys.The transcript expression was converted to gene-level expression using the R package tximport [v1.26.1] (Soneson, et al., 2016) using the edgeR scaling methods provided in the documentation.Genes with very low expression (fewer than 10 reads found in more than 7 samples) were removed (Bourgon et al., 2010).The likelihood ratio method in edgeR [v3.40.2] ( McCarthy et al., 2012) was used to test for significant differential expression (DE).The fry method from edgeR was used to do gene set enrichment analysis (GSEA) for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms.All code used for DE and GSEA can be found in https: / / github .com/IGBB/ slick .rna-seq.

Untargeted Metabolomics Analysis
At the time of mammary gland biopsies, blood was collected via coccygeal venipuncture using K2 EDTA (EDTA) blood tubes (BD Vacutainer, Franklin Lakes, NJ, USA) and immediately placed on ice.Blood samples were centrifuged at 2,000 x g at 4°C for 12 min and plasma was collected and transferred into 1.5 mL tubes and stored at −80°C until metabolome extraction.The untargeted metabolomics was performed by Creative Proteomics located in Shirley, NY, USA.Based on the developed protocol, the metabolites were separated on an Acquity UPLC HSS T3 column (1.8 µm,100 × 2.1 mm) using a Vanquish Flex UPLC combined with Q Exactive mass spectrometer with heated electro-spray ionization (Thermo Fisher Scientific, Waltham, MA, USA).The mobile phase was composed of solvent A (0.05% formic acid water) and solvent B (acetonitrile) with a gradient elution (0-1 min, 5% B; 1-12 min, 5-95% B; 12-13.5 min, 95% B; 13.5-13.6min, 95-5% B; 13.6-16 min, 5% B).The flow rate of the mobile phase was 0.30 mL/min.The column temperature was maintained at 40°C, and the sample manager temperature was set at 4°C.The mass range was from m/z 70 to 1050.The resolution was set at 70,000 for the full scan MS and 17,500 for the dd-MS2.The mass spectrometry parameters were as follows: heater temp, 300°C; sheath gas flow rate, 45 arbitrary units; auxiliary gas flow rate, 15 arbitrary units; sweep gas flow rate, 1 arbitrary unit; and capillary temperature, 350°C.The spray voltage in positive ion mode (ESI + ) was 3.0 kV and the S-Lens RF level was 30%, whereas the spray voltage in negative ion mode (ESI -) was −3.2 kV and the S-Lens RF level was 60%.The raw data were acquired and aligned using the Compound Discoverer (3.0, Thermo) based on the m/z value and the retention time of the ion signals.Ions from both ESI + and ESI -were merged and imported into the SIMCA-P program (version 14.1) for multivariate analysis.The metabolites were filtered and confirmed by combining the results of the variable importance projection (VIP) and from the unpaired ttest.Significant difference in metabolites were declared with VIP >1.5, │Log 2 FC│ > 1.0, and P value <0.05.The Human Metabolome Database (https: / / www .hmdb.ca/ ) was utilized to determine the metabolites chemical structures using the data of accurate masses and MS/MS fragments.Mean values of metabolite contents from biological replicates of group SLICK and group WT were used to calculate the metabolite ratio.After log-transformation of the data, median centered ratios were normalized.All significant metabolites were imported to KEGG databases and MetaboAnalyst to obtain the categorical and biological annotations, including pathways and enzyme interactions.Spearman correlations between untargeted metabolomics and physiological variables were performed in R (v4.2.0) (R Core Team 2020) using rcorr in the Hmisc package (v5.0-1)(Harrell 2023) and visualized using ggplot2 (v3.4.0) (Wickham 2016) and cowplot (v1.1.1)(Wilke 2020).Significant differences were declared at P < 0.05.

Vaginal Temperatures and Voluntary Solar Radiation Exposure
Figure 1C depicts the vaginal temperatures in each genotype over a 24-h period.The time by hair coat interaction did not influence vaginal temperatures (P = 0.14).However, vaginal temperature was affected by the hair coat type (P < 0.01) and time (P < 0.0001).The SLICK cows presented lower vaginal temperatures than their WT counterparts, with daily averages of 39.08 ± 0.09 and 39.45 ± 0.10°C, respectively.Similar to the current study, others have observed simple effects of hair coat type and time, but not an interaction between both variables, for body temperature during hot season in Puerto Rican dairy cattle (August; Sánchez-Rodríguez, 2019b) and in Holstein cows derived from Senepol x Holstein crossbreds in Florida (July and August; Dikmen et al., 2014).However, the Sánchez-Rodríguez (2019b) study did not observe a significant main effect of hair coat type or an interaction between hair coat type and time when cows were evaluated during the cool season of the year (March).He also observed, during the warm season (April), a significant interaction between hair coat type and time affecting the vaginal temperatures.Thus, it appears that the thermoregulatory advantages associated with the SLICK phenotype are magnified as the environmental conditions worsen.Perhaps, this is the reason why Carmickle et al. (2022) observed that slick-haired Holstein heifers were able to better regulate their body temperature in Florida, but not in California.Finch (1986) noted that the longer hair coats of WT Bos taurus cattle may get saturated with water vapor from sweat, limiting the gradient for heat dissipation through evaporation under hot and humid environmental conditions.Thus, having a shorter hair coat, as previously reported by Jiménez-Cabán et al. (2015) on Puerto Rican SLICK Holstein cows, may allow for a greater heat transfer into the environment.Additionally, SLICK cattle possess other adaptations to hot weather, including larger sweat glands (Contreras-Correa et al., 2017;Sánchez-Rodríguez, 2019b) and sweating capacities (Dikmen et al., 2014) than their WT counterparts.Also, in the current study a minimum overall vaginal temperature daily nadir was observed between 0700 to 0800 h, averaging 38.68 ±  Ehrenreich and Bjugstad (1966) and Lyons and Machen (2000), environmental conditions conducive to heat stress decrease or eliminate livestock's grazing time during daily hours with high solar radiation or thermal humidity index (0800-1800 h).Under such environmental conditions, cattle spend most of this time resting and ruminating (Reppert, 1960;Lyons and Machen, 2000), under the shade (Carvalho de Oliveira et al., 2021) and close to available waterers (Osei-Amponsah et al., 2020).However, tropically adapted cattle may be able to spend a greater proportion of their daily time grazing, even when exposed to direct solar radiation, than their less-adapted counterparts.For instance, Hammond and Olson (1994) reported that Senepol cows spend a greater daily portion grazing than Hereford cows.Sprinkle et al. (2000) observed that crossbreed Angus x Brahman cows spent shorter daily periods under shade, as well as longer time grazing and directly exposed to solar radiation, when compared with purebred Angus cows.Also, Cooke et al. (2020) mentioned that Bos indicus cattle graze longer and maintain a greater feed intake under hot environmental conditions.In this regard, there is previous anecdotal evidence ( Sánchez-Rodríguez et al., 2019b) reporting Puerto Rican SLICK dairy cows grazing under direct solar radiation while their WT counterparts rest under the available shade.Thus, the current data indicates that SLICK cows may have returned to direct solar radiation between 0900 and 1200 h, while their WT counterparts remained under shade.When Figures 1B and 1D are visually compared, the environmental data loggers recorded light intensity values at the grazing paddocks considerable higher than those observed on the cows' necks between 0800 and 1700 h, which indicates that cows' solar radiation exposure occurred intermittently.It is worth mentioning that such solar radiation exposure period was interrupted when all cows were brought to the barn for concentrate feeding and udder washing before the afternoon milking.Thus, it is not possible from the current study to determine voluntary solar radiation exposure between 1200 to 1600 h.

Mammary Gland Hemodynamics and Udder Skin Temperatures
As demonstrated by Götze et al. ( 2010) and Mordhosrt et al. (2017), and in agreement with the current study, there is no difference in the mammary gland hemodynamics between the left and right external pudendoepigastric arteries.Therefore, the results presented in the current study are the average between both arteries, with the exception of total MBF, which is reported as the summation of the left and right pudendoepigastric arteries.The RI was lower (P < 0.05) in SLICK compared with WT Holsteins (0.56 ± 0.02 vs. 0.64 ± 0.02, respectively, Figure 2A); whereas the PI tended (P = 0.06) to be lower in SLICK compared with WT Holsteins (0.94 ± 0.07 vs.1.15± 0.07, respectively).The RI and PI are indicators of peripheral vascular resistance and are inversely proportional to blood flow.Evaluation of mammary gland RI and PI are scarce in the literature in relation to dairy cattle studies.Using Angus and Angus x Simmental crossbreds, Kennedy et al. (2019) reported at d 44 of lactation a RI and PI of 0.61 ± 0.02 and 1.04 ± 0.06, respectively.In the current study, the average diameter of the left and right pudendoepigastric arteries was greater in the SLICK compared with their WT counterparts (1.31 ± 0.07 vs. 1.08 ± 0.07 cm, respectively; P = 0.04; Figure 2B).These values conflict with those previously reported by Götze et al. ( 2010) who observed a pudendoepigastric trunk diameter of approximately 1.62 cm in lactating Holstein cows.Total MBF was increased in SLICK compared with WT cattle (9.03 ± 1.11 vs. 5.14 ± 1.11 L/min; P = 0.03; Figure 2C).Utilizing electromagnetic blood flow probes, the udder half mammary blood flow in Holstein cattle under thermal stress was 4.5 L/min versus 5.1 L/min in similar animals under thermal comfort (Lough et al., 1990).Moreover, these authors found a tendency in mammary blood flow to be lower in heat stressed cows with ad libitum diets and in cows under thermal comfort with restricted daily feed intake compared with thermal comfort counterparts with ad libitum feed intake, suggesting that the effects of Contreras-Correa et al.: Slick-haired Holstein physiology thermal stress altering mammary blood flow might be mediated by reductions in dry matter intake (Lough et al., 1990).It is important to note that even though the current study was not designed to measure DMI, the SLICK Holsteins showed a greater voluntary direct solar radiation exposure during approximately 3 h/d while their WT counterparts sought shade (Figure 1D).In cattle, exposure to direct solar radiation during the hottest periods of the day is commonly associated with grazing, while resting under the shade is not.Therefore, the lower MBF observed in the WT cows in the current study might be related to decreased DMI.
Interestingly, Götze et al. ( 2010) and the current study employed Doppler ultrasonography to measure MBF, but they reported a range of 19.9 to 27.9 L/ min in total MBF which represents a 3-fold greater blood flow compared with our findings.Lactating dairy cows have been reported to have a total cardiac output of 45 to 55 L/min (Davis et al., 1988), which implies that in the animals used at the Götze et al. ( 2010) study approximately 50% of their cardiac output was perfusing the mammary gland.Furthermore, Davis et al., (1988) reported MBF (half-udder) range from 3.3 to 4.9 L/min and determined 14 to 19% of cardiac output was perfusing to the whole udder.These absolute blood flow values, determined via electromagnetic blood flow probes, are similar to our own results reported via Doppler ultrasonography measurements.The cardiac output and MBF values appear to have great variability and inconsistency, especially when you take into consideration if these animals were pregnant and the importance of nutrient partitioning.It is possible that the differences observed in MBF in our study compared with those in Götze et al. ( 2010) are due to differences in milk production since the SLICK and WT Holsteins were producing 23.8 and 19.8 kg/d during mid-lactation and in the Götze et al. ( 2010) study, Holstein cows were producing an average of 31.2 kg/d during peak lactation.An increase in milk yield results in a higher demand of nutrients and oxygen by the mammary gland tissue (Berger et al., 2016).In theory, these requirements are met by enhancing blood supply to the mammary gland which could be attained by increasing the diameter of the pudendoepigastric arteries and/or reducing its resistance (Berger et al., 2016).Based on our findings, it seems that the aforementioned hemodynamic mechanisms help to explain the differences in milk production observed between SLICK and WT Holsteins (Figure 2D).There was a tendency (P = 0.06) for the MBF relative to BW to be greater in the SLICK compared with the WT Holsteins (16.75 ± 2.14 vs. 10.38 ± 2.14 mL/min/kg).This trend might implicate a higher demand of blood flow by the SLICK mammary gland due to its greater metabolic activity and milk synthesis.
Previous studies have described a positive correlation between milk production and mammary gland blood flow in sheep, goats, and cattle (Davis and Collier et al., 1985;Nielsen et al., 1990;Götze et al., 2010).In the current study, a tendency for a positive correlation between MY and MBF was observed (r = 0.51, P = 0.07).Lastly, as a skin gland, the mammary gland might assist in the dissipation of body heat when cows are exposed to heat stress conditions (Lough et al., 1990).Nevertheless, no differences were observed in udder skin temperature between SLICK and WT Holsteins (P = 0.67; Figure 2E-G).In this regard several researchers from different disciplines (Hoffman et al., 2013;Church et al., 2014) have reported limitations of infrared thermography for the assessment of surface temperature due to the considerable influence that the surrounding environmental conditions exert on the obtained values.Thus, precautions should be taken when interpreting the obtained results.

Differential Expression and Gene Set Enrichment analysis
One sample from the WT group was removed from the differential gene expression analysis since it did not meet the RIN criteria.The remaining 13 samples (WT, n = 6; SLICK, n = 7) were sequenced to an average depth of 31.4M reads (ranging from 29.4M in WT_2273 to 40.4M in WT_2226).Salmon mapped the raw reads to the bovine transcriptome with an average rate of 87.62% (81.12% in SLICK_2171 to 90.13% in SLICK_2287).The 76,278-reference transcript-level expressions for each sample were reduced to 20,946 gene-level expressions.Of these, 13,392 passed the low expression filter.Figure 3 shows the multidimensional scaling (MDS) plot of the samples.
A single gene (GeneID: 513329, major allergen Equ c 1) showed significant (B.H. adjusted P-value ≤0.05 and │Log 2 FC│ > 2.0) differential expression (Table 1).We aligned the protein sequence for this gene using the BLASTp against the database NR.Interestingly, the top hits that were not equine allergen, were found to be salivary lipocalin 2 (SAL2; UFI08139.1)and salivary lipocalin-like (XP_043309444.1).The major allergen Equ C1, SAL2, and salivary lipocalin-like are known to be part of the lipocalin family; which are still poorly understood, but have been found to act as allergens, pheromone carriers, and be involved in enzymatic synthesis of prostaglandins (Flower, 1996).The SAL2 and salivary lipocalin-like proteins are pheromones found in Ovis aries and Cervus canadienses, respectively, and are known to affect sexual behavior of females.The aforementioned 2 proteins are highly similar to the major allergen Equ c 1, including the exact length of 182 amino acids.It is possible that the major allergen Equ c 1 found to be differentially expressed in the mammary gland of SLICK cattle is wrongly annotated and its function could be more similar to a pheromone that facilitates the mother-young interaction.However, since lipocalins are also involved in the enzymatic synthesis of prostaglandins (Gerena et al., 1998;Inoue et al., 2008), it is possible that the major allergen Equ C1 function is related prostaglandin D2 synthase (PTGS) which was also differentially expressed in the SLICK animals using a less stringent criteria (P < = 0.01 and │Log 2 FC│ > 0.6) similar to the one used by Sosa et al. (2022a).
Both the KEGG pathway and GO term GSEA showed no statistically significant enrichment at B.H. adjusted P-value ≤0.05.Nevertheless, using a less stringent criteria with a normal P-value ≤0.05, multiple relevant pathways that deserve further discussion were identified.For example, after enrichment, the KEGG pathway analysis revealed that the SLICK Holsteins had 44 genes enriched in arachidonic acid metabolism (bta00590) and 146 genes depleted in tight junction function (bta04530) (Table 2).Integrating the results observed from increased MBF and milk production in SLICK Holsteins, we can speculate that upregulation of genes involved in the arachidonic acid metabolism could increase the prostaglandin synthesis leading to vasoconstriction and vasodilation of blood vessels supplying and regulating nutrients to the mammary gland.Even though lactose is the primarily osmolar component of milk (Lin et al., 2016), the downregulation of genes involved in tight junctions could result in increased paracellular transfer of ions into the lumen, increasing the osmotic gradient and driving more water into milk.Based on our findings, future studies should evaluate milk composition between WT and SLICK Holsteins since this was a limitation of the current study.Lastly, among the enriched GO analysis, interesting genes were found to be downregulated in growth hormone receptor signaling pathway via JAK-STAT (GO: 0060397), hair follicle development (GO: 0001942), and vasoconstriction (GO: 0042310) in the SLICK Holsteins (Figure 4).

Untargeted Metabolomics and Physiological Correlations
A total of 471 metabolites were found in SLICK and WT Holstein samples using total ions chromatograph in positive and negative ion modes.According to the screening criteria of VIP >1.5, │Log 2 FC│ > 1.0, and P value <0.05, 35 metabolites were found to be differentially abundant in SLICK plasma relative to their WT counterparts (Table 3).The differential sequential metabolites were further imported into the KEGG pathway database (https: / / www .kegg.jp) to obtain the categorical annotations.After metabolic pathways enrichment, it was found that the SLICK Holsteins exhibited perturbed metabolites related to aminoacyl-tRNA biosynthesis (P = 0.01), tryptophan metabolism (P = 0.01), and arachidonic acid metabolism (P < 0.001).It is important to note that the arachidonic acid metabolism pathway was also enriched in the SLICK mammary tissue when we analyzed its transcriptome.Most notable were decreased prostaglandin A2 in SLICK versus WT, while prostaglandin D2, E2, and E3, as well as thromboxane B2 were all increased in SLICK versus WT cows (Table 3).Interestingly, when the mammary gland transcriptomic was evaluated in the current study, a concomitant upregulation of the prostaglandin D2 synthase in the SLICK cattle was observed.Additionally, plasma oxylipin profiles have been previously reported in the transitioning dairy cow (Putman et al., 2019) and in bovine milk and plasma at different stages of lactation (Kuhn et al., 2017).Putman et al. (2019) observed maximal concentrations of arachidonic acid 2 d after cessation of lactation, while prostaglandin E2 concentrations reached a nadir 12 d after dry-off.Kuhn et al. (2017) concluded during different stages of lactation that plasma oxylipin concentrations differed from milk oxylipin profiles, therefore indicating local mammary gland changes independent of systemic responses.Nonetheless, the current study examined associations between the blood metabolome and phenotypic responses in both SLICK and WT animals.The Spearman correlations between the untargeted plasma metabolome and physiological variables revealed potential biomarkers for MY, VT, RI, diameter, and total MBF (Figure 5).Arachidic acid (r = −0.77;P < 0.01), palmitoylethanolamide (r = −0.72;P < 0.01), and 2-dodecylcyclobutanone (r = −0.71;P < 0.01) showed strong negative correlations with MY.Moreover, 2-dodecylcyclobutanone was negatively correlated with pudendoepigastric artery diameter (P = 0.02) and total MBF (P = 0.04).Malic acid, 18-hydroxyarachidonic acid, and tranexamic acid were positively correlated with RI (P < 0.05) and negatively correlated with the pudendoepigastric artery diameter and total MBF (P < 0.05).Furthermore, malic acid and 18-hydroxyarachidonic acid were greater in the plasma from WT Holsteins (P ≤ 0.03; Table 3) which also exhibited decreased MBF.9-nitrooctadecenoic acid showed a strong positive correlation with MBF (r = 0.78; P < 0.001), while SLICK cattle exhibited an increase of this metabolite in their plasma (P = 0.01).4,5-dihydroorotic acid (r = 0.82; P < 0.001) and 1H-1,2,4-triazole-3-carboxamide (r = 0.71; P < 0.01) showed a strong positive correlation with MY.Prostaglandin A2 was negatively correlated with diameter and total MBF, while Prostaglandin D2 was positively correlated with MY (P < 0.05; Figure 5).Deoxycholylserine (r = −0.73;P < 0.01), lysoPC(20:0/0:0) (r = −0.68;P < 0.01), and N-stearoyl phenylalanine (r = −0.70;P < 0.01) exhibited strong negative correlations with VT; surprisingly, all of these metabolites were increased (P < 0.05) in the SLICK Holstein plasma compared with the WT.The associations between these metabolites and temperature regulation are unclear; however, based on our current results future targeted lipidomic research may reveal novel systemic pathways in dairy cattle experiencing heat stress.

CONCLUSION
Based on the average THI of 76.6, animals in the current study were chronically exposed to heat stress conditions throughout the day.However, the SLICK Holstein cows exhibited lower vaginal temperatures, even with a greater voluntary solar radiation exposure during the most environmentally stressful hours of the day.Such prolonged voluntary exposure to direct solar radiation may indicate a longer grazing period compared with their WT counterparts, which could partially help to explain the greater daily milk production commonly observed in SLICK cows.Interestingly, the pudendoepigastric arteries exhibited lower resistance index, as well as a greater diameter which resulted in a 76% increase in mammary gland blood flow in the SLICK compared with WT Holsteins.Such differences further help to explain the greater milk yields associated with the SLICK cows.The mammary gland transcriptome and blood metabolome revealed genes and molecular mechanisms of interest that deserve further exploration such as the arachidonic acid metabolism pathway.Future research should evaluate DMI and milk composition between both genotypes and include a greater sample size since this was a limitation of the current study.Findings on the current study will allow other researchers to explore the usage of the metabolites mentioned in this study as potential biomarkers for lactating Holstein cattle physiology.and light sensor data loggers and related materials were funded by the U.S. Department of Agriculture National Institute of Food and Agriculture, Hatch initiatives .The authors would like to extend appreciation to Dr. Katherine Domenech-Pérez for providing the animal's genotype and Dr. Michael Pesato, DVM for providing the drugs for the mammary biopsies' procedures, and Zoetis for donating the progesterone free CIDR's used to measure vaginal temperatures.  .Spearman correlations between metabolites in blood plasma and physiological variables measured in both hair coat types in the current study such as milk yield (MY), vaginal temperatures (VT), left and right pudendoepigastric arteries resistance index (RI), left and right pudendoepigastric arteries diameter (Diameter), and left and right pudendoepigastric arteries total blood flow (Total_MBF).Circle size and color intensity represent the spearman correlation between the metabolite and variable; blue is positive correlation, red is negative.A silver border highlights significant correlations (P < 0.05), gold highlights high significance (P < 0.01).A total of 14 samples were used for these correlations from WT (n = 7) and SLICK (n = 7) cattle.
Contreras-Correa et al.: Slick-haired Holstein physiology Figure1.A) Descriptive statistics for the thermal humidity index and B) environmental light intensity values averaged over 24 h for 2 consecutive days.GP represents when the cattle were located at the grazing paddock and Barn denotes when cattle were brought to a shaded facility for concentrate feeding, udder washing, and milking at 0300-0600 h and 1300-1600 h.The red horizontal line in Panel A denotes THI heat stress threshold for wild-type Holstein lactating cows.For Panel A and B the gray dotted line represents the environmental conditions at the grazing paddock versus the black continuous line represents the environmental conditions recorded at the shaded barn.Comparison of C) vaginal temperatures and D) light intensity as an indicator of voluntary solar radiation exposure for wild type-haired Holstein (black continuous line) versus SLICK (gray dash line) averaged over 24 h for 2 consecutive days.Vaginal temperatures and light intensities were recorded every 15 s for 48 h and averaged by hour for statistical analysis.Cows received concentrate feed and udders were washed at 0300-0400 h and 1300-1400 h at a shaded facility (Feeding, Washing; blue bars).Cows were milked twice a day at 0400-0600 h and 1400-1600 h (Milking; gray bars).Data for panels C and D are presented as LSMEANS ± SEM.Asterisks (*) represent significant differences at P < 0.05 and dagger ( †) represent tendency at 0.05 < P < 0.1.
Contreras-Correa et al.: Slick-haired Holstein physiology Contreras-Correa et al.: Slick-haired Holstein physiology 0.08°C.Conversely, an overall maximum daily vaginal temperature peak of 39.58 ± 0.09°C was observed between 1700 to 2200 h.Similar to the current study, both Sánchez-Rodríguez (2019b) and Dikmen et al. (2014) observed minimum and maximum daily peaks in vaginal temperature during the morning and afternoon, respectively.Figure 1D depicts the light intensity values recorded by the data loggers attached on top of cows' necks as indicators of voluntary solar radiation exposure.There was a time by hair coat type interaction (P = 0.03) affecting these light intensity values.The SLICK cows recorded greater light intensity values than their WT counterparts at 1000 h (15,467.00 ± 1,584.92 vs. 10,057.00± 787.35 lx; P < 0.0001) and1100 h (12,028 ± 1,660 vs. 8,571.11± 888.13 lx; P < 0.01).A tendency toward greater light intensities in SLICK Holsteins was observed at 1200 h (13,191.00 ± 2,188.86 vs. 10,923.00± 1,904.31lx; P = 0.08).No other differences in light intensity were observed.The data loggers used in the current study were previously validated by Sánchez and Domenech (2021) as a viable estimator of voluntary solar radiation exposure in dairy cattle.According to Brumby (1955), Contreras-Correa et al.: Slick-haired Holstein physiology

Figure 2 .
Figure 2. Mammary gland hemodynamics and infrared thermography of the rear udder from WT-haired (n = 7) and SLICK (n = 7) Holsteins at 160 ± 3 DIM.A) Average resistance index, B) average diameter, and C) total blood flow of the left and right external pudendoepigastric arteries.D) Representative images of the color Doppler ultrasonography performed to collect the mammary gland hemodynamics in wild-type (WT) and slick-haired (SLICK) cows.E) Visual image and corresponding F) infrared thermograph of the rear udder from one of the evaluated cows.G) Udder skin temperature data.Black bars represent WT cows and white bar with diagonal stripes represents SLICK cows.Data are presented as LSMEANS ± SEM.Asterisks (*) represent significant differences at P < 0.05.
Figure 3. Multidimensional scaling plot of the samples by hair coat type.

Figure 4 .
Figure 4. Heat map of genes of interest found to be upregulated or downregulated in the enriched GO analysis for slick-haired (SLICK) and wild-type (WT) Holstein.Each column represents an individual animal.The numbers inside each box are the counts per million reads (CPM)for each sample.The color intensity is tied to the CPM, the more intense a color, the higher the CPM for each GO gene set.The Up/Down next to each go term is the direction of change reported by the fry method.Up GO terms are enriched in SLICK, Down terms are depleted.

Figure 5
Figure5.Spearman correlations between metabolites in blood plasma and physiological variables measured in both hair coat types in the current study such as milk yield (MY), vaginal temperatures (VT), left and right pudendoepigastric arteries resistance index (RI), left and right pudendoepigastric arteries diameter (Diameter), and left and right pudendoepigastric arteries total blood flow (Total_MBF).Circle size and color intensity represent the spearman correlation between the metabolite and variable; blue is positive correlation, red is negative.A silver border highlights significant correlations (P < 0.05), gold highlights high significance (P < 0.01).A total of 14 samples were used for these correlations from WT (n = 7) and SLICK (n = 7) cattle.

Table 1 .
Contreras-Correa et al.: Slick-haired Holstein physiology List of differentially expressed genes in SLICK relative to WT Holsteins 1Log 2 FC, log 2 fold-change relative to wild-type haired.2 logCPM, log counts per million.3 LR, fold change ratio.4 FDR, false discovery rate.

Table 2 .
List of significantly enriched KEGG pathways in SLICK HolsteinsUp pathways are enriched in SLICK, down are depleted. 3

Table 3 .
List of differentially expressed metabolites in the SLICK Holsteins relative to WT