# The Identification of Runs of Homozygosity Gives a Focus on the Genetic Diversity and Adaptation of the “Charolais de Cuba” Cattle

Yoel Rodríguez-Valera, Dominique Rocha, Michel Naves, Gilles Renand, Eliecer Pérez-Pineda, Yuliaxis Ramayo-Caldas, Sebastian E Ramos-Onsins

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yoel Rodríguez-Valera, Dominique Rocha, Michel Naves, Gilles Renand, Eliecer Pérez-Pineda, et al. The Identification of Runs of Homozygosity Gives a Focus on the Genetic Diversity and Adaptation of the “Charolais de Cuba” Cattle. Animals, 2020, 10 (12), pp.2233. 10.3390/ani10122233. hal-03237098

HAL Id: hal-03237098

https://hal.inrae.fr/hal-03237098

Submitted on 26 May 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Identification of Runs of Homozygosity Gives a Focus on the Genetic Diversity and Adaptation of the "Charolais de Cuba" Cattle

Article in Animals · November 2020

DOI: 10.3390/ani10122233

CITATIONS

0

READS

74

7 authors, including:

<figure><img src="image_3.png" /></figure>

Yoel rodriguez valera

University of Granma

35 PUBLICATIONS

40 CITATIONS

SEE PROFILE

<figure><img src="image_4.png" /></figure>

Gilles Renand

French National Institute for Agriculture, Food, and Environment (INRAE)

208 PUBLICATIONS

3,526 CITATIONS

SEE PROFILE

<figure><img src="image_5.png" /></figure>

Michel Naves

French National Institute for Agriculture, Food, and Environment (INRAE)

9 PUBLICATIONS

10 CITATIONS

SEE PROFILE

<figure><img src="image_6.png" /></figure>

Yuliaxis Ramayo-Caldas

French National Institute for Agriculture, Food, and Environment (INRAE)

180 PUBLICATIONS

1,543 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

<figure><img src="image_7.png" /></figure>

Feed-a-Gene View project

<figure><img src="image_8.png" /></figure>

Relationship between digestive microbiota and methane emission in cattle View project
<figure><img src="image_9.png" /></figure>

Article

# The Identification of Runs of Homozygosity Gives a Focus on the Genetic Diversity and Adaptation of the "Charolais de Cuba" Cattle

Yoel Rodríguez-Valera<fnref n="1" />, Dominique Rocha<fnref n="2" />, Michel Naves<fnref n="3" />, Gilles Renand<fnref n="2" />,
Eliecer Pérez-Pineda<fnref n="1" />, Yuliaxis Ramayo-Caldas<fnref n="2" />,*, and Sebastian E. Ramos-Onsins<fnref n="5" />,\*

Received: 9 November 2020; Accepted: 23 November 2020; Published: 27 November 2020

**Simple Summary:** The Charolais de Cuba cattle is a tropical adapted breed founded in Cuba around 120 years ago from Charolais French specimens. Nowadays, it is still a closed breed and remains as a small population. In this work, we analyzed the inbreeding and diversity patterns, as well as the population size, of this recent adapted breed via a run of homozygosity (ROH) analysis. We found that the genomic inbreeding levels are higher in the Charolais de Cuba breed compared to French and British Charolais populations. Nevertheless, we detected that the effective population size experienced a very similar decline during the last century in the three Charolais populations studied. Finally, a number of regions with exceptional patterns of long homozygosity were identified in this breed, and these could be related to processes of adaptation to tropical conditions.

**Abstract:** Inbreeding and effective population size (Ne) are fundamental indicators for the management and conservation of genetic diversity in populations. Genomic inbreeding gives accurate estimates of inbreeding, and the Ne determines the rate of the loss of genetic variation. The objective of this work was to study the distribution of runs of homozygosity (ROHs) in order to estimate genomic inbreeding (FROH) and an effective population size using 38,789 Single Nucleotide Polymorphisms (SNPs) from the Illumina Bovine 50K BeadChip in 86 samples from populations of Charolais de Cuba (n = 40) cattle and to compare this information with French (n= 20) and British Charolais (n = 26) populations. In the Cuban, French, and British Charolais populations, the average estimated genomic inbreeding values using the FROH statistics were 5.7%, 3.4%, and 4%, respectively. The dispersion measured by variation coefficient was high at 43.9%, 37.0%, and 54.2%, respectively. The effective population size experienced a very similar decline during the last century in Charolais de Cuba (from 139 to 23 individuals), in French Charolais (from 142 to 12), and in British Charolais (from 145 to 14) for the ~20 last generations. However, the high variability found in the ROH indicators and FROH reveals an opportunity for maintaining the genetic diversity of this breed with an adequate mating strategy, which can be favored with the use of molecular markers. Moreover, the detected ROH were compared to previous results obtained on the detection of signatures of selection in the same breed. Some of the observed signatures were confirmed by the ROHs, emphasizing the process of adaptation to tropical climate experienced by the Charolais de Cuba population.

<footnote n="1">¹ Faculty of Agricultural Sciences, University of Granma, Bayamo 95100, Cuba;
yrodriguezvalera@udg.co.cu (Y.R.-V.); eperezp@udg.co.cu (E.P.-P.)</footnote>

<footnote n="2">² GABI, INRAE, AgroParisTech, University Paris-Saclay, F-78350 Jouy-en-Josas, France;
dominique.rocha@inrae.fr (D.R.); gilles.renand@inrae.fr (G.R.)</footnote>

<footnote n="3">³ INRAE, URZ, 97170 Petit Bourg, Guadeloupe, France; michel.naves@inrae.fr</footnote>

<footnote n="4">⁴ Animal Breeding and Genetics Program, Institute for Research and Technology in Food and Agriculture (IRTA), Torre Marimon, 08140 Caldes de Montbui, Spain</footnote>

<footnote n="5">5 Plant and Animal Genomics, Centre of Research in Agricultural Genomics (CRAG) Consortium CSIC-IRTA-UAB-UB, Campus UAB, 08193 Bellaterra, Spain</footnote>

<footnote>* Correspondence: yuliaxis.ramayo@irta.cat (Y.R.-C.); sebastian.ramos@cragenomica.es (S.E.R.-O.)</footnote>
**Keywords:** genomic inbreeding; effective population size; runs of homozygosity; candidate genes; tropical climate; adaptation

# 1. Introduction

As a result of importing superior breeds for meat production, which was the main objective before 1959 for the Cuban livestock industry, Charolais sires and dams were imported directly from France to Cuba from 1900 until 1923. These animals were the ancestors of the current Charolais de Cuba breed. The main destination of these imported animals was the east of the country, exactly in Jiguani (today Granma province) on the Farm San José del Retiro with a current census of 700 animals. The Charolais de Cuba (CHCU) breed was then created and became adapted to tropical conditions. Other subsequent imports of Charolais animals were made in 1944 from Mexico, in 1957 from Texas, and finally in 1969 from Canada [1]. Since then, no new animals have been introduced to the CHCU population, which has evolved in a closed and isolated way, which implies the possibility that there has been a decrease in genetic diversity in this breed, as well as an increase in consanguinity.

Inbreeding is defined as the probability that two alleles at the same locus are identical by descent from a common ancestor [2]. The increase in inbreeding results in an increase in the frequency of homozygous genotypes, so this indicator can be used as an estimator of inbreeding. However, alleles that are identical by descent (IBD) and by state (IBS) cannot be differentiated, and both could be included in the same indicator. Runs of homozygosity (ROH), which are continuous and uninterrupted regions of DNA sequences without heterozygosity in the diploid state [3], have been proposed to solve this problem. ROHs allow for the distinguishing of these differences in many cases, and have been analyzed in several species [4–6].

In finite-size populations, inbreeding produces changes in genotypic frequencies, increasing homozygosity at expenses of heterozygosity, while allelic frequencies remain constant [7]. The presence of autozygosity is a consequence of individuals with a common ancestor that transmitted an identical segment of chromosome to their progeny, resulting in homozygous regions and an increase in the number of ROHs [8].

