Genomic regions associated with muscularity in beef cattle differ in five contrasting cattle breeds

Background Linear type traits, which reflect the muscular characteristics of an animal, could provide insight into how, in some cases, morphologically very different animals can yield the same carcass weight. Such variability may contribute to differences in the overall value of the carcass since primal cuts vary greatly in price; such variability may also hinder successful genome-based association studies. Therefore, the objective of our study was to identify genomic regions that are associated with five muscularity linear type traits and to determine if these significant regions are common across five different breeds. Analyses were carried out using linear mixed models on imputed whole-genome sequence data in each of the five breeds, separately. Then, the results of the within-breed analyses were used to conduct an across-breed meta-analysis per trait. Results We identified many quantitative trait loci (QTL) that are located across the whole genome and associated with each trait in each breed. The only commonality among the breeds and traits was a large-effect pleiotropic QTL on BTA2 that contained the MSTN gene, which was associated with all traits in the Charolais and Limousin breeds. Other plausible candidate genes were identified for muscularity traits including PDE1A, PPP1R1C and multiple collagen and HOXD genes. In addition, associated (gene ontology) GO terms and KEGG pathways tended to differ between breeds and between traits especially in the numerically smaller populations of Angus, Hereford, and Simmental breeds. Most of the SNPs that were associated with any of the traits were intergenic or intronic SNPs located within regulatory regions of the genome. Conclusions The commonality between the Charolais and Limousin breeds indicates that the genetic architecture of the muscularity traits may be similar in these breeds due to their similar origins. Conversely, there were vast differences in the QTL associated with muscularity in Angus, Hereford, and Simmental. Knowledge of these differences in genetic architecture between breeds is useful to develop accurate genomic prediction equations that can operate effectively across breeds. Overall, the associated QTL differed according to trait, which suggests that breeding for a morphologically different (e.g. longer and wider versus shorter and smaller) more efficient animal may become possible in the future.


Background
Linear type traits have been used extensively to characterize conformation in both dairy [1][2][3] and beef cattle [4,5]. Muscularity linear type traits have previously been documented as moderate to highly heritable traits in beef cattle [5][6][7] and are known to be genetically associated with carcass merit [8,9] and with both animal live weight and price [4]. Therefore, the genetic merit of a young animal for these traits may be a good representation of its merit for carcass traits. While both carcass value and conformation have been reported to be correlated with linear type traits [9], the correlation with any one type trait is not equal to 1 which implies that the same carcass value can be achieved with morphologically different animals; by extension then, this implies that, for example, an animal with a better developed loin and a shallow chest may have the same yield as an animal with a lesser developed loin and a deep chest. Such morphological differences could contribute, in turn, to differences in individual carcass retail cut weights, and thus overall carcass value.
Many previous genomic studies in cattle have focused on live weight and carcass traits as the phenotypes of interest [10][11][12], but only a few have been published on the underlying features that contribute to differences in linear type traits in either beef cattle [13] or dairy cattle [14]. While previous studies have attempted to compare and contrast putative mutations, genes, and associated biological pathways across multiple breeds of beef cattle for carcass traits [15], no study has attempted to do this using linear type traits. Knowledge of any kind of similarities or differences between breeds could enable the introduction of more accurate multi-breed genomic evaluations for both pure and crossbred animals. Therefore, the objective of the present study was to identify genomic regions associated with five muscularity linear type traits and to determine if these associated regions are common across multiple beef cattle breeds.

Phenotypic data
As part of the Irish national beef breeding program, routine scoring of linear type traits is carried out on both registered and commercial beef herds by trained classifiers who are employed by the Irish Cattle Breeding Federation [4,16], with each classifier scoring animals from a range of different breeds. The muscularity type traits used in the present study describe the development of the hind quarter (DHQ), inner thigh (DIT), and loin (DL), and the width of the thigh (TW) and withers (WOW). Each trait was scored on a scale from 1 to 15 where 1 = low and 15 = high for DHQ, DIT and DL, and 1 = narrow and 15 = wide for TW and WOW (see Additional file 1: Table S1). Data on these five linear type traits were available for 147,704 purebred Angus (AA), Charolais (CH), Hereford (HE), Limousin (LM), or Simmental (SI) beef cattle scored between the age of 6 and 16 months from 2000 to 2016 [7]. Animals were discarded from the dataset if the sire, dam, herd, or classifier was unknown, or if the parity of the dam was not recorded. Parity of the dam was recoded as 1, 2, 3, 4, and ≥ 5. Contemporary group was defined as herd-by-scoring date generated separately per breed. Each contemporary group had to have at least five records. Following these edits, data were available on 81,200 animals: 3356 AA, 31,049 CH, 3004 HE, 35,159 LM and 8632 SI.

