Genomic analysis reveals population structure and selection signatures in plateau dairy cattle

Abstract

Background

To solve the problem of an insufficient supply of dairy products in Tibet, work has been carried out to improve native dairy cattle and introduce purebred dairy cattle from low-altitude areas. The harsh environment of the plateau not only severely limits the production performance of high-yielding dairy cattle, such as Holstein and Jersey cattle, but also challenges their survival. The population structure and plateau adaptation mechanism of plateau dairy cattle are rarely reported. In this study, key genes and pathways affecting plateau purebred and crossbred dairy cattle were explored using genetic chip information.

Results

The results showed that the genetic diversity of the Tibet dairy cattle population was higher than that of the native cattle and plains dairy cattle. Purebred Holstein and Jersey cattle in Tibet were genetically closer to dairy cattle in the plains, and crossbred dairy cattle were admixed with more Tibet cattle and Apaijiza cattle. Based on the fixation index (F_ST), integrated haplotype score (iHS), and cross-population extend haplotype homozygosity (XP-EHH) approaches, 60 and 40 genes were identified in plateau Holstein and Jersey cattle, respectively. A total of 78 and 70 genes were identified in crossbred cattle compared to Holstein and Tibet cattle respectively. These genes are related to cardiac health and development, neuronal development and function, angiogenesis and hematopoietic, pigmentation, growth and development, and immune response.

Conclusions

Our results provide a glimpse into diverse selection signatures in plateau dairy cattle, which can be used to enhance our understanding of the genomic basis of plateau adaptation in dairy cattle. These results support further research on breeding strategies such as marker-assisted selection and gene editing in plateau dairy cattle populations.

Introduction

The Qinghai‒Tibet Plateau is located in western China and has an average elevation of more than 4,000 m above sea level. This high-altitude area is characterized by extreme environmental conditions, including low oxygen concentrations, alpine cold temperatures and high intensities of ultraviolet (UV) radiation, which restrict the survival of domestic animals [1]. Milk is indispensable for the Tibetan people, and due to the insufficient lactation performance of yaks and local cattle, a considerable portion of dairy products need to be transferred from the mainland [2]. The Tibet Autonomous Region has been utilizing Holstein cows for local cattle improvement since the early 1980s, and the limited stock of crossbred cows is still unable to meet the large demand for dairy products in Tibet [3, 4]. A large number of good-breed dairy cattle have been introduced to Tibet since 2014, the mortality of Holstein is 26.5% and the mortality of Jersey is 8%. They have not demonstrated long-term adaptability to the plateau environment in terms of either disease resistance or production performance [5], to the extent that the original breeding objectives and appointment requirements could not be met. Therefore, it is urgent to make it clear that how plateau dairy cattle adapt to plateau environment and improve their performance through genetic breeding.

Plateau adaptation is considered a complex trait that is influenced by multiple genes [6, 7]. A series of genes with functions related to energy metabolism, oxidative response, stress response, hypoxic response, temperature acclimatization, cardiovascular system and body morphology were considered as candidate genes in studies on selection signatures for high-altitude acclimatization [8,9,10]. Several studies have used genomic data from different breeds of cattle and yaks in different regions to analyze genetic variation resulting from selection at high altitudes, to explore the mechanisms of plateau adaptation and to search for genetic markers [10,11,12,13,14]. Since plateau dairy cattle are still undergoing plateau habituation and such groups are relatively rare, the analysis of their population structure and selection signatures will differ from previous studies.

In this study, the SNP chip data was used to analyze the genetic diversity and population structure of the populations in two large-scale dairy farms in Lhasa, to determine the genetic status of dairy cattle populations and to provide a research basis for genetic improvement and breeding planning for plateau dairy cattle. Selection signatures analysis was also used to detect genomic differences between purebred dairy cattle at different altitudes, as well as genomic differences between hybrid cattle and possible parent breeds. The aim was to reveal the genes that were selected after several generations of introduction of dairy cattle to the plateau, as well as the genes carried in the crossbred cows that differed from those of the parents. The present results can be useful for future plateau adaptation studies and genetic improvement of dairy cattle in Tibet.

Materials and methods

Animal information and genomic data

The plateau dairy cattle came from Gaba Ecological Farm (Tibet, China) and Zhizhao Industrial Park’s high-standard dairy farm (Tibet, China). Two plains dairy cattle populations used for the comparison were from Shounong Animal Husbandry Jingyindao Farm (Beijng, China) and Hebei Shounong Modern Agricultural Science and Technology Co. (Hebei, China). Three plateau native cattle populations used for comparison came from herdsmen’s homes in different parts of Tibet. About 2 mL of blood was collected from the tail-heel vein of each cattle and Blood Cards were submitted to Neogen Bio-Scientific Technology (Shanghai) Co., Ltd. (Shanghai, China) for processing. This study did not involve the euthanasia or sacrifice of animals.

DNA from each blood card sample was extracted using the Sbeadex Livestock Nucleic Acid Livestock Extraction Kit (LGC Biosearch Technologies, Petaluma, CA). DNA samples were quantified using NanoDrop™ 8000 (Thermo Scientific, Waltham, MA) and genotyped according to the Illumina Infinium XT protocol (Illumina, San Diego, CA). SNP markers were scored on the Bovine Geneseek Genomic Profiler –100 K Beadchip or Bovine Geneseek Genomic Profiler –150 K Beadchip (Neogen Inc, Lincoln, NE) at Neogen. The beadchips were scanned using the Illumina IScan (Illumina) and the raw data was processed and exported using Genome Studio 2.0 software according to the Genome Studio Framework User Guide (Illumina) based on selection criteria of at least 90% animal call rate.

