Whole-Genome Analysis Deciphers Population Structure and Genetic Introgression Among Bovine Species

It is known that throughout history and presently, taurine (Bos taurus) and indicine/zebu (Bos indicus) cattle were crossed with other bovine species (e.g., gayal, gaur, banteng, yak, wisent, and bison). Information on the role of interspecific hybridization to facilitate faster adaptation of the newly arrived domestic species to new environments is poorly known. Herein, we collected 266 samples of bovine species of the taurine, zebu, yak, and gaur from West Europe, Southwest China, Indian subcontinent, and Southeast Asia to conduct the principal component analysis (PCA), admixture, gene flow, and selection signature analyses by using SNPs distributed across the bovine autosomes. The results showed that the genetic relationships between the zebu, yak, and gaur mirrored their geographical origins. Three ancestral components of the European taurine, East Asian taurine, and Indian zebu were found in domestic cattle, and the bidirectional genetic introgression between the Diqing cattle and Zhongdian yak was also detected. Simultaneously, the introgressed genes from the Zhongdian yak to the Diqing cattle were mainly enriched with immune-related pathways, and the ENPEP, FLT1, and PIK3CA genes related to the adaptation of high-altitude hypoxia were detected. Additionally, we found the genetic components of the Zhongdian yak had introgressed into Tibetan cattle. The 30 selected genes were detected in Tibetan cattle, which were significantly enriched in the chemokine signaling pathway. Interestingly, some genes (CDC42, SLC39A2, and EPAS1) associated with hypoxia response were discovered, in which CDC42 and SLC39A2 played important roles in angiogenesis and erythropoiesis, and heart function, respectively. This result showed that genetic introgression was one of the important ways for the environmental adaptation of domestic cattle.

Introduction

The genus Bos belongs to the subfamily Bovidae and includes several species, such as Bos taurus (taurine), Bos indicus (indicine/zebu), Bos frontalis (gayal), Bos gaurus (gaur), Bos javanicus (banteng), and Bos grunniens (yak) (Wu et al., 2018), which play significant roles in human history. The high-throughput genome sequencing effort has revealed new and helpful insights on the domestication and environmental adaptation mechanisms operating on these Bos species (Hu et al., 2012; Qiu et al., 2012; Gautier et al., 2016; Chen et al., 2018; Wu et al., 2018).

Information from high-density genome-wide SNP panels has been used to classify the ancestral components and reconstruct the historical migration routes of domestic cattle in the world (Decker et al., 2014). For example, a previous study found that all Eurasian cattle stem from five ancestral components: European taurine, Eurasian taurine, East Asian taurine, Chinese indicine, and Indian indicine (Chen et al., 2018). Yet, little is known about the contribution of other Bos species into this genome architecture.

The domestication of B. taurus and B. indicus occurred in the Near East and Indus Valley, respectively. There is no reproductive isolation between the taurine and zebu, the two types of cattle can be hybridized and their offspring are completely fertile (Wei et al., 2016). Studies showed that domestic cattle crossed with other Bos species broke reproductive isolation by hybridization or backcrossing, resulting in genetic introgression among Bos species (Chen et al., 2018; Wu et al., 2018). It is expected that introgression will help Bos species to adapt faster to the local environment. Gao et al. used a model-based clustering and f4-statistics analyses to quantify banteng and gayal introgression into southern Chinese cattle (Gao et al., 2017). Bidirectional gene introgression between the yak and cattle was also detected from whole-genome sequencing analysis (Wu et al., 2018). Finally, the latter study also found signs of introgression between the yak and Tibetan cattle at genes related with hypoxia adaptation (EGLN1, EGLN2, and HIF3a) and other study and quantified the proportion of taurine cattle introgression (∼1.3%) in the yak genome (Medugorac et al., 2017). As demonstrated in other studies, high-density BovineHD BeadChip genotyping arrays can be useful tools to estimate for population structure, genetic introgression, and selection signatures of Bos species (Medugorac et al., 2017; Mukherjee et al., 2018). In this study, we analyzed a total of 266 individuals of 10 populations of four Bos species (taurine, zebu, yak, and gaur) from West Europe, Southwest China, Indian subcontinent, and Southeast Asia by using the Illumina BovineHD BeadChip. The main purposes of this study were to 1) examine whether the genetic relationships of these Bos species reflect their geographical origins; 2) detect genetic admixture and genetic introgression among these cattle populations; and 3) screen out the candidate genes related to high-altitude adaptation of Tibetan cattle.

Materials and Methods

Ethics Approval Statement

Ear skin tissue samples were collected by veterinary practitioners with the permission and in presence of the owners. Veterinarians followed standard procedures and relevant national guidelines to ensure appropriate animal care. This study was approved by the Committee on Animal Research and Ethics of Yunnan University (ynucae20190011).

Sample Collection and Genomic DNA Extraction

We collected a total of 266 samples of 10 populations of four Bos species (Table 1), including Holstein (n = 34), Diqing cattle (n = 16), Tibetan cattle (n = 39), Indian zebu (n = 25), Bangladesh zebu (n = 21), Myanmar cattle (n = 38), Zhongdian yak (n = 13), Tibetan yak (n = 58), Maiwa yak (n = 8), and Gaur (n = 14). Diqing cattle and Zhongdian yak live in Shangri-La, Diqing Prefecture, which is at an average altitude of 3,300 m. All samples were divided into four geographic regions: West Europe, Southwest China, Indian subcontinent, and Southeast Asia. Genomic DNA was extracted from ear skin tissues using the DNeasy Blood & Tissue kit (Qiagen).

TABLE 1

TABLE 1. Sample information of 266 samples from 10 bovine populations.

Quality Control Procedures

