Genomic Breed Composition of Selection Signatures in Brangus Beef Cattle

Cattle breeding routinely uses crossbreeding between subspecies (Bos taurus taurus and Bos taurus indicus) to form composite breeds, such as Brangus. These composite breeds provide an opportunity to identify recent selection signatures formed in the new population and evaluate the genomic composition of these regions of the genome. Using high-density genotyping, we first identified runs of homozygosity (ROH) and calculated genomic inbreeding. Then, we evaluated the genomic composition of the regions identified as selected (selective sweeps) using a chromosome painting approach. The genomic inbreeding increased at approximately 1% per generation after composite breed formation, showing the need of inbreeding control even in composite breeds. Three selected regions in Brangus were also identified as Angus selection signatures. Two regions (chromosomes 14 and 21) were identified as signatures of selection in Brangus and both founder breeds. Five of the 10 homozygous regions in Brangus were predominantly Angus in origin (probability >80%), and the other five regions had a mixed origin but always with Brahman contributing less than 50%. Therefore, genetic events, such as drift, selection, and complementarity, are likely shaping the genetic composition of founder breeds in specific genomic regions. Such findings highlight a variety of opportunities to better control the selection process and explore heterosis and complementarity at the genomic level in composite breeds.


INTRODUCTION
Breeding methods that exploit heterosis are common in livestock production. In cattle, the challenge for adopting terminal crossbreeding systems is consistent genetic composition of replacement heifers from the maternal breed (Lightner and Williams, 2018). Composite breeds (also referred to as synthetic breeds) allow for consistency in heterosis retention and heifer replacement. Brangus, developed in the United States, are an example of a composite breed, defined as 62.5% Angus and 37.5% Brahman (International Brangus Breeders Association 1 ). The breed represents the complementarity between the tropically adapted Bos taurus indicus and the temperate high-valued carcass of Bos taurus taurus cattle (Gregory and Cundiff, 1980;Buzanskas et al., 2017). Brangus registration by the International Brangus Breeders Association started in 1949.
After a composite breed is formed, a genetic improvement program can be applied, selecting animals across generations based on expected progeny differences (EPD) for specific traits. The United States Brangus breeder association (IBBA) has developed EPDs for birth weight, weaning weight, yearling weight, milk production, total maternal, calving ease direct, calving ease maternal, scrotal circumference, ribeye area, and intramuscular fat (International Brangus Breeders Association; 1 ). Therefore, artificial selection pressure, at varying levels of intensity, would have been employed on these traits. In this process, the inbreeding level can increase due the selection of few parents, especially bulls, and, consequently, decreasing heterosis. In this scenario, genomic selection signatures may arise after composite formation (Goszczynski et al., 2017). These genomic regions with selective sweeps may have different genomic breed composition than expected due to selective advantages of genes coming from one of the founder breeds.
The evaluation of selection signatures and genomic breed composition in composite breeds can contribute to a better understanding of the genetic effects associated with traits under selection and the inheritance of loci in crossbreeding systems (Grigoletto et al., 2019). Concerning the dynamics of composite breed development (Paim et al., 2020), we can gain new insights for crossbreeding systems based on a genomic perspective. The aim of this work was to expand our knowledge of composite breed genomics by identifying genomic inbreeding and selection signatures in Brangus. Further, we aimed to evaluate the genomic breed composition of these selected regions, identifying differential founder (Angus or Brahman) contributions to that region.

Pedigree Evaluation
The pedigree file was evaluated using the optiSel package (Wellmann, 2017) in R 3.4.2 software (R Core Team, 2017). The Angus, Brahman, and crossbred animals (with pedigree breed composition other than the 5/8 Angus, 3/8 Brangus) were considered as ancestors, totaling 75,449 ancestors in the pedigree file. The number of equivalent generations for each animal (hereinafter called generations) was calculated by the equation: where g is the equivalent generation number and n is the number of generations separating the individual from each known ancestor. The method used is similar to the equation described by Welsh et al. (2010).
A summary of the pedigree analysis of the Brangus bulls used is shown in the Supplementary Material (Supplementary Figure 1). The index of pedigree completeness (PCI) was 0.94 (±0.143 SD), computed following the MacCluer et al. (1983) algorithm. PCI is the harmonic mean of the pedigree completenesses of the parents, summarizing the proportion of known ancestors in each ascending generation. Pedigree inbreeding of 0.04 (±0.035) was found for the breed. The average pedigree relationship was 0.086 (±0.081), and only 1.69% of the pairs had a pedigree relationship higher than 0.3.