Because accurate pedigrees were not available and crossbred dairy cattle were indistinguishable from Holsteins in coat color, we obtained the bloodline composition calculated from chip data at Neogen. The percentage of DNA contributed to the animal by each of 5 evaluated breeds: Holstein, Jersey, Brown Swiss, Ayrshire, and Guernsey was estimated and termed as Breed base representation (BBR) [15]. BBR values represent the breed composition of an individual that is more accurate and much easier to interpret than breed-check markers. A BBR value of 94 to 99% is typically set to 100%, which representing animals with 100% purebred ancestry. Therefore, the plateau dairy cattle were categorized to 5 subpopulations according to the BBR values: 119 crossbred dairy cattle with 40% ~ 74% Holstein bloodlines (GH40), 92 dairy cattle with 75% ~ 84% Holstein bloodlines (GH75), 92 dairy cattle with 85% ~ 93% Holstein bloodlines (GH85), 129 Holstein cattle with 94% ~ 100% Holstein bloodlines (GH100) and 28 Jersey cattle with 76% ~ 100% Jersey bloodlines (GJ). The genomic data used in this study were shown in Table 1.

Table 1 Sources of experimental animal and genetic data

The liftover tool [16] was used to update the reference genome version from Bos_taurus_UMD_3.1 to ARS-UCD1.2 for the 150 K chip. Data from chips of the same density were merged and subjected to quality control, including (1) excluding individuals with a genotyping detection rate less than 0.9; (2) excluding SNPs with a genotyping detection rate less than 0.9; (3) excluding SNPs with a minor allele frequency (MAF) less than 0.01; (4) excluding sites severely deviating from Hardy‒Weinberg equilibrium (P < 10^–6); and (5) retaining SNPs located on autosomes and the X chromosome. The Beagle (v5.2) [17] tool was used to impute missing genotypes in the quality-controlled data. After retaining common sites between the two versions of the chips, a total of 752 cattle with 62,422 SNPs were obtained for subsequent analysis.

Genetic diversity and linkage disequilibrium

The average MAF and observed heterozygosity rate (HET) were estimated using PLINK (v1.90) [18]. The parameters for detecting Runs of homozygosity (ROH) were set as follows: for each 70 kb segment of ROH, there must be one SNP; the interval between consecutive SNPs in the same ROH must be less than 500 kb; one heterozygous locus is allowed; only ROH with a length exceeding 30 SNPs and greater than 500 kb are detected. The ratio of ROH length to chromosome length was taken as inbreeding coefficient (F_ROH).To estimate and compare decay patterns of linkage disequilibrium (LD) at the genome-wide level among populations and subpopulations, PopLDdecay v.3.40 [19] was used to calculate the square correlation coefficient (r²) between SNP pairs up to 500 kb. The effective population size 13 generations ago (N_E13) was estimated from LD by using the SNeP tool [20].

Population structure

The genetic distance based on identity-by-state (IBS) was calculated using PLINK (v1.90), and multidimensional scaling (MDS) was carried out to show the populations and subpopulations structure. The genetic distance matrix between individuals, computed through VCF2Dis (v1.47) [21], was uploaded to FastME (v2.0) [22] for the construction of a neighbor-joining evolutionary tree, which was ultimately visualized using iTOL (v6) [23]. Base on the LD decay result, PLINK (v1.90) was submitted to use the command “– indep-pairwise 50 10 0.2” to remove linkage sites, followed by the use of ADMIXTURE (v1.3.0) [24] to conduct population structure analysis under the assumption of ancestral numbers k = 1 ~ 20.

Selection signature detection

Detection of selection signatures in the plateau dairy cattle genome was based either on the evaluation of allele frequency differences among populations and on properties of haplotypes segregating in the populations. The applied methods included within-population (integrated haplotype score, iHS) and between-population (fixation index, F_ST and cross-population extend haplotype homozygosity, XP-EHH) tests. In comparisons between the same breeds at different altitudes, the within population signatures of selection in GH100 and GJ were computed using iHS, and between-population tests were applied to identify potential sweeps that occurred in the GH100 compared to PH and GJ compared to PJ. In comparisons between crossbred dairy cattle and their possible parent breeds, the within population signatures of selection in GH40 were computed using iHS, and between-population tests were applied to identify potential sweeps that occurred in the GH40 compared to PH and XZ.

F_ST detection for each SNP was performed using vcftools (v0.1.16) [25]. After phasing with Beagle (v5.2) [17], iHS and XP-EHH was calculated using selscan (v1.3.0) [26], and the XP-EHH values were standardized using the norm in selscan (v1.3.0). Large positive or negative iHS values indicate unusually long haplotypes carrying the ancestral or derived alleles, respectively. SNP iHS scores were further transformed into piHS = -log[1–2|Φ(iHS)-0.5|] whereΦ(x) represent the Gaussian cumulative distribution function [27]. Assuming iHS are normally distributed (under neutrality), piHS might thus be interpreted as log10(1⁄P) where P is the two-sided P-value associated to the neutral hypothesis (no selection). Negative XP-EHH scores suggest selection happened in the reference population, otherwise happened in the observed population. The XP-EHH is transformed into piXP-EHH = -log[1–2|Φ(XP-RHH)-0.5|] where Φ(x) represents the Gaussian cumulative distribution function [27]. Assuming XP-EHH are normally distributed (under neutrality), piXP-EHH may be interpreted as a two-sided P-value in a -log10 scale. In this study, sites within the top 0.1% of the F_ST detection results and P < 0.001 (0.1%) of the iHS and XP-EHH detection results were defined as selected or differential sites [28,29,30,31].