Genotyping was carried out using the Illumina BovineHD BeadChip (777,962 SNPs) (Rincon et al., 2011). Markers were filtered to exclude the unmapped UMD3.1 bovine genome assembly (Zimin et al., 2009), X, Y, and mitochondrial chromosomes. Finally, 735,293 autosomal SNPs were used for quality control analysis. During quality control, markers with a call rate less than 90% and minor allele frequency (MAF) lower than 0.05 and significantly deviated from Hardy Weinberg equilibrium (p < 0.001) were excluded for all analysis based on previous research (Uzzaman et al., 2014). Finally, 638,034 SNPs were excluded, leaving 97,259 SNPs for downstream analyses.

Genomic Diversity Analysis

The observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (f) for each population of four Bos species were calculated using PLINK v1.90 (Purcell et al., 2007; Chang et al., 2015). Moreover, heterozygosity estimates are sensitive to various ascertainment biases when SNPs discovered in one breed are used to genotype other breeds (Nielsen, 2004; Edea et al., 2015). Analyses of diversity indexes based on ascertained SNP data would produce false positive results (Lachance and Tishkoff, 2013). Previous studies demonstrated that pruning of SNPs in high linkage disequilibrium (LD) would reduce the influence of ascertainment biases (López Herráez et al., 2009; Kijas et al., 2012). Therefore, we used SNP data sets after quality control and LD-pruning to calculate the genomic diversity indexes.

Principal Component and Admixture Analyses

Before performing the principal component analysis (PCA) and admixture analyses, the autosomal SNPs data were further pruned for linkage disequilibrium (LD) higher than 0.5 using a sliding window approach of 50 SNPs and a step size of 5 SNPs (50 5 0.5) in PLINK v1.9 (Purcell et al., 2007; Chang et al., 2015). The PCA was executed using the SMARTPCA program in EIGENSOFT v6.1.4 (Price et al., 2006) to discern genetic relationships among three data sets: 1) all 266 samples from 10 populations; 2) all 79 yak samples, including 58 Tibetan yak, 13 Zhongdian yak, and 8 Maiwa yak from Southwest China; 3) all 84 zebu samples, including 38 Myanmar cattle from Southeast Asia and 25 Indian zebu and 21 Bangladesh zebu from the Indian subcontinent. Furthermore, an identity-by-state (IBS) matrix of genetic distance containing each pairwise combination of all individuals was generated by using PLINK v1.9 referred to a previously study (Zhou et al., 2016). We constructed a neighbor-joining (NJ) tree based on the IBS matrix in MEGA v7.0 (Kumar et al., 2016).

To further determine the relative contribution of different potential ancestors on the genomic structure of these 10 bovine populations, population admixture analysis was carried out by using ADMIXTURE v1.30 (Alexander et al., 2009). Additionally, the corresponding cross-validation error value of cluster (K = 1–10) was calculated in ADMIXTURE v1.30. These results of PCA and admixture analyses were visualized using the R package ggplot2 (Wickham, 2016).

Gene Flow Detection

To detect gene flow among populations, we used the three-population test (f3 statistics) and calculated their corresponding normalized value (z-scores) in the “qp3Pop” program implemented in the ADMIXTOOLS v6.0 (Patterson et al., 2012). The f3 statistics considers triplet of the populations (A, B, and C), where C is the test (target) population with A and B as reference (source) populations. If the z-score (z ≤−3.0) is significantly negative, showing the test population C has admixture from both the reference populations A and B.

To assess the direction of gene flow, the D statistics were implemented in the ADMIXTOOLS v6.0. The D statistics method considers the tree topology [[[W, X], Y], Z], where Z represents outgroup. The gaur population (14 individuals) was selected as outgroup based on previous similar research (Barbato et al., 2020). Y is the admixture population, and W and X are the test populations. The D statistics method counts the “ABBA” sites, where W and Z share the outgroup allele (A) and X and Y share the derived allele (B) as well as the “BABA” sites, where W and Y shares the derived allele and X and Z shares the outgroup allele. Admixture between Y and either of the test populations creates a significant difference between the ABBA and BABA counts, with a z-score >3.0 (gene flow between W and Y) or ≤ −3.0 (gene flow between X and Y).

To observe migration events among 10 populations, we constructed a maximum likelihood (ML) tree by using the TREEMIX v1.12 (Pickrell et al., 2012). We set gaur as the root and grouped together 1,000 SNPs to account for the linkage disequilibrium and used 5,000 SNPs to calculate standard errors of migration proportions based on a previous study (Wu et al., 2018). The migration edges were added to the ML tree until 99.8% of the variance of the ancestry between the populations was explained. After building the ML tree, we respectively allowed 1–6 migration events in the tree to draw the ML graphs and corresponding residuals.

Genetic Introgression Detection and Enrichment Analyses

