Effects of the bovine SLICK1 mutation in PRLR on sweat gland area, FOXA1 abundance, and global gene expression in skin

The SLICK1 mutation in the prolactin receptor ( PRLR ) results in a short-hair coat and increased ability to regulate body temperature during heat stress. It is unclear whether the mutation affects capacity for sweating. The objective of this observational study was to evaluate whether the SLICK1 mutation in PRLR alters characteristics of skin related to sweat gland abundance or function. Skin biopsies from 31 Holstein heifers, including 14 wild-type (SL −/− ) and 17 heterozygous slick (SL +/− ), were subjected to histological analysis to determine the percent of the surface area of skin sections that are occupied by sweat glands. We detected no effect of genotype on this variable. Immunohistochemical analysis of the forkhead transcription factor A1 (FOXA1), a protein essential for sweating in mice, from 6 SL −/− and 6 SL +/− heifers indicated twice as much FOXA1 in sweat glandular epithelia of SL +/− heifers as in SL −/− heifers. Results from RNA sequencing of skin biopsies from 5 SL −/− and 7 SL +/− heifers revealed few genes that were differentially expressed and none that have been associated with sweat gland development or function. In conclusion, results do not support the idea that the SLICK1 mutation changes the abundance of sweat glands in skin, but do show that functional properties of sweat glands, as indicated by increased abundance of immunoreactive FOXA1, are modified by inheritance of the mutation in PRLR .