Gene annotation and functional enrichment

GALLO (v1.1) [32] was used for gene annotation within 200 kb upstream and downstream of selected sites. To gain a better understanding of the biological processes (BP) shared by these candidate genes, we also performed functional enrichment analysis using DAVID [33] and the ‘clusterProfiler’ package [34].

Results

Genetic diversity and linkage disequilibrium

In this study, the genetic diversity indicators were calculated for 10 populations and subpopulations, and the results are shown in Table 2. Compared with those of other breeds, the Holstein cattle populations exhibited higher heterozygosity and MAF, with heterozygosity consistently above 0.4. Jersey cattle showed the highest level of inbreeding, followed by plains Holstein cattle. GH100 and GJ showed large differences in indicators of genetic diversity. The heterozygosity, MAF and F_ROH of GJ were closer to those of PJ. As the proportion of Holstein ancestry decreased, the MAF increased, the F_ROH decreased and the N_E13 increased. GH40 and GH75 showed higher genetic diversity than plateau native cattle. Among the plateau native cattle, TF had the lowest heterozygosity, MAF, and the highest F_ROH. XZ had a larger N_E13 compared to AP and TF.

Table 2 Statistics for the studied populations and subpopulations calculated on the basis of the SNP genotypes

Figure 1 shows the LD decay diagram within 500 kb of each population and subpopulations. The LD levels of Jersey cattle in the plains and plateaus were higher than that of other populations and subpopulations. The LD level of PH is higher than that of native cattle breeds in Tibet and Holstein cattle with different bloodline purities in Tibet. The LD decay rates of XZ and GH40 were similar and the fastest. The LD decay rate of different bloodline purities Holstein cattle subpopulations in Tibet increased as the proportion of Holstein ancestry decreased.

Fig. 1

LD decay plot for each population and subpopulation. GH100 is Holstein cattle with 94% ~ 100% Holstein bloodlines; GH85 is dairy cattle with 85% ~ 93% Holstein bloodlines; GH75 is dairy cattle with 75% ~ 84% Holstein bloodlines; GH40 is crossbred dairy cattle with 40% ~ 74% Holstein bloodlines; GJ is Jersey cattle with 76% ~ 100% Jersey bloodlines; PH is Holstein cattle in Beijing; PJ is Jersey cattle in Hebei; XZ is Tibet cattle; AP is Apeijia cattle; TF is Shigatse hump cattle

Population genetic structure

The MDS analysis based on the IBS distance is presented in Fig. 2(a). The first dimension separates Tibetan native cattle breeds from other Holstein and Jersey cattle populations, with XZ being genetically closest to other Holstein cattle populations, while the TF exhibits the farthest genetic distance from other Holstein cattle populations. Holstein cattle of different bloodline purities in Tibet are distributed along the first dimension according to the proportion of Holstein ancestry, with GH100 having the closest genetic distance to PH and GH40 having the closest genetic distance to the native Tibetan cattle breeds. The second dimension separates Jersey cattle, while PJ and GJ do not exhibit distinct clustering. Five individuals from the GH40 subpopulation are positioned between the Holstein and Jersey cattle populations.

Fig. 2

Population structure, as revealed by Multidimensional scaling analysis (a) and Neighborhood phylogenetic tree (b). GH100 is Holstein cattle with 94% ~ 100% Holstein bloodlines; GH85 is dairy cattle with 85% ~ 93% Holstein bloodlines; GH75 is dairy cattle with 75% ~ 84% Holstein bloodlines; GH40 is crossbred dairy cattle with 40% ~ 74% Holstein bloodlines; GJ is Jersey cattle with 76% ~ 100% Jersey bloodlines; PH is Holstein cattle in Beijing; PJ is Jersey cattle in Hebei; XZ is Tibet cattle; AP is Apeijia cattle; TF is Shigatse hump cattle

The neighbor-joining phylogenetic tree constructed based on genetic distance (Fig. 2(b)) divides 752 cattle into three distinct branches. PJ and GJ clustered in the second branch. PH, most of the GH100 strains, a few of the GH85 strains and a few of the GH75 strains clustered in the first branch. Other cattle were distributed in the third branch, with TF being more distant from the other populations, and XZ and AP clustered together and closer to TF. Most of the GH40s were closer to Tibetan native cattle, and GH85 and GH75 did not show significant clustering.

In the population structure analysis (Fig. 3), three native Tibetan cattle breeds were separated from Holstein and Jersey cattle at k = 2, and Jersey cattle were clearly demarcated at k = 3. TF cattle were clearly demarcated at k = 4, and Holstein cattle of different bloodline purities in Tibet were admixed with more XZ and AP bloodlines as the proportion of Holstein bloodlines decreased. Ancestry of Holstein cattle was added at k = 5, and at k = 6 TF cattle appeared to have ancestry different from AP and XZ. XZ and AP showed a clear difference in ancestry by k = 7. GJ and PJ both had good consistency in terms of pedigree composition from k = 2 to 7. The CV error continued to decrease from k = 1 to k = 20 (Additional file 1:Table S1).

Fig. 3

