Genome-Wide Association Analysis for Liveweight Traits in Braunvieh Cattle

A genome-wide association study (GWAS) for liveweight traits of Braunvieh cattle was performed. The study included 300 genotyped animals by the GeneSeek® Genomic Proler Bovine LDv.4 panel; after quality control, 22,734 SNP and 276 animals were maintained in the analysis. The examined phenotypic data considered birth, weaning, and yearling weights. The association analysis was performed using the principal components method via the egscore function of the GenABEL version 1.8-0 package in the R environment. The marker rs133262280 located in BTA 22 was associated with birth weight, and two SNPs were associated with weaning weight, rs43668789 (BTA 11) and rs136155567 (BTA 27). New QTL associated with these liveweight traits and four positional and functional candidate genes potentially involved in variations of the analyzed traits were identied. The most important genes in these genomic regions were MCM2 (minichromosome maintenance complex component 2), TPRA1 (transmembrane protein adipocyte associated 1), GALM (galactose mutarotase), and NRG1 (neuregulin 1), with relationships with embryonic cleavage, bone and tissue growth, cell adhesion, and organic development. This study is the rst to present a GWAS conducted in Braunvieh cattle in Mexico and represents a basis for future research. Further analyses of found associated regions will clarify its contribution to the genetic basis of growth-related traits.