To further detect the genomic regions of potential introgression among these 10 populations of Bos species, the genomic local ancestry was inferred using PCAdmix v1.0 (Brisbin et al., 2012). This software scans the target genome using a sliding window approach to estimate the relative ancestry proportions of the reference populations for each chromosome and haploid individual. Before running the PCAdmix software, genotypes of target population and reference population were phased using default parameters in BEAGLE v4.1 (Browning and Browning, 2013). Due to the absence of genetic distance, we assumed that the genetic distance was equal to physical distance (1 Mb ≈ 1 cM) (Rothammer et al., 2013). Gene annotation was performed for the introgression regions among populations based on the UMD3.1 reference genome. Additionally, to understand biological functions of introgressed genes, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were executed using the DAVID v6.8 (https://david-d.ncifcrf.gov/). Among the GO enrichment analysis included Molecular function (MF), Cellular component (CC) and Biological process (BP).

Selection Signature Analysis

To detect genes under positive selection in Tibetan cattle, Holstein as a reference population of a lowland taurine breed was used to investigate selection signatures in Tibetan cattle by using the population differentiation index (F_ST), cross-population extended haplotype homozygosity (XP-EHH), and cross-population composite likelihood ratio (XP-CLR) methods. The VCFTOOLS v0.1.15 (Danecek et al., 2011) was used to estimate the F_ST values using non-overlapping 50 kb windows sliding in 10 kb increments, and the result showed that the negative and missing F_ST values were discarded because these values have no biological interpretation (Akey et al., 2002).

Before running the XP-EHH and XP-CLR methods, the haplotypes were reconstructed for every chromosome using default parameters in BEAGLE v4.1 (Browning and Browning, 2013). The XP-EHH values was computed in SELSCAN v1.1.0 (Szpiech and Hernandez, 2014), then all XP-EHH values were normalized using the norm program implemented in SELSCAN v1.1.0b. Finally, the standardized XP-EHH values were sorted from largest to smallest; the SNP corresponding to the top 1% XP-EHH values was selected as the core site, and the upstream and downstream extensions were each extended by 50 kb as the selected regions. For the XP-CLR scores were computed in XP-CLR v1.0 (Chen H. et al., 2010). We used the non-overlapping sliding windows of 50 kb and a maximum number of 600 SNPs within each window to calculate the XP-CLR values. Furthermore, we downweighted the correlation level from which the SNPs contribute to XP-CLR results to 0.95. The regions with XP-CLR values in the top 1% of the empirical distribution were considered as candidate signatures, and the genes that span the window regions were defined as candidate genes (Taye et al., 2018).

All markers with top 1% F_ST, XP-EHH, and XP-CLR values were considered to represent a selection signature, and then gene annotation was performed based on the UMD3.1. We carried out the functional enrichment analysis and detected important candidate genes related to high-altitude adaptation by literature review.

Results

Genomic Diversity Index

We used both the 97,259 SNPs and 64,704 SNP data sets to calculate the genomic diversity index of 10 bovine populations (Table 2). For 97,259 SNP data sets, the results showed that the average Ho ranged from 0.029 in the Tibetan yak to 0.396 in Holstein and the average He ranged from 0.028 in the Tibetan yak/Maiwa yak to 0.362 in Holstein. Similarly, using the 64,704 SNP data sets, the lowest and highest heterozygosities were found in the Tibetan yak and Holstein, respectively. Furthermore, the minimum and maximum levels of inbreeding were found in Holstein and gaur population, respectively.

TABLE 2

TABLE 2. Genetic diversity index among Bos species.

Population Structure of Bos Species

After LD-pruning, the 64,704 SNP data sets were used for the PCA and admixture analyses. The PCA results of 10 cattle populations from four Bos species are shown in Figure 1. The first (PC1) and second (PC2) principal components accounted for 46.0% and 16.7% of the total variation, respectively. The PCA clearly separated the individuals into five clusters: one group consisting of Holstein cattle, second cluster formed by Diqing cattle and Tibetan cattle, third including all zebu cattle, fourth gathering all yak populations, and fifth groups constituted by gaur individuals (Figure 1B). Furthermore, the PC2 and third (PC3) principal components accounted for 16.7% and 8.3% of the total variation, which was similar to the results of PC1 and PC2 (Figure 1C).

FIGURE 1

FIGURE 1. Geographical distribution and principal component analysis (PCA) results of 10 bovine populations. (A) Region information of 10 populations; (B) PCA results of PC1 vs. PC2; (C) PCA results of PC2 vs. PC3; (D) PCA results of 84 samples from all zebu populations (Indian zebu, Bangladesh zebu, and Myanmar cattle); (E) PCA results of 79 samples from all yak populations (Tibetan yak, Zhongdian yak, and Maiwa yak).

Additionally, the genetic relationships between the zebu, yak, and gaur mirrored their geographical origins (Figure 1A). For zebu population, the PC1 accounted for 2.2% and PC2 explained 2.0% of the total variation. The PCA divided the zebu population into a distinct group formed by Myanmar individuals and a combined cluster including Indian and Bangladesh samples (Figure 1D). For the yak population, we found that PC1 and PC2 accounted for 2.7% and 2.1% of the total variation, respectively. The PC1 clearly separated the Zhongdian yak from the Tibetan yak and Maiwa yak (Figure 1E). These PCA results (Figures 1B,C) were strongly supported by the population clustering observed in a neighbor-joining tree (Figure 2).

FIGURE 2

FIGURE 2. Neighbor-joining (NJ) tree of 10 populations.

We performed the admixture analysis among populations based on 64,704 SNPs and obtained patterns of admixture for each population in our study (Figure 3). When K = 2, the yak population (“sky blue”) from southwestern China were first separated from other populations. Subsequently, the zebu population (“red”) was separated from other populations at K = 3. When K = 4, the Holstein population (“purple”) were separated from Diqing cattle and Tibetan cattle (“green”). The lowest cross-validation error value (0.35968) was detected at K = 5 (Supplementary Figure S1), all individuals formed five clusters, including Holstein as a cluster, Diqing cattle and Tibetan cattle from Southwest China as another cluster, the remaining zebu population, yak population, and gaur population (“grey”) each clustered together. These results were consistent with those of the PCAs (Figures 1B,C). Interestingly, we observed genetic components of the yak population may have introgressed into Diqing cattle and Tibetan cattle at K = 4–5, while the Zhongdian yak has genetic components of Diqing cattle/Tibetan cattle. Additionally, we found that three ancestral components of the European taurine (e.g., Holstein), East Asian taurine (e.g., Diqing cattle and Tibetan cattle), Indian-type zebu (Indian zebu, Bangladesh zebu, and Myanmar cattle) in domestic cattle by combining the PCA and NJ tree results. It should be noted that Indian zebu, Bangladesh zebu, and Myanmar cattle from different geographical regions formed a distinct cluster (Figure 3).

FIGURE 3

FIGURE 3. ADMIXTURE analyses of 10 populations (K = 2–5). The bovine population abbreviations are given in Table 1.

Diqing Cattle and Zhongdian Yak Have a Relevant Gene Flow

In order to further test the gene flow between Diqing cattle and yak. The f3 statistics analysis was performed by using Diqing cattle as a target and Holstein (or Tibetan cattle) and other Bovine species as the source population. When using the Zhongdian yak and Holstein as the sources has a significant z-score (−5.754), showing that there is gene flow between Diqing cattle and Zhongdian yak (Supplementary Table S1). Additionally, the D statistics was performed following the tree topology [[[Holstein, Diqing cattle], Zhongdian yak], Gaur] and found a significant z-score (−13.164) (Supplementary Table S2), displaying that Zhongdian yak introgression into Diqing cattle. TreeMix analysis (m = 5) also confirmed the gene flow from the Zhongdian yak to Diqing cattle (Supplementary Figure S2).

The same methods (f3 statistics, D statistics, and TreeMix analyses) were used to detect the genetic introgression from Diqing cattle to the Zhongdian yak. When using the Tibetan yak (or Maiwa yak) and Diqing cattle as the source populations, it had a significant z-score (−17.987/−14.108) for f3 statistics (Supplementary Table S3). Furthermore, the D statistics was executed following the tree topology [[[Tibetan yak, Zhongdian yak], Diqing cattle], Gaur] and had a significant z-score (−15.981) (Supplementary Table S4), showing that Diqing cattle possibly had introgressed into the Zhongdian yak. The TreeMix results (m = 4) also confirmed the gene flow from the Diqing cattle to the Zhongdian yak (Supplementary Figure S2). The bidirectional genetic introgression between the Diqing cattle and Zhongdian yak was detected based on aforementioned research results.

Diqing Cattle Has the Genetic Components of Zhongdian Yak

We selected Diqing cattle as the target population, and the Zhongdian yak and Holstein as the reference populations to detect the introgressed proportion from the Zhongdian yak into Diqing cattle using PCAdmix v1.0. We found that the average introgressed proportion was 2.37% (0.84%–7.68%) from the Zhongdian yak into Diqing cattle (Supplementary Table S5). The average introgressed proportion was 5.04% (0.83%–21.25%) from Diqing cattle into the Zhongdian yak (Supplementary Table S6). To identify the biological function of the introgression regions from the Zhongdian yak into Diqing cattle, we exploited the gene of 402 introgressed fragments that were shared by at least two haplotypes in Diqing cattle. The introgressed genes were enriched for GO terms and KEGG pathways (Supplementary Table S7, p < 0.01), the results displayed that the GO terms mainly include CCR chemokine receptor binding (GO:0048020), chemokine activity (GO:0008009), inflammatory response (GO:0006954), and lymphocyte chemotaxis (GO:0048247). The KEGG pathways were mainly involved in the chemokine signaling pathway (bta04062) and leukocyte transendothelial migration (bta04670), and these terms and pathways are associated with immunity. Furthermore, to test whether the genetic introgression of the Zhongdian yak influences Diqing cattle to high-altitude and hypoxia adaptations, we found that some key genes related to high-altitude and hypoxia adaptations, such as ENPEP (Eichstaedt et al., 2015; Wang et al., 2016; Wei et al., 2016), FLT1 (Kraft et al., 2016), and PIK3CA (Shimomura et al., 2013).

Tibetan Cattle Have the Genetic Components of Zhongdian Yak

Tibetan cattle mainly live in Tibetan areas with an altitude of 2,300–3,800 m. Studies have shown that Tibetan cattle have obtained gene alleles from the yak through genetic introgression (Chen et al., 2018; Wu et al., 2018). To further detect whether there was genetic introgression from the yak into Tibetan cattle, we performed the f3 statistics, D statistics, and TreeMix analyses. When using the Zhongdian yak and Holstein as the source populations has a significant z-score (−6.877) for f3 statistics, it shows that there is gene flow between Tibetan cattle and Zhongdian yak (Supplementary Table S8). The D statistics was performed following the tree topology [[[Holstein, Tibetan cattle], Zhongdian yak], Gaur] and found a significant z-score (−10.250) (Supplementary Table S9), displaying that Zhongdian yak had introgressed into Tibetan cattle. Furthermore, the TreeMix analysis (m = 5) also confirmed the gene flow from the Zhongdian yak to Tibetan cattle (Supplementary Figure S2).

In this study, we selected the Tibetan cattle (39 individuals) as target population and Holstein (34 individuals) of lowland as reference population. The F_ST, XP-EHH, and XP-CLR methods were used to detect selection signatures on the Tibetan cattle genome. Manhattan plots of the genome-wide F_ST, XP-EHH, and XP-CLR analyses are depicted in Figure 4. Finally, 270, 162, and 296 genes were identified from the F_ST, XP-EHH, and XP-CLR methods, respectively. The 30 common genes were obtained in at least two methods, of which CDC42 and SLC39A2 genes were found to be associated with high-altitude adaptation. The EPAS1 related to high-altitude adaptation was detected in the separate XP-EHH method. Additionally, we executed the GO and KEGG analyses for the 30 common genes (Table 3) and found that the GO terms were mainly enriched in chemotaxis (GO:0006935) and small GTPase-mediated signal transduction (GO:0007264). The KEGG pathway mainly enriched in chemokine signaling pathway (bta04062) (p < 0.05).

FIGURE 4

FIGURE 4. Genome-wide distribution of F_ST, XP-EHH, and XP-CLR of Tibetan cattle vs. Holstein on autosomes. The horizontal black line in the graph indicates the values of the top 1%.

TABLE 3

TABLE 3. Enrichment pathways of the 30 common genes in Tibetan cattle vs. Holstein.

Discussion

In this study, a total of 266 individuals of 10 populations of four Bos species were used to explore their population structure and genetic introgression. The PCA showed that all individuals formed five clusters (Figures 1B,C), and the genetic relationship among other three Bos species mirrored their geographic origins except taurine based on the 64,704 SNP data sets (Figure 1A). Furthermore, the PCA and a neighbor-joining (NJ) tree showed that the Indian zebu and Bangladesh zebu and Tibetan yak and Maiwa yak were clustered into a mixed cluster, respectively. It is possible that their geographical connection leads to frequent gene flow. The admixture analysis displayed that all samples formed five clusters, which were consistent with the PCA results, illustrating the reliability of clustering results. We observed European taurine, East Asian taurine, and Indian zebu ancestral components in domestic cattle by performing the admixture analysis; the result was similar to a previous study based on the whole-genome re-sequencing analysis (Chen et al., 2018), and it further showed that the genome-wide SNP chip is reliable for the classification of ancestral components. Studies have shown that zebu may have dispersed from the Indian subcontinent to East Asia about 3,500–2,500 years before present (YBP) (Chen S. et al., 2010), passed Southeast Asia (e.g., Myanmar), and entered into southwestern China (Higham, 1996). Therefore, it can be used to explain that Indian zebu, Bangladesh zebu, and Myanmar cattle from different geographical regions gathered in a separate cluster.

Bos species can provide production and living resources for human beings, including meat, milk, leather, and fuel. Thus, human migration also drove the continuous migration of different Bos species in history, leading to frequent gene flow among these bovine populations. For example, genetic introgression from the yak into Tibetan cattle (Chen et al., 2018; Wu et al., 2018), Bali cattle introgressed into Chinese Hainan and Luxi cattle breeds (Decker et al., 2014), and domestic cattle (taurine and zebu) introgressed into the yak (Medugorac et al., 2017). These studies demonstrated that genetic introgression can help domestic cattle better adapt to extreme environments (e.g., tropical climate and high-altitude environments). Previous studies have shown that genetic component of the yak introgressed into Diqing cattle based on mtDNA data (Nie et al., 1999; Yu et al., 1999). In this study, both Diqing cattle and Zhongdian yak were collected from Shangri-La, Diqing Prefecture. The bidirectional genetic introgression between Diqing cattle and Zhongdian yak was detected, and this result was similar to previously reported results based on the genetic introgression between Tibetan cattle and yak (Chen et al., 2018; Wu et al., 2018). Simultaneously, the enrichment analysis showed that the introgressed genes from the Zhongdian yak to Diqing cattle were mainly enriched in immune-related pathways, being different from those of previous studies (Chen et al., 2018). Their studies also detected some pathways related to diseases, speculating that the genes related to the disease were not subject to strong selection. Three introgressed genes (ENPEP, FLT1, and PIK3CA) related to high-altitude and hypoxia adaptations were detected in Diqing cattle, of which the ENPEP gene has been reported in high-altitude adaptability studies of the Andean (Eichstaedt et al., 2015) and Tibetan goats (Wang et al., 2016) and Tibetan sheep (Wei et al., 2016). The previously reported common EGLN1 gene (Chen et al., 2018; Wu et al., 2018) was not detected in our study. Here, the PIK3CA gene was an important component of the hypoxia-inducible factor 1 (HIF-1) pathway (Shimomura et al., 2013). These genes detected may have contributed to the adaptation of Diqing cattle to the high-altitude and hypoxia environment.

The introgression event was also found between the Zhongdian yak and Tibetan cattle in our study, being similar to the previous study on genetic components from the yak introgressed into Tibetan cattle (Chen et al., 2018; Wu et al., 2018). The 30 common genes under positive selection were detected in Tibetan cattle, which were enriched in the chemokine signaling pathway, indicating that immune-related pathways probably play an important role in resisting diseases and maintaining energy balance in high-altitude environments. Interestingly, the EPAS1 gene was detected by using the XP-EHH method, which was congruent with the previous study using F_ST, XP-EHH, and ΔDAF methods based on re-sequencing technology (Wu et al., 2019). Studies showed that the EPAS1 gene was under positive selection and had undergone convergent evolution in Tibetan dogs, Tibetans, and Tibetan pigs (Li et al., 2014; Wang et al., 2014; Wu et al., 2019). Additionally, the CDC42 and SLC39A2 genes were also found in Tibetan cattle. Among CDC42 related to the morphology, adhesion, migration, and erythropoiesis of hematopoietic stem cells (Wang et al., 2006; Yang et al., 2007), this gene can mediate bmp-induced sprouting angiogenesis (Wakayama et al., 2015). CDC42 was under positive selection and enriched with vascular endothelial growth factor (VEGF) for Tibetan goats (Jin et al., 2020). SLC39A2 and SLC9A6 genes were related to heart functions (Du et al., 2019), which have also been reported in Yunnan golden monkeys (Zhou et al., 2016). It is speculated that the regulation of genes related to hypoxia response, vascular development and erythropoiesis, and heart functions may have extremely important significance for the high-altitude adaptation of Tibetan cattle by the aforementioned comparative studies.

In conclusion, we used the high-density BovineHD genotyping arrays to assess the population structure and genetic introgression of Bos species. The genetic relationships between the zebu, yak, and gaur mirrored their geographical origins. Furthermore, it confirmed that there were at least three ancestral components of the European taurine, East Asian taurine, and Indian indicine in domestic cattle. We found the bidirectional genetic introgression between Diqing cattle and Zhongdian yak. Especially, the Zhongdian yak introgression into Diqing cattle will help the latter to quickly adapt to the high-altitude and hypoxia environment. Furthermore, we also found that the genetic components of the Zhongdian yak had introgressed into Tibetan cattle and detected some genes (CDC42, SLC39A2, and EPAS1) associated with high-altitude adaptation by using three methods of selection signatures (F_ST, XP-EHH, and XP-CLR). The functions of these genes might have great significance for high-altitude adaptation of Tibetan cattle. These results showed that genetic introgression was one of the important ways for the environmental adaptation of domestic cattle. Future studies are needed to validate these results by using whole-genome re-sequencing and functional verification studies.

Data Availability Statement

The data sets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://datadryad.org/stash/dataset/doi:10.5061/dryad.h9w0vt4k8.

Ethics Statement

The animal study was reviewed and approved by the Committee on Animal Research and Ethics of Yunnan University (ynucae20190011). Written informed consent for participation was not obtained from the owners because ear skin tissue samples were collected by veterinary practitioners with the permission and in presence of the owners. Veterinarians followed standard procedures and relevant national guidelines to ensure appropriate animal care.

Author Contributions

SC and AB-P conceived and designed the research. SC, CL, MSB, and MB collected samples. RL and VC extracted genomic DNA. RL and SC performed bioinformatics analyses. RL, SC, HX, and AB-P wrote and revised the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by National Key R&D Program of China (No. 2021YFD1200903), National Natural Science Foundation of China (No. 32160146 and No. 31460281), and Top Young Talents Program of Ten Thousand Plan of Yunnan Province (YNWRQNBJ-2018-024).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank local farmers and veterinarians for their kind supports in collecting samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.847492/full#supplementary-material

References

Akey, J. M., Zhang, G., Zhang, K., Jin, L., and Shriver, M. D. (2002). Interrogating a High-Density SNP Map for Signatures of Natural Selection. Genome Res. 12, 1805–1814. doi:10.1101/gr.631202

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast Model-Based Estimation of Ancestry in Unrelated Individuals. Genome Res. 19, 1655–1664. doi:10.1101/gr.094052.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbato, M., Reichel, M. P., Passamonti, M., Low, W. Y., Colli, L., Tearle, R., et al. (2020). A Genetically Unique Chinese Cattle Population Shows Evidence of Common Ancestry with Wild Species when Analysed with a Reduced Ascertainment Bias SNP Panel. PloS One 15, e0231162. doi:10.1371/journal.pone.0231162

PubMed Abstract | CrossRef Full Text | Google Scholar

Brisbin, A., Bryc, K., Byrnes, J., Zakharia, F., Omberg, L., Degenhardt, J., et al. (2012). PCAdmix: Principal Components-Based Assignment of Ancestry along Each Chromosome in Individuals with Admixed Ancestry from Two or More Populations. Hum. Biol. 84, 343–364. doi:10.3378/027.084.0401

PubMed Abstract | CrossRef Full Text | Google Scholar

Browning, B. L., and Browning, S. R. (2013). Improving the Accuracy and Efficiency of Identity-By-Descent Detection in Population Data. Genetics 194, 459–471. doi:10.1534/genetics.113.150029

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: Rising to the challenge of Larger and Richer Datasets. GigaSci 4, 7. doi:10.1186/s13742-015-0047-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Patterson, N., and Reich, D. (2010). Population Differentiation as a Test for Selective Sweeps. Genome Res. 20, 393–402. doi:10.1101/gr.100545.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, N., Cai, Y., Chen, Q., Li, R., Wang, K., Huang, Y., et al. (2018). Whole-genome Resequencing Reveals World-wide Ancestry and Adaptive Introgression Events of Domesticated Cattle in East Asia. Nat. Commun. 9, 2337. doi:10.1038/s41467-018-04737-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Lin, B.-Z., Baig, M., Mitra, B., Lopes, R. J., Santos, A. M., et al. (2010). Zebu Cattle Are an Exclusive Legacy of the South Asia Neolithic. Mol. Biol. Evol. 27, 1–6. doi:10.1093/molbev/msp213

PubMed Abstract | CrossRef Full Text | Google Scholar

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., et al. (2011). The Variant Call Format and VCFtools. Bioinformatics 27, 2156–2158. doi:10.1093/bioinformatics/btr330

PubMed Abstract | CrossRef Full Text | Google Scholar

Decker, J. E., McKay, S. D., Rolf, M. M., Kim, J., Molina Alcalá, A., Sonstegard, T. S., et al. (2014). Worldwide Patterns of Ancestry, Divergence, and Admixture in Domesticated Cattle. Plos Genet. 10, e1004254. doi:10.1371/journal.pgen.1004254

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, L., Zhang, H., Zhao, H., Cheng, X., Qin, J., Teng, T., et al. (2019). The Critical Role of the Zinc Transporter Zip2 (SLC39A2) in Ischemia/reperfusion Injury in Mouse Hearts. J. Mol. Cell Cardiol. 132, 136–145. doi:10.1016/j.yjmcc.2019.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Edea, Z., Bhuiyan, M. S. A., Dessie, T., RothschildKim, M. F., Dadi, H., and Kim, K. S. (2015). Genome-wide Genetic Diversity, Population Structure and Admixture Analysis in African and Asian Cattle Breeds. Animal 9, 218–226. doi:10.1017/s1751731114002560

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichstaedt, C. A., Antão, T., Cardona, A., Pagani, L., Kivisild, T., and Mormina, M. (2015). Genetic and Phenotypic Differentiation of an Andean Intermediate Altitude Population. Physiol. Rep. 3, e12376. doi:10.14814/phy2.12376

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Gautier, M., Ding, X., Zhang, H., Wang, Y., Wang, X., et al. (2017). Species Composition and Environmental Adaptation of Indigenous Chinese Cattle. Sci. Rep. 7, 16196. doi:10.1038/s41598-017-16438-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gautier, M., Moazami-Goudarzi, K., Levéziel, H., Parinello, H., Grohs, C., Rialle, S., et al. (2016). Deciphering the Wisent Demographic and Adaptive Histories from Individual Whole-Genome Sequences. Mol. Biol. Evol. 33, 2801–2814. doi:10.1093/molbev/msw144

PubMed Abstract | CrossRef Full Text | Google Scholar

Higham, C. (1996). The Bronze Age of Southeast Asia. Cambridge: Cambridge University Press.

Google Scholar

Hu, Q., Ma, T., Wang, K., Xu, T., Liu, J., and Qiu, Q. (2012). The Yak Genome Database: an Integrative Database for Studying Yak Biology and High-Altitude Adaption. BMC Genomics 13, 600. doi:10.1186/1471-2164-13-600

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, M., Lu, J., Fei, X., Lu, Z., Quan, K., Liu, Y., et al. (2020). Selection Signatures Analysis Reveals Genes Associated with High-Altitude Adaptation in Tibetan Goats from Nagqu, Tibet. Animals 10, 1599. doi:10.3390/ani10091599

PubMed Abstract | CrossRef Full Text | Google Scholar

Kijas, J. W., Lenstra, J. A., Hayes, B., Boitard, S., Porto Neto, L. R., San Cristobal, M., et al. (2012). Genome-wide Analysis of the World's Sheep Breeds Reveals High Levels of Historic Mixture and strong Recent Selection. Plos Biol. 10, e1001258. doi:10.1371/journal.pbio.1001258

PubMed Abstract | CrossRef Full Text | Google Scholar

Kraft, B. D., Suliman, H. B., Colman, E. C., Mahmood, K., Hartwig, M. G., Piantadosi, C. A., et al. (2016). Hypoxic Gene Expression of Donor Bronchi Linked to Airway Complications after Lung Transplantation. Am. J. Respir. Crit. Care Med. 193, 552–560. doi:10.1164/rccm.201508-1634oc

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Stecher, G., and Tamura, K. (2016). MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol. Biol. Evol. 33, 1870–1874. doi:10.1093/molbev/msw054

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachance, J., and Tishkoff, S. A. (2013). SNP Ascertainment Bias in Population Genetic Analyses: Why it Is Important, and How to Correct it. Bioessays 35, 780–786. doi:10.1002/bies.201300014

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Wu, D.-D., Boyko, A. R., Wang, G.-D., Wu, S.-F., Irwin, D. M., et al. (2014). Population Variation Revealed High-Altitude Adaptation of Tibetan Mastiffs. Mol. Biol. Evol. 31, 1200–1205. doi:10.1093/molbev/msu070

PubMed Abstract | CrossRef Full Text | Google Scholar

López Herráez, D., Bauchet, M., Tang, K., Theunert, C., Pugach, I., Li, J., et al. (2009). Genetic Variation and Recent Positive Selection in Worldwide Human Populations: Evidence from Nearly 1 Million SNPs. PloS one 4, e7888. doi:10.1371/journal.pone.0007888

PubMed Abstract | CrossRef Full Text | Google Scholar

Medugorac, I., Graf, A., Grohs, C., Rothammer, S., Zagdsuren, Y., Gladyr, E., et al. (2017). Whole-genome Analysis of Introgressive Hybridization and Characterization of the Bovine Legacy of Mongolian Yaks. Nat. Genet. 49, 470–475. doi:10.1038/ng.3775

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukherjee, A., Mukherjee, S., Dhakal, R., Mech, M., Longkumer, I., Haque, N., et al. (2018). High-density Genotyping Reveals Genomic Characterization, Population Structure and Genetic Diversity of Indian Mithun (Bos frontalis). Sci. Rep. 8, 10316–10410. doi:10.1038/s41598-018-28718-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nie, L., Yu, Y., Zhang, X. Q., Yang, G. F., Wen, J. K., and Zhang, Y. P. (1999). Genetic Diversity of Cattle in South China as Revealed by Blood Protein Electrophoresis. Biochem. Genet. 37, 257–265. doi:10.1023/a:1018798924778

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, R. (2004). Population Genetic Analysis of Ascertained SNP Data. Hum. Genomics 1, 218–224. doi:10.1186/1479-7364-1-3-218

PubMed Abstract | CrossRef Full Text | Google Scholar

Patterson, N., Moorjani, P., Luo, Y., Mallick, S., Rohland, N., Zhan, Y., et al. (2012). Ancient Admixture in Human History. Genetics 192, 1065–1093. doi:10.1534/genetics.112.145037

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickrell, J. K., Pritchard, J. K., and Tang, H. (2012). Inference of Population Splits and Mixtures from Genome-wide Allele Frequency Data. Plos Genet. 8, e1002967. doi:10.1371/journal.pgen.1002967

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006). Principal Components Analysis Corrects for Stratification in Genome-wide Association Studies. Nat. Genet. 38, 904–909. doi:10.1038/ng1847