Admixture plot represents the cluster structure of the studied populations and subpopulations for the number of clusters (K) 2 to 7. GH100 is Holstein cattle with 94% ~ 100% Holstein bloodlines; GH85 is dairy cattle with 85% ~ 93% Holstein bloodlines; GH75 is dairy cattle with 75% ~ 84% Holstein bloodlines; GH40 is crossbred dairy cattle with 40% ~ 74% Holstein bloodlines; GJ is Jersey cattle with 76% ~ 100% Jersey bloodlines; PH is Holstein cattle in Beijing; PJ is Jersey cattle in Hebei; XZ is Tibet cattle; AP is Apeijia cattle; TF is Shigatse hump cattle

Selection signatures based on dairy cattle in plains and plateaus

The genome-wide distribution of |iHS| values in the GH100 is shown in Fig. 4(a). The Manhattan plots of pairwise XP-EHH and F_ST analysis between GH100 and PH are presented in Fig. 4(b-c). The mean F_ST value of GH100 in comparison to PH was 0.0123. A total of 79 significant SNPs (p < 0.001) distributed across 23 chromosomes were identified in iHS analysis (Additional file 2: Table S2). BTA3 and BTA9 were the chromosomes harboring the largest number of significant |iHS| markers (10 and 8, respectivel). The XP-EHH selection signals were detected on 20 chromosomes (Additional file 2: Table S3), and positive XP-EHH scores (indicating a selection on the GH100 genome) are observed on BTA2, BTA4, BTA6, BTA7 and BTA20. F_ST signals were identified on 19 chromosomes (Additional file 2: Table S4). The regions on the chromosomes having highest F_ST values were located on BTA13 (0.30), BTA1 (0.27), BTA30 (0.26) and BTA15 (0.26), respectively.

Fig. 4

a iHS plot for the GH100. b XP-EHH plot for GH100 comparison with PH. c F_ST plot for GH100 comparison with PH. d iHS plot for the GJ. e XP-EHH plot for GJ comparison with PJ. f F_ST plot for GJ comparison with PJ. GH100 is Holstein cattle in Tibet; PH is Holstein cattle in Beijing; GJ is Jersey cattle in Tibet; PJ is Jersey cattle in Hebei

The genome-wide distribution of |iHS| values in the GJ is shown in Fig. 4(d). The Manhattan plots of pairwise XP-EHH and F_ST analysis between GJ and PJ are presented in Fig. 4(e–f). The mean F_ST value of GJ in comparison to PJ was 0.0117. A total of 107 significant SNPs (p < 0.001) distributed across all chromosomes except BTA25, BTA27 and BTA29 were identified in iHS analysis (Additional file 2: Table S5). The XP-EHH selection signals were detected on 19 chromosomes (Additional file 2: Table S6), and positive XP-EHH scores (indicating a selection on the GJ genome) are observed on BTA8, BTA13, BTA16, BTA20, BTA23 and BTA29. BTA5 and BTA7 were the chromosomes harboring the largest number of these SNPs (83 and 61, respectively), followed by BTA10 (36 SNPs). F_ST signals were identified on 20 chromosomes (Additional file 2: Table S7). The regions on the chromosomes having highest F_ST values were located on BTA16 (0.35), BTA18 (0.34), and BTA27 (0.32), respectively.

The top 10 SNPs with their iHS values or positive XP-EHH values and annotated genes are shown in Table 3. The enriched biological processes were shown in Additional file 3: Table S8—S9. Among the top 10 |iHS| markers in GH100, 3 are located on BTA3, in which annotated STYXL2. STYXL2 (BTA3, iHS = 4.10), ALPL(BTA2, XP-EHH = 3.32) and PPM1K(BTA6, XP-EHH = 3.29) enriched the dephosphorylation BP (GO:0016311, p < 0.05). Among the top 10 positive XP-EHH markers in GJ, 8 are located on BTA8, in which annotated TYRP1. TYRP1 (BTA8, XP-EHH = 4.16) was contained in the cellular pigmentation BP (GO:0033059, p < 0.05). GRM4 (BTA23, XP-EHH = 3.89) was contained in the G protein-coupled glutamate receptor signaling pathway BP (GO:0007216, p < 0.05).The strongest iHS signal (5.67) in GJ was found on BTA1 (BovineHD0100035549), in which annotated DIPK2A. DIPK2A (BTA1, iHS = 5.67) was contained in the cardiac muscle tissue development BP (GO:0048738, p < 0.05). BPI (BTA13, iHS = 5.03) and LBP (BTA13, iHS = 5.03) enriched the regulation of macrophage activation (GO:0043030, p < 0.01).

Table 3 List of the top 10 |iHS| and positive XP-EHH for the SNP markers in GH100 and GJ, and their closest genes, with information on the Bos taurus chromosome (BTA) position

A total of 29 genes on BTA3, BTA5, BTA11, BTA22 and BTA23 were determined by more than one approach in GH100 in comparison to PH (Table 4). SCN5A, SCN10A and SCN11A on BTA22 enriched to the neuronal action potential BP (GO:0019228, p < 0.01). SCN5A and SCN10A enriched to the AV node cell action potential BP (GO:0086016, p < 0.01). SCN5A and SCN11A enriched to the cardiac muscle cell action potential involved in contraction BP (GO:0086002, p < 0.05). RAP1GAP (BTA2, XP-EHH = 3.32), SNX13 (BTA4, XP-EHH = 3.68), STAC (BTA22) and EPM2AIP1 (BTA22) enriched to the positive regulation of molecular function BP (GO:0044093, p < 0.05). TTLL2 (BTA9, iHS = 3.96), PRMT8 (BTA5, iHS = 4.45), HIPK2 (BTA4, XP-EHH = 3.42) and STK38L (BTA5) enriched to the peptidyl-amino acid modification BP (GO:0018193, p < 0.05).

