Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 https://doi.org/10.1186/s12711-020-00597-9 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Whole‑genome sequencing reveals insights into the adaptation of French Charolais cattle to Cuban tropical conditions Lino C. Ramírez‑Ayala1, Dominique Rocha2, Sebas E. Ramos‑Onsins1, Jordi Leno‑Colorado1, Mathieu Charles2,3, Olivier Bouchez4, Yoel Rodríguez‑Valera5, Miguel Pérez‑Enciso1,6† and Yuliaxis Ramayo‑Caldas2,7*† Abstract Background: In the early 20th century, Cuban farmers imported Charolais cattle (CHFR) directly from France. These animals are now known as Chacuba (CHCU) and have become adapted to the rough environmental tropical condi‑ tions in Cuba. These conditions include long periods of drought and food shortage with extreme temperatures that European taurine cattle have difficulty coping with. Results: In this study, we used whole-genome sequence data from 12 CHCU individuals together with 60 whole- genome sequences from six additional taurine, indicus and crossed breeds to estimate the genetic diversity, structure and accurate ancestral origin of the CHCU animals. Although CHCU animals are assumed to form a closed population, the results of our admixture analysis indicate a limited introgression of Bos indicus. We used the extended haplotype homozygosity (EHH) approach to identify regions in the genome that may have had an important role in the adapta‑ tion of CHCU to tropical conditions. Putative selection events occurred in genomic regions with a high proportion of Bos indicus, but they were not sufficient to explain adaptation of CHCU to tropical conditions by Bos indicus introgres‑ sion only. EHH suggested signals of potential adaptation in genomic windows that include genes of taurine origin involved in thermogenesis (ATP9A, GABBR1, PGR, PTPN1 and UCP1) and hair development (CCHCR1 and CDSN). Within these genes, we identified single nucleotide polymorphisms (SNPs) that may have a functional impact and contribute to some of the observed phenotypic differences between CHCU and CHFR animals. Conclusions: Whole-genome data confirm that CHCU cattle are closely related to Charolais from France (CHFR) and Canada, but also reveal a limited introgression of Bos indicus genes in CHCU. We observed possible signals of recent adaptation to tropical conditions between CHCU and CHFR founder populations, which were largely independent of the Bos indicus introgression. Finally, we report candidate genes and variants that may have a functional impact and explain some of the phenotypic differences observed between CHCU and CHFR cattle. Background Climate change and global warming are among the main challenges currently faced by Agriculture and Livestock husbandry. In this scenario, it is fundamen- *Correspondence: yuliaxis@gmail.com tal to investigate the mechanisms that allow animals † Miguel Pérez-Enciso and Yuliaxis Ramayo-Caldas were equally to adapt to high temperatures. The hot temperature responsible for this work 2 conditions of tropical climates today might resemble Université Paris-Saclay, INRAE, Jouy‑En‑Josas, AgroParisTech, GABI 78350, France those that animals raised outdoors under temperate Full list of author information is available at the end of the article climates, such as most beef cattle and small ruminants, © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat​iveco​ mmons​.org/licen​ses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creat​iveco​mmons​.org/publi​cdoma​in/ zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 2 of 11 will face in the future. Therefore, animals that currently changes produced by genetic drift from those due to live in hot climates and that have a European or tem- adaptation (e.g., [9–11]). perate climate origin can provide clues into the genetic Because SNP array genotype data are biased and have a mechanisms underlying adaptation to increasing tem- low resolution, for our study, we obtained whole-genome peratures [1]. sequence data from 72 animals, including 12 CHCU ani- Cattle breeds can be divided into temperate taurine mals, to provide an unbiased estimation of the popula- breeds Bos taurus, of European origin, and Indian zebu tion structure and to fine-map regions that could have breeds Bos indicus, which diverged ~ 250,000 years ago played a role in the adaptation to tropical conditions of [2]. B. indicus breeds produce less meat and of lower the Chacuba population. quality but are more adapted to heat and parasites than taurine breeds. For that reason, they were imported Methods to tropical American regions starting in the mid-19th Samples century [3, 4]. Under these tropical climates, they Seventy-two whole-genome sequences from taurine, have largely replaced the primigenious cattle that were indicine and crossbred cattle were analyzed. CHCU cat- imported by the first Spanish and Portuguese settlers tle have been maintained under pedigree control at the [5]. Numerous hybrid populations between B. taurus “Manuel Fajardo” genetic center that is located in Jiguani and B. indicus, such as the Brangus, Texas Longhorn, (Granma Province) with a current census of ~ 700 ani- Santa Gertrudis, among others, also coexist with pure mals. Therefore, the genetic relationship between animals B. indicus breeds. These mixed breeds exhibit a good can be accurately tracked and we used this information to resistance to parasites and heat, and also produce car- select 12 unrelated CHCU animals that were sequenced casses of much higher quality than pure B. indicus in this work. We also used 15 French Charolais (CHFR) breeds. [12], six Limousine (LIMS) from France [13] and At the beginning of the 20th century, Cuban farmers sequences from 39 additional individuals that were imported Charolais animals from France, which have downloaded from the sequence read archive (SRA) resulted in a population of cattle now known as ‘Cha- database [see Additional file 1: Table S1]: 15 Canadian cuba’ (CHCU). This population has adapted to the breed- Charolais (CHCA), five Limousine from Canada (LIMS), ing conditions in the Cuban tropical environment in ~ 20 five Brangus (BRG), 10 Brahman (BRM) and four Texas generations of breeding. In this short period of time, Longhorn (TXL). Brahman is a pure B. indicus breed several clear phenotypic differences between the origi- whereas BRG and TXL are admixed breeds between nal French Charolais (CHFR) and its Cuban counterpart B. indicus and B. taurus cattle. Data on the French and have appeared. CHCU cattle are smaller than CHFR, Canadian LIMS individuals were merged in the analyses with weights of 34 vs. 46 kg at birth and 290 vs. 493 kg reported here. for heifer’s weight at 18 months, respectively [6, 7]. Also, CHCU cattle are hairless and their carcasses have a lower Bioinformatic analysis grade and higher fat content than CHFR [7]. Although All the sequences were mapped against the bovine refer- the CHCU cattle are thought to form a closed population ence assembly (UMD3.1.1) using the Burrows-Wheeler with no records of interbreeding, Ribas [8] reported the aligner (BWA) v. 0.7.12-r1039 software [14]. PCR dupli- presence of a specific B. indicus blood group allele (U’1), cates were removed using the Picard MarkDuplicates at a very low frequency. More recently, Rodriguez-Valera (v2.18.9) program and realigned around InDels with et al. [6] used the Illumina Bovine 50 k single nucleotide the GATK IndelRealigner tool [15]. For each individual, polymorphism (SNP) BeadChip to investigate the genetic SNP calling was done with the SAMtools mpileup and structure and putative ancestral origin of this popula- bcftools call (v. 0.1.19-96b5f2294a) tools with the follow- tion and showed that CHCU clusters with the taurine ing parameters: minimum and maximum depths between breeds. Nevertheless, in spite of the short period of time 5 × and twice the average sample’s depth; a minimum since the importation of Charolais animals from France, SNP quality of 10; and a minimum mapping quality and a marked differentiation (fixation index (FST) = 0.049) is minimum base quality of 20. Next, we merged individ- observed between CHFR and CHCU cattle. Genetic and ual gVCF files into a multi-individual VCF file, with all phenotypic differentiation can be caused by genetic drift, the SNPs from the 72 samples. For this purpose, we fol- i.e., as a consequence of the importation of a small num- lowed a two-step approach as detailed in [16], using a ber of individuals and by selection of beneficial variants pipeline available at https​://githu​b.com/migue​lpere​zenci​ that promote adaptation to the tropical environment. A so/NGSpi​pelin​e. In brief, to identify whether a position number of statistics have been proposed to distinguish is equal to that in the reference genome, polymorphic or missing, first we generated a fasta file from the gVCF file Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 3 of 11 for each individual and generated a multi-individual VCF ascertaining B. indicus introgression in CHCU. Never- file by using the individual file. Once the multiple sample theless, K = 2 was also the value that resulted in the low- file was obtained, SNPs with more than 20% missing data est cross-validation error. The program was run either across samples and populations were removed. Finally, by including all the SNPs or with pruned data from we imputed the missing genotypes and inferred phases which SNPs in strong disequilibrium were removed, but with the Beagle 4.1 software [17]. the results were identical. To get a more precise map of Genetic variants that alter transcription factor binding potential admixture, we used the ELAI software [27], sites (TFBS) were predicted with a custom script using which is a partially supervised algorithm that requires the TFBS models from the JASPAR (JASPAR CORE 2018 data from the putative founder populations (CHFR and collection, [18]), HOCOMOCO (version v10, [19] and BRM) and the potentially admixed population (CHCU). TRANSFAC (version v3.2 public, [20] databases. These Therefore, ELAI was run using CHCU, CHFR and BRM databases contain a curated set of TFBS models repre- genotypes only. ELAI implements a two-layer hidden sented as position weight matrices (PWM), which are Markov model and was run by setting the recommended derived from published collections of experimentally default parameters, which included removal of SNPs with defined eukaryotic TFBS. Only vertebrate PWM were a minimum allele frequency (MAF) lower than 0.01. downloaded and used in our study. Finally, we identified microRNA binding sites by using the TargetScan (release Results 7.2) software [21]. Population structure and impact of B. indicus introgression on CHCU Population genomics Average read depth across breeds varied between 8.7 We estimated Watterson’s nucleotide variability (θ) [22] (LIMS) and 12.5 (CHCU) (see Table 1). We found and differentiation values (FST, [23] between popula- 42,144,809 SNPs among which 14,929,949 were specific tions with the mstatspop (v.0.1beta, https​://githu​b.com/ to the pure B. indicus breed (BRM), 6,839,436 to the tau- CRAGE​NOMIC​A/mstat​spop) software in consecutive rine breeds and 1,176,249 to CHCU. As expected [28], B. non-overlapping 30-kb windows. This software imple- indicus samples were more variable than B. taurus sam- ments algorithms that allow for missing data [24]. For the ples, i.e. the nucleotide variability was equal to 0.0035 remaining analyses, we imputed missing genotypes with per bp for BRM but was two times lower for LIMS and Beagle 4.1 [17]. The EHH-derived statistics (Rsb, and iHs) European Charolais (Table 1). The hybrid breeds TXL [25] were computed between CHCU and CHFR for each and BRG had an intermediate level of genome diversity, SNP (https​://githu​b.com/CRAGE​NOMIC​A/Tang_Rsb). i.e. between those of the B. indicus and B. taurus breeds, As putative selection events, we retained the windows and CHCU had a level of genome diversity similar to with a ‘permutation p-value’ lower than 0.05 among the that of TXL (Table 1). Thus, it is interesting to note that 2000 windows that contained the SNPs with the largest although the number of founders of the CHCU popula- Rsb value. The ‘permutation p-value’ was obtained by tion is small, it has a higher nucleotide variability than its randomly shuffling the CHCU and CHFR samples and ancestral CHFR population. running the Rsb algorithm along the whole genome. The In line with previous observations [29, 30], the prin- process was repeated 100 times and we computed the cipal component analysis (PCA) plot shows a clear number of times the observed Rsb statistics was larger separation between indicine and taurine breeds, with than the values obtained from permutation, and obtained the latter breeds being tightly clustered (Fig. 1a). As a ‘permutation p-value’ for each SNP. This process aims expected, animals sampled from the admixed breeds at correcting for different levels of disequilibrium along BRG and TXL, are positioned towards the cluster of the genome that may locally inflate Rsb values. In a sec- indicine breeds but much closer to the taurine than to ond step, we performed the same procedure but this time the indicine clusters, because their proportion of indi- between CHFR and CHCA, to exclude the Rsb intervals cine genome is less than 50% or even 11% in the case of that were common between CHFR and CHCA and not TXL breed [31]. Regarding the CHCU individuals, they specific to CHCU. Therefore, we focused only on the Rsb are positioned near the LIMS and Canadian Charolais, intervals with signals exclusively from CHCU. but separated from the original French Charolais. Since one CHCU individual appeared to be an out- Admixture lier, we inspected its genotype heterozygosity patterns The software ADMIXTURE v1.3.0 [26] was run in an (11%), but we could not find any anomalous deviation. unsupervised manner with a number of clusters K = 2 Moreover, this individual does not appear as an outlier using CHCU, BRM, CHFR and CHCA genotypes. when only BRM, CHCU and CHFR are represented We chose K = 2 because we were only interested in in a separate PCA plot (Fig. 1b). In terms of FST, the Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Table 1 Description of the breeds analyzed, shared and breed-exclusive SNPs Breed Country Number Wattersons % missing Mean depth CHCU-exclusive SNPs Shared SNPs Shared SNPs Shared SNPs Shared SNPs Shared Shared of samples variability with CHCU with CHFR with CHCA with LIMS SNPs with SNPs with TXL BRG CHCU CU 12 0.0021 0.03 12.5 1,176,249 – CHFR FR 15 0.0015 0.06 11.4 582,972 10,110,342 – CHCA CA 15 0.0018 0.13 10.4 737,160 10,667,212 10,847,599 – LIMS FR/CA 11 0.0016 0.21 8.7 695,594 9,429,688 9,796,472 9,961,440 – TXL US 4 0.0023 0.08 10.3 576,849 8,243,412 7,639,145 8,085,077 7,370,881 – BRG US 5 0.0028 0.05 11.5 842,771 9,514,147 6,854,966 9,033,938 7,934,515 7,655,123 – BRM US/AU 10 0.0035 0.11 10.8 14,929,949 11,442,083 8,653,635 10,134,007 8,267,282 9,169,187 12,519,875 Page 4 of 11 Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 5 of 11 a b BRG BRM 0.00 BRM CHCU 0.10 CHCA CHFR CHCU −0.10 PC2= 0.17 PC2= 0.24 CHFR 0.00 LIMS TXL −0.20 −0.10 −0.30 −0.20 −0.13 −0.11 −0.09 −0.07 −0.25 −0.20 −0.15 −0.10 PC1= 0.34 PC1= 0.31 c d 1.0 CHCU − BRG 15 CHCU − BRM 0.8 CHCU − CHCA CHCU − CHFR 10 0.6 Ancestry CHCU − LIMS Density CHCU − TXL 0.4 5 0.2 0.0 0 −0.2 0.0 0.2 0.4 0.6 0.8 BRM CHCU CHCA CHFR Fst Samples Fig. 1 Population structure. (a) Principal component analysis using all samples. Individuals are grouped into Bos taurus, Bos indicus and Hybrid clusters. Black: Brangus, red: Brahman, blue: Canadian Charolais, cyan: Cuban Charolais, green: French Charolais, magenta: Limousin, purple: Texas Longhorn (b) Principal component analysis using all samples. Red: Brahman, cyan: Cuban Charolais, green: French Charolais (c) FST between CHCU and the other breeds. (d) Results of admixture analyses with two ancestral populations (K = 2) Red: Brahman (Bos indicus) and green: French and Canadian Charolais (Bos taurus) closest populations to CHCU were the Canadian (0.04) and CHCU breeds are represented in the PCA plot and French Charolais (0.05) (Table 1). The plot of FST (Fig. 1b). values across 1-Mb windows between CHCU and the To evaluate more precisely the extent and impact of B. taurine breeds shows a modal value near zero, whereas indicus introgression into CHCU, we ran the ELAI soft- that between CHCU and the indicine breeds shows that ware [27], which provides a map, for each hybrid indi- they are clearly distinct (Fig. 1c). vidual, showing the probability for each SNP to descend In previous work based on the Illumina 50 k SNP Chip from one of the two putative founder populations. The genotype data, we reported a putative introgression of B. results show that the B. indicus introgression is not indicus into CHCU and although it was small [6], it could homogeneous, neither across individuals nor across explain why CHCU and CHFR are clearly separated, in chromosomes (Fig. 2) and [see Additional file 2: Fig- spite of the short period since the original importation of ure S1]. For instance, on average, B. taurus (BTA) auto- Charolais from France. Here, we confirm this introgres- somes BTA12, 13 and 23 had the highest proportions of sion with a better resolution using whole-genome data. B. indicus introgression (i.e. 24, 19 and 18%, respectively) As shown in Fig. 1d, an unsupervised Admixture analysis across all animals. On a per individual basis, the high- with K = 2 clusters clearly separates the taurine CHCA est proportion of B. indicus introgression was found for and CHFR animals from the pure B. indicus BRM breed, the CHCU11 individual, i.e. 47% on BTA12, followed whereas it reveals a small introgression of B. indicus in by CHCU10, i.e. 45% on BTA13. In contrast, on average CHCU. This is also reflected when only the BRM, CHFR the lowest proportion of B. indicus introgression was Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 6 of 11 Fig. 2 Frequency of indicus component per CHCU sample and chromosome inferred from analysis with the ELAI software found for individuals CHCU12 (6.5%), CHCU4 (7.4%) and CHCU2 (7.5%) and for chromosomes BTA5 (2.4%), BTA15 (2.7%) and BTA28 (3.9%). Additional file 3: Figure S2 [see Additional file 3: Figure S2] shows the correlation between the proportion of B. indicus genome in CHCU and the FST value between CHCU and CHFR across win- dows. A higher proportion of B. indicus was moderately correlated (0.19) with a greater differentiation between the Charolais populations. Putative selective sweeps between Cuban and French Charolais Among the 2000 windows with the largest positive Rsb value, we selected the 407 windows with a ‘permutation p-value’ lower than 0.05 and exclusive of CHCU (Fig. 3), i.e., we excluded the significant intervals observed in the comparison between CHFR and CHCA, since introgres- sion and adaptation to a hot climate are not expected to have played a role in any of the founder populations Fig. 3 Extended haplotype homozygosity (EHH) intervals per chromosome. The x-axis represents the chromosome positions and [see Additional file 4: Table S2]. In addition, we focused the y-axis the chromosome number on positive Rsb values only, since we are interested in regions where the putative selective pressure is specific to the Cuban population, i.e., where the disequilibrium is larger in CHCU than in CHFR. Figure 4 shows the genome. The average percentages of Brahman introgres- distribution of the percentage of indicus introgression sion were 0.22 and 0.17% in the selected and control for the putative selective windows and the rest of the Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 7 of 11 within the GABBR1 (gamma-aminobutyric acid B recep- 0.5 Selected regions tor 1) gene, which, if inactivated, induces hypothermia in Rest of the genome mice [33], 19 of these SNPs are within the 3′UTR and two SNPs (rs110080552 and rs210443447) have the poten- 0.4 tial to impact splicing. In addition, 227 of the 229 SNPs located in the upstream region of GABBR1 can poten- tially modify the binding site of transcription factors (see 0.3 Frequency Additional file 5: Table S3). Interestingly, nine of these potential regulatory SNPs are fixed for the alternative allele in CHCU. We found 1182 SNPs in the PGR gene, 0.2 which encodes the progesterone receptor and is associ- ated with thermotolerance in cattle [34]. Among these 1182 SNPs, two can impact splicing and 16 are located 0.1 in the upstream region and can potentially modify TFBS (see Additional file 5: Table S3). 0.0 Within the PTPN1 (protein-tyrosine phosphatase 1B) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 gene, we detected 382 SNPs including three located in Proportion BRM the 3′UTR and seven located in the upstream region Fig. 4 Percentage of BRM in windows containing putative selection events vs. the other windows, inferred from analysis with the ELAI that can potentially alter TFBS. Inactivation of PTPN1 software in mice results in an increase in adaptative thermogen- esis [35] and therefore some of these candidate regula- tory SNPs might explain the difference in heat tolerance between CHFR and CHCU. Finally, among the candidate windows, respectively. The difference, although small, genes involved in thermogenesis, we found two deleteri- was significant (P < 5e-09) according to a Wilcoxon rank ous missense variations (rs443726914 and rs715309385) test. within the UCP1 gene and 77 upstream SNPs, includ- We identified 243 genes within the 407 selected win- ing four that are fixed (rs438305189, rs211174809, dows, including five (ATP9A, GABBR1, PGR, PTPN1 and rs209939359 and rs211622720) for the alternative alleles UCP1) involved in thermotolerance and two (CCHCR1 in CHCU (see Additional file 5: Table S3). Interestingly, and CDSN) involved in hair development (Table 2). The these four upstream SNPs can potentially alter the bind- use of whole-genome data, allowed us to also pinpoint ing sites of several transcription factors. genetic variants within these genes that may have a func- Regarding hair development, 133 SNPs were found tional impact and thus may explain some of the pheno- within the CCHCR1 gene, which is involved in hair typic differences that exist between CHCU and CHFR loss [36]. Three SNPs (rs110552603, rs207611773 and animals. For example, 1253 SNPs were detected within rs381805999) are located within the 3′UTR and two the ATP9A gene, which is associated with heat tolerance SNPs (rs17871433 and rs377855638) can potentially in pigs [32]. Among these SNPs, one (rs207874965) is impact splicing. In addition, we found ten missense located within the 3′UTR, five SNPs have the potential SNPs and, in the upstream region, seven SNPs, which to alter splicing, and 22 SNPs are in the upstream region, all change TFBS and one of them (rs109692098) is fixed including 20 that modify TFBS. We found 251 SNPs Table 2 Selected candidate genes located within the putative selective intervals Genome position Gene name Size (bp) Number of SNPs Function Reference 13:79263518-79324185 PTPN1 60,667 382 Heat tolerance Klaman et al. [35]) 13:80167668-80262763 ATP9A 95,095 1253 Heat tolerance Kim et al. [32]) 15:8104485-8222755 PGR 118,270 1182 Heat tolerance Tsubota et al. [41] 17:17467450-17473822 UCP1 6372 118 Heat tolerance Charkoudian and Stachenfeld [38] 23:27779849-27791241 CCHCR1 11,392 133 Hair development Oka et al. [36] 23:27808496-27812710 CDSN 4214 466 Hair development Leclerc et al. [48] 23:28775534-28803895 GABBR1 28,361 251 Heat tolerance Haller et al. [33] Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 8 of 11 for the alternative allele in CHCU. Another candidate the main causes of the whole genomic differentiation at gene involved in hair development is the corneodesmo- the non-functional positions of the genome, together sin gene (CDSN), which contained 466 SNPs within our with B. indicus introgression, since an increase of exclu- bovine samples, among which one SNP (rs462034580) sive variability also increases population differentiation. can potentially impact splicing of this gene, 194 SNPs Nevertheless, selection of genetically beneficial variants are located in the upstream region, and three are mis- can also shape the local genetic differentiation at specific sense deleterious SNPs (rs209222317, rs434552200 and regions, as well as the phenotypic patterns of this breed. rs479537418). Interestingly, one of the CDSN deleterious The statistics that we used should correct for demo- mutations (rs434552200) is nearly fixed for the alterna- graphic effects, and should allow us to detect unusual tive allele in CHCU, with a frequency of 0.994. Further patterns that are compatible with positive selection. experimental analysis of the impact of these variants is The B. indicus component had a measurable impact needed to determine whether they are involved in some on increasing nucleotide variability in CHCU compared of the phenotypic differences between CHCU and CHFR. to that in the European breed, more than offsetting the effect of the founder bottleneck (Table 1). It is tempting Discussion to hypothesize that the detected B. indicus introgres- Creole cattle refer to the descendants of the first Euro- sion is related to the adaptation of CHCU. If this was pean animals that have adapted to local tropical condi- the case, genomic regions with a high proportion of B. tions on the American continent. However, as Burgos-Paz indicus should be enriched in signatures of selection. et al. [37] showed for ‘creole’ pigs, the origin of these ani- We did find a significant excess of B. indicus component mals is usually mixed and, usually there is little trace of within selective windows (Fig. 4), but it was too small to the original founders However, in very few cases, the ori- explain all the adaptive events detected. Although the gin of extant animals can be tracked accurately and, for B. indicus genomic component is associated with adap- that reason, CHCU is a unique population for which ped- tation (Fig. 4), the alleles of B. indicus origin cannot be igree records have been maintained over most of the time considered as the main drivers of selection. This leads us and isolation has been, in principle, guaranteed. to hypothesize that, likely, most of the adaptation events In spite of this assumed isolation, our findings confirm that have occurred in CHCU are due to changes in allele that CHCU was crossed with B. indicus animals, as previ- frequencies that were already present in the French ously suggested by Rodriguez-Valera et al. [6] based on Charolais, i.e., soft sweeps. SNP array data. In our study, we estimated that the per- Our results confirmed 17 of the 104 genomic regions centage of the CHCU genome originating from B. indicus reported by Rodriguez-Valera et al. [6]. The difference is in the order of 4 to 8%, depending on the method used in the number of identified genomic intervals between or the individual considered (Figs. 1d and 2, and (see these two studies can be partly explained by the fact Additional file 2: Figure S1)). The proportion of B. indi- that Rodriguez-Valera et al. [6] annotated both posi- cus component varied largely, both across chromosomes tive and negative extreme Rsb values and did not apply and individuals (Fig. 2). If the B. indicus introgression is a second filter based on the permutated p value. In our assumed to be a recent and sporadic event, it explains study, we focused on positive and significant (permutated this imbalance across chromosomes and individuals, p-value < 0.05) genomic intervals, since our aim was to which will be smoothed out in the future generations as detect regions in which the putative selective pressure is recombination events increase, and the correlation of specific to the Cuban population. Moreover, the detec- the proportion of the B. indicus alleles between individu- tion of putative signatures of selection in [6] was based als across windows is expected to increase with time. We on SNP chip data and, consequently, the sizes of the found a positive low correlation between proportion of B. genomic intervals were much larger (~ 700 kb vs. 30 kb indicus and FST (see Additional file 3: Figure S2), which here). Therefore, our results provide a better resolution suggests that the divergence between founder (CHFR) and accurate identification of putative signatures of selec- and derived (CHCU) populations is probably not due to tion between CHCU, CHFR and CHCA, which in turn B. indicus introgression only. facilitates the identification of genetic variants within The genetic differentiation in CHCU may also result candidate genes related to the adaptation to tropical from the strong bottleneck that occurred during the conditions. importation process, since it is likely that only a few hap- Within the selective windows, we found seven genes lotypes were introduced in Cuba. This differentiation related to thermotolerance and hair development might have increased due to genetic drift if the effective (Table 2), which are both key traits in the adaptation to size of this population remained low during the follow- tropical conditions and might explain some of the pheno- ing generations. Founder effect and genetic drift are likely typic differences between the CHFR and CHCU breeds. Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 9 of 11 Regarding thermotolerance, the genes ATP9A, GABBR1, had large differences in allele frequency between CHCU PGR, PTPN1 and UCP1 were putatively under selec- and CHFR. For example, two of these SNPs can poten- tion. Previously, Dikmen et al. [34] reported an associa- tially disrupt the binding sites of heat shock factors. As tion between a SNP in PGR and rectal temperature in US previously mentioned, the corneodesmosin (CDSN) gene Holstein lactating cows exposed to heat stress. This SNP is also related to hair development. Mutations in the (rs109506766) is an intronic G/C SNP and, interestingly, human CDSN are associated with hypotrichosis simplex, the frequency of the G allele differs between CHCU and a scalp-specific hair loss [46] and mice having undergone CHFR (0.38 versus 0.53). We also found several SNPs targeted inactivation of Cdsn showed rapid hair loss, in PGR that had large differences in allele frequency which confirms the essential role of Cdsn for maintaining between CHCU and CHFR and that might have a poten- the architecture of the hair follicle [47, 48]. tial functional impact. For example, 16 of these SNPs can alter the binding sites of transcription factors, and Conclusions interestingly, seven of these candidate regulatory variants In conclusion, analysis of whole-genome data confirms modify the binding sites of heat shock factors. that the CHCU animals included in this study are closely We also found one deleterious missense variant related to Charolais cattle from France and Canada, but (r42676011) in the codon for an arginine amino-acid also reveal a limited introgression of B. indicus in CHCU. within the ligand-binding domain of the progesterone We observed signals of recent adaptation to tropical receptor (PGR) gene, which is conserved among 70 euthe- conditions between CHCU and CHFR founder popula- rian mammals. Several studies have shown that proges- tions, which are largely independent of the B. indicus terone has a vasoconstrictive effect, which reduces the introgression, which suggests that most of the selection cutaneous blood flow when the temperature of the skin events are caused by soft sweeps. Some of the identified increases, and ultimately decreases heat dissipation (e.g. regions harbor genes that are involved in thermogenesis [38]. Moreover, glucocorticoids induce heat resistance (e.g., ATP9A, GABBR1, PGR, PTPN1 and UCP1) and hair in mammalian cells, whereas progesterone, a glucocorti- development (CCHCR1 and CDSN). We also report the coid antagonist, inhibits the development of this resist- presence of SNPs within these genes that can have a func- ance [39]. These findings suggest that PGR might play an tional impact and might explain some of the phenotypic important physiological role in reducing body heat loss. differences observed between CHCU and CHFR animals. The mitochondrial uncoupling protein 1 (UCP1) gene Future experimental work is needed to evaluate the role is predominantly expressed in brown adipose tissue and of these a priori relevant genes in CHCU adaptation. plays major roles in regulating body temperature, meta- bolic rate and controlling energy expenditure via both Supplementary Information non-shivering thermogenesis and diet-induced ther- The online version contains supplementary material available at https​://doi. mogenesis [40]. UCP1 is essential for maintaining body org/10.1186/s1271​1-020-00597​-9. temperature in non-cold conditions [41]. Variants in the bovine UCP1 have been found to be associated with Additional file 1: Table S1. Description of the samples included in this study. milk yield, milk fat percentage and milk protein percent- Additional file 2: Figure S1. Proportion of indicus per CHCU samples. age [42], but to our knowledge not to thermotolerance. We can hypothesize that genetic variations within UCP1 Additional file 3: Figure S2. Relationship between proportion of indicus in CHCU and FST CHCU-CHFR, each dot corresponds to the FST value of a might partly explain the heat tolerance of the CHCU 30-kb window. cattle. These variants should have a negative impact on Additional file 4: Table S2. Description of the intervals identified as UCP1 by reducing its activity or expression, and thus selective sweeps between Cuban and French Charolais. it would be interesting to perform an experimental Additional file 5: Table S3. Putative regulatory SNPs (rSNPs) located confirmation. upstream of interesting genes. CCHCR1 is involved in hair development and encodes the coiled-coil alpha-helical rod protein 1, which is Acknowledgements expressed in basal keratinocytes [43]. Its exact function We thank the Empresa Genética y Cría Manuel Fajardo and farmers of San José del Retiro for logistic support and help during the sampling. remains unknown, but it has been shown to play a wide variety of roles in steroidogenesis, proliferation, differ- Authors’ contributions entiation and cytoskeletal organization [44, 45]. Interest- YRC, MPE and DR designed the study. YRC and YRV performed animal sam‑ pling. LCRA, SERO, JLA, MPE and DR analyzed the data and interpreted the ingly, Oka et al. [36] have identified a missense variant in results. YRC, MPE and DR drafted the manuscript with the help of the rest of CCHCR1 that is associated with alopecia areata, an auto- coauthors. All authors read and approved the final manuscript. immune disease affecting the hair follicle. We found sev- eral functional candidate variants within CCHCR1 that Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 10 of 11 Funding the milestones of past selection in the bovine genome. Front Genet. MPE and SERO were funded by MINECO grant AGL2016-78709-R and from the 2015;6:36. EU through the BFU2016-77236-P (MINECO/AEI/FEDER, EU) and the “Centro de 12. Letaief R, Rebours E, Grohs C, Meersseman C, Fritz S, Trouilh L, et al. Iden‑ Excelencia Severo Ochoa 2016-2019″ award SEV-2015-0533. LCRA is funded tification of copy number variation in French dairy and beef breeds using by “Don Carlos Antonio López” Graduate program (BECAL) from Paraguay. next-generation sequencing. Genet Sel Evol. 2017;49:77. YRC was funded by the European Union H2020 Research and Innovation 13. Guillocheau GM, El Hou A, Meersseman C, Esquerré D, Rebours E, Letaief programme under Marie Skłodowska-Curie grant (P-Sphere) agreement No R, et al. Survey of allele specific expression in bovine muscle. Sci Rep. 6655919. This work was partially funded by the Animal Genetics Division of 2019;9:4297. INRAE as part of the CHACUBA project. Y. Rodriguez-Valera was funded by the 14. Li H, Durbin R. Fast and accurate short read alignment with Burrows- France Embassy in Havana (Cuba) and INRAE. Wheeler transform. Bioinformatics. 2009;25:1754–60. 15. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Ethics approval and consent to participate et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing Not applicable. next-generation DNA sequencing data. Genome Res. 2010;20:1297–303. 16. Pérez-Enciso M, de los Campos G, Hudson N, Kijas J, Reverter A. The Consent for publication “heritability” of domestication and its functional partitioning in the pig. Not applicable. Heredity. 2017;118:160–8. 17. Browning BL, Browning SR. Genotype imputation with millions of refer‑ Competing interests ence samples. Am J Hum Genet. 2016;98:116–26. The authors declare that they have no competing interests. 18. Sandeli A, Alkema W, Engström P, Wasserman WW, Lenhard B. JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Author details Nucleic Acids Res. 2004;32:D91–4. 1 Plant and Animal Genomics, Centre de Recerca en Agrigenòmica (CRAG) 19. Khamis AM, Motwalli O, Oliva R, Jankovic BR, Medvedeva YA, Ashoor H, CSIC-IRTA-UAB-UB, Campus UAB, Bellaterra, Spain. 2 Université Paris-Saclay, et al. A novel method for improved accuracy of transcription factor bind‑ INRAE, Jouy‑En‑Josas, AgroParisTech, GABI 78350, France. 3 INRAE, SIGENAE, ing site prediction. Nucleic Acids Res. 2018;46:e72. Jouy‑En‑Josas 78350, France. 4 INRAE, GeT-PlaGe, Genotoul, Castanet‑Tolosan, 20. Knüppel R, Dietze P, Lehnberg W, Frech K, Wingender E. TRANSFAC US 1426, France. 5 Facultad de Ciencias Agropecuarias, Universidad de retrieval program: a network model database of eukaryotic transcription Granma, Granma, Cuba. 6 Institut Català de Recerca I Estudis Avançats regulating sequences and proteins. J Comput Biol. 1994;1:191–8. (ICREA), Barcelona, Spain. 7 Animal Breeding and Genetics Program, Institute 21. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA for Research and Technology in Food and Agriculture (IRTA), Torre Marimon, target sites in mammalian mRNAs. eLife. 2015;4:e05005. Caldes De Montbui 08140, Spain. 22. Watterson GA. On the number of segregating Sstes in genetical models without recombination. Theor Pop Biol. 1975;7:256–76. Received: 3 September 2020 Accepted: 11 December 2020 23. Hudson RR, Slatkin M, Maddison WP. Estimation of levels of gene flow from DNA sequence data. Genetics. 1992;132:583–9. 24. Ferretti L, Raineri E, Ramos-Onsins S. Neutrality tests for sequences with missing data. Genetics. 2012;191:1397–401. 25. Tang K, Thornton KR, Stoneking M. A new approach for using genome References scans to detect recent positive selection in the human genome. PLoS 1. Naves M, Flori L, Thevenon S, Gautier M. Adaptation of tropical cattle Biol. 2007;5:e171. breeds to their environment, in the perspective of climatic change. 26. Alexander DH, Novembre J, Lange K. Fast model-based estimation of In Proceedings of the Climate-Smart Agriculture 2015 Global Science ancestry in unrelated individuals. Genome Res. 2009;19:1655–64. Conference: 16-18 March 2015;Montpellier; 2015. 27. Guan Y. Detecting structure of haplotypes and local ancestry. Genetics. 2. Upadhyay M, Bortoluzzi C, Barbato M, Ajmone-Marsan P, Colli L, Ginja 2014;196:625–42. C, et al. Deciphering the patterns of genetic admixture and diversity 28. The Bovine HapMap Consortium, Gibbs RA, Taylor JF, Van Tassell in southern European cattle using genome-wide SNPs. Evol Appl. CP, Barendse W, Eversole KA, et al. Genome-wide survey of SNP 2019;12:951–63. variation uncovers the genetic structure of cattle breeds. Science. 3. Ajmone-Marsan P, Lenstra JA, Fernando Garcia J, The Globaldiv Con‑ 2009;324:528–32. sortium. On the origin of cattle: how aurochs became domestic and 29. Kasarapu P, Porto-Neto LR, Fortes MRS, Lehnert SA, Mudadu MA, colonized the world. Evol Anthropol Issues News Rev. 2010;19:148–57. Coutinho L, et al. The Bos taurus-Bos indicus balance in fertility and milk 4. Koufariotis L, Hayes BJ, Kelly M, Burns BM, Lyons R, Stothard P, et al. related genes. PLoS One. 2017;12:e0181930. Sequencing the mosaic genome of Brahman cattle identifies historic and 30. Lin BZ, Sasazaki S, Mannen H. Genetic diversity and structure in bos recent introgression including polled. Sci Rep. 2018;8:17761. taurus and bos indicus populations analyzed by snp markers. Anim Sci J. 5. Ginja C, Gama LT, Cortés O, Burriel IM, Vega-Pla JL, Penedo C, et al. The 2010;81:281–9. genetic ancestry of American Creole cattle inferred from uniparental and 31. McTavish EJ, Decker JE, Schnabel RD, Taylor JF, Hillis DM. New World cattle autosomal genetic markers. Sci Rep. 2019;9:11486. show ancestry from multiple independent domestication events. Proc 6. Rodriguez-Valera Y, Renand G, Naves M, Fonseca-Jiménez Y, Moreno- Natl Acad Sci USA. 2013;110:E1398–406. Probance TI, Ramos-Onsins S, et al. Genetic diversity and selection 32. Kim KS, Seibert JT, Edea Z, Graves KL, Kim ES, Keating AF, et al. Char‑ signatures of the beef ‘Charolais de Cuba’ breed. Sci Rep. 2018;8:11005. acterization of the acute heat stress response in gilts: III. Genome- 7. Renand G, Menendez-Buxadera A, Krauss D, Menissier F. Experimental wide association studies of thermotolerance traits in pigs. J Anim Sci. comparison of the Cuban Charolais cattle selected in tropical environ‑ 2018;96:2074–85. ment with the French Charolais cattle. In Proceedings of the 4­ th Rencon‑ 33. Haller C, Casanova E, Müller M, Vacher CM, Vigot R, Doll T, et al. Floxed tres Recherches Ruminants: 4-5 December 1997; Paris. 1997. allele for conditional inactivation of the GABAB(1) gene. Genesis. 8. Ribas M. Gene frequencies in the blood group systems of the Cuban 2004;40:125–30. Charolais. Ann Genet Sel Anim. 1981;13:293–300. 34. Dikmen S, Wang XZ, Ortega MS, Cole JB, Null DJ, Hansen PJ. Single 9. Peripolli E, Metzger J, de Lemos MVA, Stafuzza NB, Kluska S, Olivieri BF, nucleotide polymorphisms associated with thermoregulation in lactating et al. Autozygosity islands and ROH patterns in Nellore lineages: evidence dairy cows exposed to heat stress. J Anim Breed Genet. 2015;132:409–19. of selection for functionally important traits. BMC Genomics. 2018;19:680. 35. Klaman LD, Boss O, Peroni OD, Kim JK, Martino JL, Zabolotny JM, et al. 10. Alshawi A, Essa A, Al-Bayatti S, Hanotte O. Genome analysis reveals Increased energy expenditure, decreased adiposity, and tissue-specific genetic admixture and signature of selection for productivity and envi‑ insulin sensitivity in protein-tyrosine phosphatase 1B-deficient mice. Mol ronmental traits in Iraqi cattle. Front Genet. 2019;10:609. Cell Biol. 2000;20:5479–89. 11. Utsunomiya YT, Pérez O’Brien AMP, Sonstegard TS, Sölkner J, Garcia JF. 36. Oka A, Takagi A, Komiyama E, Yoshihara N, Mano S, Hosomichi K, et al. Genomic data as the “hitchhiker’s guide” to cattle adaptation: tracking Alopecia areata susceptibility variant in MHC region impacts expressions Ramírez‑Ayala et al. Genet Sel Evol (2021) 53:3 Page 11 of 11 of genes contributing to hair keratinization and is involved in hair loss. 45. Tervaniemi MH, Siitonen HA, Söderhäll C, Minhas G, Vuola J, Tiala I, et al. EBioMed. 2020;57:102810. Centrosomal localization of the psoriasis candidate gene product, 37. Burgos-Paz W, Souza CA, Megens HJ, Ramayo-Caldas Y, Melo M, Lemús- CCHCR1, supports a role in cytoskeletal organization. PLoS One. Flores C, et al. Porcine colonization of the Americas: a 60 k SNP story. 2012;7:e49920. Heredity. 2013;110:321–30. 46. Levy-Nissenbaum E, Betz RC, Frydman M, Simon M, Lahat H, Bakhan 38. Charkoudian N, Stachenfeld N. Sex hormone effects on autonomic mech‑ T, et al. Hypotrichosis simplex of the scalp is associated with non‑ anisms of thermoregulation in humans. Auton Neurosci. 2016;196:75–80. sense mutations in CDSN encoding corneodesmosin. Nat Genet. 39. Fisher GA, Anderson RL, Hahn GM. Glucocorticoid-induced heat resist‑ 2003;34:151–3. ance in mammalian cells. J Cell Physiol. 1986;128:127–32. 47. Matsumoto M, Zhou Y, Matsuo S, Nakanishi H, Hirose K, Oura H, et al. 40. Chouchani ET, Kazak L, Spiegelman BM. New advances in adaptive ther‑ Targeted deletion of the murine corneodesmosin gene delineates mogenesis: UCP1 and beyond. Cell Metab. 2019;29:27–37. its essential role in skin and hair physiology. Proc Natl Acad Sci USA. 41. Tsubota A, Okamatsu-Ogura Y, Bariuan JV, Mae J, Matsuoka S, Nio-Kob‑ 2008;105:6720–4. ayashi J, et al. Role of brown adipose tissue in body temperature control 48. Leclerc EA, Huchenq A, Mattiuzzo NR, Metzger D, Chambon P, Ghyselinck during the early postnatal period in syrian hamsters and mice. J Vet Med NB, et al. Corneodesmosin gene ablation induces lethal skinbarrier Sci. 2019;81:1461–7. disruption and hair-follicle degeneration related to desmosome dysfunc‑ 42. Zhou H, Cheng L, Byun SO, Aizimu W, Hodge MC, Edwards GR, et al. Hap‑ tion. J Cell Sci. 2009;122:2699–709. lotypic variation in the UCP1 gene is associated with milk traits in dairy cows. J Dairy Res. 2017;84:68–75. 43. Tiala I, Suomela S, Huuhtanen J, Wakkinen J, Hölttä-Vuori M, Kainu K, et al. Publisher’s Note The CCHCR1 (HCR) gene is relevant for skin steroidogenesis and downreg‑ Springer Nature remains neutral with regard to jurisdictional claims in pub‑ ulated in cultured psoriatic keratinocytes. J Mol Med. 2007;85:589–601. lished maps and institutional affiliations. 44. Tiala I, Wakkinen J, Suomela S, Puolakkainen P, Tammi R, Forsberg S, et al. The PSORS1 locus gene CCHCR1 affects keratinocyte proliferation in transgenic mice. Hum Mol Genet. 2008;17:1043–51. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions