Creole cattle were brought to the American continent by Spanish conquerors in 1493. They were initially introduced to tropical environments in the Caribbean islands and then spread across the continent (Primo, 1992). Towards the south, Lima (Peru) was the main focus of dispersal and they crossed to western Bolivia, a region characterised for having one of the three main plateaus across the world, the Andean Plateau, located at ~3,000 m above the sea level (Primo, 1992; Friedrich and Wiener, 2020). These high-altitude environments have lower barometric pressure and ~40% lower atmospheric oxygen pressure than the sea level, colder temperatures, and increased UV radiation. As a result of long-term evolution, Creole cattle from the highland developed several strategies and mechanisms to evolve variations in the cardio-respiratory system, metabolic pathways, and morphological traits in order to adapt to their local environment at high altitude (Li et al., 2021). During the XVI century, the migration continued from the west highland to the tropical and subtropical plains in eastern Bolivia. This region required the adaptation to new environmental conditions and the development of resistance to specific diseases (e.g. tick, babesiosis, anaplasmosis) and heat stress, and changes in phenotypic features like coat, as shorter and thinner hair is beneficial in warm climates (da Silva et al., 2003; Bayssa et al., 2021). Adaptive evolutionary mechanisms may involve changes in one or few loci with major effects or in multiple loci with small effects (Storz et al., 2010; Friedrich and Wiener, 2020). These genotypic, and consequently phenotypic, variations were essential to allow them to survive and develop in these harsh habitats (Storz et al., 2010; Aguirre-Rofrio et al., 2019; Friedrich and Wiener, 2020; Rojas-Espinosa et al., 2023; Alvarez-Cecco et al., 2024).
To assess selection footprints, it is essential to propose a correct hypothesis and select the appropriate population to validate it. For this reason, firstly the population structure and purity of HBC and LBC were evaluated. The fastSTRUCTURE and PCA clearly differentiate both populations. The cluster analysis showed 3-8% of common ancestry between zebu breeds (Nelore and Bahman) and Bolivian Creole populations. This could be a consequence of a Zebu introgression on Creole cattle and/or the foundational origin of those Zebu breeds that include absorption of Creole dams. Lirón et al. (2006) reported Zebu introgression using autosomal (4-8%) and holandric (≈10%) microsatellites in Bolivian Creole cattle. Additionally, HBC presented high levels of Ho introgression which could be related to the introduction of Holstein animals into the Highland region according to oral history. This crossbreeding strategy was carried out in order to improve milk production since pure Holstein individuals exhibited high mortality rates in this highland environment. Considering this context, HBC-Ho and LBC-Ho comparisons were added in the current footprint analysis to identify regions under selection due to Ho introgression.
The analyses used for selection footprints and gene ontology were interpreted together in order to find common regions and reveal the potential mechanisms involved in divergent adaptation to highland and subtropical environments. The most prominent common region was located in the centre of BTA20. This region contains the candidate gene responsible for the slick-hair phenotype often found in Creole cattle breeds (Olson et al., 2003; Porto-Neto et al., 2018; Sosa et al., 2021; Sosa et al., 2022). Littlejohn et al. (2014) first reported a mutation in the prolactin receptor (PRLR) responsible for the slick-hair phenotype in purebred Senepol caused by a single base deletion [20:39136558 GC > G] in exon 10 that introduces a stop codon and results in a truncated protein. Additional studies in South American Creole cattle breeds from warm environments found that individuals with slick-hair phenotype were discordant with this reported frameshift mutation. Instead, three nonsense variants leading to stop codons and a SNP which produced a synonymous mutation were found, all in the same region of the PRLR sequence encoding the cytoplasmic portion of the protein receptor (Porto-Neto et al., 2018). All these collectively termed slick mutations generate truncated proteins with nearly identical effects on the protein function. Although the specific mechanism by which these mutations alter prolactin signalling is not known, they appear to enhance the inhibition of hair growth caused by the prolactin (Sosa et al., 2022). Therefore, the slick-haired animals are characterised by a short sleek hair coat and lower follicle density. This feature confers superior ability to thermoregulate under heat stress conditions, through heat loss from skin convection and conduction (Olson et al., 2003; Porto-Neto et al., 2018; Florez-Murillo et al., 2021; Sosa et al., 2021). This selection sweep was detected when LBC was compared to HBC and Ho, but was not observed in HBC-Ho. In LBC, this chromosomal region exhibited higher LD and lower haplotypic diversity than HBC, evidencing the positive selection of slick-hair phenotype in Creole cattle from eastern Bolivia. Moreover, considering the reported QTLs within this region, the matched peaks of selection footprints and Ho ancestry may indicate selection towards dairy production while maintaining the long hair necessary to adapt to western highland Bolivia. Other candidate genes found in this region were SLC45A2, HSPB3, and DNAJC21. SLC45A2 is involved in the melanin synthesis and is associated with skin and coat pigmentation variation in several species (Mariat et al., 2003; Soejima et al., 2007; Dooley et al., 2013; Wang et al., 2016; Bâlteanu et al., 2021). Ding et al. (2022) found a relationship between the different alleles of SLC45A2 and heat tolerance in indigenous Chinese cattle. Meanwhile, HSPB3 and DNAJC21 are heat shock proteins (HSPs) and have been related to thermotolerance through association, selection footprint and transcriptomic studies (Otto et al., 2019; Lemal et al., 2023; Wang et al., 2024; Alvarez-Cecco et al., 2024). In agreement, QTLs for coat texture and hair length were reported in these regions (https://www.animalgenome.org/cgi-bin/QTLdb/index). These results support the divergent adaptation related with the local environment of Bolivian Creole cattle, while short slick-hair is beneficial for heat loss through skin convection and conduction in lowland subtropical environments, long hair is desirable for highland Creole cattle which are exposed to colder temperatures.
Considering environmental microorganisms and adaptation, livestock in high or lowland areas are exposed to different pathogens. It has been reported that the diversity and composition of the skin microbiome, which is associated with animal’s health, is different when comparing high and lowland adapted individuals (Zeng et al., 2017; Sun et al., 2019; Ma et al., 2019a). In this sense, it was expected that the second main peak was detected in the Bovine Lymphocyte Antigen (BoLA) system located in BTA23. This selection sweep included Class I (BoLA Class I, BoLA-NC1, TRIM, JSP.1) and Class IIa (e.g., BoLA-DQA, BoLA-DQB, BoLA-DRB3) genes and QTLs related to immune response, infectious diseases, and parasite load. Previous works have demonstrated association between different resistance to specific diseases and the allelic variability of BoLA genes. Particularly, BoLA-DRB3 alleles were widely associated with infectious diseases, such as virus-induced lymphoma and proviral load in bovine leukemia virus (BLV) infection, somatic cell count in milk in mastitis, endo and ectoparasites, immune response traits, response to vaccination and production traits (e.g., milk yield). Furthermore, BoLA-DQA1 was associated with proviral load in BLV infection and mastitis (Reviewed in Takeshima and Aida, 2008; Aida et al., 2015). As mentioned above, this observed peak on BTA23 could be due to the differential exposure to pathogens. Cattle from tropical regions exhibit high resistance to infestation by ectoparasites (e.g., ticks; Ortega et al., 2023) and endoparasites (e.g., anaplasma, babesia; Casa et al., 2023) while highland Creole cattle are more susceptible, particularly evidenced when they are moved to tropical plains. Additionally, one of the main features of the BoLA region is the gene copy number variation. Qiu et al. (2012), studying the high altitude adapted Yak species, proposed that the presence of multicopy genes plays a relevant role in the divergent genetic architecture of adaptation, particularly in the functional categories ‘olfactory sensation’ and ‘host defence immunity’. Remarkably, this region on BTA23 contained several genes that belong to the olfactory receptor family (OR genes) and reported QTLs related to immune response, tick resistance and disease susceptibility (https://www.animalgenome.org/cgi-bin/QTLdb/index). Other peaks included genes involved in immune response, located on BTA20 (IL7R, IL6ST, IL31RA, C6, C7, OSMR, LIFR), and BTA5 (STAT6, NKG2A, IRAK4, KLR and CLEC genes). GO analysis evidenced significant enriched KEGG pathways and biological processes related to adaptive immunity and immune response to multiple diseases. These findings support the hypothesis that Bolivian Creole cattle from highland and tropical environments present divergent adaptation in response to the differential exposure to pathogens.
The adaptation to highland or tropical habitats also involves metabolic pathways and morphological traits. Reduced oxygen availability is the primary stressor of high altitude conditions and restricts the correct functioning of respiratory and cardiovascular systems (Ivy & Scott, 2015). Moreover, chronic hypoxia increases the production of reactive oxygen species (ROS; Wang et al., 2024). The reduction of the overall metabolic rate and modifications of the oxygen cascade and haematological system, such as red blood cell count and amount of haemoglobin, are necessary to cope with low oxygen levels (Hochachka et al., 1996; Weber, 2007; Stortz et al., 2010). In addition, smaller body size decreases the energy demands (Friedrich and Wiener, 2020). Experimental works have demonstrated that hypoxia upregulates the expression of ANXA2 and NDUFA4L2 genes, found in BTA10 and BTA5 respectively, through the direct action of hypoxia-inducible factor-1 (HIF-1; Huang et al., 2011; Liu et al., 2021). Remarkably, the HIF gene family was extensively associated with altitude adaptation in livestock species and other mammals (Reviewed in Friedrich and Wiener, 2020). It has also been shown that ANXA2 belongs to a common pathway relevant to fibrin homeostasis and angiogenesis, while NDUFA4L2 plays a key role in the development of pulmonary artery hypertension (PAH; Jacovina et al., 2009; Huang et al., 2011; Hajjar, 2015; Liu et al., 2021). In addition, a candidate gene for haematological parameters, CPLANE1, and two vascular endothelial growth factor (VEGF) genes, NRP1 and NRP2, were identified (Oh et al., 2002; Alghamdi et al., 2020; Yang et al., 2024). While angiogenesis helps to increase blood flow and oxygen supply under low oxygen conditions in high altitude, this physiological process of growing new blood vessels in the skin improves the heat dissipation in tropical environments. The enrichment analysis also evidenced pathways and terms related to glutathione (GSTA1-5 in BTA23) and retinol metabolism (HSD17B6, RDH16 and SDR9C7 in BTA5), which could be indicative of an oxidative and heat stress response. Previous works evidenced the selection of antioxidase-related genes in hypoxia-tolerant mammals (Wang et al., 2024). While GST gene family plays an important role in cellular detoxification to reduce the damage caused by ROS, HSD17B6 catalyses the oxido-reduction of different molecules (Deng et al., 2024). Regarding retinol metabolism, heat stress can decrease vitamin levels, including retinol, retinoyl β-glucuronide and biotin which have anti-oxidative properties removing ROS (Yang et al., 2022). Finally, in agreement with previous reports about body size and adaptation to different environments, four candidate genes related to height and carcass conformation were spotted in the sweep selection regions including PLAG1, CHCHD7, CAP2 and ARL15 (Purfield et al., 2019; Ghoreishifar et al., 2020; Zhang et al., 2022; Zhao et al., 2022).
In conclusion, the findings in the present work suggest that the divergent adaptation of Bolivian Creole cattle populations involves multiple and complex mechanisms. This includes changes in coat features and other morphological traits, immune response, and metabolic processes such as hypoxia and stress response. The inference ancestry analysis evidenced an uneven distribution of Ho introgression in the HBC genome. Except for BTA20, the enrichment regions did not match the selection footprints. Therefore, these sweeps could be a consequence of divergent adaptation of Bolivian Creole cattle to highland and tropical lowland environments.