The concept of genomic inbreeding was recently introduced [9] in order to obtain a direct and accurate estimate of inbreeding. Genomic inbreeding is measured by $F_{ROH}$, defined as the proportion of the (autosomal) genome in homozygous state while assuming that homozygous regions (ROHs) of the genome could only be the result of inbreeding, resulting in autozygosity [4]. The existence of medium to high correlation between $F_{ROH}$ and other inbreeding indices suggests that the levels of autozygosity derived from ROHs can be used as estimators of individual inbreeding [10].

ROHs are more common in regions of high linkage disequilibrium (LD) and low percentage of recombination [3]. For example, in a study at the genomic level in individuals from human European populations [11], a non-uniform distribution of ROHs was detected, with almost identical patterns in all the studied populations. This structure could not only be explained by a reduced genetic diversity (genetic diversity explained only 7% of the ROH variation), a localized heterozygosity deficit, or a regional increase in LD. Historical and demographic events significantly modulated ROH patterns, detecting a geographic decline associated with population density and mobility in European populations.

In addition to inbreeding, one of the most widely used indicators to assess the genetic diversity of a population is its effective population size (Ne) [12–16], which is defined as the size of an ideal population that could have the same amount of genetic drift or the same inbreeding as the population under consideration [17]. It is an important parameter in quantitative and population genetics and in evolutionary and conservation biology since it determines the rate of loss of genetic variation [18]. Maintaining genetic diversity within a population is accomplished by maximizing effective size (Ne) or equivalently minimizing the increase in inbreeding across generations with appropriate mating schemes [10].
A review of the main methods used to estimate Ne [19] includes those that use LD. In 1971, a simple method based on the assumption of two loci of a finite population—monoecious and diploid—with random mating, the absence of selection, and mutation was proposed [20]. In 1981, a method to estimate Ne was implemented for a population with random crossover between neutral genes at various polymorphic loci [21]. The efficiency of the method was found to increase as the sample size and the number of genes with high LD values increased. A procedure based on the multilocus LD of a homozygous chromosome segment (CSH) [21] showed less variability in this measure than using LD, and it was found that CSH over long distances reflects better estimates of recent Ne under the assumption of constant linear growth and no mutation. However, in natural populations, the assumption of a linear change in Ne is usually unrealistic, as usually occurs in a much more complex way [15]. Nevertheless, the very versatile algorithm included in the software SNeP v1.0 [14] is based on different Ne dynamics of population size change across time [13] and, thus, allows for the calculation of the expected value of the linkage disequilibrium (E (r2)) under almost all population dynamics [15].

The objective of this work was to calculate ROH values to estimate the genomic inbreeding and effective size of the CHCU breed before comparing this information with the values observed in French (CHF) and British Charolais (CHUK) populations. In addition, this study explored the relationship of the ROHs observed with previous results obtained in the same breed with signatures of selection [22].

## 2. Materials and Methods

### 2.1. Data Used and Marker Quality Controls

The analysis was performed on populations of CHCU (n = 40), CHF (n = 20), and CHUK (n = 26). The used data were obtained from previous studies using the Illumina Bovine SNP50 genotyping chip on CHUK [23], CHF [24], and CHCU [22]. A single genotype file was made with the three populations using Plink 1.9 [25]. The markers were filtered with the following parameters: missing genotypes (geno > 0.1), minor allele frequency (MAF < 0.01), and Hardy-Weinberg equilibrium (p < 0.01). SNPs that did not meet these criteria were excluded. Therefore, we retained 38,789 SNPs for the analysis of Ne estimation. LD pruning was additionally applied to eliminate SNPs with high LD values that were not really IBD, as previously recommended [25], thus leaving a final number of 33,038 SNPs. After filtering, genomic inbreeding was estimated by a procedure recommended for evaluating genomic ROHs in small populations (small in relation to the perception of the risk that this population would become extinct [26]).

### 2.2. Detection of Runs of Homozygosity (ROH)

We used the same procedure implemented by the authors of [27] with the library DetectRuns [28] in R 3.4.2. Briefly, ROH detection was determined in sliding windows of 20 SNPs (--homozyg-window-snp 20), allowing for no more than one missing SNP (--homozyg-window-missing 1) and one heterozygous SNP per window (--homozyg-window-het 1). The minimum length of an ROH segment was 4 Mb (--homozyg-kb 4000), as recommended when using the Illumina Bovine 50K BeadChip [29]. The minimum SNP density was 1 SNP per 100 kb (--homozyg-density 1/100), and the maximum gap between two consecutive SNPs was 1000 kb (--homozyg-gap 10⁶ bp). Finally, the rate in which an SNP was included in the total of sliding windows was at least 0.05 (--homozyg-window-threshold 0.05). ROHs were classified into three classes according to their size: 4–8, 8–16, and >16 Mb, as previously reported [5,26,29,30].

### 2.3. Estimation of Genomic Inbreeding

The parameter used to calculate genomic inbreeding was the $F_{ROH}$ through the formula $F_{ROH} = \Sigma L_{ROH}/L_{AUTOSOME}$ (proportion of the genome in ROHs over the total genome of 2,504,168,970 bp covered by SNPs) [4]. The estimation, summary, and graphs of $F_{ROH}$ were made using the
aforementioned DetectRuns R library. The average size of ROH (MN<sub>ROH</sub>), defined as the average of the size of each ROH segment (L<sub>ROH</sub>), and the sum of all ROH segments per animal (S<sub>ROH</sub>) were estimated. Furthermore, the F<sub>ROH</sub> by chromosomes (BTA) (F<sub>ROH_BTA</sub>) was estimated similarly, using the formula:

$$F_{ROH\_BTA} = \sum L_{ROH\_BTA}/L_{BTA} \qquad (1)$$

where ΣL<sub>ROH_BTA</sub> is the sum of the size of each ROH segment per BTA and L<sub>BTA</sub> is the size of each BTA covered by SNPs [6].

## 2.4. Estimation of Effective Population Size

The estimation of the recent demographic history of the three populations was calculated through the trend of the effective size of the population (Ne) using SneP v1.1 [14]. This software estimates the effective population size through the following equation [13]:

$$N_{T(t)} = (4f(c_t))^{-1} (E[r_{adj}^2|c_t]^{-1} - \alpha) \qquad (2)$$

where N<sub>T(t)</sub> is the effective size of the population in t past generations, equivalent to t = (2f(c<sub>t</sub>))<sup>-1</sup> [31], f(c<sub>t</sub>) is a function related to the recombination rate (which is approximately equal to c<sub>t</sub> for small recombination values and a second order function for larger values), c<sub>t</sub> is the recombination rate defined by a specific physical distance between markers, r<sup>adj</sup> is the adjusted LD value for the sample size (40 animals for CHCU), and α is a correction for the occurrence of mutations. The proposed α values were 1, 2, or 2.2 [20,32,33].

This procedure is based on the relationships between the variance of the LD at adjacent SNPs and the effective size of the population in the presence of the mutation, which allows one to estimate recent and old Ne [34]. Thus, we studied the trend of the effective population size of the three populations, firstly in ~20 generations, that corresponded to the approximate time since the import of the French Charolais in Cuba [1,35]; and secondly, estimating the population size until ~35 generations, which approximately corresponds to the creation of the French Charolais Herd Book in 1842 [36].

A minimum allele frequency (MAF) ≥ 0.05 was used to avoid low frequency errors and to improve the Ne estimate of [13,27]. Other included options were sample size (40, 20, and 26 for CHCU, CHF, and CHUK, respectively), -mindist (2,400,000 and 1,400,000 pb) and -maxdist (50,000,000 and 2,400,000 pb—20 and from 21 to 35 past generations, respectively), which are the minimum and maximum distances to be analyzed between SNP pairs in bp.

