pro-Genome-wide association study for milk production and conformation traits in Canadian Alpine and Saanen dairy goats

Increasing the productivity of Canadian dairy goats is critical to the competitiveness of the sector; however, little is known about the underlying genetic architecture of economically important traits in these populations. Consequently, the objectives of this study were as follows: (1) to perform a single-step GWAS for milk production traits (milk, protein, and fat yields, and protein and fat percentages in first and later lacta-tions) and conformation traits (body capacity, dairy character, feet and legs, fore udder, general appearance, rear udder, suspensory ligament, and teats) in the Canadian Alpine and Saanen breeds; and (2) to identify positional and functional candidate genes related to these traits. The data available for analysis included 305-d milk production records for 6,409 Alpine and 3,434 Saanen does in first lactation and 5,827 Alpine and 2,632 Saanen does in later lactations; as well as linear type conformation records for 5,158 Alpine and 2,342 Saanen does. Genotypes were available for 833 Alpine and 874 Saanen animals. Both single-breed and multiple-breed GWAS were performed using single-trait animal models. Positional and functional candidate genes were then identified in downstream analyses. The GWAS identified 189 unique SNP that were significant at the chromosomal level, corresponding to 271 unique positional candidate genes within 50 kb up-and downstream, across breeds and traits. This study provides evidence for the economic importance of several candidate genes (e.g., CSN1S1 , CSN2 , CSN1S2 , CSN3 , DGAT1 , and ZNF16 ) in the Canadian Alpine and Saa-nen populations that have been previously reported in other dairy goat populations. Moreover, several novel positional and functional candidate genes (e.g., RPL8 , DCK , and MOB1B ) were also identified. Overall, the results of this study have provided greater insight into the


