Original Research published: 20 November 2019 doi: 10.3389/fgene.2019.01170 Mapping Recombination Rate on the Autosomal Chromosomes Based on the Persistency of Linkage Disequilibrium Phase Among Autochthonous Beef Cattle Edited by: Populations in Spain Francesca Bertolini, Technical University of Denmark, Elena Flavia Mouresan 1, Aldemar González-Rodríguez 1, Jhon Jacobo Cañas-Álvarez 2, Denmark Sebastián Munilla 1,3, Juan Altarriba 1, Clara Díaz 4, Jesús A. Baró 5, Antonio Molina 6, Reviewed by: Pascual Lopez-Buesa 1, Jesús Piedrafita 2 and Luis Varona 1* Alessandro Bagnato, University of Milan, Italy 1Departamento de Anatomía, Embriología y Genética Animal, Universidad de Zaragoza, Zaragoza, Spain, 2 Departament Diercles Francisco Cardoso, de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, Barcelona, Spain, 3 Departamento de Producción São Paulo State University, Brazil Animal, Facultad de Agronomía, Universidad de Buenos Aires, CONICET, Buenos Aires, Argentina, 4 Departamento de Andressa Oliveira De Lima, Mejora Genética Animal, Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA), Madrid, Spain, 5 Instituto Federal University of São Carlos, Agroalimentario de Aragón (IA2), Zaragoza, Spain, 6 Departamento de Ciencias Agroforestales, Universidad de Valladolid, Brazil Valladolid, Spain, 7 Departamento de Genética, Universidad de Córdoba, Córdoba, Spain, 8 Departamento de Producción Animal Christos Palaiokostas, y Ciencia de los Alimentos, Universidad de Zaragoza, Zaragoza, Spain Swedish University of Agricultural Sciences, Sweden In organisms with sexual reproduction, genetic diversity, and genome evolution are *Correspondence: Luis Varona governed by meiotic recombination caused by crossing-over, which is known to lvarona@unizar.es vary within the genome. In this study, we propose a simple method to estimate the recombination rate that makes use of the persistency of linkage disequilibrium (LD) phase Specialty section: This article was submitted to among closely related populations. The biological material comprised 171 triplets (sire/ Livestock Genomics, dam/offspring) from seven populations of autochthonous beef cattle in Spain (Asturiana a section of the journal Frontiers in Genetics de los Valles, Avileña-Negra Ibérica, Bruna dels Pirineus, Morucha, Pirenaica, Retinta, Received: 16 November 2018 and Rubia Gallega), which were genotyped for 777,962 SNPs with the BovineHD Accepted: 23 October 2019 BeadChip. After standard quality filtering, we reconstructed the haplotype phases in Published: 20 November 2019 the parental individuals and calculated the LD by the correlation -r- between each pair Citation: of markers that had a genetic distance < 1 Mb. Subsequently, these correlations were Mouresan EF, González-Rodríguez A, Cañas-Álvarez JJ, Munilla S, used to calculate the persistency of LD phase between each pair of populations along Altarriba J, Díaz C, Baró JA, the autosomal genome. Therefore, the distribution of the recombination rate along the Molina A, Lopez-Buesa P, Piedrafita J and Varona L (2019) genome can be inferred since the effect of the number of generations of divergence Mapping Recombination Rate on should be equivalent throughout the genome. In our study, the recombination rate was the Autosomal Chromosomes highest in the largest chromosomes and at the distal portion of the chromosomes. Based on the Persistency of Linkage Disequilibrium Phase In addition, the persistency of LD phase was highly heterogeneous throughout the Among Autochthonous Beef Cattle genome, with a ratio of 25.4 times between the estimates of the recombination rates Populations in Spain. Front. Genet. 10:1170. from the genomic regions that had the highest (BTA18-7.1 Mb) and the lowest (BTA12- doi: 10.3389/fgene.2019.01170 42.4 Mb) estimates. Finally, an overrepresentation enrichment analysis (ORA) showed Frontiers in Genetics | www.frontiersin.org 1 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase differences in the enriched gene ontology (GO) terms between the genes located in the genomic regions with estimates of the recombination rate over (or below) the 95th (or 5th) percentile throughout the autosomal genome. Keywords: recombination rate, linkage disequilibrium, beef cattle, multiple populations, gene ontology INTRODUCTION families or sperm cells. Second, other methods are based on the identification of local patterns in linkage disequilibrium Recombination caused by crossing-over during meiosis play a (LD) with coalescent methods (McVean et al., 2002; Li and crucial role in the genetic diversity and the genome evolution Stephens, 2003; Wall and Stevison, 2016), which estimate of organisms with sexual reproduction (Arnheim et al., 2007). It the background recombination rate ρw=4Nwcw −, where Nw creates new genetic variation by generating novel combinations and cw are the indistinguishable effective population size and of grand-paternal and grand-maternal genetic information, and recombination rate for a specific window of the genome, it helps to remove deleterious mutations that might otherwise respectively. The main limitation of the last approach is that the accumulate (Tiemann-Boege et al., 2017). effective population size can vary dramatically over time. In In most studies of genome-wide association or genomic fact, the decay of LD has been used to estimate past population selection, the distribution of crossing-over events had been history in humans (Hayes et al., 2003; Tenesa et al., 2007; Park, considered uniform, although there is strong evidence that 2012) and livestock populations (Hayes et al., 2003; De Roos recombination rate is heterogeneous along the genome (Myers et al., 2008; Xu et al., 2019). et al., 2005; Stapley et al., 2017). In general, recombination is In addition, stratification of the population can severely higher in the regions of the telomeres and smaller near the distort the estimates of recombination rates because centromere (Coop and Przeworski, 2007; Ma et al., 2015). Due subdivisions of the population have a strong effect on LD to recombination, the genome is organized into haplotype estimates (Hinrichs et al., 2009). After reproductive isolation, blocks of varying lengths, as described in humans (Gabriel et al., the structure of LD tends to differ between subpopulations 2002) and other species as rat and mouse (Guryev et al., 2006) and the similarity (or persistency) of those LD patterns and cattle (Mokry et al., 2014). The reason of this structure is depends on the number of generations of divergence and the the presence of small genomic regions that have a higher rate of recombination rate plus other evolutionary events such as recombination, known as recombination hotspots (Paigen and admixture or variations on the effective size of the populations Petkov, 2010). (Hill and Robertson, 1968). For this reason, if genotypic In addition, patterns of the recombination rate throughout information is available for closely related populations, the genome vary among species, populations, or even within measures of genome-wide persistency of LD phase among individuals in different environments (Stapley et al., 2017). The populations throughout the genome can be used to infer evolution of the distribution of the recombination rate along the the distribution of recombination rate. The rationale of this genome is an active research field (Dapper and Payseur, 2017). approach is that genetic drift, admixture, or variations of the In general, it differs according to the genomic scale in which the effective size should affect the entire genome with similar recombination rate is measured (Smukowski and Noor, 2011). In intensity and the heterogeneity of the persistency of LD phase a very fine scale (few kb), a rapid divergence of the recombination is linked to variations on the recombination rate. rate between mammal populations is observed (Auton et al., Therefore, the objective of this study was to develop a 2012; Stevison et al., 2016), whereas greater correlations are procedure to infer the distribution of the recombination observed between closely related populations when they are rate from the persistency of LD phase among closely related calculated through larger chromosomal segments (Smukowski populations and to apply it to genotypic data from seven beef and Noor, 2011; Shen et al., 2018). cattle populations in Spain. Traditionally, the distribution of the crossing-overs or recombination events within the genome has been studied by counting the number of chiasmata during meiosis (Hulten MATERIALS AND METHODS et al., 1982) or from linkage maps created from a limited number of genetic markers or phenotypes (Sturtevant, 1913). The genomic data comprised the BovineHD Genotyping Beadchip In recent years, high-throughput sequencing and genotyping (777,962 SNPs, Illumina) genotypes from 171 non-related triplets technologies have provided a valuable new tool for measuring of sire, dam, and one offspring from seven breeds, being 25 recombination rates with two main group of methods. First, Asturiana de los Valles (AV), 24 Avileña - Negra Ibérica (ANI), estimates of recombination rates are based on observations 25 Bruna dels Pirineus (BP), 25 Morucha (Mo), 24 Pirenaica (Pi), of recombination events in large pedigrees between pairs of 24 Retinta (Re), and 24 Rubia Gallega (RG) triplets. This dataset parent-offspring genotypes (Kong et al., 2010; Ma et al., 2015; has been used to analyze genetic differentiation (Cañas-Álvarez Shen et al., 2018) or in sperm typing (Sarbajna et al., 2012) et al., 2015; Cañas-Álvarez et al., 2016; González-Rodríguez et al., and require genotypic information from a large number of 2017), signatures of selection (González-Rodríguez et al., 2016), Frontiers in Genetics | www.frontiersin.org 2 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase and haplotype diversity (Mouresan et al., 2017). These breeds in each bin of 20 kb (0–20, 20–40, 40–60, 960–980, 980– represent 72% of the total census of local beef breeds in Spain 1,000) within a 1 Mb window were calculated for each pair of (Ministerio de Medio Ambiente y Medio Rural y Marino, 2010) populations. We obtained 50 correlation estimates (one per bin) and their production systems are extensive or semi-extensive. The between the LD estimates of each breed pair per window. populations are reared in mountainous regions near the Pyrenees Under the assumption of constant variance of r (or effective (Pirenaica and Bruna dels Pirineus) in the humid regions in population size) in both populations, the expectation of the northwestern Spain (Rubia Gallega and Asturiana de los Valles) correlation between LD estimates from a pair of SNP markers or in pastures in semi-arid zones of the west and southwest of is e-2cT (Hill and Robertson, 1968; De Roos et al., 2008), where Spain (Retinta, Avileña Negra-Ibérica, and Morucha). The breeds c is the recombination rate between the markers and T is the differ in production and carcass traits (Piedrafita et al., 2003) and number of generations of divergence between populations. in their meat quality (Gil et al., 2001). Initially, it was assumed that the recombination rate was 1.25 The triplets were sampled by the breeders associations with cM per Mb (Arias et al., 2009). The regression of the natural the aim of capturing most of the genetic variability of each logarithm of the correlations on the genomic distance was population. We used an ad-hoc procedure that started with one calculated and the slope was equated to -2cT to estimate triplet and incorporated the new ones by minimizing the total T between each pair of populations. coancestry between them. The SNP filtering process included Once the numbers of generations of divergence (T) were the following: 1) Mendelian error < 0.05, 2) SNP and individual estimated from all available SNP markers, they were assumed as call rates > 95%, and 3) Minor allele frequency (MAF) > 0.05 known and replaced by their estimates. Subsequently, the same in pairs of populations. Only the SNPs that were located on expression (-2cT) was used to estimate c, although the analysis autosomal chromosomes were retained. The filtering process was restricted only to the SNP markers within 1 Mb, which were yielded approximately 550,000 segregating markers for each pair centered every 0.1 Mb along the autosomal genome in sliding of populations (see Table 1). windows. Therefore, 25,098 estimates of c were calculated for The genomic information of the triplets was used to reconstruct each of the 21 population pairs. the parental haplotypes with the TRIO option of the BEAGLE Afterwards, the presence of a common pattern for the software (Browning and Browning, 2007), which were used to distribution of the recombination rate was checked calculating calculate the LD in each population and between each pair of the correlation between the estimates from all population pairs. markers (i.e. with alleles A and a, and B and b, respectively) that Next, the estimates of c were averaged by chromosome, relative had a genomic distance < 1 Mb. LD was estimated as a correlation position within the chromosome, and location within the –r-, as follows: genome to infer the distribution of the recombination rate (or the persistency of LD phase) throughout the bovine autosomal D genome. The degree of inequality of the recombination rate along r= pA pa pB pb the autosomal genome was measured with the Gini index –G- (Ceriani and Verme, 2012) as: where D = PAB Pab PAbPaB (Falconer and Mackay, 1996), PAB , Pab, ng ng PAb and Pab were the haplotype frequencies, and PA, Pa, PB, and Pb were the allelic frequencies. G= ∑ ∑ i =1 j =1 ci − c j 2 To estimate the persistency of LD phase between pairs of 2ng c populations, the Pearson correlations between LD estimates Where ci and cj were the average estimates of the recombination rate for the 21 population pairs ad for the ith and jth genomic regions, ng was the total number of genomic regions (25,098) TABLE 1 | Number of co-segregating SNP markers between all possible pairs of seven beef cattle populations in Spain [Asturiana de los Valles (AV), Avileña defined along the autosomal genome, and c was the average of - Negra Ibérica (ANI), Bruna dels Pirineus (BP), Morucha (Mo), Pirenaica (P), the 25,098 estimates of the recombination rate. Retinta (Re), Rubia Gallega (RG)]. Finally, the 5th and 95th percentiles of the average recombination rates between pairs of populations in 0.1 Mb steps Pairs of N° SNP markers Pairs of N° SNP populations populations markers throughout the autosomal genome were calculated. Thereinafter, we prospected genes mapped within genomic regions falling AV-ANI 555,373 BP-Mo 543,305 out the 5–95th the percentile using the Biomart tool of AV-BP 557,588 BP-Pi 534,336 AV-Mo 555,769 BP-Re 535,997 Ensembl (Flicek et al., 2013) (www.ensembl.org). Further, we AV-Pi 540,390 BP-RG 544,350 performed an Overrepresentation Enrichment Analysis (ORA) AV-Re 547,893 Mo-Pi 529,281 to determine if the overrepresentation of gene ontology (GO) AV-RG 553,868 Mo-Re 541,225 terms differed between the two tails of the empirical distribution ANI-BP 538,327 Mo-RG 542,682 of recombination rates. We used the WEB-based Gene SeT ANI-Mo 545,324 Pi-Re 522,670 ANI-Pi 524,630 Pi-RG 529,577 AnaLysis Toolkit (www.webgestalt.org) using the Homo sapiens ANI-Re 536,595 Re-RG 535,677 and Bos taurus annotation databases and with the complete ANI-RG 537,882 genome as the reference set. Frontiers in Genetics | www.frontiersin.org 3 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase RESULTS AND DISCUSSION between them (De Roos et al., 2008) and consistent with a decrease in effective population size in cattle (Hayes et al., 2003). We have used the persistency of LD phase among seven closely In a second step, we restricted the analysis to the LD within related populations (Cañas-Álvarez et al., 2015) to infer the 250 Kb and the slope of the regression analyses were equated landscape of the recombination events in a sequential procedure to -2cT, with c set to 1.25 cM per Mb (Arias et al., 2009). We that involved several steps. As in other studies in livestock obtained 21 estimates of the number of generations of divergence populations (Brito et al., 2015; Biegelmeyer et al., 2016; Brito et al., (T) between populations (Table 2), that ranged from 132.3 2017; Grossi et al., 2017), the similitude of LD among populations (AV-BP) to 281.9 (Pi-Re). Estimates were in concordance with was very high between adjacent markers and decreased rapidly the results obtained by Cañas-Álvarez et al. (2016) for the with genomic distance. In the first step, we used all the SNP same populations and dataset, and by De Roos et al. (2008) markers to estimate the persistency of LD phase for each pair of between two dairy populations (Holstein-Friesian and Jersey). populations as the slope of the regression analysis between the However, the divergence times found in our study were lower natural logarithm of the correlations between the r measures of than the observed between a dairy (Holstein) and beef (Angus) LD in bins of 20 kb on the genomic distance. Theoretically, this populations (De Roos et al., 2008). relationship should be linear (Hill and Robertson, 1968; De Roos The estimates of T assumed that the variance of the LD (r) et al., 2008). However, the results varied according to the genomic remained constant in each population and this is probably distance evaluated (1 Mb, 500 Kb, 250 Kb, or 100 Kb). To illustrate far from the truth. The effective size of the populations has this phenomenon, the results of the regression analysis between decreased in the last generations (Cañas-Alvarez et al., 2016) Re and RG are shown in Figure 1. The regression analysis of the and, thus, the estimates of T are probably overestimated. data within a range of 250 Kb had the highest adjusted R2 (0.999) Nevertheless, the bias caused by the variations in the effective and the linear relationship between persistency of LD phase and size should be similar throughout the autosomal genome and, genomic distance was evident only in those first 250 Kb. The therefore, local variations in the persistency of LD phase are results for the other population pairs were similar (Supplementary informative for the inference of the recombination landscape Information, Figures S1 to S20) and all of them had the highest along the autosomal genome. adjusted R2 (> 0.998) with the first 250 Kb. In all population pairs, After this preliminary step, we used the same expression persistency of LD phase decayed rapidly over short distances, but (-2cT) to infer the distribution of the recombination rate c. Now, in larger genomic distances remained > 0, as observed by De Roos the numbers of generations of divergence (T) were assumed to be et al. (2008). This is probably due to the fact that the populations known and they were replaced by their estimates with all the SNP were not totally divergent or due to the presence of some migration markers. In this case, we equated the slope of the regression of the FIGURE 1 | Persistency of Linkage Disequilibrium (LD) measured as the correlation and the log-correlation between the estimates of LD in Retinta (Re) and Rubia Gallega (RG) populations, and estimates of the slope of the regression with respect to the genomic distance for different ranges [0-1 Mb – black, 0-500 Kb-green, 0-250 Kb- red, and 0-100 Kb-blue]. Frontiers in Genetics | www.frontiersin.org 4 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase TABLE 2 | Estimated number of generations of divergence between seven beef cattle populations in Spain [Asturiana de los Valles (AV), Avileña - Negra Ibérica (ANI), Bruna dels Pirineus (BP), Morucha (Mo), Pirenaica (P), Retinta (Re), Rubia Gallega (RG)] based on the architecture of Linkage Disequilibrium. ANI BP Mo Pi Re RG AV 181.2 132.3 160.7 184.1 185.9 157.2 ANI – 244.9 133.6 268.6 175.1 225.1 BP – – 232.8 168.1 258.1 176.4 Mo – – – 252.8 168.8 205.2 Pi – – – – 281.9 215.3 Re – – – – – 229.1 persistency of LD phase on the genomic distance to -2cT, but the could also reflect differences in the intensity of the mutation rate analysis was restricted to the SNP markers located within a sliding or the occurrence of selection events. window of 1 Mb (500 kb downstream and 500 kb upstream) in Once the local estimates of c for the 21 pair of populations steps of 0.1 Mb. We obtained as many as 25,098 estimates of c for were available, we calculated the correlations between them each pair of populations (Figures S21 to S41). The rationale to (Figure 2). Values ranged from moderate (0.37) to strong (0.76), this approach was that genetic drift and variations on the effective with an average of 0.55 and a standard deviation of 0.08. The size of the populations should have affected the entire autosomal results were consistent with the output of Table 2. The minimum genome with the same intensity and, therefore, regional variation correlation (0.37) was obtained between ANI-Mo and BP-Pi, in the persistency of LD phase should reflect variations in the which is consistent with the large number of generations of recombination rate. However, as in other LD-based procedures divergence between ANI and BP (244.9), ANI and Pi (268.6), for estimating the recombination rate (Li and Stephens, 2003), it Mo and BP (232.8), and Mo and Pi (252.8). In contrast, the FIGURE 2 | Correlations between the estimates of the recombination rate (c) in bins of 1 Mb between pairs of seven beef cattle populations in Spain [Asturiana de los Valles (AV), Avileña - Negra Ibérica (ANI), Bruna dels Pirineus (BP), Morucha (Mo), Pirenaica (P), Retinta (Re), Rubia Gallega (RG)]. Frontiers in Genetics | www.frontiersin.org 5 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase greatest correlation was between ANI-BP and Mo-BP, given that length (Lynch, 2006) has been reported in several species and the number of generations of divergence between ANI and Mo may be associated to the difficulties of small chromosomes to find is only 133.6 generations. Therefore, given the noise in the LD their homologues during meiosis (Tiemann-Boege et al., 2017). estimates with small sample sizes, this average correlation and In fact, Fledel-Alon et al. (2009) suggested that, in meiosis, each the consistency of estimates between pairs of populations should chromosome usually undergoes at least one crossing-over, which be considered very relevant. The similarity between estimates produces a very strong correlation between the average number of was within the same range than the reported between estimates crossovers and chromosome length. of the recombination rate in human (Graffelman et al., 2007; Next, we evaluated the relationship between the recombination Laayouni et al., 2011; Manu et al., 2018) and livestock (Petit et al., rate and the relative physical position within the chromosome by 2017; Shen et al., 2018) populations. Thus, somehow it confirms averaging the c estimates for each percentile along the length of that similarities in the distribution of the recombination rate are each chromosome (Figure 4). The results were similar to those of achieved between closely related populations (Smukowski and Sandor et al. (2012); Ma et al. (2015), and Shen et al. (2018) in the Noor, 2011), such as the analyzed in this study (Cañas-Álvarez bovine specie and to the results of haplotype diversity measured et al., 2015; González-Rodríguez et al., 2017). This similarity in the same populations (Mouresan et al., 2017). In cattle, all was even observed between pairs of populations that do not autosomal chromosomes are acrocentric (Popescu, 1990) and, in share any population (i. e. ANI-AV and BP-Pi), whose average our study, the recombination rate was lowest at the beginning correlation (0.52 ± 0.06) was only slightly lower than the average of the chromosome, near the centromere. Furthermore, a of the correlations between estimates from pairs that share (i.e. low recombination rate was evident at the middle of the ANI-AV vs ANI-Mo) (0.58 ± 0.08). It reinforces the hypothesis chromosome, although the centromere of chromosomes in cattle that the similarity between persistence of LD phase at different is not located there. Ma et al. (2015) argued that the bimodal locations of the genome and between pairs of populations is distribution of recombination rates might be caused by positive related with variations in the recombination rate throughout crossover interference. The highest recombination rate was at the the genome. distal portion of the chromosome (over the 95 % percentile of Despite divergences associated with a specific population pair the relative position within chromosomes), in agreement with and possible selection events, the distribution of the persistency of studies that have shown that recombination rate is highest at the LD phase appeared to follow a global pattern. Therefore, we used telomeres (Nachman, 2002; Coop and Przeworski, 2007). the estimates of c to describe the distribution of the recombination Additionally, a map of the recombination rates throughout the rate along the autosomal genome. Initially, we calculated genome was calculated by averaging the estimates from the 21 pairs an average of all of the estimates of c for each chromosome of populations in 0.1 Mb steps along the autosomal chromosomes (Figure 3), which ranged from 1.12 cM per Mb in BTA9 to 1.50 (Figure 5). The average estimate of the recombination rate cM per Mb in BTA25. In general, the largest chromosomes tended was 1.275 cM per Mb with a standard deviation of 0.381. As to have the lowest recombination rates. The relationship between expected, the recombination rate was very similar to the rate recombination rate and chromosome length (Kaback et al., 1992; assumed in the initial step of the study (1.25 cM per Mb). The Jensen-Seaman et al., 2004; Li and Freudenberg, 2009) or genome estimated recombination rates within the genome were highly FIGURE 3 | Average estimate per chromosome of the recombination rate (cM/Mb) in bins of 1 Mb and for all pairs of seven breeds of beef cattle populations in Spain. Frontiers in Genetics | www.frontiersin.org 6 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase FIGURE 4 | Average estimate of the recombination rate (cM/Mb) in bins of 1 Mb and for all pairs of populations and the relative physical position in the chromosome in seven breeds of beef cattle in Spain. FIGURE 5 | Average estimate of the recombination rate (cM/Mb) in bins of 1 Mb and throughout the autosomal genome for all pairs of seven beef cattle populations in Spain. heterogeneous, and the ratio between the genomic regions that Figure 6 displays the average recombination rate and all specific had the highest (BTA18-7.1 Mb) and the lowest (BTA12-42.4 persistencies of LD phase for each population pair. Mb) estimated recombination rates was 25.37. To illustrate the The heterogeneity of the recombination rate reflected the differences between those two extreme regions of the genome, presence of highly recombining genomic regions. Most of Frontiers in Genetics | www.frontiersin.org 7 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase FIGURE 6 | Average (bold lines) and specific population pair (transparent lines) persistency of LD (correlation between LD estimates in pairs of populations) in the genomic regions that had the highest (BTA12-42.4 Mb, continuous lines) or the lowest persistency (BTA18-7,1 Mb, dotted lines) among seven beef breeds in Spain. the recombination events may occur in a small portion of the suggested that there is a negative correlation between gene genome, as observed in other species. However, the Gini index density and the frequency of recombination hotspots (Myers between the cumulative distributions of the recombination rate et al., 2005; Freudenberg et al., 2009; Stapley et al., 2017), which and the genetic distance was 0.1803, which is lower than others was not detectable with the methods used in our study. reported in human (Kong et al., 2002), apes (Stevison et al., 2016), Furthermore, some authors (The International HapMap and livestock (Petit et al., 2017) populations. It is likely that the Consortium, 2007) have suggested that genes with highly Gini index was low because the method used in our study was conserved function are located surrounding regions with low only able to distinguish among the rates of recombination (as a recombination rate; on the other hand, HRR regions contains measure of the persistency in LD phase) within relatively large genes that are exposed to recurrent adaptive process to allow genomic regions (1 Mb), and recombination hotspots often are plasticity of organism to coming circumstances. In our study, restricted to 1–2 kb (Myers et al., 2005; Mancera et al., 2008). we have tried to corroborate these statements using ORA with The 5th and 95th percentiles of the recombination rate estimates the GO terms for biological processes, cellular components, and were 0.593 and 1.856 cM per Mb, respectively, and we identified molecular functions. The results are presented in Tables 3 to 5 the genes within the genomic regions that had values that were (Homo sapiens annotation database) and in Supplementary either above or below those percentiles. The number of genes Tables 1 to 3 (Bos taurus annotation database). In general, the within the regions that were above the 95th (high recombination results with the Homo sapiens database yielded results with lower rate – HRR-) or below the 5th (low recombination rate –LRR-) FDR than with the Bos taurus database, probably because the percentiles were 665 and 669, respectively. Some studies have human genome in notably more annotated than the bovine one. TABLE 3 | False discovery rate (FDR) for the top 10 enriched Gene Ontology (GO) terms for Biological processes with the Homo sapiens database for genes within the genomic regions located over the 95th (high recombination rate) and below the 5th (low recombination rate) percentiles of the average recombination rates. High Recombination Rate Low Recombination Rate GO TERM FDR GO TERM FDR Protein citrullination 2.00e-04 Homophilic cell adhesion via plasma membrane adhesion molecules 2.45e-01 Histone citrullination 2.00e-04 DNA replication initiation 2.45e-01 Extracellular vesicle biogenesis 3.09e-01 Vitamin A metabolic process 2.45e-01 Regulation of action potential 3.27-01 Cell-cell adhesion via plasma-membrane adhesion molecules 5.84e-01 Regulation of substrate adhesion- 3.27-01 Regulation of chemokine biosynthetic process 5.84e-01 dependent cell spreading Sodium ion transmembrane transport 3.27-01 Celullar response to electrical stimulus 5.84e-01 Peptidyl-arginine modification 3.27-01 Chemokine biosynthetic process 5.84e-01 phagosome acidification 3.27-01 Chemokine metabolic process 5.84e-01 Cardiac muscle cell action potential 3.27-01 Neutrophil mediated killing of bacterium 1.00e-00 Regulation of fibroblast growth factor 3.27-01 Sulfur amino acid catabolic process 1.00e-00 receptor signaling pathway Frontiers in Genetics | www.frontiersin.org 8 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase TABLE 4 | False discovery rate (FDR) for the top 10 enriched Gene Ontology (GO) terms for Cellular Components with the Homo sapiens database for genes within the genomic regions located over the 95th (high recombination rate) and below the 5th (low recombination rate) percentiles of average recombination rates. High Recombination rate Low Recombination Rate GO TERM FDR GO TERM FDR Neuron part 0.02e-04 Plasma membrane region 2.00e-04 Synapse 0.02e-04 Golgi apparatus 2.70e-03 Synapse part 1.5e-03 Intrinsic component of the plasma membrane 1.25-02 Golgi apparatus 1.5e-03 Cytoplasmic vesicle part 2.83e-02 Cell projection part 1.5e-03 Integral component of the plasma membrane 2.83e-02 Plasma membrane bounded cell projection part 1.5e-03 Perinuclear region of cytoplasm 2.83e-02 Postsynapse 3.0e-03 Cell projection part 3.05e-02 Neuron projection 3.2e-03 Plasma membrane bounded cell projection part 3.05e-02 Vacuole 3.7e-03 Golgi subcompartment 3.05e-02 Phagocytic vesicle 4.8e-03 Plasma membrane protein complex 4.6e-02 TABLE 5 | False discovery rate (FDR) for the top 10 enriched Gene Ontology (GO) terms for Molecular Functions for genes with the Homo sapiens database within the genomic regions located over the 95th (high recombination rate) and below the 5th (low recombination rate) percentiles of the average recombination rates. High Recombination Rate High Recombination Rate GO TERM FDR GO TERM FDR Protein-arginine deaminase activity 1.00e-04 GTP-dependent protein binding 8.01-0e2 Hydrolase activity, acting on carbon nitrogen (but not 1.13e-02 Aminoacyl-tRNA ligase activity 4.76e-01 peptide) bonds, in linear amidines Solute: cation antiporter activity 4.74e-02 Ligase activity, forming carbon-oxigen bonds 4.76e-01 Cation: cation antiporter activity 7.49e-02 Molecular carrier activity 4.76e-01 Potassium: proton antiporter activity 1.01e-01 Drug binding 4.76e-01 Metal ion transmembrane transported activity 1.11e-01 Phosphatidylinositol bisphosphate kinase activity 4.76e-01 Solute: Proton antiporter activity 1.16e-01 ARF guanyl-nucleotide exchange factor activity 4.76e-01 Monovalent cation: proton antiporter activity 1.16e-01 Nucleocytoplasmic carrier activity 4.76e-01 Sodium: proton antiporter activity 1.16e-01 Identical protein binding 4.76e-01 Monovalent inorganic cation transmembrane 1.19e-01 Metalloexopeptidase activity 4.76e-01 transported activity The results of the enrichment analysis for biological processes regions were enriched with genes whose expression is located at the with the human database (Table 3) only provided enriched GO extracellular space and related with neuronal interactions (neuron terms with a False Discovery Rate (FDR) lower than 0.05 with part, synapse, synapse part, postsynapse, or neuron projection). In the genes present in the HRR genomic regions. The significant contrast, the genes located at LRR regions were associated with very GO terms correspond to Protein citrullination and Histone basic intracellular (Cytoplasmic vesicle part, Perinuclear region of the citrullination. Citrullination, the conversion of the amino acid cytoplasm) or membrane components (Plasma membrane region, arginine in a protein into the amino acid citrulline, has been related Intrinsic component of the plasma membrane Integral component to an increase in antigenic diversity (Nguyen and James, 2016). of the plasma membrane, Plasma membrane protein complex). The The higher recombination rate of those genomic regions might results obtained from the Bos taurus database were significant indicate that they have evolved to have high plasticity to adapt (FDR < 0.05) only for cytosol and nuclear lumen in the LRR regions, to changing environments (Charlesworth et al., 2009; Campos confirming the results provided by the Homo sapiens database. et al., 2014). Therefore, the generation of new genetic variants by The enriched GO terms for molecular functions are presented recombination may help the antigen diversity from the perspective in Table 5 (Homo sapiens) and Supplementary Table 3 (Bos of the host. Thus, it works as a mechanism to adapt its immune Taurus). The three significantly (FDR < 0.05) enriched GO terms response to fight against the ability of the pathogen to modify its for the HRR genomic regions were coherent with the enriched antigenic targets. The results obtained from the Bos taurus database GO terms for biological processes. In fact, two of those were (Supplementary Table 1) were not significant (FDR < 0.05). clearly associated with citrullination [Protein-arginine deiminase The top 10 enriched GO terms for cellular components for the activity and Hydrolase activity, acting on carbon nitrogen (but not genomic regions that had either high or low recombination rate are peptide) bonds, in linear amidines] and the other (Solute: cation presented in Table 4 (Homo sapiens) and Supplementary Table 2 (Bos antiporter activity) was linked to transmembrane transportation of taurus). For Homo sapiens, the FDR was generally lower than it was solutes. In contrast, the only significantly enriched GO terms for for biological processes (FDR < 4.6e-02). Some cellular components LRR genomic regions was GTP-dependent protein binding (FDR (Golgi apparatus, Cell projection part, Plasma membrane bounded bell = 8.01e-02), which was confirmed with the results from the Bos projection part) occurred in both types of genomic regions, but there taurus database (FDR = 2.7e-02). The genes located in the LRR were some important differences between them. The HRR genomic genomic regions should be more conserved, since they may be Frontiers in Genetics | www.frontiersin.org 9 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase necessary for basic functions of the organism. In this sense, the it requires populations to be close enough that the recombination genes belonging to the GTP-dependent protein GO term regulate rate is well conserved. guanine nucleotide-binding proteins that play a crucial role in signal transduction and in a large number of cellular processes (Zachariou et al., 2012). AUTHOR CONTRIBUTIONS The results of this study confirm that the genomic architecture The study was conceived by LV. Preliminary data preparation of persistency of LD phase is well conserved among closely was done by EM, AG-R, JC-Á, and SL. Data analysis was done by related populations, such as the Spanish autochthonous beef cattle EM and LV. EM and LV produced the draft manuscript. JA, CD, breeds, and is heterogeneous within the autosomal genome and AM, and JP collaborated in generating the data. SL, JA, CD, JB, that this heterogeneity can be used to estimate the recombination AM, PL-B, and JP discussed the results. All authors reviewed and rate. Several studies have estimated the persistency of LD phase revised the manuscript. between populations as a measure of genetic diversity (De Roos et al., 2008; Villa-Angulo et al., 2009; Cañas-Álvarez et al., 2016) and as a mean of predicting the marker density required for multi- FUNDING breed genomic evaluation (De Roos et al., 2008; Cañas-Álvarez We thank the AGL 2010-15903 grant from the Spanish government et al., 2016). Nevertheless, the genetic architecture of persistency and the KBBE.2011.1.3-06 project from the European Union’s of phase within the genome has received limited attention. In this Seventh Framework. We also thank the Breed Societies for their study, we estimated the persistency of LD phase among seven collaboration in collecting samples. The support of FEAGAS is beef cattle populations in Spain and its distribution within the acknowledged. JC-Á acknowledges the COLCIENCIAS support genome, which has been related with the genetic architecture of by the Francisco José De Caldas Fellowship 497/2009, and AG-R the recombination rates. Even though the recombination rate acknowledges the financial support provided by the BES-2011- varies among species, sex, and populations (Dapper and Payseur, 045434 fellowship. 2017), some general patterns were described (Tiemann-Boege et al., 2017). The patterns were confirmed in our analysis by the similitude of our results with some studies in other species or ACKNOWLEDGMENTS populations (Ma et al., 2015; Shen et al., 2018). Therefore, the main conclusion of this study is that the We thank Bruce Macwirther and María Martínez-Castillero for heterogeneity of persistency of LD phase between closely their help in the language editing. Preliminary results of this related populations can be used to estimate the recombination study were presented at the World Congress of Genetics Applied rate with the procedure developed here, which is simpler and to Livestock Production, 11.285. brings similar results as more complex and coalescent dependent methods. This implies that it may help to identify regions related SUPPLEMENTARY MATERIAL to hot and cold spots when data from several populations of the same ancestral origin are available. Nevertheless, our procedure The Supplementary Material for this article can be found online at: has several limitations since differences in the mutation rate or https://www.frontiersin.org/articles/10.3389/fgene.2019.01170/ selection events can locally affect the persistency of LD phase and full#supplementary-material REFERENCES and composite breeds’ history revealed by a high-density SNP chip. BMC Genet. 18 (1), 25. doi: 10.1186/s12863-017-0492-8 Arias, J. A., Keehan, M., Fisher, P., Coppieters, W., and Spelman, R. (2009). Browning, S. R., and Browning, B. L. (2007). Rapid and accurate haplotype A high density linkage map of the bovine genome. BMC Genet. 10, 18. doi: phasing and missing-data inference for whole-genome association studies by 10.1186/1471-2156-10-18 use of localized haplotype clustering. Am. J. Hum. Genet. 81, 1084–1097. doi: Arnheim, N., Calabrese, P., and Tiemann-Boege, I. (2007). Mammalian meiotic 10.1086/521987 recombination hot spots. Annu. Rev. Genet. 41, 369–399. doi: 10.1146/annurev. Cañas-Álvarez, J. J., González-Rodríguez, A., Munilla, S., Varona, L., Díaz, C., genet.41.110306.130301 Baro, J. A., et al. (2015). Genetic diversity and divergence among spanish beef Auton, A., Fledel-Alon, A., Pfeifer, S., Venn, O., Ségurel, L., Street, T., et al. (2012). cattle breeds assessed by a bovine high-density SNP chip. J. Anim. Sci. 93 (11), A fine-scale chimpanzee genetic map from population sequencing. Science 13, 5164–5174. doi: 10.2527/jas2015-9271 193–198. doi: 10.1126/science.1216872 Cañas-Álvarez, J. J., Mouresan, E. F., Varona, L., Díaz, C., Molina, A., Baro, J. A., et al. Biegelmeyer, P., Gulias-Gomes, C. C., Caetano, A. R., Steibel, J. P., and Cardoso, F. F. (2016). Linkage disequilibrium, persistence of phase, and effective population size (2016). Linkage disequilibrium, persistence of phase and effective population in spanish local beef cattle breeds assessed through a high-density single nucleotide size estimates in hereford and braford cattle. BMC Genet. 17 (1), 32. doi: polymorphism chip. J. Anim. Sci. 94 (7), 2779–2788. doi: 10.2527/jas2016-0425 10.1186/s12863-016-0339-8 Campos, J. L., Halligan, D. L., Haddrill, P. R., and Charlesworth, B. (2014). The Brito, L. F., Jafarikia, M., . Grossi, D. A., Kijas, J. W., . Porto-Neto, L. R., Ventura, R. V., relation between recombination rate and patterns of molecular evolution and et al. (2015). Characterization of linkage disequilibrium, consistency of gametic variation in drosophila melanogaster. Mol. Biol. Evol. 31 (4), 1010–1028. doi: phase and admixture in australian and canadian goats. BMC Genet. 16 (1), 67. 10.1093/molbev/msu056 doi: 10.1186/s12863-015-0220-1 Ceriani, L., and Verme, P. (2012). The origins of the gini index: extracts from Brito, L. F., McEwan, J. C., Miller, S. P., Pickering, N. K., Bain, W. E., Dodds, K. G., variabilità e mutabilità (1912) by corrado gini. J. Economic Inequality 10 (3), et al. (2017). Genetic diversity of a new zealand multi-breed sheep population 421–443. doi: 10.1007/s10888-011-9188-x Frontiers in Genetics | www.frontiersin.org 10 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase Charlesworth, B., Betancourt, A. J., Kaiser, V. B., and Gordo, I. (2009). Genetic Kong, A., Gudbjartsson, D. F., Sainz, J., Jonsdottir, G. M., Gudjonsson, S. A., recombination and molecular evolution. Cold Spring Harb. Symp. Quant. Biol. Richardsson, B., et al. (2002). A high-resolution recombination map of the 74, 177–186. doi: 10.1101/sqb.2009.74.015 human genome. Nat. Genet. 31 (3), 241–247. doi: 10.1038/ng917 Coop, G., and Przeworski, M. (2007). An evolutionary view of human Kong, A., Thorleifsson, G., Gudbjartsson, D. F., Masson, G., Sigurdsson, A., recombination. Nat. Rev. Genet. 8 (1), 23–34. doi: 10.1038/nrg1947 Jonasdottir, A., et al. (2010). Fine-scale recombination rate differences between Dapper, A. L., and Payseur, B. A. (2017). Connecting theory and data to understand sexes, populations and individuals. Nature 467 (7319), 1099–1103. doi: recombination rate evolution. Philos. Trans. R. Soc. B: Biol. Sci. 372 (1736), 10.1038/nature09525 20160469. doi: 10.1098/rstb.2016.0469 Laayouni, H., Montanucci, L., Sikora, M., Melé, M., Dall’Olio, G. M., Lorente- De Roos, A. P. W., Hayes, B. J., Spelman, R. J., and Goddard, M. E. (2008). Linkage Galdos, B., et al. (2011). Similarity in recombination rate estimates highly disequilibrium and persistence of phase in holstein-friesian, jersey and angus correlates with genetic differentiation in humans. Edited by Carles Lalueza- cattle. Genetics 179, 1503–1512. doi: 10.1534/genetics.107.084301 Fox. PloS One 6 (3), e17913. doi: 10.1371/journal.pone.0017913. Falconer, D. S., and Mackay, T. F. C. (1996). Introduction to Quantitative Genetics. Li, W., and Freudenberg, J. (2009). Two-parameter characterization of Harlow, UK: Longman Group. chromosome-scale recombination rate. Genome Res. 19 (12), 2300–2307. doi: Fledel-Alon, A., Wilson, D. J., Broman, K., Wen, X., and Ober, C. (2009). Broad- 10.1101/gr.092676.109 scale recombination patterns underlying proper disjunction in humans. PloS Li, N., and Stephens, M. (2003). Modeling linkage disequilibrium and identifying Genet. 5 (9), 1000658. doi: 10.1371/journal.pgen.1000658 recombination hotspots using single-nucleotide polymorphism data. Genetics Flicek, P., Ahmed, I., Ridwan Amode, M., Barrell, D., Beal, K., Brent, S., et al. 165 (4), 2213–2233. (2013). Ensembl 2013. Nucleic Acids Res. 41, D48–D55. doi: 10.1093/nar/ Lynch, M. (2006). The origins of eukaryotic gene structure. Mol. Biol. Evol. 23 (2), gks1236 450–468. doi: 10.1093/molbev/msj050 Freudenberg, J., Wang, M., Yang, Y., and Li, W. (2009). Partial correlation Ma, L., O’Connell, J. R., VanRaden, P. M., Shen, B., Padhi, A., Sun, C., et al. (2015). analysis indicates causal relationships between gc-content, exon density Cattle sex-specific recombination and genetic control from a large pedigree and recombination rate in the human genome. BMC Bioinf. 10, S66. doi: analysis. PloS Genet. 11 (11), e1005387. doi: 10.1371/journal.pgen.1005387. 10.1186/1471-2105-10-S1-S66 Edited by Adam J. Auton. Gabriel, S. B., Schaffner, S. F., Nguyen, H., Moore, J. M., Roy, J., Blumenstiel, B., Mancera, E., Bourgon, R., Brozzi, A., Huber, W., and Steinmetz, L. M. (2008). et al. (2002). The Structure of Haplotype Blocks in the Human Genome. Science High-resolution mapping of meiotic crossovers and non-crossovers in yeast. 296 (5576), 2225–2229. doi: 10.1126/science.1069424 Nature 454 (7203), 479–485. doi: 10.1038/nature07135 Gil, M., Serra, X., Gispert, M., Oliver, M.A.Á., SanÄ udo, C., Panea, B., et al. Manu, S., Acharya, K. K., and Thiyagarajan, S. (2018). Systematic analyses of (2001). The effect of breed-production systems on the myosin heavy chain 1, autosomal recombination rates from the 1000 genomes project uncovers the biochemical characteristics and the colour variables of longissimus thoracis the global recombination landscape in humans. BioRxiv, 246702. doi: from seven spanish beef cattle breeds. Meat Sci. 58, 181–188. doi: 10.1016/ 10.1101/246702 s0309-1740(00)00150-9 McVean, G., Awadalla, P., and Fearnhead, P. (2002). A coalescent-based method González-Rodríguez, A., Munilla, S., Mouresan, E. F., Cañas-Álvarez, J. J., Díaz, C., for detecting and estimating recombination from gene sequences. Genetics 160 Piedrafita, J., et al. (2016). On the performance of tests for the detection of (3), 1231–1241. signatures of selection: a case study with the Spanish autochthonous beef cattle Ministerio de Medio Ambiente y Medio Rural y Marino (2010). Razas de Ganado populations. Genet. Selection Evol. 48 (1). doi: 10.1186/s12711-016-0258-1 Del Catálogo Oficial En España. Madrid, Spain: Ministerio de Medio Ambiente González-Rodríguez, A., Munilla, S., Mouresan, E. F., Cañas-Álvarez, J. J., Baro, y Medio Rural y Marino. J. A., Molina, A., et al. (2017). Genomic differentiation between asturiana de los Mokry, F., Buzanskas, M., Mudadu, M. A., Grossi, D. A., Higa, R., Ventura, R., valles, avileña-negra ibérica, bruna dels pirineus, morucha, pirenaica, retinta et al. (2014). Linkage disequilibrium and haplotype block structure and rubia gallega cattle breeds. Animal 11 (10), 1667–1679. doi: 10.1017/ in a composite beef cattle breed. BMC Genomics 15 (Suppl 7), S6. doi: S1751731117000398 10.1186/1471-2164-15-S7-S6 Graffelman, J., Balding, D. J., Gonzalez-Neira, A., and Bertranpetit, J. (2007). Mouresan, E. F., González-Rodríguez, A., Cañas-Álvarez, J. J., Díaz, C., Altarriba, Variation in estimated recombination rates across human populations. Hum. J., Baro, J. A., et al. (2017). On the haplotype diversity along the genome in Genet. 122 (3–4), 301–310. doi: 10.1007/s00439-007-0391-6 spanish beef cattle populations. Livestock Sci. 201, 33–33. doi: 10.1016/j. Grossi, D. A., Jafarikia, M., Brito, L. F., Buzanskas, M. E., Sargolzaei, M., and livsci.2017.04.015 Schenkel, F. S. (2017). Genetic diversity, extent of linkage disequilibrium and Myers, S., Bottolo, L., Freeman, C., McVean, G., and Donnelly, P. (2005). A fine- persistence of gametic phase in canadian pigs. BMC Genet. 18 (1), 6. doi: scale map of recombination rates and hotspots across the human genome. 10.1186/s12863-017-0473-y Science 310 (5746), 321–324. doi: 10.1126/science.1117196 Guryev, V., Smits, B. M. G., Belt, J. V. D., Verheul, M., Hubner, N., and Cuppen, E. Nachman, M. W. (2002). Variation in recombination rate across the genome: (2006). Haplotype block structure is conserved across mammals. PloS Genet. 7, evidence and implications. Curr. Opin. Genet. Dev. 12 (6), 657–663. doi: 1111–1118. doi: 10.1371/journal.pgen.0020121 10.1016/s0959-437x(02)00358-1 Hayes, B. J., Visscher, P. M., McPartlan, H. C., and Goddard, M. E. (2003). Novel Nguyen, H., and James, E. A. (2016). Immune recognition of citrullinated epitopes. multilocus measure of linkage disequilibrium to estimate past effective Immunology 149 (2), 131–138. doi: 10.1111/imm.12640 population size. Genome Res. 13 (4), 635–643. doi: 10.1101/gr.387103 Paigen, K., and Petkov, P. (2010). Mammalian recombination hot spots: Hill, W. G., and Robertson, A. (1968). Linkage disequilibrium in finite populations. properties, control and evolution. Nat. Rev. Genet. 11 (3), 221–233. doi: Theor. Appl. Genet. 38 (6), 226–231. doi: 10.1007/BF01245622 10.1038/nrg2712 Hinrichs, A. L., Larkin, E. K., and Suarez, B. K. (2009). Population stratification Park, L. (2012). Linkage disequilibrium decay and past population history in the and patterns of linkage disequilibrium. Genet. Epidemiol. 33, 88–92. doi: human genome. PloS One 7 (10), e46603. doi: 10.1371/journal.pone.0046603. 10.1002/gepi.20478 Edited by Thomas Mailund. Hulten, M. A., Palmer, R. W., and Laurie, D. A. (1982). Chiasma derived genetic Petit, M., Astruc, J.-M., Sarry, J., Drouilhet, L., Fabre, S., Moreno, C., et al. maps and recombination fractions: chromosome 1. Ann. Hum. Genet. 46 (2), (2017). Variation in recombination rate and its genetic determinism in sheep 167–175. doi: 10.1111/j.1469-1809.1982.tb00707.x populations. Genetics 207 (2). doi: 10.1534/genetics.117.300123 Jensen-Seaman, M. I., Furey, T. S., Payseur, B. A., Lu, Y., Roskin, K. M., Chen, C.-F., Piedrafita, J., Quintanilla, R., Sañudo, C., Olleta, J.-L., Campo, M.-M., Panea, B., et al. (2004). Comparative recombination rates in the rat, mouse, and human et al. (2003). Carcass quality of 10 beef cattle breeds of the southwest of europe genomes. Genome Res. 14, 528–538. doi: 10.1101/gr.1970304 in their typical production systems. Livestock Prod. Sci. 82 (1), 1–13. doi: Kaback, D., Guacci, V., Barber, D., and Mahon, J. (1992). Chromosome size- 10.1016/S0301-6226(03)00006-X dependent control of meiotic recombination. Science 256 (5054), 228–232. doi: Popescu, P. C. (1990). Chromosomes of the cow and the bull. Adv. Vet. Sci. Comp. 10.1126/science.1566070 Med. 34, 41–71. doi: 10.1016/B978-0-12-039234-6.50007-0 Frontiers in Genetics | www.frontiersin.org 11 November 2019 | Volume 10 | Article 1170 Mouresan et al. Mapping Recombination Rate from the LD Phase Sandor, C., Li, W., Coppieters, W., Druet, T., Charlier, C., and Georges, M. (2012). Tiemann-Boege, I., Schwarz, T., Striedner, Y., and Heissl, A. (2017). The consequences Genetic variants in REC8, RNF212, and PRDM9 influence male recombination of sequence erosion in the evolution of recombination hotspots. Philos. Trans. R. in cattle. Edited by Kenneth Paigen. PloS Genet. 8 (7), e1002854. doi: 10.1371/ Soc. B: Biol. Sci. 372 (1736), 20160462. doi: 10.1098/rstb.2016.0462 journal.pgen.1002854 Villa-Angulo, R., Matukumalli, L. K., Gill, C. A., Choi, J., Tassell, C. P. V., and Sarbajna, S., Denniff, M., Jeffreys, A. J., Neumann, R., Artigas, M. S., Veselis, A., Grefenstette, J. J. (2009). High-resolution haplotype block structure in the cattle et al. (2012). A major recombination hotspot in the xqyq pseudoautosomal genome. BMC Genet. 10 (1), 19. doi: doi.org/10.1186/1471-2156-10-19 region gives new insight into processing of human gene conversion events. Wall, J. D., and Stevison, L. S. (2016). Detecting recombination hotspots from Hum. Mol. Genet. 21 (9), 2029–2038. doi: 10.1093/hmg/dds019 patterns of linkage disequilibrium. G3 (Bethesda) 6 (8), 2265–2271. doi: Shen, B., Jiang, J., Seroussi, E., Liu, G. E., and Ma, L. (2018). Characterization of 10.1534/g3.116.029587 recombination features and the genetic basis in multiple cattle breeds. BMC Xu, L., Zhu, B., Wang, Z., Xu, L., Liu, Y., Chen, Y., et al. (2019). Evaluation of Genomics 19 (1), 1–10. doi: 10.1186/s12864-018-4705-y linkage disequilibrium, effective population size and haplotype block structure Smukowski, C. S., and Noor, M. A. F. (2011). Recombination rate variation in chinese cattle. Animals 9 (3), 83. doi: 10.3390/ani9030083 in closely related species. Heredity 107 (6), 496–508. doi: 10.1038/hdy.2011.44 Zachariou, V., Duman, R. S., and Nestler, E. (2012). “G proteins,” in Basic Stapley, J., Feulner, P. G. D., Johnston, S. E., Santure, A. W., and Smadja, C. M. neurochemistry. Eds. Brady, S. T., Siegel, G. J., Albers, R. W., and Price, D. L. (2017). Variation in recombination frequency and distribution across (Waltham, MA. USA: Academic Press), 411–422. doi: 10.1016/B978-0-12- eukaryotes: patterns and processes. Philos. Trans. R. Soc. B: Biol. Sci. 372 (1736), 374947-5.​00021-3 20160455. doi: 10.1098/rstb.2016.0455 Stevison, L. S., Woerner, A. E., Kidd, J. M., Kelley, J. L., Veeramah, K. R., Conflict of Interest: The authors declare that the research was conducted in the McManus, K. F., et al. (2016). The time scale of recombination rate evolution absence of any commercial or financial relationships that could be construed as a in great apes. Mol. Biol. Evol. 33 (4), 928–945. doi: 10.1093/molbev/msv331 potential conflict of interest. Sturtevant, A. H. (1913). The linear arrangement of six sex-linked factors in drosophila, as shown by their mode of association. J. Exp. Zool. 14 (1), 43–59. Copyright © 2019 Mouresan, González-Rodríguez, Cañas-Álvarez, Munilla, doi: 10.1002/jez.1400140104 Altarriba, Díaz, Baró, Molina, Lopez-Buesa, Piedrafita and Varona. This is an open- Tenesa, A., Navarro, P., Hayes, B. J., Duffy, D. L., Clarke, G. M., Goddard, M. E., access article distributed under the terms of the Creative Commons Attribution et al. (2007). Recent human effective population size estimated from License (CC BY). The use, distribution or reproduction in other forums is permitted, linkage disequilibrium. Genome Res. 17 (4), 520–526. doi: 10.1101/gr.6023607 provided the original author(s) and the copyright owner(s) are credited and that the The International HapMap Consortium (2007). A second generation human original publication in this journal is cited, in accordance with accepted academic haplotype map of over 3.1 million SNPs. Nature 449 (7164), 851–861. doi: practice. No use, distribution or reproduction is permitted which does not comply 10.1038/nature06258 with these terms. Frontiers in Genetics | www.frontiersin.org 12 November 2019 | Volume 10 | Article 1170