Filtering and Quality Control of Genomic Data
Markers with a call rate lower than 95% or not physically mapped to the bovine genome assembly UMD3.1 were removed from the analyses. The remaining genotypes were 698,282 SNP markers on the autosomes and 38,581 SNP on the sex chromosomes (37,538 in X and 1,043 in Y). Markers with minor allele frequency lower than 1% were removed. One Brangus sample with a call rate lower than 90% was removed.

Runs of Homozygosity and Selection Signatures
The runs of homozygosity (ROH) analyses were conducted in SNP and Variation Suite R v8.7 (Golden Helix, Inc., Bozeman, MT, 2 ). The parameters were set to a minimum run length equal to 1000 kb with minimum of 70 SNP, allowing runs to contain up to two heterozygotes and five missing genotypes with a maximum gap equal to 50 kb and minimal SNP density of 1 SNP per 50 kb.
The minimal number of SNP to constitute a ROH (l) was determined by the same method used by Purfield et al. (2012) and determined by Lencz et al. (2007): where ns is the number of SNPs per individual, ni is the number of individuals, α is the percentage of false positive ROH (set to 0.01 in this study), and het is the mean SNP heterozygosity across all SNP. The incidence of common ROH was transformed to each breed's frequency, dividing by the number of animals of each breed in the analysis. Normality tests were performed, and the frequency threshold defining the top 1% of the observations for each breed was determined. The homozygous regions above the frequency threshold of each breed (38% for Angus, 25.4% for Brahman, and 25.9% for Brangus) were considered as selected regions.
According to the length of the ROH, it is possible to estimate the number of generations traced back to the common ancestor, which generates the homozygosity in that region. We classified the ROH into 4 classes (1 = more than 10 generations, 2 = between 5 and 10, 3 = between 3 and 5, and 4 = less than 3 generations) using the equation proposed by Curik et al. (2014): is the expected length of an identical by descendent (IBD) haplotype (in centiMorgans -cM), and gcA is the number of generations from the common ancestor. The conversion from the recombination rate metric to physical distance (from cM to Mb) was performed using the average of the results of Arias et al. (2009) and Weng et al. (2014). Based on the Curik et al. (2014) equation, for example, an ROH longer than 13 Mb has most likely originated from a common ancestor less than three generations ago.
A genomic inbreeding coefficient based on ROH (F ROH ) was calculated on each animal according to McQuillan et al. (2008) with the equation where L ROHj is the length of ROH j , and L total is the total size of the autosomes (using the estimated value in the Btau5.0.1 genome assembly of 2,522,199,562 bp). For each animal, F ROH was calculated based on each of the four classes explained before and for each chromosome using the total size of each chromosome as L total (following the chromosome size estimated by the Btau5.0.1 genome assembly).

Chromosome Painting
We used the copying model, implemented in ChromoPainter (Lawson et al., 2012), to estimate the ancestry of regions across each genomic region. This copying model relates the patterns of linkage disequilibrium (LD) across chromosomes to the underlying recombination process. The method uses a hidden Markov model to reconstruct a sampled haplotype.
To reinforce chromosome-painting results, we ran Fst analyses (Weir and Cockerham, 1984) for each region comparing the pairs (Angus vs. Brangus and Brahman vs Brangus). The function-fst in the plink1.9 software 3 was used. We used the founder breeds, Angus and Brahman, as haplotype donors to the Brangus haplotypes. The ChromoPainter analyses were performed twice (allowing or not allowing selfcopying) using the linked model. The recombination files were created using the Perl scripts provided on the ChromoPainter website 4 . Beagle3.3 (Browning and Browning, 2007) was used to phase the genotypes (using 20 iterations).

Simulation Model
We performed a population genetics simulation using the online tool 5 . The initial parameters were set to an initial allele frequency of 62.5% (representing the Angus allele in the first generation of Brangus); 10 generations; effective population size of 100; no selection, mutation, migration, and inbreeding (similar to a neutral model). We performed 50 simulations for each generation. The raw data were used to calculate the summary statistics (mean and standard deviation) and to determine the expected lower and upper value (within 99% of the Gaussian distribution) of the expected founder composition for each locus. These lower and upper values were applied as a threshold in the visualization of chromosome painting results to identify regions with significant enrichment of alleles coming from one of the founders.