INTRODUCTION
The Canadian dairy goat sector has rapidly expanded in the last decade due to growing demand for goat cheese, milk powder, and other products.Canada produced an estimated 65,302 tonnes of goat milk in 2019 on 329 registered operations (Canadian Dairy Information Centre, 2021).Although this is a small amount of milk when compared with the global production of about 20 million tonnes, national production in Canada increased by 52% between 2009 and 2019, whereas global goat milk production increased by 15.2% in the same period (Food and Agriculture Organization of the United Nations, 2019; Canadian Dairy Information Centre, 2021).Nevertheless, increasing the productivity of Canadian dairy goats is critical to the profitability of individual operations and the continued growth of the sector.
Dairy goat producers are paid based on milk components; thus, milk production and composition traits have a direct effect on producer revenue.Average yearly milk production has been estimated to be between 917 and 1,003 L per doe, or about 3.0 to 3.3 L per day, based on a standard 305-d lactation (Ontario Goat, 2017; M. Paibomesai, Ontario Ministry of Agriculture Food and Rural Affairs, Elora, ON, Canada, personal communication).However, the average yearly production per doe across herds varies greatly (e.g., 600 to 1,400 L across 20 herds; M. Paibomesai, personal communication), indicating substantial opportunity to improve dairy goat productivity.Although milk pro-duction traits have a direct effect on revenue, other functionally important traits such as conformation can also have a major influence on doe productive life and welfare, veterinary expenses, and replacement rearing costs, and thus contribute to the overall cost of production.
Compared with Canadian dairy cattle, the cost of production per liter of milk is considerably higher for goats ($1.27 vs. $0.77 per liter in 2016;Canadian Dairy Commission, 2017;Ontario Goat, 2017).The weighted average producer revenue per liter of goat milk was reported to be approximately $1.00 in 2020 (M.Paibomesai, personal communication); however, goat milk production is not supply managed, and prices continuously fluctuate.In comparison, average producer revenue from the sale of cow milk in Ontario is reported to be approximately $0.81, with approximately 9 to 10 times the quantity produced (Canadian Dairy Commission and Dairy Farmers of Ontario, 2021).Overall, these statistics highlight the need to increase production efficiency on Canadian dairy goat operations.
Genetic selection is one proven method of cumulatively and permanently increasing the efficiency of livestock production.However, genetic selection has been underutilized in the Canadian dairy goat sector because low participation in phenotype recording and registration programs has hindered traditional selection based on EBV.However, recently, Massender et al. (2022a,b) demonstrated that substantial gains in selection accuracy can be expected from the implementation of genomic selection for milk production and conformation traits of Canadian dairy goats.The implementation of genomic selection would also enable herds that have not traditionally participated in phenotype or pedigree recording to use genomic evaluations as a selection tool.In addition to genomic selection, there is interest from the sector in the development of genetic tests for important genes to use in marker-assisted selection schemes.
Although milk production and conformation traits are highly polygenic, several genes (e.g., CSN1S1, DGAT1) have been reported to have major influence on milk composition (Martin et al., 2017).Genetic testing for CSN1S1 has been used as a selection tool in dairy goat breeding programs in several countries, including Canada.Accounting for the effects of major genes has also been shown to increase the accuracy of genomic predictions (Carillier-Jacquin et al., 2016;Teissier et al., 2018Teissier et al., , 2019)).Thus, understanding the regions of the genome that influence milk production and conformation traits in the Canadian populations could improve genetic evaluations and may provide additional selection tools.
Genome-wide association studies are one method of exploring the genetic architecture of economically im-portant traits.To date, GWAS have been performed for milk production traits conducted in the French (Martin et al., 2017(Martin et al., , 2018;;Talouarn et al., 2020), UK (Mucha et al., 2018), and New Zealand (Scholtens et al., 2020) dairy goat populations.A GWAS was also performed for the Canadian Alpine and Saanen breeds (Vermette et al., 2013); however, fewer genotypes were available for the analysis (n = 720) and a multiple-step approach was used for the GWAS, with EBV of genotyped animals being used as pseudo-phenotypes.Additionally, Vermette et al. (2013) did not identify candidate genes in genomic regions with significant SNP nor discuss their potential biological relevance.Single-step GWAS can be beneficial when only a fraction of animals are genotyped, because the genomic EBV are estimated using all phenotypic and genotypic information available, which are then used to more accurately estimate the marker effects (Wang et al., 2012;Aguilar et al., 2019).
Consequently, the objectives of this study were (1) to perform a single-step GWAS to identify SNP significantly associated with milk production and conformation traits in the Canadian Alpine and Saanen breeds; and (2) to identify positional and functional candidate genes related to these traits to increase our understanding of the genetic architecture of economically important traits in Canadian dairy goat populations.

Phenotypes and Pedigree
The data used in this research were obtained from industry organizations or samples collected by commercial producers; thus, institutional animal care approval was not required.Phenotypic records used in this study were provided by the Canadian Centre for Swine Improvement (www .ccsi.ca)through the Canadian Dairy Goat Genetic Improvement Program (www .goatgenetics .ca).The data consisted of 305-d lactation milk production records for 5 traits and linear conformation scores for 8 traits.The 5 milk production traits were milk yield (MY, kg), protein yield (kg), fat yield (FY, kg), protein percentage (PP, %), and fat percentage (FP, %).The first and later lactation records were considered separate traits, as described in Massender et al. (2022a), and are noted throughout the text with the number 1 for first lactation records (i.e., MY1) and 2+ for later lactation records (i.e., MY2+).The 8 linear conformation traits were body capacity, dairy character, feet and legs, fore udder, general appearance, rear udder, suspensory ligament, and teats.
The data editing procedures, descriptive analyses, and genetic parameter estimation were presented in Massender et al. (2022a) and Massender et al. (2022b) for milk production and conformation traits, respectively.The final milk production data sets included records for 6,409 Alpine and 3,434 Saanen does in first lactation, and 12,236 Alpine records (5,827 does) and 5,008 Saanen records (2,632 does), recorded in lactations 2 to 11.The conformation trait data sets included linear type classification records for 5,158 Alpine and 2,342 Saanen does.The phenotypes of Alpine and Saanen were also combined for multiple-breed analyses.
Pedigree information for registered animals was obtained from the Canadian Livestock Records Corporation (www .clrc.ca)and trimmed to include only ancestors of animals with phenotypes or genotypes.The final pedigrees for the milk production analyses of Alpine and Saanen had 13,437 and 7,379 animals, respectively, and the multiple-breed pedigree had 20,239 animals.For the conformation traits, 11,486 Alpine and 6,270 Saanen animals were included in the singlebreed pedigrees, and 17,362 animals were included in the multiple-breed pedigree.Animals with phenotypic or genotypic records had an average pedigree depth of 15.8 to 20.3 generations, depending on breed and trait (Massender et al., 2022a,b).

Genotypes
A total of 1,707 genotypes were used in this study, of which 833 were from the Alpine breed (78 bucks, 755 does) and 874 from the Saanen breed (97 bucks, 777 does).All animals were genotyped with the first version of the GoatSNP50 Bead Chip (Tosser-Klopp et al., 2014) from Illumina Inc., which contains 53,347 SNP.Genotypic quality control was performed within breed, as described in Massender et al. (2022a), retaining 45,221 SNP (84.8%) for the multiple-breed analyses, and 44,598 (83.6%) and 43,598 SNP (81.7%) for single-breed analyses of the Alpine and Saanen breeds, respectively.

GWAS
The single-step GWAS analyses were performed using the blupf90 family of programs (Aguilar et al., 2014(Aguilar et al., , 2019;;Misztal et al., 2014;Lourenco et al., 2020) and the models described in Massender et al. (2022a,b), with the exception that default scaling parameters were used for the H matrix to reduce deflation of P-values.Single-trait animal models were used for all traits, with the effects in the models differing between trait groups and lactations.The models for milk production traits included fixed effects of doe age (in months, 8 to 108+), parity (2 to 7+) in the models for later lactation traits, and a linear covariate of DIM at the final milk test for the lactation.The random effects in-cluded contemporary group (herd and year of test) and animal additive genetic effects, as well as permanent environmental effects in the models for later lactation traits.For the conformation traits only the doe's first classification record was used, and the models included a fixed effect of parity (first or later), linear covariates of doe age (days) and DIM, and random contemporary group (herd-year-classifier) and animal additive genetic effects.For multiple-breed analyses, the phenotypes recorded on Alpine and Saanen were considered a single trait, and a fixed effect of breed (Alpine or Saanen) was included.More details on the models and genetic parameters used for each trait group are available in Massender et al. (2022a) for the milk production traits and Massender et al. (2022b) for the conformation traits.
Under the single-step GWAS framework, GEBV are back-solved to estimate marker effects (Wang et al., 2012).The P-values for the marker effects were obtained following the methods of Lu et al. (2018) and Aguilar et al. (2019), as implemented in the postGSf90 program (Aguilar et al., 2014).In the preliminary analyses, normal quantile-quantile plots and genomic inflation factor values were found to vary considerably between traits and breeds (Supplemental Figures S1 to S9; https: / / doi .org/ 10 .7910/DVN/ QJOS4D; Massender, 2022).Consequently, the P-values were adjusted using the method of genomic control to reduce the likelihood of false-positive or false-negative associations (Devlin and Roeder, 1999;Devlin et al., 2001).Specifically, genomic inflation factors were estimated as the observed median of chi-squared test statistics (i.e., square of the z-test statistic obtained from postGSf90 as in Aguilar et al., 2019) for all SNP tested in an analysis divided by the expected median of the chi-squared distribution with one degree of freedom (Devlin and Roeder, 1999;Devlin et al., 2001).The chi-squared test statistics were then divided by λ, and P-values were calculated from the adjusted test statistics (Devlin and Roeder, 1999;Devlin et al., 2001).Thereafter, SNP were considered significantly associated with the trait if the adjusted −log 10 P-value was greater than chromosomal-or genome-wise significance thresholds, after a Bonferroni correction for multiple comparisons at a 5% family-wise error rate [i.e., −log 10 P > −log 10 (0.05/n), where n is the number of SNP tested either for an individual chromosome or across the genome for a given breed].The genome-wise significance threshold was found to be 5.95, whereas chromosomal-wise significance thresholds ranged from 4.18 to 4.77 (mean ± SD: 4.46 ± 0.16).

Functional Analyses
The GALLO package (Fonseca et al., 2020), available in R software version 4.0.4(R Core Team, 2021), was used to retrieve positional candidate genes located within 50 kb up-and downstream of the significant SNP using the gene annotation information for ARS1 (Bickhart et al., 2017) available in NCBI data sets (Sayers et al., 2022) release 102 (www .ncbi.nlm.nih.gov/assembly/ GCF _001704415 .1).Gene ontology (GO) was performed with PANTHER version 16.0 (www .pantherdb.org/ ) to determine biological processes, cellular components, and molecular functions significantly (P-value <0.05) overrepresented among the positional candidate genes (Mi et al., 2019(Mi et al., , 2021)).Finally, the STRING database (www .string-db .org/ ) was used to construct gene networks based on predicted protein interactions between the annotated genes (Szklarczyk et al., 2021).

Positional Candidate Genes
Genes were identified as positional candidates if they were located within 50 kb upstream or downstream of significant SNP (−log 10 P > 4.18 to 4.77).The number of unique annotated genes by breed and trait (i.e., not double-counted if multiple SNP were annotated to the same gene) is summarized in Table 1.There were 271 unique annotated positional candidate genes across all analyses, of which 82.0% were protein coding genes and 18.0% were nonprotein coding genes (i.e., pseudogenes, long noncoding RNA, or tRNA).
Thirteen annotated genes were common across Alpine, Saanen, and the multiple-breed analyses, 105 were shared between either Alpine or Saanen and the multiple-breed analyses, and 52, 73, and 26 were unique to Alpine, Saanen, and the multiple-breed analyses, respectively.Positional candidate genes by breed and trait are described in Tables 2 to 4.

DISCUSSION
The GWAS described in this study was performed with a single-step approach, allowing phenotypic records measured on nongenotyped animals to contribute to the estimation of SNP effects (Wang et al., 2012).This approach was beneficial due to the limited number of genotypes available for analysis, which would have reduced the power to detect significant associations under a multiple-step approach.Although the Bonferroni correction used for multiple comparisons was stringent, and therefore only loci of major effect were expected to be identified, significant SNP were identified for all traits.The significant regions varied by trait and breed and were distributed across the genome.For some traits, the genomic inflation factor (Figure 4) varied considerably from one (Supplemental Figures S1 to S9; https: / / doi .org/ 10 .7910/DVN/ QJOS4D), which often indicates population stratification or "cryptic relatedness" among the genotypes (Devlin and Roeder, 1999;Devlin et al., 2001).This was unexpected, given that the models used accounted for pedigree structure, breed, and herd, similarly to the models used in routine genetic evaluations, and may warrant further investigation.The correction applied to the P-values aimed to reduce the likelihood of reporting false-positive or false-negative associations, but it should be acknowledged that this correction affects the significant regions and potential positional and functional candidate genes reported.It is important to note that genomic regions of importance to a trait are expected to result in clusters of significant SNP in linkage disequilibrium; thus, individual significant SNP may indicate spurious associations and should be validated in future research (Misztal et al., 2021).
To our best knowledge, this is the first study to identify potential positional candidate genes and perform functional analyses in Canadian dairy goat populations.Vermette et al. (2013) performed the first GWAS in this population using a subset of the genotyped animals in the present study (n = 720) and a multiple-step genomic BLUP approach.They identified significant SNP (false discovery rate <0.10) for MY (CHI17 and CHI24) and fore udder (CHI1 and CHI4) for Alpine, and significant SNP for fore udder (CHI1, CHI4, and CHI13) and teats (CHI4, CHI8, and CHI9) for Saanen.None of the significant regions identified by Vermette et al. (2013) were in common with those identified in the present study.
The most significant regions (−log 10 P > 8) were identified for PP and FP, which are highly heritable and known to be influenced by several major genes in dairy goats (Martin et al., 2017;Scholtens et al., 2020).In contrast, relatively few significant SNP were identified for the conformation traits.In general, the conformation traits were less heritable than the milk production traits and there were also fewer phenotypic records available for their analysis, both of which may have contributed to a lower power to detect significant associations.It is also worth noting that 63 unique significant SNP were not within 50 kb of any annotated genes (Supplemental Table S1; https: / / doi .org/ 10 .7910/DVN/ QIZ1F8; Massender, 2022).These regions should be studied further as the annotation of the goat genome improves.For the sake of brevity, this discussion will first highlight some key differences between breeds and lactations and then discuss some of the most functionally relevant positional candidate genes.

Comparison of Trait Groups
The physical structure of dairy animals can affect production.For example, McLaren et al. (2016) reported weak to moderate genetic correlations between several udder conformation traits and MY in UK dairy goats.Because milk production and conformation traits can be genetically correlated, there could be genomic regions and genes associated with both trait groups.However, no significant regions were found in common between the milk production and conformation traits in the present study.Similarly, Vermette et al. (2013) did not report any SNP associated with both trait groups in the population.In international research, a notable pleiotropic region on CHI19 has been found to be associated with milk production, udder conformation, and somatic cell score traits in French, UK, and New Zealand dairy goat populations (Martin et al., 2018;Mucha et al., 2018;Scholtens et al., 2020;Talouarn et al., 2020).The limited number of phenotypes available for the conformation traits is likely to have reduced the power to detect association, and these analyses should be repeated as more phenotypic and genotypic records become available.

Breed Differences
A total of 15 positional candidate genes were shared between the Alpine and Saanen breeds for milk production traits (ADAMTS3, CLOCK, CSN1S1, CSN1S2, CSN2, CSN3, DCK, LOC102178810, LOC102179276, LOC106502214, MOB1B, SHROOM3, SLC4A4, TRNAC-ACA, and TRNAC-GCA).Additionally, many of the positional candidate genes identified for milk production traits in the multiple-breed analyses were shared with at least one of the single-breed analyses.However, 26 genes were identified only in the multiple-breed analyses.Combining breeds can increase the power to detect significant associations in populations with a limited number of animals genotyped; however, this is only beneficial if the SNP are segregating and in the same gametic phase in both   breeds.Brito et al. (2015) reported that the consistency of gametic phase (i.e., Pearson correlation of signed r values between breeds) for the Canadian Alpine and Saanen breeds was 0.93 and 0.69 at distances of less than 0.02 Mb and between 0.02 and 0.03 Mb, respectively.Although both breeds have a common origin, they are substantially different at the genetic level (Brito et al., 2015(Brito et al., , 2017)).However, the use of multiple-breed GWAS could provide insight into genomic regions shared between the breeds.
Perhaps the most interesting difference between the breeds was the peak on CHI14 for PP1, FP1, FY2+, and FP2+ that was observed for Saanen and in the multiple-breed analyses, but not for Alpine.This chromosome contains DGAT1 and ZNF16, both of which have been previously identified as candidate genes in dairy goat populations (Martin et al., 2017;Scholtens et al., 2020).Martin et al. (2017) reported that the region containing DGAT1 was significantly associated with FP in both breeds.Additionally, the minor allele frequencies for the SNP in the region of DGAT1 were similar for both breeds (0.30 and 0.37 in Saanen vs. 0.43 and 0.49 in Alpine, as reported in Supplemental Table S2; https: / / doi .org/ 10 .7910/DVN/ MPHK5T; Massender, 2022) so the significant SNP in Saanen were also segregating in Alpine.These results could suggest differences in the genetic architecture for milk composition traits between breeds in the Canadian dairy goat populations.However, this finding should be validated as more animals are genotyped.

Comparison of Lactations
Overall, more significant SNP (−log 10 P > 4.18 to 4.77) were identified for later-lactation traits than for first-lactation traits (Tables 2 and 3).This might be because more phenotypic records existed for the laterlactation traits, which could lead to greater power to detect significant associations.Many of the positional candidate genes identified in first-lactation traits were also identified in later lactations.However, we found 30 and 75 unique positional candidate genes corresponding to SNP that were only significant in first or later lactations, respectively.These results highlight that many significant regions are shared between first and later lactations, but also suggest the possibility of some genetic differences between them that should be validated in the future.

Functional Candidate Genes
Genome-wide association studies identify regions of the genome associated with traits of interest in a specific population and can be a cost-effective method of exploring the genetic architecture underlying economically important traits.However, additional downstream analyses are needed to understand the functional relevance of the regions identified.Functional analyses have been limited in dairy goat populations because the goat reference genomes are not as well annotated as other major livestock species.The present study is one of the first to perform GO and network analyses in dairy goat populations; however, the functional relevance of many of the genes described have been inferred based on research in cattle and sheep.The GO analyses identified biological processes of broad importance, whereas the network analyses demonstrated the complex relationships and polygenic nature of the milk production traits.In the future, the use of other "-omics" technologies (e.g., transcriptomics, proteomics, metabolomics) could further our understanding of the biology of these traits.
Casein Genes.It is well known that the casein gene cluster (CSN1S1, CSN2, CSN1S2, CSN3) on CHI6 has a major influence on goat milk composition (Selvaggi et al., 2014).The genes in the casein cluster each encode for 1 of the 4 casein proteins, which together account for an estimated 64.5% of the total protein in goat milk (Selvaggi et al., 2014).This study identified SNP associated with CSN1S1 for PP1 in both the Alpine and multiple-breed analyses and PP2+ in all analyses; CSN2 for PP1 and PP2+ for all analyses, FP1 in the Saanen and multiple-breed analyses, and FP2+ in the multiple-breed analyses; CSN1S2 for PP1 and FP1 in Saanen and the multiple-breed analyses, PP2+ in all analyses, and FP2+ in the multiple-breed analyses; and CSN3 for PP1 and FP1 in Saanen and the multiplebreed analyses, and PP2+ for all analyses (Tables 2  and 3).This region has been previously reported in GWAS in French dairy goat populations (Martin et al., 2017).

Massender et al.: GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS
Given the large effect of the CSN1S1 genotypes on goat milk production identified in other countries (Martin et al., 2017), and the results of the present research identifying that this region is significantly associated with milk composition traits in Canadian dairy goats, the frequencies of the CSN1S1 genotypes should be investigated in the Canadian goat populations.
DGAT1.Significant SNP (−log 10 P > 6) for FP1 and FP2+ were identified in the region of the DGAT1 gene on CHI14 for Saanen and the multiple-breed analyses, explaining 0.07% to 0.21% of additive genetic variation.The DGAT1 gene is conserved across ruminant species and has a major influence on milk composition (Khan et al., 2021).It encodes for diacylglycerol O-acyltransferase 1 enzymes, which are rate-limiting and catalyze the synthesis of triglycerides (Khan et al., 2021;Mu et al., 2021).In goats, the DGAT1 gene has been found to be significantly associated with FP in the French  Alpine and Saanen populations (Martin et al., 2017).Mutations in the DGAT1 gene have been reported to explain between 6 and 46% of the phenotypic variation for FP (Martin et al., 2017).In dairy cattle, DGAT1 also has an influence on MY and PP (Nayeri and Stothard, 2016;Jiang et al., 2019) However, the results of  et al. (2017), that is, that the region containing DGAT1 only influenced FP.Several other genes in the same region as DGAT1 on CHI14 are also potential functional candidate genes for FP.The genes in this region include BOP1 (block of proliferation 1), CPSF1 (cleavage and polyadenylation specific factor 1), ADCK5 (aarF domain-containing protein kinase 5), FBXL6 (F-box and leucine rich repeat protein 6), HSF1 (heat shock transcription factor 1), SCRT1 (scratch family transcriptional repressor 1), SCX (scleraxis bHLH transcription factor), SLC52A2 (solute carrier family 52 member 2), and TMEM249 (transmembrane protein 249).Many of these genes have been significantly associated with milk production traits in dairy cattle (Nayeri and Stothard, 2016;Nayeri et al., 2017;Oliveira et al., 2019;Pedrosa et al., 2021) and displayed complex interactions in the gene network (Figure 10).However, given their proximity to DGAT1, further research is needed to determine their potential functional relevance in goats.
RPL8, ZNF16, ZNF34, and ZNF250.An additional SNP downstream of DGAT1 on CHI14 also displayed a significant association (−log 10 P > 4.5) with several milk composition traits (PP1, FP1, FY2+, FP2+).Candidate genes in this region included RPL8, which encodes for 60S ribosomal protein L8; sev-eral zinc finger protein coding genes (ZNF16, ZNF34, ZNF250); C14H8orf33, the chromosome 14 open reading frame 33 ortholog gene, and an uncharacterized locus (LOC106502817).It has been proposed that RPL8 may play an important role in bovine milk fat synthesis (Mu et al., 2021).Studies in Chinese Holstein cattle have identified a SNP in the promoter region of RPL8 that is associated with FP, independent of the effects of DGAT1 (Zheng et al., 2019;Mu et al., 2021).Furthermore, the functional importance of this gene to FP was confirmed, as silencing RPL8 resulted in significantly reduced expression of several enzymes known to be important to milk fat synthesis (Zheng et al., 2019).To the best of our knowledge, this gene has not been reported as a functional candidate gene in dairy goat populations.However, Scholtens et al. (2020) identified ZNF16 as the closest gene to a SNP significantly associated with FY in a multiple-breed New Zealand dairy goat population.The significant SNP identified by Scholtens et al. (2020) was 60 bp away from the identified SNP in the present study, supporting ZNF16 as a potential candidate gene.The ZNF16 gene is thought to have several biological functions, including cell proliferation (Li et al., 2011).
DCK and MOB1B.The DCK and MOB1B genes on CHI6 were identified as positional candidate genes for both Alpine (PP1 and PP2+) and Saanen (PP1,  (Slot Christiansen et al., 2015).In humans, DCK is also a medically important gene, as it is associated with drug resistance and sensitivity (Slot Christiansen et al., 2015).The MOB1B gene encodes MOB kinase activator 1B and, in human cells, has been reported to be involved in cell proliferation and apoptosis (Chow et al., 2010).The DCK gene has been previously identified as a positional candidate gene for MY (Pedrosa et al., 2021) and subclinical ketosis (Nayeri et al., 2019), and both genes have been associated with udder health traits (Abdel-Shafy et al., 2014;Wu et al., 2015;Moretti et al., 2021) in various dairy cattle populations, but are novel candidate genes for milk production traits in goats.Given the association of DCK and MOB1B with traits related to udder health in cattle, it would be interesting to assess whether these genes are also significantly associated with somatic cell count, as an indicator of mastitis.Unfortunately, this trait is not routinely recorded in the populations analyzed in this study.It has been previously suggested that mastitis may phenotypically influence milk composition traits (Paixão et al., 2017); thus, it is worth exploring potential pleiotropic effects for both milk composition and udder health traits of these genes in the population.

Future Research
This study confirmed that several regions associated with milk composition traits (e.g., CSN1S1, CSN2, CSN1S2, CSN3, DGAT1, and ZNF16) previously identified in other dairy goat populations (Martin et al., 2017;Scholtens et al., 2020) are also important in the Canadian Alpine and Saanen breeds.Additionally, we have proposed several novel candidate genes, based on their functional roles and importance in other dairy species (e.g., DCK, MOB1B, RPL8).
These analyses should be repeated as more animals are genotyped in the population.In particular, the breed-specific effect of regions on CHI14 for FP in Saanen and the relevance of DCK and MOB1B to udder health traits are 2 possible areas that would be interesting to further explore.It may be worthwhile to explore other approaches for GWAS in the future, such as the weighted single-step approach (Wang et al., 2012), which may have more power to detect significant associations, although it is still constrained because of its reliance on linkage disequilibrium.Performing GWAS based on whole-genome sequence data could also help identify potential causal mutations and additional genomic regions associated with the traits of interest.
Overall, the results of this study provide further evidence that the implementation of genomic selection will be beneficial to improving milk production traits in the Canadian dairy goat populations, as these traits are highly polygenic with only a few loci of major effect (Cole et al., 2009).The few major genes identified as positional and functional candidates could also be useful in marker-assisted selection programs or accounted for to improve the accuracy of genomic evaluations (Carillier-Jacquin et al., 2016;Teissier et al., 2018Teissier et al., , 2019)).Genetic testing for CSN1S1 is currently used as a selection tool by some Canadian dairy goat breeders, but the results of this research suggest that genetic testing for mutations in DGAT1 (Martin et al., 2017) could be another useful tool, at least for FP in Saanen goats.If available, genetic tests could be a useful tool, given the limited participation in industry genetic improvement programs by Canadian dairy goat breeders.

CONCLUSIONS
This study identified positional and functional candidate genes for milk production and conformation traits in the Canadian Alpine and Saanen dairy goat populations.The GWAS identified 189 unique SNP that were significant at the chromosomal level, corresponding to 271 unique positional candidate genes within 50 kb up-and downstream.We found some differences in the regions associated with the traits identified between breeds; for example, a peak on CHI14 for milk composition traits was identified only for the Saanen breed.Additionally, no overlapping regions were identified between the milk production and conformation traits.This study provides evidence that several candidate genes (e.g., CSN1S1, CSN2, CSN1S2, CSN3, DGAT1, ZNF16) are related to milk composition traits in these populations.In addition, several novel candidate genes (e.g., DCK, MOB1B, and RPL8) were proposed.Overall, this study provides insights into the genetic architecture underpinning economically important traits in Canadian dairy goat populations.
Massender et al.: GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS Figure 1.Manhattan plots for milk yield for Canadian Alpine and Saanen goats in first and later lactations, by breed.
Figure 2. Manhattan plots for protein yield for Canadian Alpine and Saanen goats in first and later lactations, by breed.

Figure 3 .
Figure 3. Manhattan plots for fat yield for Canadian Alpine and Saanen goats in first and later lactations, by breed.

Figure 4 .
Figure 4. Manhattan plots for protein percentage for Canadian Alpine and Saanen goats in first and later lactations, by breed.

Figure 5 .
Figure 5. Manhattan plots for fat percentage for Canadian Alpine and Saanen goats in first and later lactations, by breed.

Figure 6 .Figure 7 .
Figure 6.Manhattan plots for udder conformation traits for Canadian Alpine and Saanen goats, by breed.

Figure 8 .
Figure 8. Gene interaction network for the central genes associated with milk production traits for both Canadian Alpine and Saanen goats in first and later lactations.

Figure 9 .
Figure 9. Gene interaction network for the central genes associated with protein percentage for both Canadian Alpine and Saanen goats in first and later lactations.

Figure 10 .
Figure 10.Gene interaction network for the central genes associated with fat percentage for both Canadian Alpine and Saanen goats in first and later lactations.

Table 1 .
Number of chromosomally (genomically) significant SNP and unique annotated genes within 50 kb up-and downstream by breed and trait

Table 2 .
Chromosome-and genome-wise significant SNP and positional candidate genes within 50 kb up-and downstream by breed, for firstlactation milk production traits Massender et al.: GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS 1Traits: milk yield (MY), protein yield (PY), fat yield (FY), protein percentage (PP), and fat percentage (FP) in first lactation (1).2Position of significant SNP recorded as Capra hircus autosome: base pair.3Dashes indicate the SNP was not significant for a given analysis.

Table 2 (
Continued).Chromosome-and genome-wise significant SNP and positional candidate genes within 50 kb up-and downstream by breed, for first-lactation milk production traitsMassender et al.: GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS

Table 3 .
Chromosome-and genome-wise significant SNP and positional candidate genes within 50 kb up-and downstream by breed, for laterlactation milk production traits

Table 3 (
Massender et al.:GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS Continued).Chromosome-and genome-wise significant SNP and positional candidate genes within 50 kb up-and downstream by breed, for later-lactation milk production traits the present study are in line with the findings of Martin

Table 4 .
Massender et al.:GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS Chromosome-and genome-wise significant SNP and positional candidate genes within 50 kb up-and downstream by breed, for conformation traits

Table 5 .
Massender et al.: GWAS OF MILK PRODUCTION TRAITS IN DAIRY GOATS Significantly (P < 0.05) overrepresented gene ontology terms for positional candidate genes