Genetic Architecture and Signatures of Selection in the Caqueteño Creole (Colombian Native Cattle)

Abstract

Evolutionary mechanisms have shaped the genomic architecture of Colombian Creole cattle breeds. The mating and selection processes have impacted several traits, promoting differences within and between populations. Studies of population structure and selection signatures in Colombian Creole breeds are scarce, and need more attention to better understand genetic differentiation, gene flow, and genetic distance. This study aimed to analyze the population structure and identify selection imprints in the Criollo Caqueteño (CAQ) population. It used 127 CAQ animals genotyped with Chip HD 777,000 SNPs. The population structure analyses used discriminant principal component analysis (DAPC), integrated haplotype scoring (iHS), and index-fixing (Fst) methodologies to detect selection signals. We can highlight SNP regions on the genes TMPRSS15, PGAM2, and EGFR, identified by the Fst method. Additionally, the iHS regions for cluster 1 identified candidate genes on BTA 3 (CMPK1 and FOXD2), BTA 11 (RCAN1), and BTA 22 (ARPP21). In group 2, we can highlight the genes on BTA 4 (SLC13A4, BRAF), BTA 9 (ULBP), BTA 14 (CSMD3) and BTA 19 (KRTAP9-2). These candidate genes have been associated with fertility traits, precocity, growth, and environmental and disease resistance, indicating a genetic potential in CAQ animals. All this promotes a better understanding of the diversity and genetic structure in the CAQ population. Based on that, our study can significantly assist the sustainable development and conservation of the breed in the Colombian Amazon.

1. Introduction

Colombian Creole cattle breeds are descendants of cattle brought by Christopher Columbus in 1493 to the islands known today as the Dominican Republic and Haiti [1]. These animals spread throughout the South American territory and settled in different regions, and in some cases, without human intervention. For example, in Colombia, eight breeds formed that today are bred throughout the country: Romosinuano and Costeño con Cuernos on the Atlantic Coast; Blanco Orejinegro (BON) and Chino Santandereano in the mountainous areas; Hartón del Valle in the Cauca River Valley; Casanareño and San Martinero in the Orinoquía, and Criollo Caqueteño (CAQ) in the Amazon [2,3,4].

The CAQ breed’s history is associated with the colonization of the Amazon territory in the Caquetá state, carried out by Franciscan friars [5]. The Criollo Caqueteño cattle have survived due to robustness and adaptation traits, which resulted from natural selection against the conditions of the Colombian Amazon, becoming an alternative for livestock in the region. In addition, replacing native vegetation (prairies) with plating forages sustained the breed’s development in that region. Currently, the Criollo Caqueteño cattle have a low number of females with reproductive capacity [6,7,8]. With the development of next-generation sequencing (NGS) technologies, and single-nucleotide polymorphism (SNP) annotation, studies on population history, genetic diversity, and genome structure in domestic animals became possible [9,10,11,12,13]. This technology has been applied to cattle to study aspects such as the evolutionary history and genetic structure of different breed populations [14]. Genetic diversity studies help us understand breed evolution, genetic progress, and the level of differentiation between breeds. Additionally, it is fundamental to develop effective strategies to improve, manage and conserve the genetic resources of populations [11,15].

Several metrics are employed to assess genetic dissimilarity, with the Fst fixation index being a notable example [16]. Fst is used in studies to reveal genetic variations between populations and identify selection traces within a population [17]. Another approach involves haplotype-based methods in selective sweep models to identify genomic regions where mutations progress rapidly toward fixation, reducing diversity around the locus. In this context, integrated haplotype scoring (iHS) is a widely used method that is particularly effective in livestock studies [18,19,20]. The use of the Fst or iHS methodology depends on the temporal aspect of the selection events, with Fst being suitable for detecting past traces and iHS more suitable for identifying recent upbeat selections.

Few studies have been conducted on Colombian Creole breeds related to genetic structure and selection footprints. However, in purebreds and some crosses, several methods are reported to characterize genetic diversity within and between cattle breeds via single-nucleotide polymorphism (SNP) analyses on the entire genome [9,21,22], as well as selection signal identification [23,24,25]. In this study, we present a comprehensive analysis of the population structure and selection footprints of the CAQ population using high-density genomic information. The methods (iHS and Fst, XP-EHH) implemented scanned the whole genome, and we performed a functional enrichment of genes identified within regions harboring selection footprints. Our findings will be helpful in further studies on the conservation and genetic improvement of the CAQ population.