Table 4 Genes identified by more than one test in plateau dairy cattle

A total of 20 genes on BTA1, BTA4, BTA6, BTA10, BTA11, BTA16 and BTA29 were determined by more than one approach in GJ in comparison to PJ (Table 4). LBP (BTA13, iHS = 5.03) and THBS4 (BTA10) enriched to the positive regulation of neutrophil chemotaxis BP (GO:0090023, p < 0.05). Two genes WDR48 (BTA22), MAEL (BTA3, iHS = 3.91) from Holstein and three genes ASPM (BTA16), APOB (BTA11, iHS = 4.61), TDRD15 (BTA11, iHS = 4.61) from Jersey were enriched for spermatogenesis BP (GO:0007283, p < 0.05). MAEL (BTA3, iHS = 3.91) selected in GH100, NCAPD2 (BTA5) and ASPM (BTA16) selected in GJ were enriched for meiotic cell cycle process BP (GO:1903046, p < 0.01).

Selection signatures based on crossbred dairy cattle in Tibet

The genome-wide distribution of |iHS| values in the GH40 is shown in Fig. 5(a). A total of 88 significant SNPs (p < 0.001) distributed across 20 chromosomes were identified. BTA14 was the chromosome harboring the largest number of these SNPs (33). Results of the pairwise genome-wide XP-EHH analyses between the GH40 and PH/XZ are shown in the Manhattan plots included in Fig. 5(b, d). Against PH, the XP-EHH selection signals were detected across all chromosomes except BTA19 and BTA28. Positive XP-EHH values (indicating a selection on the GH40 genome) were identified on BTA14, BTA17 and BTA24. Against XZ, the XP-EHH selection signals were detected on 23 chromosomes, and positive XP-EHH values were identified on BTA2, BTA8, BTA9, BTA14, BTA16, BTA20 and BTA30.F_ST values in the two between-breeds comparisons, that is GH40 v.s. PH and GH40 v.s. XZ, were 0.0451 and 0.0651, respectively. Figure 5(c, e) reports the Manhattan plots of the window-based pairwise genome-wide F_ST analyses of the GH40 breed against other breeds. In the GH40 v.s. PH comparison, F_ST signals were identified on 23 chromosomes. The SNPs with the highest F_ST values were on BTA25 (0.50), BTA3 (0.49) and BTA4 (0.46). In the GH40 v.s. XZ comparison, F_ST signals were identified across all chromosomes except BTA21, BTA25 and BTA27. The SNPs with the highest F_ST values were on BTA4 (0.59), BTA11 (0.54) and BTA30 (0.53).

Fig. 5

a iHS plot for the GH40. b XP-EHH plot for GH40 comparison with PH. c F_ST plot for GH40 comparison with PH. d XP-EHH plot for GH40 comparison with XZ. e F_ST plot for GH40 comparison with XZ. GH40 is crossbred dairy cattle in Tibet; PH is Holstein cattle in Beijing; XZ is Tibet cattle

The top 10 SNPs with their iHS values or positive XP-EHH values and annotated genes are shown in Table 5. The enriched biological processes were shown in Additional file 3: Table S10 – S11. The highest iHS value on autosomes was on BTA20, annotated to NKX2-5. NKX2-5 (BTA20, iHS = 4.47) and SMAD7 (BTA24, XP-EHH_PH = 3.32) enriched the ventricular septum morphogenesis BP (GO:0060412, p < 0.05). Among the top 10 positive XP-EHH values in the GH40 and XZ comparisons, 5 are located on BTA2. PINK1 (BTA2, XP-EHH_XZ = 4.00), HP1BP3 (BTA2, XP-EHH_XZ = 4.00) and MTOR (BTA16, XP-EHH_XZ = 3.60) were all shown to be selected in GH40 compared to XZ, enriched for cellular response to hypoxia BP (GO:0071456, p = 0.01). TERT (BTA20, XP-EHH_XZ = 3.62) enriched the RNA-templated DNA biosynthetic process BP (GO:0006278, p < 0.05).

Table 5 List of the top 10 |iHS| and positive XP-EHH for the SNP markers in GH40 and their closest genes, with information on the Bos taurus chromosome (BTA) position

A total of 52 genes and 35 genes were determined by more than one approach in GH40 in comparison to PH and XZ, respectively (Table 6). In comparison with PH, RAB27B (BTA24, XP-EHH_PH = 3.35) and OCA2 (BTA2) enriched to the pigmentation BP (GO:0043473, p < 0.01). TMEM109 (BTA29) enriched to the cellular response to radiation BP (GO:0071478, p < 0.05). TOE1(BTA3) and CPSF7 (BTA29) enriched to the RNA 3'-end processing BP (GO:0031123, p < 0.05). P4HA2 (BTA7) was contained in the peptidyl-proline hydroxylation BP (GO:0019511, p < 0.05). BTBD3 (BTA13), LAMA3 (BTA24) and ZFYVE27 (BTA26) enriched to the cell morphogenesis involved in neuron differentiation BP (GO:0048667, p < 0.05). NRCAM (BTA4) and BTBD3 (BTA13) are genes that were selected in PH and enriched for brain development BP (GO:0007420, p < 0.05). MOS (BTA14, iHS = 3.86), LYN (BTA14, iHS = 3.86) and FGF10 (BTA20) enriched for positive regulation of MAPK cascade BP (GO:0043410, p < 0.05). SMAD7 (BTA24, XP-EHH_PH = 3.32), FGF10 (BTA20) and ZFYVE27 (BTA26) enriched for response to growth factor (GO:0070848, p < 0.05). FGF10 (BTA20) and SFRP5(BTA26) enriched for regulation of canonical Wnt signaling pathway (GO:0060828, p < 0.05).