PubMed Abstract | CrossRef Full Text | Google Scholar

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am. J. Hum. Genet. 81, 559–575. doi:10.1086/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, Q., Zhang, G., Ma, T., Qian, W., Wang, J., Ye, Z., et al. (2012). The Yak Genome and Adaptation to Life at High Altitude. Nat. Genet. 44, 946–949. doi:10.1038/ng.2343

PubMed Abstract | CrossRef Full Text | Google Scholar

Rincon, G., Weber, K. L., Van Eenennaam, A. L., Golden, B. L., and Medrano, J. F. (2011). Hot Topic: Performance of Bovine High-Density Genotyping Platforms in Holsteins and Jerseys. J. Dairy Sci. 94, 6116–6121. doi:10.3168/jds.2011-4764

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothammer, S., Seichter, D., Förster, M., and Medugorac, I. (2013). A Genome-wide Scan for Signatures of Differential Artificial Selection in Ten Cattle Breeds. BMC Genomics 14, 908. doi:10.1186/1471-2164-14-908

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimomura, M., Hinoi, T., Kuroda, S., Adachi, T., Kawaguchi, Y., Sasada, T., et al. (2013). Overexpression of Hypoxia Inducible Factor-1 Alpha Is an Independent Risk Factor for Recurrence after Curative Resection of Colorectal Liver Metastases. Ann. Surg. Oncol. 20, 527–536. doi:10.1245/s10434-013-2945-2