2. Materials and Methods

2.1. Ethics Statement

The experiment was carried out under article 2 of agreement No. 027 of 2020 of the Faculty of Agricultural Sciences, with the guidelines of the Ethics, Bioethics and Animal Welfare Committee of the University of the Amazon, code CEBBA 129, 2021.

2.2. Animals and Genotype

The genotypic database included 127 individuals (109 Females and 18 Males) genotyped with Chip HD 777,000 SNP by CAQ cattle, belonging to three farms in the department of Caquetá, Colombia. The total of individuals was selected after the quality control, which was performed at the sample level, considering the following criteria: call rate > 97, average call rate for samples passed > 98.5, and at the SNP level using the call rate threshold > 97.

2.3. Genetic Structure of the Population

In aiming to assess the structure of the CAQ cattle and understand the relationship within and between populations at the genomic level, principal component analysis (PCA) and discriminant analysis of principal component (DAPC) were performed using the adegenet v.2.1.7 package for R software [26], a similar approach utilized in our previous study with the Gir cattle breed [17]. The PCA approach allowed the classifying of individuals based on the reduced number of important orthogonal principal components (PC). The DAPC, based on the K-means clustering algorithm and Bayesian information criterion (BIC), defined the ideal number of PCs to retain in the analysis [26]. Additionally, the DAPC estimated α-score in the range of PCs chosen from 1 to 40, and the function considered the existence of three discriminants.

The groups identified during this step (Group 1—n = 24; Group 2—n = 96) were used in the following analysis.

2.4. Integrated Haplotype Scoring (iHS)

The iHS analysis was conducted using unpolarized alleles. The rehh v3.2.2 package allows the function “FALSE”, which is adequate for non-model organisms. The iHH (integrated HHE) values were calculated for the major and minor alleles according to Santos et al. [27]. Genomic regions in the iHS positive tail represent the state of the ancestral allele, while regions in the negative tail represent the state of the derived allele. For both tails, iHS values ≥ 5 (p < 0.01) were considered as signatures of selection in groups 1 and 2.

2.5. The Fixation Indice (Fst) and XP-EHH Analysis

The fixation index (Fst) is one of the most commonly used parameters to study the genetic differentiation between populations. This index is initially defined in terms of the two-gamete union correlations described by Wright [28]. The Fst was performed via VCFtools software using the --weir-fst-pop parameter to specify each population. The SNP cutoff was defined in three standard deviations above the mean for each autosome.

Additionally, we used the XP-EHH analysis, a robust cross-population statistic, to determine the selection process. It uses extended haplotype homozygosity (EHH) statistics based on specific fixed alleles inside each population [29], and we used the rehh package [30] to calculate that measure. The parameter used to consider the signature of selection was XP-EHH value ≥ 5 (p < 0.01) in the groups.

2.6. Candidate Genes and Functional Analysis