Introduction
The identi cation of causal genetic variability is one of the main goals in the genetic improvement of cattle. Commonly, liveweight traits are used as the primary selection criterion in cow-calf production systems in Mexico (Jahuey- Martínez et al., 2016). Usually, farms use these traits as e ciency and meat potential production indicators, and they are used for genetic evaluations in most of the registered breeds (Parra-Bracamonte et al., 2015).
Braunvieh is a worldwide cattle breed used in the beef industry that has been used in both specialized beef and dualpurpose production herds (Orantes-Zebadúa et al., 2014;Phillips et al., 2009). Due to its initial dual-purpose origin, most of the available information about the Braunvieh deals with dairy production traits. However, during the last 15 years, Braunvieh cattle have been selected and genetically improved for beef production traits (Chin-Colli et al., 2016;Phillips et al., 2009). In Mexico, Braunvieh is one of the breeds most used for the beef production industry either as purebred or in crosses with Bos indicus cattle (AMCGSR, 2017, Orantes-Zebadúa et al., 2014. However, despite its extensive use, there is scarce information on the breed's productive performance, and the available information is mainly related to growth traits coming from genetic evaluations or isolated studies (AMCGSR, 2017;Chin-Colli et al., 2016;Silva et al., 2002). The selection, management, and genetic improvement programs of the Braunvieh cattle could be enhanced using high-throughput genotyping technologies.
The use of microarrays of thousands of SNP markers in genome-wide association studies (GWAS) has allowed discovering the genetic basis of complex traits and diseases by detecting genotype-phenotype associations in a group of individuals (Tian et al., 2020). GWAS approaches have con rmed many QTL for growth traits in beef and crossbred cattle (Jahuey- Martínez et al., 2016;Lu et al. 2013;Martínez et al., 2016), some of which have been used as the basis for the search for speci c causal variation (Takasuga, 2016) and a better understanding of the genetic architecture of these complex traits. Many of these QTL, genomic regions and genes, affecting production traits in beef cattle has been reported (Jahuey et al., 2016;Londoño-Gil et al., 2021;Pur eld et al., 2015;Rolf et al., 2012), but most of the association studies focused on specialized beef breeds and only a few researches have been implemented in breeds such as Braunvieh (Guo et al., 2012;Maxa et al., 2012). The present study is aimed to perform a GWAS to identify QTLs and candidate genes related to liveweight traits in a Braunvieh cattle population.

Materials And Methods
Approval from the ethical committee for animal care and use was unnecessary because the samples used in this study consisted of hair follicles, and all analyses were performed using pre-existing databases.

Population and phenotypic data
Hair follicle samples from 236 females and 64 males registered in the Mexican Braunvieh Cattle Association database was collected. The cattle were born between 2000 and 2015. This population came from herds located in the east, west, and central highlands of Mexico. Herds from west and east were raised under extensive production systems, while central highlands herds were under intensive regimens. The sampled population's genetic background included Austrian, Swiss, Canadian, American, and Mexican animals. Phenotypic data were provided by the breeding association and included records of birth weight (BWT, kg), weaning weight (WWT, kg), yearling weight (YWT, kg).
Weaning and yearling weights were adjusted to perform the GWAS analysis. Table 1 shows the descriptive statistics for each trait. Genotyping and quality control The animals were genotyped using 30,125 SNP markers from the GeneSeek® Genomic Pro ler Bovine LDv.4 panel (Neogen Corp. Lincoln, NE, USA). Before association analysis, the genotypic data quality was veri ed using the SNPQC program (Gondro et al., 2014). The genotypes were considered successful if they presented a GenCall value greater than 0.50, and all SNPs with lower values were discarded (n = 1623). Those SNPs that were monomorphic (n = 3604), presented call rates of less than 90% (n = 1290) or minor allele frequencies < 0.01 (n = 1325), or deviated from Hardy-Weinberg equilibrium according to Fisher's exact test and exhibited P-values > 1 × 10 -15 (n = 0) were also eliminated. Besides, SNPs with unknown coordinates in the assembly of the bovine genome UMD v3.1 (Zimin et al., 2009) (n = 1484) and SNPs that were not located on autosomal chromosomes (n = 1820) were discarded.
Samples were also eliminated if they exhibited call rates of less than 80% (n = 0) or levels of heterozygosity (HE) above 3 SD (n = 1), considering that the mean and SD of the observed HE were 0.32 and 0.019, respectively. A Pearson correlation was computed for detecting potentially duplicate samples, considering a maximum of r = 0.98, according to their genotype information obtaining an average of r = 0.817 and minimum and maximum values of 0.66 and 0.90, respectively. A total of 22,734 SNPs and 276 samples passed the quality control procedure and were retained for further analysis. Quality control and subsequent analyses were performed in the R environment.
Population structure and association analysis Population structure was analyzed, calculating rst a genomic relationships matrix using the information on genotypes according to Van Raden (2008), besides performing a singular value decomposition and a principal components (PC) analysis.
The PC analysis indicated that the rst two PCs explained 28.6% of the variance in the data. Therefore, the genomewide association analysis was performed using the PC method proposed by Price et al. (2006). For this analysis, the egscore function from the GenABEL package (Aulchenko et al., 2007) was employed. This function accounts for population strati cation and uses the genomic kinship matrix to derive axes of genetic variation, and then both the phenotypes and genotypes are adjusted onto these axes.
A linear model for each trait was tted, including the rst two PC as covariates. For the analysis of BWT, the model also included the contemporary group (CG) and the linear and quadratic effects of cow age at the birth and weaning of her calf. The CG included herd, sex, year, and calving season. The statistical model used to adjust the other traits only included the CG and the PCs as covariates; cow age was excluded because it was not signi cant in the previous analysis. Finally, the association between corrected genotypes and phenotypes was assessed via correlation. Pvalues were obtained by calculating the square of the correlation multiplied by (N-K-1), where N was the number of genotyped individuals, and K was the number of PCs.
Minimum allele frequencies, allele substitution effect (β), and percentage of phenotypic variance explained by the SNP were estimated. SNP with P-values < 5 × 10 − 5 were considered signi cantly associated with studied traits. The proportion of phenotypic variance explained by the SNPs was estimated by dividing the X 2 value for a df by the number of individuals used to analyze each SNP marker, followed by multiplication by 100. All described analyses and estimations were performed using the GenABEL package (Aulchenko et al., 2007).

Analysis of genomic regions with signi cant SNPs
The closest genes to signi cant markers and those located within a 250-kb window on both SNP location sides were identi ed. The list of genes was obtained using the snp2gene.LD function from the Postgwas package (Hiersche et al., 2013). Distance between SNPs and genes was calculated as the difference between the marker position and the beginning or end of the gene, according to coordinates from bovine genome assembly UMD v3.1. Gene functions were investigated in the UniProt database (The UniProt Consortium, 2021).
Annotations from humans or mice were used when there was no information on the genes in cattle. Genes were considered functional and positional candidates if they were biologically related to the trait under study, supported by experimental evidence in the literature. Finally, we determined whether signi cant SNPs mapped against QTLs previously associated with growth-related traits such as BW, carcass, and reproduction traits, deposited in the cattle AnimalQTLdb (Hu et al., 2013). For this purpose, SNP positions according to the Btau4.6 genome sequence were used because many of the previously reported QTLs had no well-de ned positions in the bovine genome assembly UMD v3. According to the signi cance threshold considered (P < 5 × 10 − 5 ), 3 SNP were associated with the live weight traits (  Fig. 2). The rs133262280 located in BTA22 was associated with BW, showing an allelic substitution effect of 0.320 ± 0.02 kg. The rs43668789 and rs136155567, located in BTA11 and 27, respectively, were associated with WW.
These markers showed allelic substitution effects of -9.590 ± 0.25 and 1.110 ± 0.72 kg, respectively (Table 2).  QTLs previously associated with growth-related traits are shown in Table 3. Tables 4 to 6 show complete descriptions, including the identi er number and exact location identi ed in this study.

Discussion
The inclusion of the population's genetic structure into the analysis model allowed the better tting of the GWAS model for all traits, as showed by quantile-quantile plots (Fig. 1). This genetic structure was expected because tested herds presented different selection criteria, and perhaps, ancestors from the imported genetic material (i.e., semen, sires). Strati cation results could include extensive use of sires or semen that breeders usually choose in their genetic improvement programs. Some studies (Erbe et al., 2012;Plieschke et al., 2015) have used subdivisions to estimate QTLs using genome-wide association studies (GWAS). Smitz et al. (2014) concluded that the strati cation in the studied populations needs to be considered in genetic improvement programs to conserve those populations' "genetic health". Jemaa et al. (2015) indicated that some QTLs found in GWAS could not be present in all the studied animals due to the population's strati cation.
Birth weight in Braunvieh cattle represents an important trait to consider in the genetic improvement programs due to its association with calving di culty in young heifers, especially when the Braunvieh is used as a sire for smaller-size breeds (Hagger & Hofer, 1990). In the present study, the rs133262280 was identi ed as the only marker associated with BW, located at 60.7 Mb of BTA 22. This SNP showed an allelic substitution effect of 0.320 kg, explaining 0.1% of the phenotypic variance of BW. Genes located closer to this SNP included CHCHD6 (coiled-coil-helix-coiled-coil-helix domain containing 6), MCM2 (mini-chromosome maintenance complex component 2), PLXNA1 (plexin A1), PODXL2 (podocalyxin like 2), TPRA1 (transmembrane protein adipocyte associated 1), and uncharacterized LOC10105309 (Table 4).
The most important genes identi ed in this region were MCM2 and TPRA1. MCM2 gene is located at 177.6 kb and TPRA1 at 160.1 kb; both genes are upstream of the rs133262280 SNP. MCM2 acts as a component of the MCM2-7 complex (MCM complex) which is the putative replicative helicase essential for 'once per cell cycle' DNA replication initiation and elongation in eukaryotic cells (Todorov et al., 1994). Additionally, it plays a role in cell division and apoptosis (Gao et al., 2015). Gao et al. (2015) reported MCM2 protein expression in the cochlea of rats and guinea pigs slightly increase the apoptosis rate of the cells without any changes in proliferation or cell cycle. Recently, Khan et al. (2020) found by a transcriptomic analysis that supplementation with folic acid in perinatal Holstein cows signi cantly increases the expression of the MCM2 gene.
The other associated gene with biological importance was the TPRA1 gene belongs to the G protein-coupled receptor (GPCR) family. Functions related to this gene include regulating early embryonic cleavage and enhancing the hedgehog signaling pathway (Aki et al., 2008;Singh et al., 2015). Several studies have highlighted its importance in pre-and perinatal tissue development in mice. Aki et al. (2015) determined that the TPRA1 gene in uenced the hedgehog signaling pathway, which plays an essential role in vertebrate embryonic tissue patterning of many developing organs, showing differences around 50% in the signaling levels comparing homozygotes and heterozygotes animals.
This evidence suggests that MCM2 and TPRA1 could participate in the early stages of cattle development and, therefore, in uence BW. There were no quantitative trait loci previously located in this region, which could be a speci c QTL of the studied population.
The present study identi ed two regions (Table 3,  Besides, Boichard et al. (2003) and Buitenhuis et al. (2007) reported associations between the identi ed regions in this study and conformation traits, explaining between 5.9 to 8.9 % of the structural soundness in ten European dairy cattle breeds. On the other hand, Sherman et al. (2009) and Rolf et al. (2012) reported associations with allele substitution effects between − 0.319 to 2.199 kg for feeding traits like average daily gain and residual feed intake in Angus, Charolais, and Canadian beef hybrid cattle.
The most important gene identi ed in this region was GALM. This gene is located 217.4 kb upstream of the rs43668789 and belongs to the proteins that convert the α-aldose to β-anomer. GALM is involved in the pathway hexose metabolism, which is part of carbohydrate metabolism (Thoden et al., 2004). McClure et al. (2010) reported a positive association of GALM with the weaning weight in Angus cattle. Shin et al. (2014) mentioned that the association between GALM and the weaning weight in Holstein and Hanwoo cattle lies in quantity and the quality of the calves' milk consumption. Quantitative trait loci located in this region have been previously associated with weaning weight in Angus (McClure et al., 2010), conformation in dairy cattle breed (Boichard et al., 2003;Buitenhuis et al., 2007), and residual feed intake in Canadian beef synthetic cattle (Sherman et al., 2009).
The second marker associated with WW was rs136155567, located at 27.0 Mb of BTA 27, and its allele substitution effect was 1.110 kg which explains 1.1% of the phenotypic variance. Genes located closer to this SNP (± 600 kb) included LOC104976093 (uncharacterized LOC104976093) and NRG1 (neuregulin 1) ( Table 5). NRG1 was the most important gene identi ed. This gene is located at 567.1 kb downstream of the rs136155567. It is considered the direct ligand for ERBB3 and ERBB4 tyrosine kinase receptors. The multiple isoforms perform diverse functions such as inducing growth and differentiation of epithelial, glial, neuronal, and skeletal muscle cells and in uence motor and sensory neuron development (Ieguchi et al., 2010;Plowman et al., 1993). In cattle, NRG1 has been highly associated with organ development (Sweeney et al., 2001). Zhao (2013) mentioned that this gene could in uence the weaning weight as an emerging regulator of prolactin secretion.
In general, the phenotypic variance explained by the SNPs identi ed by this study was marginal (1.39 % on average).
In growth trait studies, it is expected that most SNP markers will explain only a tiny proportion of the observed phenotypic variance due to the polygenic control over such traits and because individual genes only slightly in uence a phenotype. However, consideration of SNPs' sets that are signi cantly associated with each trait may allow a greater proportion of phenotypic variance to be explained. For example, the two SNP associated with WW could explain 4.08 % of the variance in that trait. However, the present outcomes increase knowledge of the genetic architecture of live weight traits important in beef cattle production.
In conclusion, in the present study, three SNP were associated with live weight traits of Braunvieh cattle. Two SNPs were located in intergenic regions, and one was located in an intronic region of the ARHGEF33 gene. Evidence shows that some of the genes closer to the three identi ed SNPs markers are functionally related to growth through embryonic cleavage, bone and tissue growth, cell adhesion, and organ development. There were four candidate genes with potential associations with assessed live weight traits in Braunvieh cattle, including MCM2, TPRA1, GALM, and NRG1. Subsequent studies examining these genomic regions could lead to the identi cation of polymorphisms with potential uses in the marker-assisted selection, providing a deeper understanding of the genetic basis of growth traits in cattle. This study represents the rst study to describe a GWAS conducted in Braunvieh cattle in Mexico. Further analysis using the present information would allow conducting assessments on the ontogeny and speci c search of causative mutations for live weight traits. Furthermore, examining particular and general genic effects would indicate the possibility of including genomic information into current genetic evaluations.

ACKNOWLEDGMENTS
Authors thank the National Council of Science and Technology (CONACYT) for supporting the rst author's doctorate studies grant and to Universidad Autónoma Chapingo for the research grant DGIP-166701012. The authors also thank the Asociación Mexicana de Criadores de Ganado Suizo A. C. for their collaborative support.

Funding
This study was funded by research grant DGIP-166701012 by the Universidad Autónoma Chapingo.

Con icts of interest/Competing interests
The authors have no con ict of interest to declare Availability of data and material (data transparency) Data could be available upon reasonable request to the authors.
Code availability (software application or custom code) Code used for analysis could be available upon reasonable request to the authors.

Ethics approval
No live animals were used in this study, then was not necessary the ethics approval.

Consent to participate
Not applicable

Consent for publication
All the authors have read the content and consented to submit the manuscript Figure 1 Quantile-quantile (QQ) plots for the genome wide association study of live weight traits in Braunvieh cattle. The straight line in the QQ plots indicates the distribution of SNP markers under the null hypothesis, and the skew at the edge indicates that these markers are more strongly associated with the traits than would be expected by chance. BW = birth weight; WW = weaning weight; YW = yearling weight.

Figure 2
Manhattan plots of the P-values for the genome-wide association study of growth traits in Braunvieh cattle. The horizontal line indicates the signi cance threshold for signi cant associations (P < 5 × 10−5). BW = birth weight; WW = weaning weight; YW = yearling weight.