## 2.5. Candidate Regions Responsible of the CHCU Adaptation to Tropical Conditions

The ROHs reported in this study were compared with regions under selection previously reported in the same CHCU population [22]. In addition, genes within the ROH intervals were annotated using the Map Viewer of National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/mapview/) of the bovine genome assembly release UMD3.1.1. Gene ontology and pathway analyses were carried out using the Panther version 13.1 software tool [37] (http://pantherdb.org).

# 3. Results

## 3.1. Estimation of Runs of Homozygosity (ROH) and Genomic Inbreeding

After merging and performing the quality control of the data, we finally retained 33,038 SNPs for data analyses. In concordance with [29], we considered the definition of Regions of Homozygosity (ROHs) as the segments of homozygous genotypes greater than 4 Mb in length, with a minimum number of 20 SNPs; we allowed for at least one heterozygous SNP to estimate genomic inbreeding [26] in order to avoid overestimating the number of ROHs [29]. Note that the utilization of genotype markers does not allow for the detection of all existent variants, so the number of ROHs was possibly
overestimated. We assessed the impact of Single Nucleotide Polymorphism (SNP) density and genotyping errors using the Illumina Bovine 50 SNP genotyping chip on ROH identification, and we revealed an abundance of small segments (overestimated the numbers of segments 1–4 Mb long), suggesting that it was not sensitive enough for the precise determination of small segments. However, it proved suitable for detecting segments longer than 4 Mb. In the current study, the values of FROH > 4 Mb are shown in Table 1. These values were found to be higher on average in CHCU (0.057), intermediate in CHUK (0.04), and lower in CHF (0.034), and they had medium-high levels of dispersion, with coefficients of variation of 0.44, 0.54, and 0.37 for Cuban, French and British Charolais, respectively.

Table 1 shows a decrease in FROH in the three populations as the segment size increased. In all cases, estimates of recent inbreeding values (FROH at larger sizes) were lower compared to estimates of ancestral inbreeding (FROH at smaller sizes). The CHCU breed showed the highest FROH values for the 4–8, 8–16, and >16 Mb classes compared to CHF and CHUK. An inversely proportional relationship between the FROH values and the size of the segments of each size interval was also observed. The distributions of the number of ROHs according to the different categories are shown in Table 2.

<table><caption>**Table 1.** Genomic inbreeding (FROH) according to three intervals (Mb) in three Charolais populations.
CHCU: Charolais de Cuba; CHF: French Charolais; CHUK: British Charolais.</caption><thead><tr><th>Breeds</th><th>4–8 Mb<br>FROH (±SD)</th><th>8–16 Mb<br>FROH (±SD)</th><th>>16 Mb<br>FROH (±SD)</th></tr></thead><tbody><tr><td>CHF</td><td>0.034 (0.012)</td><td>0.017 (0.011)</td><td>0.007 (0.007)</td></tr><tr><td>CHCU</td><td>0.057 (0.025)</td><td>0.037 (0.022)</td><td>0.016 (0.016)</td></tr><tr><td>CHUK</td><td>0.040 (0.022)</td><td>0.023 (0.019)</td><td>0.011 (0.014)</td></tr></tbody></table>

FROH: Genomic inbreeding obtained from runs of homozygosity (ROHs) of 4–8, 8–16, and >16 Mb.
±SD: standard deviation.

<table><caption>**Table 2.** Number of runs of homozygosity (ROHs) by breed and interval in Mb.</caption><thead><tr><th rowspan="2">Interval (Mb)</th><th colspan="3">CHF</th><th colspan="3">CHCU</th><th colspan="3">CHUK</th></tr><tr><th>Number</th><th>Freq (%)</th><th>MN<sub>ROH</sub> (Mb)</th><th>Number</th><th>Freq (%)</th><th>MN<sub>ROH</sub> (Mb)</th><th>Number</th><th>Freq (%)</th><th>MN<sub>ROH</sub> (Mb)</th></tr></thead><tbody><tr><td>4–8</td><td>158</td><td>71.5</td><td>5.35</td><td>354</td><td>58.3</td><td>5.549</td><td>193</td><td>65.6</td><td>5.598</td></tr><tr><td>8–16</td><td>47</td><td>21.3</td><td>11.072</td><td>182</td><td>29.9</td><td>11.24</td><td>73</td><td>24.9</td><td>10.633</td></tr><tr><td>>16</td><td>16</td><td>7.2</td><td>22.195</td><td>72</td><td>11.8</td><td>23.634</td><td>28</td><td>9.5</td><td>27.377</td></tr></tbody></table>

Number: ROH number found in each class. Freq: percentage of ROHs in each class. MN<sub>ROH</sub>: average ROH size in Mb in each class.

For each chromosome, the FROH_BTA was determined as the proportion of BTA in ROHs over the length of the chromosome (BTA) covered by the studied SNPs (Figure 1 and Supplementary Materials Table S1). We observed a high variability within and between chromosomes (BTA), as well as high variability in the size and location of ROHs. Some animals showed high inbreeding values per chromosome; in CHF, FROH_BTA was >0.4 in BTA15 and 16; in CHCU, FROH_BTA exceeded 0.8 in BTA 25 and 0.5 in six other BTAs (8, 12, 14, 17, 21, and 26; see Figure 1); and in CHUK, three chromosomes (BTAs 4, 15, and 24) presented values > 0.5. In addition, there were overlaps between breeds with FROH_BTA values higher than 0.04: CHF and CHUK in BTA15 and CHCU and CHUK in BTA8 and BTA25. The distribution of the number, average size, and frequency of ROHs by chromosomes are shown in Table S2.

The distribution of the number, average size, and frequency of ROHs per chromosome are shown in Table S2. We also estimated the average size of ROHs (Table 3) and the relationship between
the total homozygous regions higher than 4 Mb with the size of each ROH in the three Charolais populations (Figure 2).

<figure><img src="image_10.png" /><figcaption>**Figure 1.** Distribution of genomic inbreeding estimates for each chromosome (FROH_BTA), calculated as the proportion of BTA (FROH by chromosomes) in ROHs over the length of the BTA covered by the involved SNPs. The x-axis represents the chromosomes, and the y-axis represents FROH values in boxplot; black points are maximum values. ROH, run of homozygosity; CHF, French Charolais; CHCU, Charolais de Cuba; and CHUK, British Charolais.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 2.** Relationship between the total number of runs of homozygosity (ROHs) >4 Mb and the total length (Mbps) of genome in such ROHs for individuals from the three analyzed Charolais populations. Each dot represents an individual.</figcaption></figure>

<table><caption>**Table 3.** Average size of ROH (MN<sub>ROH</sub>) in three Charolais populations.</caption><thead><tr><th>Populations</th><th>Class Mb</th><th>N<sub>ROH</sub></th><th>MN<sub>ROH</sub> (Mb)</th><th>±SD</th><th>CV %</th><th>IC -95%</th><th>IC +95%</th></tr></thead><tbody><tr><td rowspan="3">CHF</td><td>4-8</td><td>158</td><td>5.34</td><td>1.08</td><td>20.28</td><td>5.18</td><td>5.52</td></tr><tr><td>8-16</td><td>47</td><td>11.07</td><td>2.27</td><td>20.49</td><td>10.40</td><td>11.73</td></tr><tr><td>>16</td><td>16</td><td>22.19</td><td>5.65</td><td>25.49</td><td>19.18</td><td>25.21</td></tr><tr><td rowspan="3">CHCU</td><td>4-8</td><td>354</td><td>5.54</td><td>1.14</td><td>20.69</td><td>5.42</td><td>5.67</td></tr><tr><td>8-16</td><td>182</td><td>11.24</td><td>2.19</td><td>19.51</td><td>10.91</td><td>11.56</td></tr><tr><td>>16</td><td>72</td><td>23.63</td><td>9.35</td><td>39.58</td><td>21.43</td><td>25.83</td></tr><tr><td rowspan="3">CHUK</td><td>4-8</td><td>193</td><td>5.60</td><td>1.14</td><td>20.37</td><td>5.44</td><td>5.76</td></tr><tr><td>8-16</td><td>73</td><td>10.63</td><td>2.07</td><td>19.47</td><td>10.15</td><td>11.11</td></tr><tr><td>>16</td><td>28</td><td>27.37</td><td>13.22</td><td>48.29</td><td>22.25</td><td>32.50</td></tr></tbody></table>

N<sub>ROH</sub>: Numbers of ROH. MN<sub>ROH</sub> (Mb): average ROH size in Mb. ±SD: Standard deviation. CV (%): Coefficient of variation. IC -95%: Confidence interval; IC +95%: Confidence interval.

## 3.2. Estimation of Effective Populations Size

During the recent period (Figure 3a), the observed trend showed a similar decrease in the three populations (CHCU: from 139 to 23; CHF: from 142 to 12; and CHUK: from 145 to 14). However, from generation 11 to the most recent generation, the CHCU population was found to have experienced a lower decrease than the others. Before this period (Figure 3b), the observed trend of the Ne was decreasing in a very similar way between the CHCU (from 216 to 143), CHF (from 220 to 146), and CHUK (from 225 to 147) as expected, because CHCU originated from ancestral individuals of these breeds.
<figure><img src="image_12.png" /><figcaption>**Figure 3.** Effective size of three populations of the Charolais breed: (a) approximately 20 generations and (b) from 21 to 35 generations. Blue line, Charolais de Cuba; red line, French Charolais; and yellow line, British Charolais. x-axis: past generations; y-axis: effective population size.</figcaption></figure>

### 3.3. *Candidate Regions Responsible of the CHCU Adaptation to Tropical Conditions*

We compared the patterns of homozygosity between CHCU and CHF with the Extended Haplotype Homozygosity (EHH) intervals reported in a previous study [22]. We observed overlaps in the regions of the genome that were previously classified as candidates for being affected by
adaptation processes [22] by using the estimates of FROH: 102 ROHs (12.30% of the total ROHs identified in the CHF and CHCU samples) distributed in the 29 BTAs that contained previously reported 50 regions (48.07% of the total). The FROH estimate fluctuated in these regions between a minimum of 0.002 and a maximum of 0.022.

# 4. Discussion

## 4.1. Limitations and Trade-Offs of the Analysis Based on Genotype Datasets and Small Sample Sizes

The use of genotypes obtained from genotyping SNP chips that are not specifically designed for the target population can bias the results because many exclusive variants would not be detected. This effect, named ascertainment bias, is a potential problem for the analysis of variability. Here, the effect of ascertainment bias (the SNP50 chip was designed mostly for the Holstein-Friesian breed) is apparently inexistent for the following reasons: (i) Charolais is a taurine breed that is closely related to the Holstein-Friesian breed in terms of the origin and evolution of the mutations [22]; (ii) the number of SNPs detected (38,789 before filtering for LD; see the Materials and Methods section) was on the order of found in other taurine populations, (e.g., Brown Swiss at 38,710, Pinzganer at 38,198, and Tyrol Grey at 42,997 for the same SNP chip [29]); (iii) here, small ROH regions that could be more affected by ascertainment bias (<4 Mb; see [29]) were excluded from the analysis; (iv) the comparative analysis of the distribution of FROH between populations was expected to be weakly affected by ascertainment bias because this effect would be very similar for all related Charolais populations. Nevertheless, the identified ROHs must be interpreted with caution, as they were only candidate regions for selection.

From a strict point of view, homozygosity is only accurately calculated with genome sequencing data, because it is the only way to unbiasedly estimate the levels of heterozygosity and, thus, homozygosity. The use of SNP chips, whatever considered density, eliminate private SNPs that only appear in one population. The use of medium- or high-density SNPs are approaches that detect a number of candidate regions effected by positive selection. Indeed, an optimal genotyping chip would satisfy the presence of separated SNPs with a good enough recombination rate to allow for separated segregation between them. This optimal separation is difficult to estimate and depends mainly on the population size, demography, and number of generations that mutations create. In the case of the Cuban Charolais population, which is a very recent population with small number of individuals, the chances for recombination are relatively low, suggesting that a medium SNP genotyping chip may be enough for such analyses. That is, in case of acting disparate evolutionary forces between sister populations (such as selective pressure), differences between populations would be observed because such forces would affect large genomic regions.

A large sample size would increase the number of detected variants but would also increase the effort required to search for samples in the field, as well as the genotyping cost. In this work, the optimal sample size necessary for such analysis was that which allowed for the detection of the effect of differential patterns in homozygosity. Here, we analyzed 40 samples of our target population (CHCU) plus 20 and 26 samples from the CHF and CHUK populations, respectively. Though this was a small sample size, the number of SNPs expected to be detected for a sample size of n = 700 diploid individuals (i.e., the entire CHCU population, which is around 200 times in relation to n = 40) is only 1.6 times larger than that expected to be detected using n = 40 (under standard assumptions [38]). In other words, although a small sample size was used and the results were limited, the study of a large sample in homozygosity analysis would not greatly increase the number of detected SNPs. Note that an analysis of homozygosity needs a very different experimental design from that in Genome Wide Association Studies (GWAS) studies, where a large number of samples is needed to achieve meaningful association results of the genome variants.

## 4.2. Estimation of Runs of Homozygosity (ROHs) and Genomic Inbreeding
Inbreeding is a highly variable measure between different cattle populations that largely depends on the isolation of the population, size, and genetic improvement programs used, among other factors [39]. Additionally, the number of SNPs, the density of the SNP chip, and the selection criteria for the used SNPs to determine the ROHs can have a huge effect on FROH values [40]. It was expected that CHCU would have higher LD as a result of their small size and isolation. This was observed (Table 1) using the global statistics of FROH and MNROH, although we also observed that the LD determined by MNROH at long fragments (Table 2, MNROH > 16 Kb, affected by recent events) was higher in CHCU in relation to CHF but surprisingly lower compared to CHUK. This surprising result may be explained by an event of recent strong selection pressure in CHUK or by the contribution (although limited) of few animals from other breeds to CHCU [41].

A comparison of our results to the bibliography is presented in Table S3. Similar results in FROH values have been found in local cattle breeds in Italy [26,30], Spain [10], dual-purpose cattle in Austria [5], and US Holsteins [40]. These authors also pointed out the importance of mating strategies in maintaining genetic diversity in local breeds. Lower values have been found in local breeds in Africa [42], while higher values have been observed in Blanc Bleu Belge [43], Spanish Lidia [44], and the Gyr breed in Brazil, which was explained by a small Ne and a high LD value [45,46].

The estimates of FROH not only allowed for the measurement of individual autozygosity but also the distinguishing between ancient and recent inbreeding values. This explains why recombination events interrupt large chromosome segments. ROHs (~10 Mb) appear as a result of recent inbreeding (up to five past generations), while small ROHs (~1 Mb) may indicate a greater ancestral distance (up to 50 generations past) [4]. For example, an inversely proportional relationship of FROH and class segments was found when estimating FROH in Fleckvieh cattle from Austria [5]. In addition, similar results were observed when evaluating the autozygosity, effects of genotyping errors, and density of SNPs using high and medium density panels in the three following bovine breeds: Brown Swiss, Pinzgauer, and Tyrol Grey [29]. Other authors have also observed the same patterns for different populations [47,48].

We observed a high variability within and between chromosomes (BTA) in the size and location of ROHs, which was in agreement with previous studies performed in bovines [10,29,49].

### 4.3. *Estimation of Effective Populations Size*

The method used to estimate Ne uses the LD for this purpose. Several authors have found an inverse relationship between these two indicators when comparing them between various breeds of cattle [16,50,51]. In this same way, it has been reported [22] that the LD model between CHF and CHCU is lower in CHCU compared to CHF, with this decrease in LD observed between both breeds seen as a result of strong selection pressure in CHF (a consequence of the improvement programs of this breed in France) compared to the CHCU.

Some limitations in the estimation of Ne have been attributed to the used methodology [13]. Estimating the Ne variable at different times found results that caused certain doubts, especially for recent generations ($T < 10$) even for a corrected sample size. However, they indicated that methodology is highly sensitive to a value of $\alpha$ (correction for the occurrence of mutations). Finally, these authors advised being extremely careful, especially in terms of conservation, with decisions made based on estimates of Ne when it is assumed to vary. However, these limitations have been addressed in SneP according to its developers [14], although they clarified that the quantitative estimates are still highly dependent on the sample size, the type of estimation of the LD, and discrete interval processes.

The severe reduction of the estimated population size of CHCU can be explained by the founder effect related to the reduced number of sires that were exported to Cuba, as well as to the inbreeding created by the repeated use of these sires, which has been even more dramatic due to the selection scheme and artificial insemination that was massively introduced in Cuba in the 1960s [52]. The occurrence of a bottleneck in French Charolais cattle was previously detected by pedigree analysis by the authors of [53], who argued that it occurred due to the increase in the use of artificial insemination (AI) in the 1990s and the intensity of selection.
Estimates of an effective size smaller or similar to those found in CHCU were obtained in Chinese cattle breeds [16] in the north (Ne = 34), south (Ne = 27), southwest (Ne = 31), and the Wagyu breeds (Ne = 14), with the difference being that in those studies, the Ne values were determined for 1000 generations up to the most recent generation. These authors suggested that it was possible to use genomic selection programs adapted to breeds with a small effective size in these breeds. There is no evidence that new Charolais animals have been imported since 1969 to Cuba [1], so CHCU has a small degree of differentiation with its ancestor CHF (Fst = 0.0497) [22]. In addition, CHCU is a breed adapted to tropical conditions, so CHCU could generate an economic benefit. This represents a real possibility of extending its use through the flow of CHCU alleles to other latitudes, especially in the current context of global warming.

The existence of a high variance in the FROH indicators determined in CHCU (Figure 2) could also have been the result of possible crossbreeding events with animals of other breeds of Charolais [41], as explained by imports made from other countries [1], or the possible introgression of alleles of zebu origin [54]. To avoid the accumulation of deleterious variants that would seriously affect the fitness of a breed, it is recommended to have at least Ne = 50 sires per generation [55], which corresponds to a consanguinity rate of 1%. According to this, all populations in conservation plans should achieve this minimum as a first objective. In each conservation program, despite the existence of a bottleneck, as manifested in CHCU, attempts should be made to conserve genetic diversity. If the population is very small, as in the investigated case, each animal is extremely important because its loss has a great impact on the total size of the population [55].

### 4.4. Candidate Regions Responsible of the CHCU Adaptation to Tropical Conditions

CHCU animals have adapted well to the tropical conditions and have developed a tolerance to high temperatures [22]. The analysis of ROHs allowed us to detect large regions that may be affected by positive selection. These regions should be considered with caution, because there are many reasons to find autozygosity (e.g., a low density of the SNP genotyping chip, genetic drift caused by population reduction or other demographic causes, and regions containing low recombination or low mutation rates, in addition to selective sweeps). Furthermore, the detected candidate regions contain a large number of genes, so it was mostly expected to find genes with functions related to the causal process of interest for the researcher, although they are not involved in such a process [56].

In this work, we found several ROHs overlapping with the regions under selection, in concordance with [22]. A total of 1195 genes in CHF and 4086 genes in CHCU were in the ROHs. BTA5:17.363:33.013, BTA5:6.285:36.84, and BTA5:21.708:33.687 were found to overlap with BTA5:28.660:30.114 [22], and they were considered together as a unique selection signal. Many genes located within the selected regions in CHCU have already been linked to heat stress. For instance, *AQP5* is located in three ROHs, as also reported as putative selective sweeps by the authors of [22]. *AQP5* encodes a channel protein that selectively transports water through the plasma membrane of secretory and absorptive cells found, for example, in salivary or sweating glands [57].

Other genes, such as heat shock proteins (HSPs), play central role in cytoprotection during heat stress by protecting the animal against hyperthermia, circulatory shock, and cerebral ischemia via the overexpression of HSPs [58]. The activation of these systems appears to be initiated at skin surface temperatures exceeding 35 °C as animals begin to store heat and rapidly increase evaporative heat loss mechanisms [59]. One potential cytoprotection system that can be present in CHCU is the activation of heat shock transcription factor 1 (*HSF1*). This gene coordinates thermal tolerance, and it was found in CHCU ROHs BTA14:1.464–14.362. In addition, we found that the heat shock protein-related gene *HSPB1* is located in ROH BTA25:34.446:42.851 in CHCU, and it encodes a small heat shock protein that functions as a molecular chaperone [60]. This ROH overlaps with another interesting region reported as a selective sweep by the authors of [22].

The importance of the *SLICK* locus in tropical adaptation has been previously shown [61]. A study (GWAS) of the Limonero animals using Illumina BovineHD genotypes from 20 *SLICK* and 53 non-*SLICK* individuals was made by the authors of [62]; they identified significant associations to markers in and around *PRLR* (prolactin receptor). Two SNP markers intragenic to *PRLR* presented
genome-wide significance, namely rs42551770 (20:39104658, *p* = 2.51 × 10⁻⁸bp) and rs137009256 (BTA 20:39110968, *p* = 1.67 × 10⁻⁸ bp). Later, the authors of [63] found PRLR genetic variants associated with both the undercoat and the topcoat that have the potential to be used for selection in order to breed for animals with shorter hair, thus leading to greater thermotolerance and potentially increasing production in hot and humid climates. Interestingly, we found PRLR annotated in ROHs BTA20:38.341–71.794 in CHCU.

ACHE is another gene potentially involved in thermotolerance because it has been found to be associated, for instance, with an acute heat stress response in chickens [64]. Interestingly, we also found *ITGA9*, which was found to be associated with heat tolerance in summer in a genome-wide association study for body temperature performed in US beef cattle [65]. This gene is located within ROH BTA22:4.113:12.47 and was also identified as a selective sweep by the authors of [22] on BTA22:10.477:11.159.

One study on the genetic diversity of 21 autochthonous cattle breeds from the whole Mediterranean basin performed genome-wide association analyses with covariables discriminating the different Mediterranean climate subtypes [66]; the authors of this work found various genes associated with climatic covariables that overlapped within the ROHs found in the current study: *PTPRF* (protein tyrosine phosphatase receptor Type F) on BTA3:90.368–117.379 in CHF; *NRG1* (neuregulin 1) on BTA27:15.885–28.268 in CHF, *OTUD7A* (OTU deubiquitinase 7A) on BTA21:2.75–53.659 in CHCU, *C5H12orf54* and *RACGAP1* (Rac GTPase activating protein 1) on BTA5:17.363–33.013, BTA5:6.285–36.84, and BTA5:21.708–33.687 in CHCU are associated with climatic covariables—particularly *ANTXR2* (anthrax toxin receptor 2) on BTA6:90.666–119.217 in CHCU, which is associated with both the annual mean temperature and THI (temperature humidity index), and *TRAM2* (translocation associated membrane protein 2) on BTA23:2.651–25.55, BTA23:19.069–25.553, and BTA23:20.348–25.553 in CHCU, which is associated with annual mean radiation and annual precipitation [66].

**5. Conclusions**

In this work, we estimated, for the first time, the genomic inbreeding value and effective size of the CHCU breed using a medium-sized SNP panel, and then we compared this information to that of two closely related European Charolais populations.

Genomic inbreeding levels were found to be higher in the CHCU breed compared to the CHF and CHUK populations. The genomic inbreeding patterns were in agreement with the results obtained with the estimates of the effective population size, with low estimated Ne values and a tendency to decrease along generations. However, the high variance found in the indicators of ROHs and FROH suggest an opportunity to maintain the genetic diversity of this population with an adequate mating strategy that, accompanied with the benefits of molecular tools, could guarantee a better management of this important genetic resource.

Finally, we found several ROHs in CHCU on BTA2, 5, 14, 25, 20, 22, 21, and 23 that might be associated with environmental adaptation traits.

**Supplementary Materials:** The following are available online at www.mdpi.com/2076-2615/10/12/2233/s1, Table S1: Statistic summary the FROH_BTA per chromosome in three Charolais cattle populations. Table S2: Distribution of the number and size of runs of homozygosity (ROH) by populations and chromosomes (BTA). Table S3: Value estimated of FROH in other cattle breeds.

**Author Contributions:** Conceptualization, Y.R.-V. and Y.R.-C.; methodology, Y.R.-V., Y.R.-C. and S.E.R.-O.; software, Y.R.-V.; validation, S.E.R.-O., Y.R.-C. and D.R.; formal analysis, Y.R.-V., Y.R.-C. and S.E.R.-O.; investigation, Y.R.-V.; resources, Y.R.C, D.R., S.E.R.-O. and M.N.; data curation, Y.R.-V. and Y.R.-C.; writing—original draft preparation, Y.R.-V.; writing—review and editing, S.E.R.-O., Y.R.-C., M.N., D.R., G.R. and E.P.-P.; visualization, Y.R.-V.; supervision, S.E.R.-O., Y.R.-C. and M.N.; project administration, Y.R.C, D.R and S.E.R.-O.; funding acquisition, Y.R.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by Ministerio de Economía y Competividad grants AGL2013-41834-R (MEC, Spain), AGL2016-78709-R (MEC, Spain) and by the CERCA Programme/Generalitat de Catalunya. We acknowledge financial support from the Spanish Ministry of Economy and Competitiveness, through the
"Severo Ochoa Programme for Centres of Excellence in R&D" 2016-2019 (SEV-2015-0533). Genotyping of the CHCU samples was funded by the Animal Genetics Division of INRAE, as part of the CHACUBA project. Y. Rodriguez-Valera was funded by the French Embassy in Havana (Cuba) and INRAE. Y. Ramayo-Caldas was funded by the European Union H2020 Research and Innovation programme under Marie Skłodowska-Curie grant (P-Sphere) agreement No 6655919.

**Acknowledgments:** We thank the Empresa Genética y Cría Manuel Fajardo, farmers of San José del Retiro as well as the Departamento de Genética Nacional del Ministerio de Agricultura de Cuba.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. López: D.; Miriam, R.; Díaz, J.; Menéndez, A. Introducción y desarrollo del Charolais en Cuba. In *En El Charolais Cubano*.; Editorial Científico-Técnica: La Habana, Cuba, 1977; pp. 13–17. (In Spanish).

2. Wright, S. Coefficients of Inbreeding and Relationship. *Am. Nat.* **1922**, *56*, 330–338, doi:10.1086/279872.

3. Gibson, J.; Morton, N.E.; Collins, A. Extended tracts of homozygosity in outbred human populations. *Hum. Mol. Genet.* **2006**, *15*, 789–795, doi:10.1093/hmg/ddi493.

4. McQuillan, R.; Leutenegger, A.L.; Abdel-Rahman, R.; Franklin, C.S.; Pericic, M.; Barac-Lauc, L.; Smolej-Narancic, N.; Janicijevic, B.; Polasek, O.; Tenesa, A.; et al. Runs of Homozygosity in European Populations. *Am. J. Hum. Genet.* **2008**, *83*, 359–372, doi:10.1016/j.ajhg.2008.08.007.

5. Ferenčaković, M.; Hamzić, E.; Gredler, B.; Solberg, T.R.; Klemetsdal, G.; Curik, I.; Sölkner, J. Estimates of autozygosity derived from runs of homozygosity: Empirical evidence from selected cattle populations. *J. Anim. Breed. Genet.* **2013**, *130*, 286–360 doi:10.1111/jbg.12012.

6. Silió, L.; Rodríguez, M.C.; Fernández, A.; Barragán, C.; Benítez, R.; Óvilo, C.; Fernández, A.I. Measuring inbreeding and inbreeding depression on pig growth from pedigree or SNP-derived metrics. *J. Anim. Breed. Genet.* **2013**, *130*, 349–360, doi:10.1111/jbg.12031.

7. Fernández, A.; Angel Toro, M.; López-Fanjul, C. The effect of inbreeding on the redistribution of genetic variance of fecundity and viability in *Tribolium castaneum*. *Heredity* **1995**, *75*, 376–381, doi:10.1038/hdy.1995.149.

8. Doekes, H.P.; Veerkamp, R.F.; Bijma, P.; De Jong, G.; Hiemstra, S.J.; Windig, J.J. Inbreeding depression due to recent and ancient inbreeding in Dutch Holstein-Friesian dairy cattle. *Genet. Sel. Evol.* **2019**, *51*, 54, doi:10.1186/s12711-019-0497-z.

9. Curik, I.; Ferenčakovic, M.; Sölkner, J. Genomic dissection of inbreeding depression: A gate to new opportunities. *Rev. Bras. Zootec.* **2017**, *46*, 773–782, doi:10.1590/S1806-92902017000900010.

10. Goszczynski, D.; Molina, A.; Terán, E.; Morales-Durand, H.; Ross, P.; Cheng, H.; Giovambattista, G.; Demyda-Peyrás, S. Runs of homozygosity in a selected cattle population with extremely inbred bulls: Descriptive and functional analyses revealed highly variable patterns. *PLoS ONE* **2018**, *13*, e0200069, doi:10.1371/journal.pone.0200069.

11. Nothnagel, M.; Tehua, T.L.; Kayser, M.; Krawczak, M. Genomic and geographic distribution of SNP-defined runs of homozygosity in Europeans. *Hum. Mol. Genet.* **2010**, *19*, 2927–2935, doi:10.1093/hmg/ddq198.

12. Flury, C.; Tapio, M.; Sonstegard, T.; Drögemüller, C.; Leeb, T.; Simianer, H.; Hanotte, O.; Rieder, S. Effective population size of an indigenous Swiss cattle breed estimated from linkage disequilibrium. *J. Anim. Breed. Genet.* **2010**, *127*, 339–347, doi:10.1111/j.1439-0388.2010.00862.x.

13. Corbin, L.J.; Liu, A.Y.H.; Bishop, S.C.; Woolliams, J.A. Estimation of historical effective population size using linkage disequilibria with marker data. *J. Anim. Breed. Genet.* **2012**, *129*, 257–270, doi:10.1111/j.1439-0388.2012.01003.x.

14. Barbato, M.; Orozco-terWengel, P.; Tapio, M.; Bruford, M.W. SNeP: A tool to estimate trends in recent effective population size trajectories using genome-wide SNP data. *Front. Genet.* **2015**, *6*, 109, doi:10.3389/fgene.2015.00109.

15. Hui, T.J. Estimating Effective Population Size from Genetic Data: The Past, Present, and the Future. Ph.D. Thesis, Imperial College London, London, UK, January 2017.

16. Xu, L.; Zhu, B.; Wang, Z.; Xu, L.; Liu, Y.; Chen, Y.; Zhang, L.; Gao, X.; Gao, H.; Zhang, S.; et al. Evaluation of linkage disequilibrium, effective population size and haplotype block structure in Chinese cattle. *Animals* **2019**, *9*, 83, doi:10.3390/ani9030083.
17. Falconer, D.S.; Mackay, T.F.C. Effective populations size. In *Introduction into Quantitative Genetics*; 3rd ed.; Longman Scientific & Technical: Harlow, UK, 1989; pp. 70–76.

18. Saura, M.; Tenesa, A.; Woolliams, J.A.; Fernández, A.; Villanueva, B. Evaluation of the linkage-disequilibrium method for the estimation of effective population size when generations overlap: An empirical case. *BMC Genom.* **2015**, *16*, 922, doi:10.1186/s12864-015-2167-z.

19. Wang, J.; Santiago, E.; Caballero, A. Prediction and estimation of effective population size. *Heredity* **2016**, *117*, 193–206, doi:10.1038/hdy.2016.43.

20. Sved, J.A. Linkage disequilibrium and homozygosity of chromosome segments in finite populations. *Theor. Popul. Biol.* **1971**, *2*, 125–141, doi:10.1016/0040-5809(71)90011-6.

21. Hill, W.G. Estimation of effective population size from data on linkage disequilibrium. *Genet. Res.* **1981**, *38*, 209–216.

22. Rodriguez-Valera, Y.; Renand, G.; Naves, M.; Fonseca-Jiménez, Y.; Moreno-Probance, T.I.; Ramos-Onsins, S.; Rocha, D.; Ramayo-Caldas, Y. Genetic diversity and selection signatures of the beef "Charolais de Cuba" breed. *Sci. Rep.* **2018**, *8*, 1–9, doi:10.1038/s41598-018-29453-z.

23. Matukumalli, L.K.; Lawley, C.T.; Schnabel, R.D.; Taylor, J.F.; Allan, M.F.; Heaton, M.P.; O'Connell, J.; Moore, S.S.; Smith, T.P.L.; Sonstegard, T.S.; et al. Development and Characterization of a High Density SNP Genotyping Assay for Cattle. *PLoS ONE* **2009**, *4*, e5350, doi:10.1371/journal.pone.0005350.

24. Gautier, M.; Laloe, D.; Moazami-Goudarzi, K. Insights into the Genetic History of French Cattle from Dense SNP Data on 47 Worldwide Breeds. *PLoS ONE* **2010**, *5*, e13038, doi:10.1371/journal.pone.0013038.

25. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.R.; Bender, D.; Maller, J.; Sklar, P.; de Bakker, P.I.W.; Daly, M.J.; et al. Plink: A toolset for whole-genome association and population-based linkage analysis. *Am. J. Hum. Genet.* **2007**, *81*, 559–575, doi:10.1086/519795.

26. Mastrangelo, S.; Tolone, M.; Di Gerlando, R.; Fontanesi, L.; Sardina, M.T.; Portolano, B. Genomic inbreeding estimation in small populations: Evaluation of runs of homozygosity in three local dairy cattle breeds. *Animal* **2016**, *10*, 746–754, doi:10.1017/S1751731115002943.

27. Ben, J.S.; Boussaha, M.; Thamri, N.; Mnara, S.; Rebours, E.; Rocha, D. Linkage disequilibrium and past effective population size in native Tunisian cattle. *Genet. Mol. Biol. Genet.* **2019**, *42*, 52–61, doi:10.1590/1678-4685-GMB-2017-0342.

28. Biscarini, F.; Cozzi, P.; Gaspa, G.; Marras, G. detectRUNS: An R Package to Detect Runs of Homozygosity and Heterozygosity in Diploid Genomes. Sample Data 2018. Available online: https://cran.r-project.org/web/packages/detectRUNS/index.html (accessed on 25 January 2020).

29. Ferenčaković, M.; Sölkner, J.; Curik, I. Estimating autozygosity from high-throughput information: Effects of SNP density and genotyping errors. *Genet. Sel. Evol.* **2013**, *45*, 1–9, doi:10.1186/1297-9686-45-42.

30. Marras, G.; Gaspa, G.; Sorbolini, S.; Dimauro, C.; Ajmone-Marsan, P.; Valentini, A.; Williams, J.L.; MacCiotta, N.P.P. Analysis of runs of homozygosity and their relationship with inbreeding in five cattle breeds farmed in Italy. *Anim. Genet.* **2015**, *46*, 110–121, doi:10.1111/age.12259.

31. Hayes, B.J.; Visscher, P.M.; McPartlan, H.C.; Goddard, M.E. Novel multilocus measure of linkage disequilibrium to estimate past effective population size. *Genome Res.* **2003**, *13*, 635–643, doi:10.1101/gr.387103.

32. Ohta, T.; Kimura, M. Linkage disequilibrium between two segregating nucleotide sites under the steady flux of mutations in a finite population. *Genetics* **1971**, *68*, 571–580.

33. Tenesa, A.; Navarro, P.; Hayes, B.J.; Duffy, D.L.; Clarke, G.M.; Goddard, M.E.; Visscher, P.M. Recent human effective population size estimated from linkage disequilibrium. *Genome Res.* **2007**, *17*, 520–526, doi:10.1101/gr.6023607.8.

34. Purfield, C.D.; Berry, P.D.; McParland, S.; Bradley, G.D.; Burgess, H.A.; Martinez, S.; Reiner, O. Runs of homozygosity and population history in cattle. *BMC Genet.* **2012**, *58*, 70, doi:10.1186/1471-2156-13-70.

35. Menéndez-Buxadera, A. Variación genética de la composición de la canal en la raza Charolais. Ph.D. Thesis, Instituto Superior de Ciencias Agropecuarias de La Habana, La Habana, Cuba, November 1978; pp. 2–6. (In Spanish)

36. Felius, M.; Beerling, M.L.; Buchanan, D.S.; Theunissen, B.; Koolmees, P.A.; Lenstra, J.A. On the history of cattle genetic resources. *Diversity* **2014**, *6*, 705–750, doi:10.3390/d6040705.

37. Mi, H.; Thomas, P. Panther pathway: An ontology-based pathway database coupled with data analysis tools. *Methods Mol. Biol.* **2009**, *563*, 123–140.
38. Watterson, G. On the number of segregating sites in genetical models without recombination. *Theor Popul Biol* **1975**, *7*, 256.

39. Weigel, K.A. Controlling Inbreeding in Modern Breeding Programs. *J. Dairy Sci.* **2001**, *84*, E177-E184, doi:10.3168/jds.s0022-0302(01)70213-5.

40. Bjelland, D.W.; Weigel, K.A.; Vukasinovic, N.; Nkrumah, J.D. Evaluation of inbreeding depression in Holstein cattle using whole-genome SNP markers and alternative measures of genomic inbreeding. *J. Dairy Sci.* **2013**, *96*, 4697-4706, doi:10.3168/jds.2012-6435.

41. François, L.; Wijnrocx, K.; Colinet, F.G.; Gengler, N.; Hulsegge, B.; Windig, J.J.; Buys, N.; Janssens, S. Genomics of a revived breed: Case study of the Belgian campine cattle. *PLoS ONE* **2017**, *12*, e0175916, doi:10.1371/journal.pone.0175916.

42. Msalya, G.; Kim, E.S.; Laisser, E.L.K.K.; Kipanyula, M.J.; Karimuribo, E.D.; Kusiluka, L.J.M.M.; Chenyambuga, S.W.; Rothschild, M.F. Determination of genetic structure and signatures of selection in three strains of Tanzania Shorthorn Zebu, Boran and Friesian cattle by genome-wide SNP analyses. *PLoS ONE* **2017**, *12*, e171088, doi:10.1371/journal.pone.0171088.

43. Solé, M.; Gori, A.; Faux, P.; Bertrand, A.; Farnir, F.; Gautier, M.; Druet, T. Age-based partitioning of individual genomic inbreeding levels in Belgian Blue cattle. *Genet. Sel. Evol.* **2017**, *49*, 92, doi:10.1186/s12711-017-0370-x.

44. Cortés, O.; Eusebi, P.; Dunner, S.; Sevane, N.; Cañón, J. Comparison of diversity parameters from SNP, microsatellites and pedigree records in the Lidia cattle breed. *Livestig. Sci.* **2019**, *21*, 80-85, doi:10.1016/j.livsci.2018.11.006.

45. Toro-Ospina, M.A.; da Silva-Faria, R.A.; Vercesi-Filho, A.E.; dos Santos-Gonçalves, J.N.C.; Zerlotti-Mercadante, M.E.; Abdallah-Curi, R.; de Vasconcelos-Silva, J.A.I. Genome-wide identification of runs of homozygosity islands in the Gyr breed (Bos indicus). *Reprod. Domest. Anim.* **2020**, *55*, 333-342, doi:10.1111/rda.13639.

46. Toro-Ospina, A.M.; Maiorano, A.M.; Curi, R.A.; Pereira, G.L.; Zerlotti-Mercadante, M.E.; dos Santos Gonçalves Cyrillo, J.N.; Aspilcueta-Borquis, R.R.; Josineudson, J.A.I.I. Linkage disequilibrium and effective population size in Gir cattle selected for yearling weight. *Reprod. Domest. Anim.* **2019**, *54*, 1524-1531, doi:10.1111/rda.13559.

47. Šidlová, V.; Kasarda, R.; Moravčíková, N.; Trakovická, A.; Curik, I.; Ferenčaković, M. Genomic variability among cattle populations based on runs of homozygosity. *Poljoprivreda* **2015**, *21*, 44-47, doi:10.18047/poljo.21.1.sup.9.

48. Kasarda, R.; Kadlečík, O.; Trakovická, A.; Moravčíková, N.; Biology, B. Genomic and pedigree-based inbreeding in Slovak Spotted Cattle. *AGROFOR Int. J.* **2019**, *4*, 102-110, doi:10.7251/AGRENG1901102K.

49. Szmatoła, T.; Gurgul, A.; Ropka-Molik, K.; Jasielczuk, I.; Zabek, T.; Bugno-Poniewierska, M. Characteristics of runs of homozygosity in selected cattle breeds maintained in Poland. *Livestig. Sci.* **2016**, *188*, 72-80, doi:10.1016/j.livsci.2016.04.006.

50. Lu, D.; Sargolzaei, M.; Kelly, M.; Li, C.; Voort, G.V.; Wang, Z.; Plastow, G.; Moore, S.; Miller, S.P. Linkage disequilibrium in Angus, Charolais, and Crossbred beef cattle. *Front. Genet.* **2012**, *3*, 1-10, doi:10.3389/fgene.2012.00152.

51. Kim, S.; Cheong, H.S.; Shin, H.D.; Lee, S.S.; Roh, H.J.; Jeon, D.Y.; Cho, C.Y. Genetic diversity and divergence among Korean cattle breeds assessed using a BovineHD single-nucleotide polymorphism chip. *Asian-Australas. J. Anim. Sci.* **2018**, *31*, 1691-1699, doi:10.5713/ajas.17.0419.

52. Aguilar, R.; Bu, A.; Dresdner, J.; Fernández, P.; González, A.; Polanco, C.; Tansini, R. Capítulo 1. Antecedentes de la ganadería vacuna en Cuba. In *La ganadería en Cuba: Desempeño y Desafíos*.; ASDI: INIE: Udelar. FCS-DE: Montevideo, Uruguay, 2004; p. 15; ISBN 959-7160-16-1.

53. Bouquet, A.; Venot, E.; Laloë, D.; Forabosco, F.; Fogh, A.; Pabiou, T.; Moore, K.; Eriksson, J.; Renand, G.; Phocas, F. Genetic structure of the European Charolais and Limousin cattle metapopulations using pedigree analyses. *J. Anim. Sci.* **2011**, *89*, 1719-1730, doi:10.2527/jas.2010-3469.

54. Ribas, M. Gene frequencies in the blood group systems of the Cuban Charolais. *Ann. Génét. Sél. Anim.* **1981**, *13*, 293-300.

55. FAO. *Secondary Guidelines for Development of National Farm Animal Genetic Resources Management Plans: Management of Small Populations at Risk*; FAO: Rome, Italy, 1998. Available online: http://www.fao.org/3/a-w9361e.pdf (accessed on 20 June 2019).
56. Pavlidis, P.; Jensen, J.D.; Stephan, W.; Stamatakis, A. A critical assessment of storytelling: Gene ontology categories and the importance of validating genomic scans. *Mol. Biol. Evol.* **2012**, *29*, 3237–3248, doi:10.1093/molbev/mss136.

57. Direito, I.; Madeira, A.; Brito, M.A.; Soveral, G. Aquaporin-5: From structure to function and dysfunction in cancer. *Cell. Mol. Life Sci.* **2016**, *73*, 1623–1640, doi:10.1007/s00018-016-2142-0.

58. Rout, M.; Panigrahi, S.; Pradhan, S.; Swain, K. Genetic basis of heat tolerance in cattle. *Pharma Innov. J.* **2018**, *7*, 183–185.

59. Collier, R.J.; Collier, J.L.; Rhoads, R.P.; Baumgard, L.H. Invited review: Genes involved in the bovine heat stress response. *J. Dairy Sci.* **2008**, *91*, 445–454, doi:10.3168/jds.2007-0540.

60. Mymrikov, E.V.; Seit-Nebi, A.S.; Gusev, N.B. Large Potentials of Small Heat Shock Proteins. *Physiol. Rev.* **2011**, *91*, 1123–1159, doi:10.1152/physrev.00023.2010.

61. Flori, L.; Gonzatti, M.I.; Thevenon, S.; Chantal, I.; Pinto, J.; Berthier, D.; Aso, P.M.; Gautier, M. A quasi-exclusive European ancestry in the Senepol tropical cattle breed highlights the importance of the slick locus in tropical adaptation. *PLoS ONE* **2012**, *7*, e36133, doi:10.1371/journal.pone.0036133.

62. Porto-Neto, L.R.; Bickhart, D.M.; Landaeta-hernandez, A.J.; Utsunomiya, Y.T.; Pagan, M.; Jimenez, E.; Hansen, P.J.; Dikmen, S.; Schroeder, S.G.; Kim, E.; et al. Convergent Evolution of Slick Coat in Cattle through Truncation Mutations in the Prolactin Receptor. *Front. Genet.* **2018**, *9*, 1–8, doi:10.3389/fgene.2018.00057.

63. Davila, K.M.S.; Howell, A.; Nunez, A.; Orelien, A.; Roe, V.; Rodriguez, E.; Dikmen, S.; Mateescu, R.G. Genome-wide association study identifies variants associated with hair length in Brangus cattle. *Anim. Genet.* **2020**, *1–4*, doi:10.1111/age.12970.

64. Cheng, C.-Y.; Tu, W.-L.; Chen, C.-J.; Chan, H.-L.; Chen, C.-F.; Chen, H.-H.; Tang, P.-C.; Lee, Y.-P.; Chen, S.-E.; Huang, S.-Y. Functional genomics study of acute heat stress response in the small yellow follicles of layer-type chickens. *Sci. Rep.* **2018**, *8*, 1320, doi:10.1038/s41598-017-18335-5.

65. Howard, J.T.; Kachman, S.D.; Snelling, W.M.; Pollak, E.J.; Ciobanu, D.C.; Kuehn, L.A.; Spangler, M.L. Beef cattle body temperature during climatic stress: A genome-wide association study. *Int. J. Biometeorol.* **2014**, *58*, 1665–1672, doi:10.1007/s00484-013-0773-5.

66. Flori, L.; Moazami, K.; Véronique, G.; Araba, A.; Boujenane, I.; Boushaba, N.; Casabianca, F.; Casu, S.; Ciampolini, R.; Coeur, A.; et al. A genomic map of climate adaptation in Mediterranean cattle breeds. *Mol. Ecol.* **2019**, *28*, 1009–1029, doi:10.1111/mec.15004.

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_13.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).