High-throughput sequencing technologies have increased the ability to detect sequence variations for complex trait improvement. A high throughput genome wide genotyping-by-sequencing (GBS) method was used to generate 515,787 single nucleotide polymorphisms (SNPs), from which 76,355 SNPs with call rates >85% and minor allele frequency ≥1.5% were used in genome wide association study (GWAS) of 44 milk traits in 1,246 Canadian Holstein cows. GWAS was accomplished with a mixed linear model procedure implementing the additive and dominant models. A strong signal within the centromeric region of bovine chromosome 14 was associated with test day fat percentage. Several SNPs were associated with eicosapentaenoic acid, docosapentaenoic acid, arachidonic acid, CLA:9c11t and gamma linolenic acid. Most of the significant SNPs for 44 traits studied are novel and located in intergenic regions or introns of genes. Novel potential candidate genes for milk traits or mammary gland functions include ERCC6, TONSL, NPAS2, ACER3, ITGB4, GGT6, ACOX3, MECR, ADAM12, ACHE, LRRC14, FUK, NPRL3, EVL, SLCO3A1, PSMA4, FTO, ADCK5, PP1R16A and TEP1. Our study further demonstrates the utility of the GBS approach for identifying population-specific SNPs for use in improvement of complex dairy traits.