Generation of adjusted phenotypes
Prior to inclusion in the analysis, all phenotypes were first adjusted within-breed in ASREML [17] using the model: where y is the linear type trait, HSD is the fixed effect of herd by scoring date (11,130 levels), Sex is the fixed effect of the sex of the animal (male or female), AM is the fixed effect of the age in months of the animal (11 classes from 6 to 16 months), DP is the fixed effect of the parity of the dam (1, 2, 3, 4 and ≥ 5), Animal is the random additive effect of the animal, and e is the random residual effect. The adjusted phenotype was the raw phenotype minus the fixed effect solutions of HSD, Sex, AM and DP.

Genotype data
Of the 81,200 animals with linear type trait information, 19,449 animals from five beef breeds (1444 AA, 6433 CH, 1129 HE, 8745 LM, and 1698 SI) were imputed to whole-genome sequence as part of a larger dataset of 638,662 multi-breed genotyped animals. All 638,662 animals were genotyped using the Bovine Illumina SNP50 panel [n = 5808; 54,001single nucleotide polymorphisms (SNPs)], the Illumina High Density (HD) panel (HD; n = 5504; 777,972 SNPs), the Illumina 3k panel (n = 2256; 2900 SNPs), the Illumina low-density (LD) genotyping panel (n = 15,107; 6909 SNPs) or a bespoke genotype panel (IDB) developed in Ireland [18] with three versions, i.e. version 1 (n = 28,288; 17,137 SNPs), version 2 (n = 147,235; 18,004 SNPs) and version 3 (n = 434,464; 53,450 SNPs). Each animal had a call rate higher than 90% and only autosomal SNPs, SNPs with a known chromosome and position on UMD 3.1, and SNPs with a call rate higher than 90% within a panel were retained for imputation.
All genotyped animals were imputed to HD using a two-step approach in FImpute2 with pedigree information [19]; this involved imputing the 3 k, LD and IDB genotyped animals to the Bovine SNP50 density, and consequently imputing all resulting genotypes (including the Bovine SNP50 genotypes) to HD using a multi-breed reference population of 5504 influential sires genotyped on the HD panel. Imputation to whole-genome sequence y = HSD + Sex + AM + DP + Animal + e, (WGS) was then undertaken using a reference population of 2333 Bos taurus animals from multiple breeds from Run6.0 of the 1000 Bull Genomes Project [20]. All variants in the sequence reference population were called using SAMtools and genotype calls were improved using the Beagle software to provide a consensus SNP density across all animals. Details of the alignment to UMD 3.1 bovine reference genome, variant calling and quality controls completed within the multi-breed reference population are described in Daetwyler et al. [20]. In total, 41.39 million SNPs were identified across the genome and the average coverage was 12.85X. Imputation of the HD genotypes to WGS was completed by first phasing all 638,662 imputed HD genotypes using Eagle (version 2.3.2) [21], and subsequently imputing to WGS using minimac3 [22]. The average genotype concordance of imputation to WGS, defined as the proportion of correctly called SNPs versus all SNPs using a validation set of 175 Irish animals, was estimated to be 0.98 [23].
Quality control edits were imposed on the imputed sequence genotypes within each breed, separately. Regions of poor WGS imputation accuracy, which could be due to local mis-assemblies or mis-orientated contigs, were removed. These regions were identified using an additional dataset of 147,309 verified parent progeny relations as described by [23], which removed 687,352 SNPs from each breed. Then, all SNPs with a minor allele frequency (MAF) lower than 0.002 were removed. Following all SNP edits, 16,342,970, 17,733,147, 16,638,022, 17,803,135 and 17,762,681 autosomal SNPs remained for the analysis of the AA, CH, HE, LM, and SI populations, respectively.

Association analyses
The association analyses were performed within each breed separately using a linear mixed model in the GCTA software [24]. Autosomal SNPs from the original HD panel (i.e., 734,159 SNPs) were used to construct the genomic relationship matrix (GRM). The model used for the within-breed analysis was the following: where y is a vector of preadjusted phenotypes, μ is the overall mean, x is the vector of imputed genotypes, b is the vector of additive fixed effects of the candidate SNP to be tested for association, u ∼ N 0, Gσ 2 u is the vector of additive genetic effects, where G is the genomic relationship matrix calculated from the HD SNP genotypes and σ 2 u is the additive genetic variance, and e ∼ N 0, Iσ 2 e is the vector of random residual effects and σ 2 e is the residual variance. Manhattan plots were created for each trait within each breed separately by using the QQman package [25] in R.