INTRODUCTION
The prolactin receptor gene (PRLR) has undergone a variety of mutations in mammals (Bogorad et al., 2008;Iso-Touru et al., 2009). One type of mutation in cattle occurs in exon 11 of the gene and is characterized by a single base-pair deletion to cause a frameshift mutation, resulting in a truncation of the C-terminal region of the protein (Littlejohn et al., 2014;Porto-Neto et al., 2018;Flórez Murillo et al., 2021). A total of 6 such mutations have been described (Flórez Murillo et al., 2021). Collectively, these mutations are termed as slick mutations because affected animals display a "slick phenotype" characterized by a short-hair coat. Slick mutations are dominant; moreover, homozygotes and heterozygotes are similar phenotypically (Olson et al., 2003;. The mutations apparently arose either in criollo breeds of cattle in the Americas or in an ancestor of those cattle. The first slick mutation described, SLICK1 (c. 1382del;rs517047387;Littlejohn et al., 2014), has been identified in the criollo breeds Carora, Romosinuano, Caracu, Blanco Orejinegro, Costeño con Cuernos, and Hartón del Valle (Porto-Neto et al., 2018;Flórez Murillo et al., 2021). The gene is also frequent in the Senepol breed located on the Caribbean island of St. Croix; presumably, some founder animals possessed the gene. Crossbreeding has also been used to deliberately introduce the SLICK1 mutation into other breeds such as Holstein by crossbreeding (Olson et al., 2003;Dikmen et al., 2008Dikmen et al., , 2014Carmickle et al., 2022).
One role of prolactin is to inhibit hair growth (Foizik et al., 2009;Littlejohn et al., 2014). Therefore, slick mutations represent gain-in-function mutations, at Effects of the bovine SLICK1 mutation in PRLR on sweat gland area, FOXA1 abundance, and global gene expression in skin least with respect to regulation of hair growth. The mechanism involved is not known because slick mutations disrupt the region of the receptor involved in activation of JAK/STAT signaling (Porto-Neto et al., 2018).
Cattle inheriting slick mutations in PRLR exhibit increased thermotolerance with respect to regulation of body temperature during heat stress (Olson et al., 2003;Dikmen et al., 2008Dikmen et al., , 2014Eisemann et al., 2020;Landaeta-Hernández et al., 2021;Carmickle et al., 2022) and production of milk in the summer (Olson et al., 2003;Dikmen et al., 2014). Much of the increased thermotolerance of slick cattle can be ascribed to the short-hair coat. Such hair coats reduce the impedance of heat flow from the surface (Berman, 2004). Moreover, evaporation from the skin is enhanced when hair length is reduced by clipping (Dikmen et al., 2008(Dikmen et al., , 2014. Short hair length in cattle was associated with superior regulation of body temperature and high growth rates in hot climates (Hamblen et al., 2018;Evangelista Façanha et al., 2020). Regulation of hair length is likely to be a physiological mechanism for adaptation to hot climates. Cows exposed to chronic heat stress experience a decline in hair weight per unit of surface area (Kibler et al., 1965).
It is not clear whether slick cattle have greater capacity for sweating than animals without the slick mutation and, if so, whether greater capacity for sweating is a consequence of the shorter hair only or because of changes in the structure or function of the sweat gland.
In some experiments during heat stress, slick Holsteins had greater sweating rates than wild-type Holsteins (Dikmen et al., 2008(Dikmen et al., , 2014, whereas no difference was found between genotypes in another experiment (Carmickle et al., 2022). Angus × Senepol heifers inheriting the slick mutation experienced greater transepidermal water loss than wild-type animals (Eisemann et al., 2020). It is difficult to interpret differences in sweating rate between animals exposed to the same heat stress. Animals with a higher sweating rate could have a greater capacity for sweating and, therefore, experience lower body temperatures (Da Silva et al., 2015;Hooper et al., 2019). Alternatively, animals exhibiting greater sweating may do so because they are experiencing a greater degree of hyperthermia (Pereira et al., 2014).
An alternative method for studying capacity for sweating is to examine histological or functional characteristics of sweat glands. Qualitative evaluation of histology of skin from a limited number of animals heterozygous for the SLICK1 mutation did not reveal any differences associated with the mutation (Littlejohn et al., 2014). In contrast, it was reported that sweat gland size was greater for slick Criollo Limonero cattle (typically possessing the SLICK2 or SLICK3 mutation of PRLR) than animals with normal hair (Landaeta-Hernández et al., 2011). One potential marker of sweat gland function is the forkhead transcription factor A1 (FOXA1). Localized in secretory cells of eccrine sweat glands (Cui et al., 2012;Li et al., 2017Li et al., , 2018, this protein is required for secretory function of the sweat gland by regulating transcription of genes, encoding 2 ion transport proteins called bestrophin 2 (BEST2) and solute carrier family 12 member 1 (SLC12A1; Cui et al., 2012).
For the present study, 3 objectives related to effect of the SLICK1 mutation on histology and function of bovine skin were pursued. The first objective was to determine whether Holstein heifers heterozygous for the SLICK1 mutation (termed SL +/− ) have a higher proportion of the skin dermis occupied by sweat glands than wild-type heifers without the mutation (termed SL −/− ). A second objective was to evaluate whether amounts of immunoreactive FOXA1 in sweat glands differed between SL −/− and SL +/− heifers. Finally, RNA sequencing was used to identify differences in skin transcriptome associated with inheritance of the SLICK1 allele.

Ethics Statement
All animal procedures were approved by the University of Florida and University of California Davis Institutional Animal Care and Use Committees.

Genotyping
Heifers were genotyped for the SLICK1 mutation in PRLR using either the CLARIFIDE Plus platform (Zoetis) or the kompetitive allele specific-PCR procedure previously described (Sosa et al., 2021). For the purposes of this paper, animals were referred to as SL −/− if they did not inherit a SLICK1 allele (i.e., were wild-type) and SL +/− if they were heterozygous for the SLICK1 allele.

Collection of Skin Biopsies for Histology
Skin biopsies were taken from 31 Holstein heifers. A total of 16 heifers were from a private dairy located in Okeechobee, Florida (27°16′N 80°46′W, including 8 SL −/− and 8 SL +/− ). Another 15 heifers were from a private dairy located in Escalon, California (37° 47′ N/120° 59′ W), including 6 SL −/− and 9 SL +/−. All heifers were sired by one of 2 SL +/− Holstein bulls 9208 [Slick-Gator Blanco (551HO03574) or Slick-Gator Lone Ranger (047HO01029)] as described elsewhere . Skin biopsies for heifers in Florida were collected on July 31, 2020, between 0900 and 1000 h, when heifers were 34 to 37 wk of age. Skin biopsies for heifers in California were collected on October 9, 2020, between 1300 and 1600 h, when the heifers were 40 to 47 wk of age.
Skin biopsy from the shoulder area was performed while the animal was restrained in a head chute or head locks. The biopsy area was clipped, wiped down with 70% (vol/vol) ethanol-soaked gauzes, and then sprayed with 5% (wt/vol) dry-spray lidocaine (Aspercreme Lidocaine Dry Spray, Chattem, Inc.). The skin sample was collected by rotating a 5-mm skin punch (Integra Miltex Standard Biopsy Punch, Fisher Scientific) while pressing downwards with mild force, and then the skin sample was taken from the skin punch using a tweezer and placed into a cassette between sponges. A sterile gauze was placed on the biopsy site to remove excess blood or stop bleeding, and subsequently the wound site was closed with super glue (Gorilla Super Glue, Gorilla Glue Co.). A tick spray consisting of 0.5% (wt/ vol) permethrin (Neogen) was applied around the biopsy site.

Histological Analysis of Skin
Cassettes containing skin were labeled and placed into a container with a fixative [10% (vol/vol) neutral buffered formalin] at 4°C for 16 h. After fixation, skin samples were cut in half and rinsed in Dulbecco's phosphate buffered saline (DPBS) before embedding in paraffin at the Molecular Pathology Core of the University of Florida Department of Pathology, Immunology, and Laboratory Medicine (Gainesville). A microtome was used to prepare 6-μm sections, which were then placed on slides and either stained using hematoxylin and eosin or remained unstained for immunohistochemistry.

Quantification of Histological Sections
Slides stained with hematoxylin and eosin were digitized using the Leica Aperio VERSA brightfield scanner. The images were saved in SDF format for further analysis using QuPath-0.2.3 software (https: / / qupath .github .io/ ). Scanned images were opened in QuPath, and brightness and contrast were adjusted to clearly distinguish the anatomic structures of the skin. The wand tool was selected to calculate the area of the epidermis, dermis, and all observable sweat glands ( Figure  1). Area was expressed as percent of total area of the image.

Immunohistochemical Localization of FOXA1
A subset of histological preparation of skin from 12 heifers (2 SL −/− and 3 SL +/− from Florida and 4 SL −/− and 3 SL +/− from California) were subjected to immunohistochemical analysis of FOXA1 in sweat glands. The antibody used was a rabbit polyclonal antibody (from antibodies-online.com) raised against the C-terminus of FOXA1 (AA sequence = SFNHPFSINN LMSSSEQQHK LDFKAYEQAL QYSPYGSTLP ASL-PLGSASV) with 100% predicted reactivity with bovine FOXA1. The antibody was used at a dilution of 5 μg/ mL in antibody dilution buffer [DPBS containing 1% (wt/vol) BSA and 0.5% (vol/vol) Triton X-100]. The negative control was prepared by adding a blocking peptide (from antibodies-online.com) at a concentration of 100 μg/mL to a working solution of primary antibody.
All incubations were conducted in Coplin jars. Slides were subjected to deparaffinization and rehydration processes by serial immersion 2 times each in xylene, 100% (vol/vol) ethanol, and DPBS. After removing excess buffer, each slide was dried at room temperature. Antigen retrieval was performed by placing slides in prewarmed 10-mM citrate buffer, pH 6.0) for 30 min while in a steamer (IHC-Tek Epitope Retrieval Steamer, Fisher Scientific). Slides were then submerged in a dish with tap water for 10 min. The tissue was circled with a pap pen to confine immunohistochemistry buffers to that defined area, and each slide was dried at room temperature.
Immunohistochemistry was performed using the ImmPRESS HRP Horse Anti-Rabbit IgG PLUS Polymer Kit, Peroxidase (LSBio). Tissue was incubated with BLOXALL blocking solution provided by the kit for 1 h, and then each slide was submerged in wash buffer [DPBS + 0.1% (vol/vol) Tween 20] for 5 min. Subsequently, tissue was incubated with normal horse serum (provided by the kit) for 45 min and, after removal of serum, each slide was submerged in a wash buffer for 5 min. Tissue was then incubated with either anti-FOXA1 antibody or the control consisting of antibody and blocking peptide as described above. Incubation proceeded at 4°C overnight. Each slide was then incubated with wash buffer for 5 min, dried out at room temperature, and rehydrated with ImmPRESS Horse Reagent provided by the kit. Incubation was for 30 min at room temperature. Each slide was then incubated twice for 5 min each in wash buffer and once for 5 min in a working solution of horseradish peroxidase substrate (prepared by mixing equal volumes of ImmPACT DAB EqV Reagent 1 and Reagent 2 in the kit). After incubation in wash buffer, each section was counterstained with hematoxylin for 5 min. Slides were rinsed with tap water, dried at room temperature, and then coverslips mounted.

Analysis of Intensity of Immunoreactive FOXA1
Slides were digitalized using a Zeiss AxioScan.Z1 scanner equipped with a Hitachi HV-F203 camera and Zen software (version 3.1; Zeiss). The images were saved in CZI format for further analysis in QuPath-0.2.3 software. Digitized images were opened in QuPath, and the image was set by selecting the brightfield H-DAB option, which is exclusively for analysis of horseradish peroxidase-based immunohistochemistry. The DAB color option (brown) was selected for performance of measurements. The wand tool was used to select the region between the outer and inner edges of the epithelium of selected sweat glands. A range of 8 to 31 separate structures were analyzed for each section; specifically, one section was analyzed for each animal. Subsequently, selected regions were imported to ImageJ (version 1.53a; Wayne Rasband, NIH; https: / / imagej .nih .gov/ ij/ ) to calculate the intensity of horseradish peroxidase labeling. Before importing the regions to ImageJ, a negative image of the selected area was created in QuPath by selecting and applying color transformation and including overlay of the region, and this negative image was opened automatedly in ImageJ for quantification. In this way, increasing amounts of label was quantified as positive pixel intensity. A representative image of the procedure is shown in Figure 2. All analyses of labeling intensity were performed without knowledge of genotype by one person.

RNA Sequencing of Skin
Skin samples were collected as described above from 5 SL −/− (65 ± 13 wk of age) and 7 SL +/− (63 ± 12 wk of age) Holstein heifers from the University of Florida Dairy (29° 46′ 14″ N, 82° 25′ 16″ W) on November 7, 2021, between 0800 and 1000 h. All heifers were sired by Slick-Gator Lone Ranger. Each biopsy was placed in a cryovial and immediately immersed in liquid nitrogen. Skin samples were stored at −80°C until RNA-Seq analysis.
The RNA extraction was performed with the Qiagen RNeasy Mini kit according to the user's guide. Quantity of RNA was determined using the QUBIT fluorescent method (Invitrogen) and Agilent Bioanalyzer. An amount of 100 ng of high-quality total RNA with RNA integrity number of 7 or higher was used for library construction, using the reagents provided in the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs) and the NEBNext Ultra II Directional RNA Library Prep Kit (New England Biolabs) according to the manufacturer's instructions. Briefly, 100 ng total RNA was used for mRNA isolation using the NEBNext Poly(A) mRNA Magnetic Isolation Module. Then, the poly(A)-enriched RNA was fragmented in NEBNext First Strand Synthesis Buffer via incubation at 94°C for the desired time. This step was followed by first-strand cDNA synthesis using reverse transcriptase and random hexamer primer. Synthesis of double-stranded cDNA was done using the second strand master mix provided in the kit, followed by endrepair and dA-tailing. At this point, Illumina adaptors were ligated to the sample. Finally, the library was amplified and then purified with AMPure beads (Beckman Coulter). Library size and mass were assessed by analysis in the Agilent TapeStation using a DNA5000 Screen Tape. Quantitative PCR was used to validate the library's functionality using the KAPA library quantification kit (Kapa Biosystems) and monitoring on the BioRad CFX 96 real-time PCR system. Twelve barcoded libraries were pooled equimolarly for sequencing simultaneously for 0.23 lanes using the Illumina NovaSeq6000. Normalized libraries were submitted to the Illumina Free Adapter Blocking Reagent protocol to minimize presence of adaptor-dimers and index hopping rates. The library pool was diluted to 0.8 nM and sequenced on one S4 flow cell lane (2 × 150 cycles) of the Illumina NovaSeq6000. The instrument's computer used the NovaSeq Control Software (version 1.6). Cluster and SBS consumables were version 1.5. The final loading concentration of the library was 120 pM with 1% PhiX spike-in control. One lane generated 2.5 to 3 billion paired-end reads (~950 Gb) with an average Q30% ≥92.5% and Cluster passing filter = 85.4%. Fastq files were generated using the BCL2fastQ function in the Illumina BaseSpace portal. The Illumina NovaSeq 6000 was used to sequence the libraries for 2 × 150 cycles. Sequencing was performed at the ICBR NextGen Sequencing (https: / / biotech .ufl .edu/ next -gen -dna/ ; RRID: SCR _019152).

Bioinformatics
Reads acquired from the Illumina NovaSeq 6000 platform were cleaned up with the cutadapt program (Martin, 2011) to trim off sequencing adaptors and low-quality bases with a quality Phred-like score <20. Reads <40 bases were excluded from RNA-Seq analysis. The genome of Bos taurus (ARS-UCD1.2) was from Joint Genome Institute. The cleaned reads of each sample were mapped individually to the reference genome using the read mapper of the STAR package (Spliced Transcripts Alignment to a Reference, version 2.7.9a; Dobin et al., 2013). The mapping results were processed with HTSeq (High-Throughput Sequence Analysis in Python, version 0.11.2; Anders et al., 2015), samtools, and with scripts developed in-house at the University of Florida Interdisciplinary Center for Biotechnology Research (Gainesville) to remove potential PCR duplicates and choose and count uniquely mapped reads for gene expression analysis. Principal component analysis (for detecting outlier samples) and volcano plot analysis based on all identified genes in each analysis were performed with the R package (version 4.1.3). The counted reads of each gene were analyzed by the DEB application (Yao and Yu, 2011), which is a DESeq2-based R pipeline. Genotype effects independent of age were determined using a multifactor analysis with effects of treatment (e.g., SL −/− vs. SL +/− ) and age (<62 wk of age vs. >62 wk of age). Significant up-and downregulated genes were selected using the P-value, log 2 -fold change, or both for downstream analysis. Volcano plots (0.6 < log 2 -fold change < −0.6 and P-value ≤ 0.01) were constructed to visualize differential gene expression for each treatment with in-house R scripts.
The RNA sequencing results were deposited to the Gene Expression Omnibus of the National Center for Biotechnology Information (https: / / www .ncbi .nlm .nih .gov/ geo/ ), and the accession number is GSE200138.

Statistical Analysis
Treatment effects on histological characteristics of skin were analyzed by ANOVA using the GLM procedure of SAS 9.4 (SAS Institute Inc.). The model included effects of genotype, sire, location (Florida vs. California), and the genotype × location interaction. Treatment effects on FOXA1 labeling were analyzed by the GLIMMIX procedure of SAS. Multiple sweat glands (range = 8 to 31) were analyzed for each animal. The mathematical model included effects of genotype, sire, genotype × sire, and the labeling procedure replicate. Heifer nested within genotype × sire × replicate was considered random. Results are expressed as least squares means ± standard error of the mean.

Skin Histology
Representative examples of skin from SL −/− and SL +/− heifers are shown in Figure 1A and Figure 1B, respectively. Each section of skin was examined to quantify the proportion of the total area that was epidermis and dermis ( Figure 1C). We found no difference in proportions between SL −/− and SL +/− (P = 0.7065). The percent area that was dermis was 94.0 ± 0.7% for SL −/− and 93.7 ± 0.5% for SL +/− . The percent area that was epidermis was 6.0 ± 0.7% for SL −/− and 6.3 ± 0.5% for SL +/− . Additionally, the proportion of skin that was occupied by the lumens of sweat glands was determined. No differences between genotypes were detected (P = 0.4346; Figure 1D).

Immunoreactive FOXA1 in Sweat Glands
The FOXA1 was localized to epidermis, dermal stroma, hair shafts, and epithelia of sweat glands ( Figure  2A). Labeling was specific because it could be eliminated by absorbing antibody with a blocking peptide corresponding to the original antigen (compare Figure  2B and D with Figure 2C and E).
The amount of labeling for FOXA1 in sweat glands was quantified after creating a negative image of the glands using ImageJ ( Figure 2F, G). Intensity of labeling was affected by genotype (P = 0.0012) and was lower in SL −/− (0.39 ± 0.06 units) heifers compared with SL +/− (0.82 ± 0.06 units) heifers ( Figure 2H).

Global Gene Expression in the Skin
A volcano plot of the log 2 fold change (SL +/− /SL −/− ) and −log 10 P-value of comparisons of each of 22,373 genes analyzed is shown in Figure 3. We detected a total of 40 genes that were differentially expressed based on the criteria of P < 0.01 and a 2-fold difference in gene expression. However, no differentially expressed gene was identified when the criterion was a false discovery rate of <0.05. Among the genes that were differentially expressed using the less stringent criteria were tryptophan 2,3-dioxygenase (TDO2, 23.6fold increase, P < 0.0001), found to be upregulated in thermotolerant cattle (Morenikeji et al., 2020), interleukin 10 (IL10; 5.2-fold increase, P = 0.0006), known to be upregulated by prolactin (Matalka, 2003), and choline-O-acetyltransferase (CHAT; −27.9-fold; P = 0.0076), which encodes for the enzyme synthesizing the neurotransmitter acetylcholine.
The list of differentially expressed genes (P < 0.01, 2-fold) was queried to evaluate possible differential expression of genes identified by Ingenuity Pathway Analysis as being regulated by prolactin (579 genes), involved in growth of hair (62 genes), or involved in sweat gland development, number, and function (6 genes). None of these genes were in the group of differentially expressed genes. Expression of FOXA1 was similar between genotypes (fold change = 1.04). Bestro-phin 2 (regulated by FOXA1 and required for sweating in mice; Cui et al., 2012) was not in the data set, but related genes BEST1, BEST3, and BEST4 were expressed but not affected by genotype. The SLC12A1 (regulated by FOXA1 in the mouse) was 4.4-fold higher in SL +/− , but the difference was not significant (P = 0.2334) and the expression of the gene was low, with number of reads ranging from 0 to 16. A key gene involved in sweat production, aquaporin 5 (AQP5), was not affected by genotype (fold change = −1.34).

DISCUSSION
The overall objective of the present study was to evaluate whether the SLICK1 mutation in PRLR causes changes in sweat gland abundance or function that could improve regulation of body temperature during heat stress. Specific objectives were to test whether SL +/− heifers have a higher proportion of the dermis of the skin occupied by sweat glands compared with SL −/− heifers, whether amounts of immunoreactive FOXA1 in sweat glands differed between SL −/− and SL +/− heifers, and whether inheritance of the SLICK1 allele causes changes in the skin transcriptome. Results were that Sosa et al.: EFFECTS OF THE SLICK1 MUTATION ON SKIN Figure 3. Volcano plot showing differences in gene expression in skin tissue between heterozygote slick (SL +/− ) and wild-type (SL −/− ) Holstein heifers. Fold change is expressed as heterozygote/wild-type. Blue (downregulated) and red (upregulated) circles represent genes that were P < 0.01 and >2-fold, whereas black circles represent genes that did not meet these criteria. Letters next to each circle are gene symbols. no difference was found in the total area of sections of skin that were occupied by sweat glands, but that the SLICK1 mutation could affect function of the sweat gland because of increased amounts of immunoreactive FOXA1 localized in the sweat gland epithelium. In the mouse, Foxa1 is required for transcription of Best2 and Slc12a1 that are in turn necessary for secretory function of the sweat gland (Cui et al., 2012).
No effort was made to count the number of sweat glands or determine the maximum diameter of each gland. Doing so accurately is very difficult because of the tortuous shape of sweat glands and because the microtome sections each gland at a variety of different angles. The total area of sweat glands is proportional to the number and size of sweat glands. Failure to find a difference in the sweat gland area between SL −/− and SL +/− heifers is consistent with observations of Littlejohn et al. (2014), but in contrast to those of Landaeta-Hernández et al. (2011) who found that the size of sweat glands in Criollo Limonero adult females was greater for slick animals than non-slick animals. Results may differ between studies for many reasons, including anatomical site evaluated and environment. Littlejohn et al. (2014) evaluated sweat gland morphology in the ear, but only 3 SLICK1 Senepol × Holstein animals were examined. In the study of Landaeta-Hernández et al. (2011), 16 slick animals were compared with 14 animals with normal hair and skin biopsies were collected from the neck. Sweat glands are not distributed uniformly but are more abundant in some regions than others. For example, there are more sweat glands in the shoulder than rump (Findlay and Yang, 1950). Production of sweat in response to heat stress is also dependent on anatomical site (Scharf et al., 2008;Dikmen et al., 2014). Structure and function of the sweat gland can also be affected by environmental conditions as shown in the patas monkey (Sato et al., 1990), cattle (Nagib Nascimento et al., 2019), and buffalo (Debbarma et al., 2020). Thus, it is possible that differences between slick-and wild-type cattle in sweat gland structure and function may vary with anatomical site or environment. Nonetheless, the data here do not support the idea that development of sweat glands is affected by the SLICK1 mutation.
Studies have reported that animals inheriting the SLICK1 mutation produce more sweat than wild-type animals under conditions of thermal stress (Dikmen et al., 2008(Dikmen et al., , 2014Eisemann et al., 2020). One possibility is that slick animals have reduced insulative layer of hair with respect to evaporation. Clipping hair increases sweating rate (Dikmen et al., 2008(Dikmen et al., , 2014. Data presented here support the idea that altered prolactin signaling in SL +/− heifers could also result in changes in function of the sweat gland. The amount of immu-noreactive FOXA1 was twice as high in sweat gland epithelia of SL +/− heifers as in epithelia of sweat glands from SL −/− heifers. Knockout of Foxa1 in mice leads to anhidrosis (Cui et al., 2012). Prolactin receptors are present in the sweat gland, at least for the sheep (Choy et al., 1997), and evidence in humans suggests that prolactin can change the chloride ion concentration in sweat (Robertson et al., 1986). Additional experiments are warranted to determine whether the amount of FOXA1 in the sweat gland epithelium affects sweat gland function and whether, as speculated by Foitzik et al. (2009), the sweat gland is one of the physiological targets of prolactin.
Few changes in the skin transcriptome were associated with inheritance of the SLICK1 mutation. Important genes related to hair growth and sweat gland development and function were not affected by the SLICK1 mutation. We detected no effect of genotype on expression of FOXA1, for example, or of AQP5, which encodes for a water channel protein expressed in the eccrine sweat gland and required for sweating (Nejsum et al., 2002). We also found no effect of genotype on expression of other related aquaporin genes. The gene encoding for the enzyme that synthesizes the neurotransmitter acetylcholine (CHAT) was downregulated in SL +/− heifers. Acetylcholine can induce sweating in cattle (Joshi et al., 1968). Interestingly, TDO2 was upregulated in animals with the SLICK1 mutation. This gene was reported to be downregulated in expression in skin of Angus and Holstein cattle as compared with cattle of 2 thermotolerant breeds: the N'Dama and White Fulani (Morenikeji et al., 2020).
We also found a lack of effect of genotype on genes whose expression is known to be regulated by prolactin. The one exception (not in the list of prolactin-regulated genes compiled by Ingenuity Pathway Analysis) was for the cytokine gene IL10, which was more expressed in SL +/− than in SL −/− heifers. Prolactin stimulates IL10 secretion in immune and mammary epithelial cells (Kim et al., 2003;Matalka, 2003;Brand et al., 2004).
The lack of effect of genotype on gene expression of skin was unexpected given the effect of the SLICK1 mutation on hair length and FOXA1 abundance. It is likely that examination of the transcriptome of the entire skin, which includes many different cell types, obscured effects of the SLICK1 mutation on cells in the hair bulb or epithelium of the sweat gland. More cell-type specific methods for evaluating gene expression of specific subpopulations of skin cells should be considered in the future, including in situ hybridization and single cell RNA sequencing or isolation of specific cell types such as sweat glands (Hamzaoui et al., 2018). Failure to find more genes associated with sweating to be differentially expressed may be because expression of some genes involved in sweating are transcribed in response to sweating and not expressed in the absence of active sweating.
In conclusion, results do not support the idea that the SLICK1 mutation changes the abundance of sweat glands in skin but do show that functional properties of sweat glands, as indicated by increased abundance of immunoreactive FOXA1, are modified by inheritance of the mutation in PRLR.