Table 6 Genes identified by more than one test in crossbred dairy cattle in Tibet

In comparison with XZ, CREBRF (BTA20) was detected to be selected in GH40 by 3 methods, contained in response to unfolded protein BP (GO:0006986). NCAPD3 (BTA15, iHS = 4.03) and NCAPG (BTA6) enriched to chromosome condensation BP (GO:0030261, p < 0.01). IRX4 (BTA20, XP-EHH_XZ = 3.62), ASTN2 (BTA8, XP-EHH_XZ = 3.59), TENM4 (BTA29), ZFYVE27 (BTA26) and CRTAC1 (BTA26) enriched to generation of neurons (GO:0048699, p < 0.01).

XKR4 (BTA14, iHS = 4.31) was detected in selection signals in both comparison groups and was annotated by 27 SNPs under selection in the iHS analysis, contained in plasma membrane phospholipid scrambling (GO:0017121, p < 0.05). ASTN2 (BTA8) was detected in XP-EHH in both comparison groups, differing in that selection occurred in GH40 compared to XZ (GH40v.s.xz) and in PH compared to PH (GH40v.s.PH). ZFYVE27 and SFRP5 on BTA26 were detected in the selection signals of both comparison groups. The generation of neurons (GO:0048699) containing ASTN2 and ZFYVE27 was significantly enriched in both comparison groups.

Discussion

Tibet has not historically had a specialized milk-producing dairy herd. Therefore, the dairy cattle in the present intensive farms are partly introduced from the plains and partly the offspring of grading up of native cattle. Based on genetic diversity analysis and LD decay, dairy cattle with different Holstein bloodline purity in Tibet exhibited low inbreeding levels, high levels of population polymorphism, while the LD levels were similar to those of native cattle breeds in Tibet. After the Holstein cattle were grouped by bloodline, significant differences were observed in the LD levels in each subpopulation. Studies have shown that artificial selection promotes the increase in LD levels within the population [35, 36]. GJ showed similar LD decay as PH and PJ that have undergone a long period of artificial selection and mating. GH40 and XZ showed a faster LD decay relative to the other populations.

Through the MDS we found that GH40 and XZ were more genetically connected compared to AP and TF. Therefore, we chose XZ as a representative of local breeds to participate in the analysis of selection signatures. Combining MDS, the NJ phylogenetic tree, and Admixture population structure analysis, GH85, GH75, and GH40 exhibited clear signs of historical hybridization, whereas GH100 and PH were not completely separate. This suggested lower genomic consistency among different-purity Holstein cattle subpopulations in Tibet. Consequently, in subsequent research, these populations should not be treated as a single group. In contrast, GJ, which has not undergone hybridization, exhibited complete genomic admixture with PJ.

It is through selection and mating systems that people in Tibet are trying to solve the problem of insufficient milk production in dairy cattle living on the plateau. This includes the use of semen from dairy cattle to mate with native cattle to improve milk production performance. At the same time, some purebred dairy cattle that could not adapt to the plateau environment were culled due to disease or death. According to the MDS, the first dimension component separated the native Tibetan breeds from the other populations, indicating that the high-altitude environment led to distinctive genomic characteristics in native Tibetan cattle. In contrast to commercial high-bred dairy cattle, Tibetan dairy cattle retain greater genomic variability, which may be crucial for addressing current and future challenges related to environmental adaptability and meeting breeding goals.

In Admixture GH100 appears in the same areas as native cattle in the Plateau. Due to the lack of pedigree, we were unable to determine the mating history of GH100 and the origin of these regions. The size of the overlap region between the first Holstein cattle introduced to the plateau and the plateau native cattle is unknown. We were also unable to determine what percentage of these regions were recognized as Holstein bloodlines in the calculations due to the absence of BBR values for plateau native cattle.

When characterizing the genetic background of strategies for surviving at high altitude, one has to differentiate between the evolutionary genetic adaptations of native highland species and populations, e.g. yaks and blue sheep, and the acclimatization of lowland species brought to high altitude short term, for which phenotypic plasticity is the decisive factor [37, 38]. Strategies and mechanisms of both acclimatization and adaptation have a genetic background and are thus linked to each other. In this study, by comparing dairy cattle living at different altitudes, we are able to access the genetic differences caused by short-term environmental changes over several generations. By comparing crossbred cattle with Holstein and native cattle, it is possible to identify genes that differed between breeds in crossbred cattle after crossbreeding. The F_ST method is best suited for detection in events occurring in the more distant past [39] whereas the iHS test is best suited for detection in recent positive selection [40]. The XP-EHH is useful for detection of entirely or approximately fixed loci [41]. As the F_ST, iHS, and XP-EHH approaches are complementary, incorporating such measures in the study of genetic differentiation and selection would strongly contribute to understanding gene flow and genetic makeup of plateau dairy cattle (purebred and hybrid) populations.