QTL detection, gene annotation and variance explained
A genome-wide SNP significance threshold of p ≤ 1 × 10 −8 and a suggestive threshold of p ≤ 1 × 10 −5 were applied to each trait. SNPs in close proximity to each other (< 500 kb) were classified as being located within the same QTL. Genes within 500 kb of the most significant SNP in a peak above the genome-wide threshold were identified using Ensembl 94 [26] on the UMD 3.1 bovine genome assembly. Moreover, the functional consequence of all significantly associated SNPs was predicted using the Variant Effect Predictor tool [27] from Ensembl. The Cattle QTLdb (https ://www.anima lgeno me.org/cgi-bin/QTLdb /BT/index ) was used to identify QTL that were known to be associated with other traits in cattle. To identify QTL regions that were suggestive in more than one breed, each chromosome was split into 1-kb genomic windows, and windows containing suggestive SNPs (p ≤ 1 × 10 −5 ) were compared across the breeds.
The proportion of genetic variance of a trait explained by a SNP was calculated as: where p is the frequency of the minor allele, a is the allele substitution effect and σ 2 g is the genetic variance of the trait in question.

Meta-analysis
Following the within-breed association analyses, metaanalyses were conducted for all traits across all five beef breeds using the weighted Z-score method in METAL [28]; only SNPs that were included in the analyses of all of the individual breeds were considered here. METAL combines the p-values and the direction of SNP effects from individual analyses, and weights the individual studies based on the sample size to compute an overall Z-score: where w i is the square root of the sample size of breed i, and z i is the Z-score for breed i calculated as where φ is the cumulative distribution function, and P i and Δ i are the p-value and direction of effect for breed i, respectively.

Conditional analyses
The summary statistics from the individual analyses for the CH population were further used to conduct conditional analyses on BTA2 based on the Q204X mutation, which was previously reported to be associated with muscularity traits in cattle [29]. These analyses were undertaken for each trait in the CH population using the conditional and joint association analysis (COJO) method in GCTA [30]. The Q204X mutation was included as a fixed effect in the association analysis model and the allele substitution effect of all remaining SNPs were re-estimated.

Pathway and enrichment analyses
Pathway analysis was conducted on all plausible candidate genes within a 500-kb region up-and downstream of SNPs that were discovered to be suggestively or significantly associated with each trait in each breed. For each gene list, DAVID 6.8 [31] was used to identify gene ontology (GO) terms and KEGG pathways which were significantly overrepresented (p < 0.05) by the set of genes. Enrichment analyses among the suggestive and significant SNPs were performed to estimate if the number of SNPs in each annotation class was greater than that expected by chance for each trait per breed [32]; this was done separately per trait and per breed and was calculated as: where a is the number of suggestive and/or significant SNPs in the annotation class of interest, b is the total number of suggestive and/or significant SNPs that were associated with the trait of interest, c is the total number of SNPs in the annotation class in the association analysis, and d is the overall number of SNPs included in the association analysis.

Results
Summary statistics of the five linear type traits for each breed are in Additional file 1: Table S1. Significant (p ≤ 1 × 10 −8 ) and/or suggestive (p ≤ 1 × 10 −5 ) SNPs were detected in all traits for the five breeds but the exact locations of these SNPs and the direction of the effects of these SNPs differed by breed. Manhattan plots for all the analyses are available in Additional file 2: Figures S1-S5.

Within-breed analyses Angus
Whereas no significant SNPs were detected for any of the muscularity linear type traits in the AA population, suggestive SNPs (p ≤ 1 × 10 −5 ) were identified for all five traits. No genomic region was common to all five type traits (see Additional file 3: Figure S6). However, there was some overlap in suggestive 1-kb windows between the traits DIT and TW; 11 windows contained SNPs of suggestive significance and the gene EMILIN22 on BTA24 was identified within those windows for both traits. Nine genomic windows were associated with both the DL and WOW traits, i.e. on BTA6 (n = 2), BTA15 (n = 6), and BTA22 (n = 1). The windows on BTA15 contained suggestive SNPs that were located within the UCP3 and CHRDL2 genes. Eighty-four SNPs within nine QTL were suggestively associated with the DHQ trait. Among these, the most strongly associated (p = 3.34 × 10 −7 ) SNP was rs433492843 on BTA23 located in an intron of the PTCHD4 gene (Table 1); it accounted for 0.002% of the genetic variance in this trait. A QTL on BTA1 was also strongly associated with DL with the most strongly associated SNP being rs465472414 (p = 1.06 × 10 −6 ), which accounted for 0.08% of the genetic variance in this trait ( Table 2). Other SNPs suggestively associated with DL were also identified within the TMEM178A gene on BTA11 and within the UCP3 and CHRDL2 genes on BTA15.
An intergenic SNP located on BTA29, rs109229230, was the most strongly associated (p = 1.82 × 10 −7 ) with DIT (Table 3). Ninety-eight SNPs were suggestively associated with TW. The strongest QTL association with TW was on BTA13, on which 10 SNPs of suggestive significance were identified in a 1-Mb region (Table 4); rs137458299 displayed the strongest association (p = 2.99 × 10 −7 ) and explained 0.9% of the genetic variation in TW. One hundred and seventy-three SNPs were associated with WOW in the AA population; among these 29.4% were located on BTA14 (Table 5) and the most strongly associated SNP, rs468048676, (p = 2.34 × 0 −9 ), was an intergenic variant on BTA6.

Hereford
No significant SNPs were detected for any of the muscularity linear type traits in the HE population, although suggestive SNPs were identified for all five traits. However, no genomic window was common to all five type traits (see Additional file 3: Figure S6); six 1-kb windows i.e. on BTA5 (n = 1), BTA7 (n = 4), and BTA25 (n = 1) were shared between DHQ and DIT with three 1-kb regions on BTA20 shared between DIT and TW.
Three hundred and eleven SNPs were suggestively associated with DHQ. The strongest association with DHQ was located within a 1-Mb QTL on BTA7 where 26 SNPs of suggestive significance were identified ( Table 1). The intergenic SNP, rs446625612 (p = 1.16 × 10 −7 ) was the most strongly associated with DL and located within a Table 1 Location of the most significant QTL, limited to the top five per breed, which were associated with development of hind quarter and the genes located within these QTL within each breed     QTL on BTA4 encompassing the ENSBTAG00000044810 gene. Most interestingly, the strongest association within the QTL on BTA2 with DL was an intronic variant, which explained 0.7% of the genetic variance and was located within the muscle related gene MYO1B.
In total, 155 SNPs were suggestively or significantly associated with DIT, and 43% of these were located within a 1-Mb QTL on BTA7 (Table 3) where a number of significant SNPs were located within the EBF1 gene. For TW, four putative candidate genes were identified (Table 4): GABRA6 on BTA7, TTLL5 on BTA10, and both ADAMTS12 and GDNF on BTA20. The SNP, rs380761563, which displayed the strongest association with WOW, explained 1% of the genetic variance and was located in an intron of the gene TNIP1 on BTA7 (Table 5).

Charolais
There were 483 1-kb suggestive genomic windows common to all five type traits in the CH population (see Additional file 3: Figure S6), among which the vast majority (n = 482) were located on BTA2 in a region encompassing the MSTN gene. The final region that was shared between all five traits was on BTA11. More overlaps were found for DHQ and DIT with 904 windows being common to just these two traits, 146 windows common to DHQ, DIT, and DL, 304 windows common to DHQ, DIT, DL, and TW, and 178 windows common to DHQ, DIT, and TW. The majority of all these windows were also located on BTA2.
For each of the muscularity linear traits, we identified a QTL on BTA2 in the CH population. DHQ had the largest number of associated SNPs, i.e. 3707 suggestive and 1851 significant SNPs (Table 1), all of which were located on BTA2 within a single QTL between positions 0.35 and 9.79 Mb. In total, 41 genes including MFSD6, MSTN, and MYO7B were located in this QTL. For DIT, a 10-Mb QTL on BTA2 was identified that contained 5075 SNPs, of which 1796 had a p-value that met the significance threshold (Table 3), whereas 178 SNPs on BTA2 in the region between 54.1 and 86.1 Mb were significantly associated with TW ( Table 4). The same SNP, an intergenic variant rs799943285, showed the strongest association with all traits. The well-known Q204X mutation within the MSTN gene was significantly associated with DHQ, DIT and TW, and this SNP explained 4.9, 0.05, and 0.01% of the genetic variation of each trait, respectively.
In the conditional analyses within the CH population, where the Q204X mutation was included as a fixed effect in the model, the most significant SNPs from the original analyses of each trait generally decreased in significance. The most significant SNP for all traits in the original analyses was rs799943285 (p-value ranging from 9.07 × 10 −49 for DIT and DHQ to 2.02 × 10 −21 for WOW). In the conditional analyses, this SNP was non-significant for DL, TW, and WOW but remained suggestive for both DIT (p = 4.02 × 10 −6 ) and DHQ (p = 4.62 × 10 −6 ). The most significant SNP in the conditional analyses of DHQ, DL, DIT, and TW was rs41638272, which is an intergenic SNP located 10 kb from the SLC40A1 gene; this SNP was significant in the original analyses but its significance actually increased when the Q204X mutation was included as a fixed effect. The most significant SNP in the conditional analysis of WOW was an intergenic variant, rs457456302 (p = 4.78 × 10 −10 ) that was located 0.1 Mb from the MSTN gene.

Limousin
There were 164 1-kb suggestive genomic regions that were common across all muscularity traits in the LM population (see Additional file 3: Figure S6); another 232 regions were common to the three traits DHQ, DIT, and TW, while 326 were common to just DHQ and DIT. All five traits had significant QTL located on BTA2, with four genes common to all traits located within these QTL, namely ASNSD1, GULP1, SLC40A1, and ANKAR.
For DHQ, there were 2983 SNPs above the suggestive threshold and most of these (n = 2610) were located in a single QTL on BTA2. The most significant SNP, rs211140207 (p = 3.22 × 10 −30 ), was located within an 8-Mb QTL on BTA2 that contains 20 genes ( Table 1). The Q204X stop-gain mutation (rs110344317) located within this QTL was significantly associated with DHQ and accounted for 2.4% of the genetic variation in this trait, although the allele frequency of the favourable mutation was only 0.02% in the LM population. The wellknown MSTN mutation in the Limousin breed, F94L (MAF = 0.3798), did not meet the suggestive threshold for association with any of the traits. Similar to DHQ, a QTL located between 4.9 and 11 Mb on BTA2 was associated with both DIT (Table 3) and TW (Table 4). In total, 2441 and 1526 SNPs were above the suggestive threshold within this QTL on BTA2, and the variant rs110344317, which was significantly associated with DHQ, was also significantly associated with both DIT and TW. For the DL trait, 748 SNPs were suggestively associated and located between 55.4 and 82.8 Mb on BTA2. The most significant SNP associated with DL (rs379791493; p = 6.69 × 10 −10 ) was also the most significantly associated SNP with DIT (p = 2.20 × 10 −28 ). The most significant SNP associated with WOW, rs211140207, (p = 8.77 × 10 −12 ), was an intergenic SNP that accounted for 0.4% of the genetic variance in this trait and was located in a QTL (between 5.9 and 8.4 Mb) that included 724 other significantly-associated SNPs (Table 5).
Suggestive QTL were also detected on autosomes other than BTA2 for all traits in the LM population except for DIT. A small QTL on BTA11 containing seven suggestive SNPs was associated with DHQ. The SNP with the strongest association, rs43666945 (p = 1.56 × 10 −6 ), was an intergenic SNP located 2.2 Mb from the DYSF gene. Both DHQ and DL had suggestively associated QTL on BTA5. The most strongly associated SNP for DHQ (p = 1.58 × 10 −7 ) was an intergenic SNP, rs718375830, located within a QTL between positions 59.6 and 60.6 Mb, whereas the most strongly associated SNP with DL (p = 2.70 × 10 −6 ) was also an intergenic SNP, rs109909829, but was located within a QTL between 71.7 and 72.8 Mb.

Simmental
For the SI breed, only a few suggestive 1-kb genomic regions overlapped for more than two traits. Sixteen 1-kb windows were suggestively associated with both DHQ and DL, eight of which were located on BTA6, seven on BTA22, and one on BTA18 (see Additional file 3: Figure  S6). Five 1-kb windows on BTA23 and one on BTA4 were common to both DHQ and DIT, while another 15 suggestive windows were associated with DHQ and WOW, 12 of which were located on BTA22.
The intergenic SNP, rs437686690 on BTA25, was the most strongly associated (p = 1.00 × 10 −7 ) with DHQ in the SI population and accounted for 0.6% of the genetic variance in DHQ (Table 1). In total, 199 SNPs were associated with DL in the SI population, among which four met the significance threshold. The most significant SNP, rs482545354 (p = 9.77 × 10 −9 ), was located in an intronic region of the SUCGL2 gene (Table 2) on BTA22. Although 194 SNPs were suggestively associated with DIT, only one, i.e., rs798946118 (p = 5.30 × 10 −8 ), achieved the significance threshold which was located on BTA21 within a 1-Mb block containing 17 other suggestive SNPs (Table 3) and accounted for 0.6% of the genetic variance of DIT. The largest 1-Mb QTL associated with TW was located on BTA29 and contained 30 suggestive SNPs (Table 4). QTL putatively associated with WOW were located on BTA1, 4, 9, 12, and 20 (Table 5) where the most significant SNP, rs801295753 (p = 5.67 × 10 −8 ), was an intronic SNP on BTA9 located within both the ROS1 and ENSBTAG000000039574 genes.

Meta-analyses
Within each of the five meta-analyses (see Additional file 4), a strong association peak on BTA2 around the MSTN gene was detected, which is consistent with the individual association results identified in the CH and LM populations. For DIT, TW, and WOW, the most significantly associated SNP was the intergenic SNP, rs799943285 (p = 5.51 × 10 −24 ), which was previously identified as the most strongly associated SNP in the CH population for each of these traits. This variant, rs799943285, was also the most significantly associated with DL in the meta-analysis, whereas the most significantly associated SNP with DHQ, rs482419628 (p = 2.06 × 10 −47 ), was located further downstream on BTA2 within 5 kb of the ASNSD1 gene.
Although the QTL on BTA2 was the most strongly associated with each of the traits analysed, we also identified several other QTL associated with muscularity. In the meta-analysis of DHQ, the most strongly associated SNP on BTA11, rs43666945 (p = 1.93 × 10 −7 ), was previously identified as being associated with DHQ in the LM population, but the level of significance increased in the meta-analysis and the QTL contained three times the number of suggestive SNPs compared to that found for the LM breed only. A 1-Mb QTL on BTA7 containing the SPRY4 and FGF1 genes was associated with both DL and WOW in the meta-analysis; the most significant SNPs in this QTL, however, differed according to trait (see Additional file 4).

Enrichment of SNPs
With the exception of WOW in the AA population, intergenic SNPs were the most common annotation class of SNPs that were significantly associated with all traits in all breeds. The 3′ UTR class was enriched for all traits in the CH and LM populations, whereas there were more downstream gene variants significantly associated with DHQ and DL in the AA, CH and HE populations, and with TW in the CH, HE, and SI populations than expected by chance ( Table 6). The intronic class of SNPs was enriched for all five traits in HE, for four traits (DHQ, DL, TW, and DIT) in SI, three traits in both AA (DHQ, DL, and WOW) and CH (DL, TW, and WOW) and two traits in LM (DHQ and DIT).

Gene ontology and KEGG pathways
Several GO terms and KEGG pathways were over-represented by the genes identified in each analysis, although this tended to differ per breed and per trait especially in the smaller AA, HE, and SI populations. In CH and LM, five GO terms were associated with each trait: skin development (GO:0043588), collagen fibril organisation (GO:0030199), extracellular matrix structural constituent (GO:0005201), cellular response to amino acid stimulus (GO:0071230), transforming growth factor beta receptor signalling pathway (GO:0007179). One KEGG pathway, i.e. protein digestion and absorption (KEGG:map04974), was also significantly associated with all traits in CH and LM. Apart from this overlap, only a limited number of terms and pathways were over-represented across  breeds. The GO term mitochondrial inner membrane (GO:0005743) was significantly over-represented for the DL trait in AA and the WOW trait in HE, although none of the same genes were significantly associated with both traits. Another GO term collagen trimer (GO:0005581) was over-represented for DIT in AA and DL in LM.

Discussion
Whereas a number of across-breed and breed-specific pleiotropic QTL have been documented for carcass traits, birth weight, weaning weight, and mature weight in beef cattle [15], as well as for dry matter intake and growth and feed efficiency [33], no study has attempted to detect across-breed or breed-specific pleiotropic QTL for muscularity linear type traits. Previous studies have been conducted on the genetic correlations between the linear type traits themselves [7] and between both meat yield and carcass cuts with the muscularity linear type traits [34]. While these genetic correlations are moderate to strong, none is equal to 1, which implies that two animals that yield a carcass of similar merit could be morphologically different. In fact, a shorter and more muscular animal or a taller and less muscular animal could have the same total carcass weight. In turn, these animals could yield very different carcass values owing to their distribution of primal cuts. For example, the loin of an animal harbours generally the most valuable cuts [35,36]. Therefore, selection for a better-developed loin could lead to a more valuable carcass in comparison to a carcass with a lesser-developed loin if that carcass was still within the factory specification for weight and conformation. Here, we have detected several genomic regions that are strongly associated with each of the muscularity traits analysed. However, most of these regions were unique to each trait or each breed, which indicates the existence of trait-specific and breed-specific QTL for muscularity traits. Thus, it is plausible to hypothesise that through more precise (i.e., targeting individual QTL) genome-based evaluations and selection, the morphology of an animal could be targeted to increase the output of high-quality carcass cuts and consequently improve the profitability of the farm system and the value to the meat processor [36]. While a similar conclusion could be achieved through traditional breeding means, exploiting the breed-and trait-specific QTL could be more efficient. This is the first published genome study on muscularity linear type traits in beef cattle using sequence data and is one of the few genome-based studies that compare multiple breeds of beef cattle. The number of animals used in our study is comparable to the number of animals used in a previous across-breed comparison that focused on carcass and birth traits in 10 cattle breeds [15] and was thought to be the largest genome based-study ever performed in beef cattle at that time. This previous across-breed study was undertaken on 12 traits including birth weight, calving ease, carcass weight, and mature weight across 10 breeds and the results were similar to what we observed here for the muscularity traits. Saatchi et al. [15] identified 159 unique QTL associated with 12 traits, but only four QTL had pleiotropic effects and segregated in more than one breed. Similar results were observed in an across-breed study on dry matter intake, growth and feed efficiency in four beef cattle breeds [33]. The QTL identified for these traits were also breed-specific with little overlap among the breeds. This is comparable to our findings that show that the majority of the QTL were also trait-specific and breed-specific.
In total, approximately 83% of all QTL that are suggestively or significantly associated with a trait in our study overlapped with previously reported QTL associated with other production traits in dairy or beef cattle in the Cattle QTLdb (accessed 08 January 2019). Approximately 36% of all QTL overlapped with other traits that were specifically related to muscle in beef cattle such as body weight, carcass weight and marbling score [31], calving traits [37], Warner-Bratzler shear force [38], and longissimus muscle area [39]. One QTL on BTA17 that was associated with DIT in the SI breed was previously associated with ribeye area in a composite beef cattle breed composed of 50% Red Angus, 25% Charolais, and 25% Tarentaise [40]. Our study is further validated by the presence of significantly associated QTL regions on BTA2, which harbours the MSTN gene, with the five muscularity traits in the CH and LM breeds, and within the meta-analysis. In a previous study on five muscularity type traits, which were combined into one singular muscular development trait in CH, a QTL on BTA2, which contained MSTN, was the only region significantly associated with these traits [13].
In general, the suggestive and significant QTL, and thus genes, associated with each trait and each breed were both trait-specific and breed-specific. The low commonality of QTL among the breeds may be due to different genetic architectures underlying the traits in these breeds, or to gene-by-environment or epistatic interactions [33], or to differences in the power to detect QTL due to the large differences in population sizes between the breeds. In many cases, the significant alleles were simply not segregating in all five breeds. The differences between breeds may also be due to limitations in the imputation process with the imputation accuracy being too low to determine strong associations between a SNP and a trait; consequently, the minor suggestive associations were interpreted with caution because of the possibility of poor imputation. Overall, the largest number of overlaps among significant genes were found between the CH and LM breeds for all traits, which is not surprising considering the relative similarities in the origins of these breeds [41] and of the selection pressures they have experienced [42].

Myostatin
MSTN was first observed as a negative regulator of skeletal muscle mass in mice [43] and since then has been identified as responsible for muscular hypertrophy in cattle [44,45] and is widely known as the causal variant for many muscularity and carcass traits in cattle [46,47]. The stop-gain mutation Q204X in MSTN was significantly associated with the muscularity traits in both the CH and LM populations in the present study. Previously published research showed that CH and LM calves carrying one copy of this mutated allele scored better for carcass traits than non-carrier animals and that young CH bulls carrying this mutation presented a carcass with less fat and more tender meat than non-carriers [47]. In the present study, the CH and LM animals carrying one copy of the minor allele scored significantly (p < 0.01) higher for muscularity type traits. The Q204X mutation was not significant in the AA population and it was removed during the data-editing step in both HE and SI as it was non-segregating. When Q204X was included as a fixed effect in the model for the CH animals, no SNPs located within the MSTN gene itself remained significant. This indicates that the significant SNPs within this gene were in tight linkage disequilibrium with Q204X, which provides evidence that this mutation may be causative for the muscularity linear type traits in the CH breed. Other genes on BTA2 that were significantly associated with some or all of the traits in CH and LM were ORMDL1, PMS1, MFSD6, and NAB1, all of which are in strong linkage disequilibrium with MSTN in mammals [48].

Other candidate genes
While the major peaks on BTA2 in the analyses on CH and LM, and all the meta-analyses contain MSTN, a known contributor to muscle development, it is also plausible that other candidate genes within the QTL on BTA2 could also contribute to muscle development. Two such genes are COL3A1 and COL5A2. Intronic variants in COL3A1 and upstream and downstream gene variants in COL5A2 were significantly associated with DHQ in both CH and LM; however, no SNPs within coding or non-coding regions of this gene were associated with any traits in AA, HE, or SI although the SNPs were indeed segregating. Collagen is abundant in muscle and the quantity and stability of these intramuscular fibres have previously been linked to eating palatability of beef [49]. The quantity and stability of muscle collagen are known to differ by breed [50], sex [51], and age [52] of cattle. Other collagen genes, COL6A1, COL6A2, and COL18A1, on BTA1 were also identified as candidate genes for DIT in the AA breed. Both type VI collagen genes have previously been linked to various muscle disorders in humans since they are known to affect muscular regeneration [53]. Type XVIII collagen has previously been proposed as a useful marker for beef marbling because it is involved in fat deposition in ruminants [54].
Another QTL on BTA2 located in the region between 13.9 and 14.9 Mb and significantly associated with four of the traits (DHQ, DIT, TW, and WOW) in the LM breed contained the PDE1A and PPP1R1C genes. The most significant SNP in this region was an intronic SNP within PDE1A. The PDE1A gene is involved in a pathway related to myofibroblast formation in smooth muscle in humans [55] while previous genome-wide studies in mice have identified the PPP1R1C gene as a possible candidate gene for muscle mass [56]. Overall, the allele frequencies of the favorable alleles in this 1-Mb region were similar in all five breeds, which support a breed-specific association with DHQ, DIT, TW, and WOW in LM rather than an imputation error.
An additional breed-specific QTL on BTA2 that contains numerous HOXD genes was associated with WOW in the LM population. The HOXD genes are documented as having a role in limb [57] and digit [58] formation, thus they probably also play a role in skeletal muscle development. The most significantly associated SNPs with WOW in this region were only segregating in the LM breed and had a very high favorable allele frequency (0.998) in this breed. These SNPs were fixed or very close to fixation in the four other breeds.
In the meta-analyses of DHQ, associated variants in all the breeds analysed were identified, which may be beneficial for across-breed genomic prediction [59]. Although the associations detected in the meta-analysis corresponded to associations identified in the CH and LM breeds, three of these QTL on BTA5, 11, and 12 increased in significance when compared to the within-breed analysis. The QTL on BTA5 which contained the AMDHD1 gene, was located close to a QTL previously associated with carcass composition [43], whereas the QTL on BTA11 contains DYSF, a gene known to be linked with muscular dystrophy in humans [60]. The QTL on BTA14 contained the PREX2 gene which was previously linked to carcass weight in Hanwoo cattle [61].
Interestingly, in the meta-analyses of DL and WOW, a 1-Mb QTL on BTA7 containing the SPRY4 and FGF1 genes became suggestively associated, although it was not associated in any breed individually. The SPRY4 gene was reported to be associated with feed intake in cattle [62], whereas FGF1, a member of the fibroblast growth factor family, is thought to be involved in embryonic muscle formation [63].
Similarly, in the meta-analysis of TW, a 3-Mb QTL on BTA6 containing the NCAPG/LCORL genes became suggestively associated, although it was not associated in any breed individually. These genes are associated with variation in body size and height in cattle [32], humans [64], and horses [65], thus they are likely plausible candidate genes associated with muscularity linear type traits describing the size of the body.

Gene ontology and KEGG pathways
Linear type traits are complex traits that are governed by many genes each with a small effect, and hence, are likely involved in many biological systems. Several GO terms were only associated with a single trait or a single breed; hence there was limited commonality among traits or breeds suggesting the absence of a central biological process that links these traits together. Over-represented GO terms in multiple traits and breeds include those related to skin development, collagen fibril organisation, and the transforming growth factor beta receptor signalling pathway. Each of these GO terms was associated with genes located in the large QTL on BTA2 that contained MSTN. Excluding the major MSTN QTL in these breeds, which is known to have a large effect on muscularity, the various GO terms and KEGG pathways represented by the genes associated with the muscularity traits suggest that the majority of genes identified as significantly associated with a trait are not only breed-specific but also trait-specific in many cases.

Regulatory regions involved in the development of muscle
Although millions of SNPs were tested for association with each trait, only 79 of the SNPs suggestively or significantly associated with a trait were located in the coding region of a gene; the vast majority of the SNPs associated with the muscularity traits in any of the breeds were located outside of the coding regions. This is consistent with previous genomic studies for complex quantitative traits in cattle using HD SNP data [66] or sequence data [32]. While the coverage of the HD study [66] may not have included the coding regions required to identify significant associations within these regions, our study and a previous study on cattle stature [32] used imputed sequence data, and thus, covered the entire genome.
Whereas many studies have previously acknowledged the importance of non-coding SNPs to genetic variability, little is actually known about the mechanisms by which these SNPs contribute to variation in complex traits [67,68]. One possibility to explain the significance of these non-coding SNPs is that the noncoding regions contain gene regulatory sequences, called enhancers, that act over long distances possibly altering the expression of a gene nearby [67]. Another possibility is that the folding of DNA into the 3-dimensional nucleus may cause distant loci, such as those in non-coding and coding regions, to become spatially close together thus enabling these regulatory regions to come into contact with genes far away or even on different chromosomes [69].
Non-coding variants such as 3′ UTR, 5′ UTR and intergenic variants were enriched for most of the traits in each breed. Downstream and upstream gene variants were also enriched in some traits. In general, the SNPs located close to and within the genes identified as candidate genes were located within non-coding or regulatory regions. For example, for DHQ in the CH breed, 60 suggestively and significantly associated SNPs were located within the MSTN gene; 10 of these were 3′UTR variants, 31 were downstream gene variants and 19 were intronic. Whereas regulatory regions may not have an effect on the coding sequence of any gene, they are thought to be particularly important for growth and development in humans [68,69] and cattle [32,70]. Thus, similar to previous observations in humans and cattle, enrichment of the non-coding classes of SNPs in our study may indicate the importance of regulatory regions for cattle muscle development.

Conclusions
Although we identified many QTL associated with muscularity in beef cattle, our results suggest that these QTL tend to be not only trait-specific but also breed-specific. Overall, the significant SNPs contained in these QTL were more likely located in regulatory regions of genes, which suggest the importance of non-coding regions that may affect gene expression for muscle development in cattle. Some shared regions associated with muscularity were found between CH and LM, with a largeeffect QTL on BTA2 containing MSTN being associated with the five traits analysed. This overlap between these breeds was somewhat expected, because they are subjected to similar selection pressures. Apart from this single QTL, extensive differences were observed between the breeds, which may be due to the much smaller sample sizes for AA, HE, and SI compared to the CH and LM populations that result in reduced power to detect QTL or they may be due to differences in genetic architecture of these traits among the populations. In many cases, the strongly associated SNPs in one breed were