Identification of Genes and QTL in Selective Signatures
Genes in the selected regions (ROH islands) were identified in the Golden Helix GenomeBrowse R visualization tool v2.1 (Golden Helix, Inc, Bozeman, MT, 6 ). The genes were identified based on the NCBI Bos taurus annotation release 105 and Btau5.0.1 genome assembly. The genes list obtained was submitted in the NetworkAnalyst online tool 7 , aiming to characterize the biological process of these genes through the Enrichment Network tool using the PANTHER database. Thereafter, a search in the literature and in the Cattle QTL database (available online at 8 ) was executed to identify traits related to genes located in each significant genomic region.

RESULTS
The runs of homozygosity (ROH) were categorized into four classes according to the expected number of prior generations to a common ancestor (>10, >5, >3, and <3 generations). The ROH classified as coming from a common ancestor within the previous FIGURE 1 | Runs of homozygosity (ROH) length observed in Brangus cattle according to the equivalent generation number of each animal. Dashed lines indicate the length threshold for ROH that relates to a common ancestor at 3, 5, and 10 prior generations (red, blue, and green, respectively) following the equation proposed by Curik et al. (2014). animals had long ROH (>13 Mb) indicative of recent inbreeding. However, chromosomes 17, 23, 26, and 28 did not have any ROH in this length range (Supplementary Figure 1).
The genomic inbreeding coefficient based on ROH (F ROH ) was significantly (p < 0.0001) higher for Angus cattle compared to Brahman and/or Brangus (Figure 2). Brangus had lower F ROH than Angus for all classes. Brahman and Brangus cattle had the same F ROH for the ROH coming from a common ancestor tracing through 10 generations (all classes with ROH > 3.9 Mb), which was not expected and suggests a high effective population size for Brahman. For ROH coming from more than 10 previous generations (ROH < 3.9 Mb), Brangus cattle had higher F ROH than Brahman cattle.
Pedigree inbreeding had a positive and significant relationship with F ROH , and a similar pattern was observed in all the classes of ROH length (Supplementary Figure 2). Brangus cattle had 8.5 ± 3.97% of genomic inbreeding and 3.9 ± 3.41% of pedigree inbreeding. Animals with no inbreeding at pedigree FIGURE 2 | Genomic inbreeding based on runs of homozygosity (F ROH ) by breed and by ROH length classes. The t-test comparison results are shown in the top (ns: not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).
had close to 6% of genomic inbreeding. Averaged across chromosomes, the rate of F ROH increased ≈1% per generation in Brangus (F ROH = 0.0196 + 0.0097 * generation, R 2 adj = 0.19, p-value = 0.0004). The increase in F ROH was not observed for all chromosomes; only chromosomes 4, 10, 13, 15, 23, 26, and 29 had a positive F ROH slope with generation number (Supplementary  Figure 3). All the aforementioned chromosomes, except for 13, had a high proportion of Angus composition.
Ten genomic regions had ROH with frequency higher than 25.9% in Brangus (the top 1% of ROH frequency). Two of the 10 regions were found to be ROH islands for both founder breeds, and three ROH islands were observed in Angus (Figure 3). ROH above a 1% threshold were identified in 10 and 21 regions for Brahman and Angus, respectively (Supplementary Figure 4).
The genes and known QTLs within homozygous Brangus regions are shown in Table 1. The main biological process observed in gene network enrichment analysis from these homozygous regions were bile acid metabolic process, fatty acid beta oxidation, pentose phosphate shunt, neuron synaptic transmission, protein folding, regulation of cell cycle, cholesterol metabolic process, and unsaturated fatty acid biosynthesis. The main traits observed in QTL analysis of these regions were body weight, milk fat, calving ease, milk production, milk protein, body weight at birth, and fat thickness at the 12th rib (Supplementary Figure 5). The breed of origin of these regions was investigated using chromosome painting (Figure 4). F ST results show Brangus had a closer relationship with Angus than with Brahman in these ROH regions ( Table 1). The haplotypes in the regions of chromosomes 1, 4, 22, 26, and 27 appear to have originated from Angus. The regions in chromosomes 8, 14, 16, 21, and 23 have a mixture of Angus and Brahman origin, falling within the range of expected ancestry based on the whole genome.