Our study did not annotate the same genes near selection signals in Holstein and Jersey cattle, but the functions associated with these genes are similar. Various within- and between-species analyses of altitude adaptation have also shown that in most cases different genetic architectures and physiological mechanisms underlie these adaptations [38]. Chronic hypoxia with severely reduced oxygen supply is a major stressor in high-altitude conditions and limits the respiratory and cardiovascular systems from functioning properly [42]. Glover and Newsom first reported in 1915 that native cattle in Colorado, USA, at elevations higher than 1524 m, suffered from altitude sickness, which was characterized by right heart failure and chest edema due to pulmonary hypertension [43]. In this study, we have identified several heart-related genes in GH100. Genetic variants of SCN5A are involved in a number of inherited cardiac channelopathies and are also found to be correlated with myocardial contractile dysfunction, dilated cardiomyopathy, and heart failure [44]. Variants in SCN10A are associated with alterations of cardiac conduction parameters and the cardiac rhythm disorder Brugada syndrome [45]. Rap1GAP plays a critical role in mediating Ang II-induced cardiomyocyte hypertrophy through its regulation on autophagy and oxidative stress [46]. Rap1GAP was verified to regulate the AMPK SIRT1/NF-κB signaling pathway and exacerbate the damage to myocardial cells caused by ischemia and hypoxia [47]. SNX13 deficiency promotes cardiomyocyte apoptosis and heart failure [48]. Kun Lian found that the cardiac function of diabetes mice overexpressing PPM1k was improved, and the myocardial infarction area and apoptosis were reduced [49]. Exercise-induced HIPK2 suppression attenuates cardiomyocytes apoptosis [50]. The serine/threonine kinase HIPK2 as novel regulatory protein participating in hypoxic gene regulation can affect apical as well as downstream events during the hypoxic response [51]. THBS4 and DIPK2A were identified in the selection signatures of GJ in this study. Studies in THBS4-knockout mice revealed that deletion of THBS4 resulted in a significantly higher heart weight to body weight ratio, severe excessive cardiac fibrosis, and aggravated cardiac dysfunction than those in the WT mice in response to pressure overload [52, 53]. DIPK2A protects cardiomyocytes from ischemia-induced cell death by activating PRKCE (protein kinase C epsilon) [54]. Many genes associated with cardiac development are subject to selection may differ in the ability of GH40 to cope with hypoxia compared to PH and XZ. NKX2-5 is involved in cardiac morphogenesis, including rightward looping and subsequent chamber specification and septation and in functional maturation and maintenance of working myocardium and conduction system [55]. SMAD7 is required for cardiac development in embryonic stage and cardiac function in adult [56]. SFRP5 diminishes cardiac inflammation and protects the heart from ischemia/reperfusion injury [57]. Serum SFRP5 levels were significantly associated with coronary artery disease in humans [58]. HP1BP3 was up-regulated to initiate a cardioprotective stress-handling response [59]. IRX4 may play a critical role in establishing chamber-specific gene expression in the developing heart [60].

In a previous study, angiogenesis is one of the key mechanisms of high-altitude acclimatization [61]. Hypoxia-stimulated vasculogenesis and angiogenesis may increase microcirculatory blood flow and capillary density, thus restoring the oxygen supply in tissues exposed to hypoxia to compensate for oxygen insufficiency [62]. ALPL and APOB were identified in the selection signatures of Gh100 and GJ in this study, respectively. ALPL might be as a distinct regulator to control exosomes secretion and subsequently regulate the pro-angiogenic potential of mesenchymal stem cells [63]. Apolipoprotein B levels can be a predictor of ischemic cardiovascular disease, and mutations and polymorphisms in APOB are associated with plasma apolipoprotein B levels, and risk of ischemic cardiovascular disease [64]. APOB and TDRD15 found to be selected in Welsh sheep breeds living on upland pasture [65]. Among genes selected in GH40, LYN has well-established functions in most haematopoietic cells, viz. progenitors via influencing c-kit signaling, through to mature cell receptor/integrin signaling, e.g. erythrocytes, platelets, mast cells and macrophages [66]. ZFYVE27 was recently found to promote endothelial cell migration and angiogenesis in humans, with its knockdown affecting these processes [67]. P53 activates P4HA2 expression exerting anti-angiogenic effects [68]. Activation of mTOR potently enhances the activity of HIF1α and vascular endothelial growth factor (VEGF)-A secretion during hypoxia [69]. TERT mutations are associated with hematopoietic disorders [70]. CREBRF was significantly upregulated by reconstituted high-density lipoproteins (rHDL) in response to both hypoxia and inflammation, and is novel target of rHDL in response to angiogenic stimuli [71].

The brain shows to be more sensitive to the fluctuation of oxygen content, brain damage caused by hypoxia impaired cognitive physiology and psychological function, and lead to learning and memory defects [72]. Among the candidate genes for GH100 and GJ, PRMT8 is an important regulatory component of membrane phospholipid composition, short-term memory function, mitochondrial function, and neuroinflammation in response to hypoxic stress [73]. GRM4 is crucial in ion channel regulation, neuronal excitation and neurotransmitter release [74]. Great concern has been attached to the mechanism of GRM4 in cancer biology in recent decades [75]. ASPM expression is critical for proper neurogenesis and neuronal migration, and ASPM is a positive regulator of Wnt signaling [76]. A large number of neuron-related genes are subject to selection among different breeds, possibly related to the different environments and life patterns in which different breeds are living. TOE1 and BTBD3 plays a crucial role in neurodevelopment [77, 78]. NRCAM plays a wide variety of roles in neural development, including cell proliferation and differentiation, axon growth and guidance, synapse formation, and the formation of the myelinated nerve structure [79]. There is extensive evidence to confirm that PINK1 functions to protect neurons against stress‐induced cell death, and PINK1 deficiency has been shown to have a profound effect on respiration and reactive oxygen species (ROS) production [80]. Aberrant expression of TENM4 may leads to lower learning ability, sleep reduction, and increased aggressiveness in animal models [81]. ASTN2 has pleiotropic effects on cardiometabolic and psychiatric traits, and associations between genetic variants in the ASTN2 locus and blood pressure, total and central obesity, neuroticism, anhedonia and mood instability were identified [82].

High-altitude environmental factors (such as UV exposure, cold and hypobaric hypoxia) can affect the immune system and increase susceptibility to cancer, various infections, and even autoimmune disease [83]. BPI and LBP were identified in the selection signatures of GJ. In mammals, bactericidal permeability/increasing protein (BPI) and lipopolysaccharide (LPS) binding protein (LBP) are two important pattern-recognition receptors involved in the innate immune response to LPS [84]. Mammalian BPI has been shown to play bactericidal, opsonic and anti-inflammatory roles [85]. The strong solar ultraviolet radiation (UVR) increased the risk of skin cancer in plateau animals. Melanin in the skin as the main component shielding UVR, prevents skin and organ damage by absorbing UV photons and free radicals to protect DNA from UV damage [86, 87]. TYRP1 is related to melanin synthesis, the biology of melanocytes and melanosomes, and the migration and survival of melanocytes during development [88]. The expression of RAB27A and RAB27B in both melanocytes and platelets makes them candidates for involvement in mouse and human disorders characterized by the combination of pigment dilution and a platelet storage pool defect [89]. OCA2 encodes the protein P, a transmembrane protein, and has been shown to play a role in pigmentation in both humans and mice [90]. In humans, it has been implicated in iris, skin, and hair pigmentation [91]. TYRP1 was identified in the selection signatures of GJ. RAB27B and OCA2 were identified in the comparisons of different breeds.

In addition, the selective signatures between breeds were annotated to a number of other growth and development-related genes. FGF10 serves as an essential regulator of lung and limb formation [92]. The LAMA3 gene belongs to the laminin family and is responsive to several epithelial-mesenchymal regulators, including keratinocyte growth factor, epidermal growth factor and skin fibrosis [93]. SNP loci within the CRTAC1 had a positive impact on milk fat, lactose, and total solids (TS) content in Gannan yak milk [94]. Variation in the XKR4 gene is significantly associated with subcutaneous rump fat thickness in indicine and composite cattle [95], and is associated with feed intake and growth phenotypes in cattle [96]. PITX1 as a key specificity factor in HIF-1α dependent responses is required for the survival and proliferation responses in hypoxia [97]. PITX1 is required for pituitary gland and hind limb development [98].

In this study, a series of genes related to cardiac health and development, neuronal development and function, angiogenesis and hematopoietic were found through genome analysis of dairy cattle at Tibet, and these genes showed significant selection pressure in the plateau environment. These findings provide potential molecular markers for the breeding and genetic improvement of dairy cows, which can help select superior dairy cattle breeds that are better adapted to the plateau environment and improve their survival, fertility and milk yield in low oxygen environment. There is still a need to find or validate variant loci in candidate genes associated with plateau acclimatization and adaptation in combination with phenotypic data, as well as other sequencing and analytical methods, and with reference to results in other species. More genotyping data of Jersey cattle can help to identify whether some of the variation is breed-specific, and provide data to support the exploration of acclimatization and adaptation mechanisms in the two breeds.

Farmers can gradually improve the plateau adaptability and production performance of dairy herds by selecting cows with favorable genes for breeding, regardless of whether the variation is acquired through crossbreeding or natural selection. Traditional breeding methods and modern genome selection techniques can be combined to more effectively select individuals with anti-hypoxia and enhanced immunity through marker-assisted selection (MAS) or whole-genome selection (GS). Scientists need to continuously monitor the performance of these genes in different environments, optimize breeding programs, and develop nutrition and management strategies suitable for the plateau environment to help farmers improve the growing environment of dairy cows. At the same time, these results also laid a foundation for other livestock adaptation research in the plateau area, with a wide range of application prospects.

Conclusion

The composition of the Tibet dairy cattle population is complex. There are significant genetic background differences between Tibet dairy cattle and Tibet native cattle populations. A number of genes associated with cardiac health and development (SCN5A, SCN10A, RAP1GAP, SNX13, PPM1k, HIPK2, THBS4, DIPK2A, NKX2-5, SMAD7, SFRP5, HP1BP3, IRX4), neuronal development and function (PRMT8, GRM4, ASPM, TOE1, BTBD3, NRCAM, PINK1, TENM4, ASTN2), angiogenesis and hematopoietic (ALPL, APOB, LYN, ZFYVE27, P4HA2, MTOR, TERT, CREBRF),and pigmentation (TYRP1, RAB27B, OCA2) are selected for in plateau dairy cattle. In addition, genes related to immune response (BPI, LBP) were selected for plateau Jersey cattle, and genes related to growth and development (FGF10, LAMA3, CRTAC1, XKR4, PITX1) were selected for in crossbred cattle. These candidate genes provide valuable information for further fine-mapping studies and functional genomic studies as a starting point to identify causal genetic variants that control plateau adaptation traits for the implementation of breeding schemes aiming genetic improvement of plateau dairy cattle populations.