CrossRef Full Text | Google Scholar

Szpiech, Z. A., and Hernandez, R. D. (2014). Selscan: An Efficient Multithreaded Program to Perform EHH-Based Scans for Positive Selection. Mol. Biol. Evol. 31, 2824–2827. doi:10.1093/molbev/msu211

PubMed Abstract | CrossRef Full Text | Google Scholar

Taye, M., Lee, W., Caetano-Anolles, K., Dessie, T., Cho, S., Jong Oh, S., et al. (2018). Exploring the Genomes of East African Indicine Cattle Breeds Reveals Signature of Selection for Tropical Environmental Adaptation Traits. Cogent Food Agric. 4, 1552552. doi:10.1080/23311932.2018.1552552

CrossRef Full Text | Google Scholar

Uzzaman, M. R., Edea, Z., Bhuiyan, M. S. A., Walker, J., Bhuiyan, A. K. F. H., and Kim, K.-S. (2014). Genome-wide Single Nucleotide Polymorphism Analyses Reveal Genetic Diversity and Structure of Wild and Domestic Cattle in Bangladesh. Asian Australas. J. Anim. Sci. 27, 1381–1386. doi:10.5713/ajas.2014.14160

PubMed Abstract | CrossRef Full Text | Google Scholar

Wakayama, Y., Fukuhara, S., Ando, K., Matsuda, M., and Mochizuki, N. (2015). Cdc42 Mediates Bmp-Induced Sprouting Angiogenesis through Fmnl3-Driven Assembly of Endothelial Filopodia in Zebrafish. Develop. Cel 32, 109–122. doi:10.1016/j.devcel.2014.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G.-D., Fan, R.-X., Zhai, W., Liu, F., Wang, L., Zhong, L., et al. (2014). Genetic Convergence in the Adaptation of Dogs and Humans to the High-Altitude Environment of the Tibetan Plateau. Genome Biol. Evol. 6, 2122–2128. doi:10.1093/gbe/evu162

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Yang, L., Filippi, M.-D., Williams, D. A., and Zheng, Y. (2006). Genetic Deletion of Cdc42GAP Reveals a Role of Cdc42 in Erythropoiesis and Hematopoietic Stem/progenitor Cell Survival, Adhesion, and Engraftment. Blood 107, 98–105. doi:10.1182/blood-2005-05-2171

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Liu, J., Zhou, G., Guo, J., Yan, H., Niu, Y., et al. (2016). Whole-genome Sequencing of Eight Goat Populations for the Detection of Selection Signatures Underlying Production and Adaptive Traits. Sci. Rep. 6, 38932. doi:10.1038/srep38932

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, C., Wang, H., Liu, G., Zhao, F., Kijas, J. W., Ma, Y., et al. (2016). Genome-wide Analysis Reveals Adaptation to High Altitudes in Tibetan Sheep. Sci. Rep. 6, 26770. doi:10.1038/srep26770

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. New York, NY: Springer.

Google Scholar

Wu, D.-D., Ding, X.-D., Wang, S., Wójcik, J. M., Zhang, Y., Tokarska, M., et al. (2018). Pervasive Introgression Facilitated Domestication and Adaptation in the Bos Species Complex. Nat. Ecol. Evol. 2, 1139–1145. doi:10.1038/s41559-018-0562-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, D.-D., Yang, C.-P., Wang, M.-S., Dong, K.-Z., Yan, D.-W., Hao, Z.-Q., et al. (2019). Convergent Genomic Signatures of High-Altitude Adaptation Among Domestic Mammals. Natl. Sci. Rev. 7, 952–963. doi:10.1093/nsr/nwz213

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Wang, L., Kalfa, T. A., Cancelas, J. A., Shang, X., Pushkaran, S., et al. (2007). Cdc42 Critically Regulates the Balance between Myelopoiesis and Erythropoiesis. Blood 110, 3853–3861. doi:10.1182/blood-2007-03-079582

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Y., Nie, L., He, Z. Q., Wen, J. K., Jian, C. S., and Zhang, Y. P. (1999). Mitochondrial DNA Variation in Cattle of south China: Origin and Introgression. Anim. Genet. 30, 245–250. doi:10.1046/j.1365-2052.1999.00483.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Meng, X., Liu, Z., Chang, J., Wang, B., Li, M., et al. (2016). Population Genomics Reveals Low Genetic Diversity and Adaptation to Hypoxia in Snub-Nosed Monkeys. Mol. Biol. Evol. 33, 2670–2681. doi:10.1093/molbev/msw150

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimin, A. V., Delcher, A. L., Florea, L., Kelley, D. R., Schatz, M. C., Puiu, D., et al. (2009). A Whole-Genome Assembly of the Domestic Cow, Bos taurus. Genome Biol. 10, R42. doi:10.1186/gb-2009-10-4-r42

PubMed Abstract | CrossRef Full Text | Google Scholar