The genomic regions with significant SNPs for Fst and iHS were considered under selection. The gene annotation used the UMD3.1 bovine genome assembly from BioMart (www.ensembl.org/biomart, accessed on 22 July 2022) and NCBI (https://www.ncbi.nlm.nih.gov, accessed on 22 July 2022) databases. In addition, the Integrated Annotation, Visualization, and Discovery (DAVID) v6.8 database was used to identify significant (p < 0.05) Gene Ontology (GO) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways.

3. Results

3.1. Genetic Structure of the Population (PCA and DAPC)

In the PC analysis, we assessed the total population of CAQ to understand its structure. Principal components explained 39.4% of the total genetic variance. In Figure 1a, we can see a clear separation of the population according to PC1, showing the genetic distances between individuals forming three groups within the population.

The individuals of cluster 1 presented a greater distance from the other groups, which indicates that they are less genetically related to the others. On the other hand, clusters 2 and 3 presented a more dispersed clustering based on PC2, indicating a small genetic distance between individuals in the population compared to cluster 1. The population segregation (clusters) possibly indicates the use of sires from different families during the breed formation. Additionally, due to the lack of prior information related to genetic architecture or animal origin, the PCA analysis would have grouped the CAQ individuals according to their different demographic histories or use of different breeders, resulting in a more dispersed pattern in the PCA.

Additionally, the DAPC analysis (Figure 1c) identified two separate groups based on the genetic distance, where group 1 (blue) is distant from group 2 (overlapping red and yellow). Therefore, only two groups were considered to explain the Caqueteño Creole cattle better, and perform the signature selection analyses: group 1 (n = 24) and group 2 (n = 96).

3.2. Selection Signature of iHS in CAQ

The iHS was used to identify regions exhibiting evidence of selection. This linkage disequilibrium-based method provides increased power for assessing selection footprints within the population using information from marker SNPs. After adjustment for within-population false discovery rate (FDR), SNPs showing piHS values ≥ 5.0 (approximately corresponding to a p-value < 0.008) were considered significant for the CAQ population.

Positive and negative iHS values were considered in our study, aiming to identify old and more recent selection signatures. In total, 292 genomic regions were detected as selection signatures in the iHS test. The selection signatures iHS by chromosome for group 1 are shown in Figure 2a. The figure shows little evidence of selective forces in different genome regions. The most significant SNPs mapped to chromosome 11 (iHS = −5.82), chromosome 22 (iHS = −5.70), and chromosome 4 (iHS = 5.10) in the IHS analysis. Group 2 shows clear evidence of selective forces in different regions of the genome (Figure 2b). The most significant SNPs mapped to chromosome 4 (iHS = −6.81), chromosomes 19 and 22 (iHS = −6.70 and −6.45), and chromosome 20 (iHS = 5.71) in the IHS analysis.

Windows with a range of 500 kb containing the significant SNPs were investigated to identify regions with strong selection signatures. As a result, 29 and 43 candidate genes overlapped with significant windows for group 1 and group 2, respectively. In addition, chromosome 22 was identified in both iHS groups. However, due to many significant signals in the iHS analysis, we considered prioritizing candidate genes within the highlighted genomic regions based on extreme iHS values. Therefore, we could highlight the candidate genes on BTA3 (CMPK1 and FOXD2), BTA11 (RCAN1), and BTA22 (ARPP21) for group 1, while in group 2, the genes located on BTA4 (SLC13A4, BRAF), BTA9 (ULBP), BTA14 (CSMD3) and BTA19 (KRTAP9-2) were pinpointed.

3.3. Gene Enrichment Analysis

The enrichment analysis for the iHS methodology was performed for the two groups of CAQ cattle. As a result, most gene enrichment Gene Ontology (GO) terms for biological processes were attributed to cellular and metabolic processes (Figure 3).

3.4. Diversity Genetics and Selection Signatures in Groups (Fst and XP-EHH)

The genetic differentiation for the two groups was 0.045, a low value. When the value of Fst is close to zero, there is no genetic differentiation between the populations. Therefore, although the DAPC evidenced two groups within the CAQ population, the lack of significant differentiation indicates that the two groups belong to the same breed.

The Fst approach also allows the detection of selection footprints based on differences in allele frequencies between populations. SNPs identified as outliers based on the Fst method may be strong evidence of traces of ancient selection. Eight SNPs had Fst values greater than 0.6, located on chromosomes 1, 10, and 22 (Figure 4). The seven most significant XP-EHH genomic regions between groups 1 and 2 of the CAQ population are shown in Figure 5. The most significant SNPs are on chromosome 3 (XP-EHH = 6.25) and chromosome 1 (XP-EHH = 6.02).

4. Discussion

The first genome-wide study on selection footprints in the CAQ cattle was conducted with high-density SNPs genomic information (777,134 autosomal SNPs) using Bovine SNPs HD BeadChip (Illumina, San Diego, CA, USA). Our findings shed light on potential candidate genes/gene clusters involved in regions under selection in this breed. In addition, it can support a better understanding of the selection footprints in other Colombian Creole breeds.

Identification of Candidate Genes within the Selection Signatures of iHS, Fst and XP-EHH in the Genome

The positive and negative values of iHS show the presence of genomic regions of ancestral and derived alleles in groups 1 and 2 of the CAQ population (Figure 2). For the most significant SNP regions in group 1, candidate genes associated with environmental resistance, immunity, growth, and meat quality were identified. We can highlight the CMPK1 gene in BTA3, which was related to response to stress and thermotolerance between zebu and taurine cattle breeds [31]. Also located in BTA3, the FOXD2 gene was pointed out. That gene was associated with metabolic processes and immune responses in native Korean cattle [32], indicating a possible selection process towards environmental resistance in the CAQ breed. A similar property can be attributed to the RCAN1 (BTA11) and ARPP21 (BTA22) genes. RCAN1 was found to play a pivotal role in muscle development/degradation, affecting carcass traits and meat quality in the Nellore breed [33]. At the same time, ARPP21 was associated with tissue development and body size in Chinese Holstein cows [34], corroborating our findings that show the CAQ breed is adapted to high-temperature and -humidity environments, developing different adaptation mechanisms to survive and reproduce.

On the other hand, for group 2, the most significant SNPs were identified close to candidate genes associated with fertility, body conformation, and immunity. In BTA4, two genes (SLC13A4, the protein of which acts as a sulfate carrier to strengthen the mammal fetus during pregnancy [35], and BRAF, associated with sexual precocity in Nellore cattle [36]), can be highlighted. Additionally, SNP markers close to the CSMD3 gene (BTA14) point to a possible selection phenomenon, since that gene was linked to body size in Swedish [37] and Chinese Wagyu cattle [38].

Finally, related to immunity, we can pinpoint ULBP (BTA9), which interacts with the NKG2D receptor to activate effector cells in the bovine immune system [39], and the KRTAP9-2 gene (BTA19), associated with tick resistance in Brahman and Angus cattle [40]. These genes, associated with fertility, precocity, and environmental resistance traits, indicate genetic potential in the CAQ breed. Additionally, some gene regions associated with growth and carcass quality could indicate a tendency for the CAQ population to express their potential to produce quality meat when exposed to adequate nutritional and sanitary management.

On the other hand, the average Fst value was close to that observed in different commercial cattle breeds [17,41,42] and native Colombian cattle, the Blanco Oreji Negro (BON) Creole breed, with an average Fst of 0.036 [43]. The level of genetic differentiation between populations was low in neutral regions of the genome or regions of balanced selection, and divergent in regions subject to directional selection. In the case of the CAQ population, there is no evidence of directed selection between the groups.

In the Fst selection signatures (Figure 4), four of the eight SNPs with significant values overlapped candidate genes previously associated with disease resistance. For example, the TMPRSS15 gene (BTA1) has been associated with mastitis resistance. Additionally, three Fst peaks were within the genomic regions of BTA22, with the PGAM2 gene associated with muscle fat deposition, body conformation, tenderness, carcass, and meat quality, and the EGFR gene associated with adaptation and reproduction in cattle populations [44,45,46,47,48,49,50,51,52,53,54].

In addition, the XP-EHH methodology identified several putative selection signature segments in the CAQ population. The identified selection signature segments could be related to production, reproduction, and adaptation traits of the breed, since the most significant SNPs were on chromosome 3 (XP-EHH = 6.25) and chromosome 1 (XP-EHH = 6.02), showing a positive XP-EHH score that indicates selection occurred in the CAQ population [55].

In summary, the Fst and iHS methodologies mainly identified gene regions associated with environmental adaptation and growth, possibly due to being a native breed without directed genetic improvement. It is worth mentioning that the approach to the selection of CAQ cattle has not been defined in relation to producing milk or meat; however, in some regions of the signatures of selection, we identified a slight trend of genes associated with growth and traits related to meat quality. For example, genomic regions related to meat production traits were identified in a recent study on the CAQ breed using the homozygosis island methodology [7].

5. Conclusions

A cattle breed’s adaptation to a changing environment is crucial nowadays. In this sense, native cattle breeds can become valuable sources of genetic variability, since they may carry specific selection signatures responsible for adapting to local environmental and nutritional conditions. Our study provides a further understanding of the genomic architecture of the CAQ breed, allowing the identification of genomic regions and genes that were affected by selection during the breed’s formation history. Additionally, our findings will be helpful in the sustainable development and conservation of the CAQ breed in the Colombian Amazon region.