DISCUSSION
Overall, Brangus had 63% of the inbreeding level of Angus, based on runs of homozygosity (ROH) (Figure 2). Brangus had higher inbreeding than Brahman only in the shortest category of ROH (<3.9 Mb), which suggests the number of generations after crossbreeding was not sufficient to break down short ROH. Brangus had a higher Angus proportion (70.4%) than expected (62.5%) in the whole genome (Paim et al., 2020). The high Angus proportion might be related with the initial crossbreeding to develop the composite associated with genetic drift and selection for specific traits (Paim et al., 2020). Therefore, the higher Angus proportion may be linked to this excess of short ROH in Brangus compared to Brahman. Moreover, inbreeding across chromosomes was not equal; this may suggest that new levels of homozygosity are starting to form as a function of selection pressure and the use of sires that are deemed superior to their contemporaries.
The length of ROH agree with the generation criteria of Curik et al. (2014). The ROH coming from a common ancestor within 3 prior generations (>13 Mb) appeared between the 4th and 5th generation and increased afterward, suggesting inbreeding has started to accumulate in this relatively new breed.
The inbreeding level increased approximately 1% per generation corresponding to an effective population size (Ne) of 51.55 (Ne = 1 2 F ) (FAO, 2013). According to FAO Guidelines for in vivo conservation of animal genetic resources (FAO, 2013), the desired inbreeding rate per generation should not exceed 1% (equal to Ne = 50). A 1% increase in inbreeding was associated with decrease of −0.23% in yearling weight and −0.64% in body condition score in a tropical composite beef cattle (Reverter et al., 2017). Therefore, selection pressure and finite population  (Weir and Cockerham, 1984). e Number of genes inside the region. f Number of QTLs identified in Cattle QTL database. g Number of traits associated with the QTLs. h Genes in the region that area associated with a trait (traits of each gene between parenthesis).
Frontiers in Genetics | www.frontiersin.org  size promotes increased inbreeding, suggesting that inbreeding management remains important for composite breeds. Three selected regions in Brangus (chromosomes 4, 16, and 23) were identified as Angus selection signatures, and two regions on chromosomes 14 and 21 were identified as selection signatures in both founder breeds (Figure 3). Chromosome painting results showed that five of the 10 homozygous regions in Brangus were predominantly Angus in origin (probability >80%), and the other five regions were of mixed origin but always with Brahman contributing less than 50% (Figure 4).
The traits associated with the predominantly Angus regions identified in the Cattle QTL database were body condition, body weight, calving ease, birth weight, fat thickness at the 12th rib, and milk traits. For example, the region on chromosome 23 (0-1090080bp) with high ROH frequency in Angus and Brangus harbored the KHDRBS2 gene, which was previously associated with calving ease (Cole et al., 2011).
One homozygous locus on chromosome 4 contained the LEP gene, which is associated with 96 traits in the Cattle QTL database. This gene is expressed in adipose tissue and codes for leptin, a hormone known to regulate feed intake and energy balance in mammals (Woronuk et al., 2012). This gene had been associated with marbling, fat thickness, rib eye area, and feed intake in several beef cattle breeds (Souza et al., 2010;Woronuk et al., 2012;Kononoff et al., 2017). Leptin is considered an extremely important gene for puberty onset (Williams et al., 2002). A high Angus contribution (90.4%) to this homozygous region was identified in Brangus (Table 1). Therefore, an allele coming from Angus was probably selected in Brangus.
Another homozygous region in chromosome 16 also was associated with first service conception in yearling Brangus heifers (Peters et al., 2013). Bos indicus-influenced heifers are known to have challenges achieving puberty early in life (Sartori et al., 2010;Fortes et al., 2012b). Therefore, high selection pressure in Brangus for early puberty since breed formation probably existed.
Another homozygous region (BTA 14) was previously identified as a QTL for weaning weight in Brangus (Weng et al., 2016). Cánovas et al. (2014) reported two genes on BTA14 at 24Mb associated with Brangus heifer fertility traits. This region harbors PLAG1 and XKR4 genes. The XKR4 was associated with subcutaneous rump fat thickness, scrotal circumference, serum concentration of prolactin, and sexual precocity (Fortes et al., 2012a;Porto Neto et al., 2012;Bastin et al., 2014;Takada et al., 2018). PLAG1 has been implicated in the regulation of stature and weight (Littlejohn et al., 2012;Pryce et al., 2012;Song et al., 2016). This gene was associated with yearling weight in Australian Tropical Composite breeds (Porto-Neto et al., 2014). The association studies of these genes used both taurine and indicine cattle, which confirms our observation of a selection signature in both founder breeds and a mixed origin of this region in Brangus.
The C allele of a putative functional mutation (rs109231213) near PLAG1 significantly increased hip height, weight, net food intake, age at puberty in males and females and decreased concentration of IGF-I in blood and fat depth (Fortes et al., 2013). These authors reported that haplotypes carrying the C allele had the same surrounding 10 SNP haplotype in B. taurus and Brahman, probably because the C allele was introgressed into Brahman from B. taurus cattle. The region with reduced heterozygosity surrounding the C allele was small in B. taurus and in Angus in this study (1.7 Mb) but 21.6 Mb long in Brahmans, here as well as in Fortes et al. (2013). Therefore, this allele represents a mutation that has been selected almost to fixation in B. taurus and, likely, introduced into Brahman cattle during crossbreeding with taurine cattle when indicine cattle were introduced into the United States (Sanders, 1980;Fortes et al., 2013).
Selection for growth and growth-related traits, such as average daily gain, feed conversion, and body size, has been conducted to improve beef productivity in both taurine and indicine breeds in the United States for several decades (Willham, 1982). Therefore, it is likely that favorable alleles for growth in genes with large phenotypic effects have also increased in frequency in both and the distribution of allele frequencies at these QTL have become similar between both populations.
The high Angus contribution for the selected genomic regions in Brangus cattle could support the use of the Brangus data for genomic selection and QTL identification (fine mapping) for Angus. This reinforces previous simulation studies that a crossbred or an admixed population can be used as training data for genomic selection and can provide reasonably accurate estimates of genomic breeding values of purebred selection candidates (Toosi et al., 2010). Marker estimates obtained from crossbred populations can be used to select purebreds looking for crossbred performance (Ibanez-Escriche et al., 2009;Toosi et al., 2010;MacNeil et al., 2011;Zeng et al., 2013;Lopes et al., 2017). Moreover, the results highlight how selection criteria can shape the genetic makeup of the composite.
The genetic composition of a composite breed is dynamic and changes across generations (Paim et al., 2020). Here, the selected regions in Brangus were mainly from Angus. The core idea of developing a composite breed is to exploit heterosis and complementarity between the breeds and, in the Brangus example, explore combining the tropical adaptation of zebu cattle and high yield and meat quality of Angus. These results and those previously reported (Paim et al., 2020) suggest Brangus is moving toward traits where Angus excel due to the selection imposed by breeders. Yield and meat quality (marbling) are measured and genetic values are available in the association's breed improvement program. The "tropical adaptation" traits, however, are not measured, and consequently, there is no genetic evaluation for their improvement. Therefore, it is important to develop and apply methods of measuring tropical adaptation and selecting for it; otherwise, this beneficial attribute of Brangus could be lost in future Brangus generations.

CONCLUSION
The majority of selection signatures in Brangus cattle came from Angus, which can be related to the traits of interest for genetic improvement and selection. These results demonstrate how quickly selection and drift can shift the genetic architecture of a population. Genomic inbreeding was found to be increasing in the composite population with advancing generations. Therefore, breeders should be aware of the need to manage inbreeding in this population. Moreover, composite cattle breeders need to be aware that selection for a set of specific traits that favor one of the progenitor breeds over the other can and will alter the original breed proportions and which, over the long term, may decrease the utility of the composite.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the all genotypic data used for this study are available in the website of The Animal-Genetic Resources Information Network (Animal-GRIN) (https://nrrc.ars.usda.gov/A-GRIN).

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because no samples were collected for this study; rather they were collected as part of other studies or program activities not associated with this study.

AUTHOR CONTRIBUTIONS
TP analyzed the data, interpreted the results, and wrote the manuscript. EH analyzed the data, discussed the results, and revised the text. CW maintained the datasets and analyzed the data. MT discussed the results and revised the text. LK was responsible for data curation, discussion of the results, and revision of the text. SP designed the study, discussed results, and revised the text. CM designed the study and revised the text. HB coordinated the study, discussed and interpreted the results, and wrote the manuscript. All authors read and approved the final manuscript.

FUNDING
Coordination for the Improvement of Higher Education Personnel (CAPES) provided the Ph.D. scholarship for the first author. USDA-ARS (National Animal Germplasm Program and the Meat Animal Research Center) provided genotypes and computing infrastructure.