Whole genome resequencing-based analysis of plateau adaptation in Meiren yak (Bos grunniens)

Abstract

The Meiren yak is an important genetic resource in Gansu Province, China. In this study, we aimed to explore the evolutionary history and population structure of the genetic resource of Meiren yak and to mine the characteristic genes of Meiren yak. We analysed a total of 93 yaks of eight yak breeds based on whole genome resequencing combined with population genomics and used θπ ratio and Fst method to screen the selected sites in the genome region. The results proved that Meiren yak can be used as a potential genetic resource in Gansu Province. The genes in Meiren yak with positive selection in selection signal analysis were subjected to the Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses, which indicated that the genes were related to the adaptability to high altitude and hypoxic environment. By analysing the genetic variation of Meiren yak at the genome-wide level, this study provided a theoretical basis for genetic improvement of Meiren yak and for the development of high-quality yak resources.

Keywords:

• Meiren yak
• whole genome resequencing
• population genetic structure
• genetic resources
• selection signals

Introduction

Yak (Bos grunniens) belonging to class Mammalia and family Bovidae, is a large mammal living in the plateau.¹ As they live in harsh environment for a long time, yaks have strong stress-resistance and high disease-resistance.² Yak is called as ‘the boat of plateau’ and ‘versatile livestock’. It is a very important gene bank of animals living on the plateau.³ The population of yaks in China is more than 16 million,⁴ accounting for more than 95% of all yaks globally.⁵ Moreover, China has the most number of yak breeds,⁶ with different regions having own representative breeds. Currently, wild yak and domestic yak are considered two subspecies of yak.^7,8 The distribution of yaks in China is mainly concentrated on the Qinghai–Tibet Plateau and extends around Pamirs to Tibet, Qinghai and Gansu.⁹ Yak is the dominant local livestock breed and source of meat and fur in the plateau area; therefore, it has high applied economic value.¹⁰ However, currently, population of yaks is decreasing.¹¹ Therefore, this should be urgently prevented by improving yak breeds and cultivating new breeds. In the long-term evolution process, yaks have developed relatively stable genetic characteristics and physiological and biochemical mechanisms.¹²

Meiren yak is an excellent local type of yak in Gannan, which is mainly distributed in the cooperative city of Gannan Tibetan Autonomous Prefecture at an altitude of approximately 3600 m with a typical plateau continental climate. The average slaughter rate of Meiren yak is 50.8%, and the average dressing percentage is 39.3%. Its meat is tender, low in fat, and high in protein. A local saying states that ‘Meiren yak is a wall’, indicating that Meiren yak has a strong physique and well-developed muscles. Indeed, Meiren yak has obvious advantages among other types of yaks in Gannan. Meiren yaks have lived long-term nomadic life; inbreeding of Meiren yak leads to serious degradation in the variety, and no improvement in meat or fur production is observed. This is causing slow growth and development, low reproduction rate, reduced production performance and weight loss in these yaks year by year.¹³ The breeding history of Meiren yak consists of a cycle of domestication–wild yak gene infiltration–domestication. Under the long-term natural selection and artificial breeding, it has developed strong disease resistance, adaptability, and unique physiological mechanisms. As an excellent subtype of Gannan yak, Meiren yak is an extremely important livestock genetic resource in China and there has a lot of development value.

Whole genome resequencing can be used to sequence the genomes of multiple samples or animal groups and compare with the published reference genomes for analysing the differences. With this, a large number of genetic variations can be identified, which has great application prospects in current research on breeding.^14,15 Whole genome resequencing is a category of population genomics and population genetics.¹⁶ Through whole genome sequencing, we can infer not only the population history of species but also the differences within and between groups of the same species..¹⁷ With the gradual development of sequencing technology and data processing, we can calculate and analyse the degree of segregation¹⁸ and speculate the population evolution (such as linkage disequilibrium (LD) analysis, inbreeding coefficient, identification of genetic relationship, domestication and migration and verification of functional genes) using population genomics.¹⁹ Liu et al.²⁰ sequenced and analysed certain locations in the genomes of 32 tigers and revealed that they had retained 1.8 million high-quality single-nucleotide polymorphisms (SNPs). This indicated that tigers had experienced at least one population bottleneck effect, and at the same time, the existing tiger population was divided into six subspecies. This provided a scientific basis for the protection of tigers. Kim et al.²¹ studied 48 native African cattle, compared their genomes with those of 4 strains of cattle from different sources, and revealed the mechanism of strong adaptability of African cattle to harsh environment. Wu et al.²² studied 98 samples of 5 species of bovine genus and reported that gene infiltration existed in different individuals and this infiltration was an important basis for adaptive variation. Guo et al.²³ sequenced the whole genome of six goat breeds, annotated the genomes, and identified several genes related to their important economic traits, such as high altitude adaptability (EPAS2) and goat coat colour (IRF2, EXOC2, etc.), enhancing the understanding of adaptation mechanisms in goats. Chen et al.²⁴ studied 49 modern breeds and eight East Asian breeds of cattle and reported for the first time that the domesticated cattle in the world can be divided into 5 continental groups;, this provided a new insight for their origin and historical migration. Fan et al.²⁵ sequenced the genomes of more than 500 sika deer, red deer, and hybrid deer and selected 1,000 specific SNPs loci from more than 13 million screened loci to synthesize a chip. This chip can accurately distinguish sika deer from hybrid deer, thus providing a basis for the protection of related germplasm resources. Li et al.²⁶ studied 70 goats from Inner Mongolia and Liaoning. They identified a large number of genetic variations, enhanced understanding of genetic diversity of different phenotypes, and provided molecular basis for livestock breeding.

As a valuable resource bank of yak genes, the genetic differentiation and classification of Chinese yak has always been a difficult problem in genetics.²⁷ Whole genome sequencing technology has been fully applied in species research, and provided theoretical data for the protection and breeding of species.²⁸ At present, there is no analysis on the genetic structure of Meiren Yak population within the whole genome resequencing.In our experiment, 93 samples from eight breeds, including Meiren Yak、 Changtai Yak、 Muli Yak、 Pamier Yak、 Sibu Yak、 Tianzhu Yak and Zhongdian Yak, were sequenced, and the population genomics was compared, so as to provide theoretical basis for genetic improvement of Meiren Yak and a good research foundation for the development and utilization of high-quality resources.

Materials and methods

Ethics statement

All animal studies conducted in this study were approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (SYXK-2014–0002).

Sample collection, blood DNA extraction and whole genome resequencing

We performed genome-wide resequencing using 93 blood samples from eight yak breeds (Meiren、 Changtai、 Jiulong、 Muli、 Pamier、 Sibu、 Tianzhu and Zhongdian yaks; n = 19、 6、 16、 15、 10、 11、 8 and 8, respectively). In total, 10 mL blood was collected from the jugular vein, placed in EDTA anticoagulation tube, and stored at −20 °C till further analysis. DNA was extracted using EasyPure Blood Genome DNA Kit (Quanjin, Beijing) and specifically adsorbed using enzymatic hydrolysis and silica gel membrane centrifugal column. The DNA integrity, concentration, and quality were determined using agarose gel electrophoresis and ultraviolet spectrophotometer NanoDrope2000. Genomic DNA was sequenced on illumina PE150 platform. Instructions for use of animals: All yaks used in this study were personally owned by local herdsmen and obtained the informed consent of the animal owner.

Whole genome sequence alignment and detection of genetic variation

Certain errors in the data generated after DNA extraction, library construction, and sequencing (e.g., sequencing errors and no detection of bases) will have corresponding effects on subsequent analysis. Therefore, to ensure the accuracy of bioinformatic analysis, we first evaluated the quality of the sequenced data and screened high-quality reads. This quality control was performed using Trimmomatic,²⁹ and high-quality reads were compared with the latest reference yak genome (GCA_005887515.1). SNP identification was performed using GATK software, and the SelectVariant and VariantFiltration parameters included in GATK software were used to set QUAL、 QD、 FS、 MQ、 HaplotypeScore and other conditions for SNP filtering.³⁰ The conditions were as follows: QUAL < 30.0、 QD < 2.0、 FS > 60.0、 MQ < 40.0 and HaplotypeScore > 13.0. To prevent the influence of sex and low allele frequency on subsequent analysis, PLINK³¹ was used to remove the following: minimum allele frequency < 0.05, the individual detection rate < 0.95, the SNP detection rate < 0.99, no chromosomal assignment, and SNP on the sex chromosome. The conditions were as follows: –mAF 0.05 –mind 0.05 –geno 0.01 –chr 1-29.

Population genomic analysis

Construction of phylogenetic tree, principal components, group structure and kinship analysis

The most important and common aspect of phylogenetic analysis is the construction of phylogenetic trees. Phylogenetic trees are also called clustering trees or evolutionary trees. Phylogenetic trees can reflect the evolutionary relationship and evolution of species, kinship of samples within the group and the classification and evolutionary relationship of gene family members. We used PLINK³¹ software to calculate the genetic distance matrix between individuals based on the genotypic information of SNP data and further used the adjacency method to construct a phylogenetic tree.

Principal component analysis (PCA) is a cluster analysis method. It is a commonly used analytical method in population genetics and widely used in many disciplines. In genetics, PCA can classify individuals according to principal components as per the differences in SNPs on the individual genome; it uses fewer variables to explain most of the variation in the original data, transforms many highly correlated variables into independent or uncorrelated variables, and helps in understanding the clustering and differences of the study population. We used PLINK³¹ for PCA of autosomal SNP data. bcftools was used to extract a subset of VCF files. First, the VCF file was converted to plink format, resulting in three files ending in ‘.map’, ‘.nosex’ and ‘.ped’. A binary bed file was generated based on the ‘.ped’ file, and two files ending in ‘.bim’ and ‘.bed’ were obtained. Further, PCA analysis was performed to obtain two files ending in ‘.eigenval’ and ‘.eigenvec’. The ‘.eigenval’ file represents the proportion of each PCA, and the ‘.eigenvec’ file records the feature vectors used for plotting the coordinates. Finally, the PCA plot was visualized using the ggplot2 package in R language.

Population structure analysis was first produced by the “STRUCTURE” software; therefore, it is also known as STRUCTURE analysis. It is an important genetic structure analysis tool, which can divide a population into several subgroups according to molecular genetics, infer the optimal number of groups of a population, and determines any gene communication among different subpopulations and the degree of hybridization of individuals within the subgroup. In this study, ADMIXTURE³² software was used to estimate the mixing ratio of each component in the population according to the expectation maximization (EM) algorithm, and finally, the ggplot2 and pophelper packages in R were used to draw the population structure stacking map.

To identify individual samples with close genetic relationship and remove the influence of these samples on subsequent analysis, in this study, we used the KING-kinship³³ software to calculate the kinship coefficient between two individuals for all SNP sites. Finally, the R-packet corrplot was used to visualize the affinity matrix.

Heterozygosity, linkage disequilibrium and continuous homozygous region analysis

Heterozygosity, also known as population mean heterozygosity, refers to the presence of one or more different alleles at a certain site. As the basis of the adaptation and evolution of a population or species, genetic diversity is an important reference index for the assessment of the status of germplasm resources. We used the mlRho³⁴ software to calculate the heterozygosity of individuals as one of the indicators to assess genetic diversity.

LD refers to the chance that alleles from two or more gene loci appear on a chromosome at the same time with frequency higher than the random frequency. Currently, it is under focus in studies related to population genetics and plays an important role in the study of population evolution. LD can be influenced by factors such as genetic drift, changes in mutation rate, and gene conversion. Therefore, studying the decay in LD can reveal the history of genetic recombination in the population. In this experiment, the PopLDdecay³⁵ software was used to calculate the degree of LD in the population, and the decay of LD was plotted using the software’s own Plot_MultiPop.pl script.

Runs of homozygosity (ROH) are continuous homozygosity fragments and refer to a class of continuous homozygosity phenomena that occur in a class of genomes, manifested by the lack of heterozygosity in a chromosomal region. Many factors affect the production of ROH, such as natural selection, gene inflection and inbreeding and gradually make it a group characteristic of different populations. At present, inbreeding is considered the most important influencing factor for ROH. In the process of biological evolution, ROH will exhibit various distribution patterns on the animal genome; therefore, the number, length, and distribution of ROH on the animal genome can reveal the rich genetic background such as the evolutionary history and inbreeding degree of animal populations. In this study, we used PLINK to detect ROH, and after calculation, files such as ‘.hom’ and ‘.hom.indiv’ were obtained. The ‘.hom’ file contains the specific information of each ROH, and the ‘.hom.indiv’ file contains a summary of ROH information for each individual, including the total length and total number. The resulting NSEG is NROH and KB is SROH. Finally, the SROH-NROH diagram is prepared according to the above two files.

Population differentiation index and genetic polymorphism analysis

The population differentiation index is a statistical test used in population genetics to measure the degree of differentiation among populations. The larger the differentiation index, the greater the difference. This index is suitable for the comparison of diversity between subgroups. The index value ranges from 0 to 1; 0 indicates that the two populations are randomly mated, and the genotypes are completely similar; 1 indicates that the populations are completely isolated and completely dissimilar. In this study, VCFtools³⁶ software was used for the sliding window analysis, and the results were visualized. We calculated the nucleotide diversity π (window-size 50 KB, window-step 25 KB) of all sample populations and other 8 populations using Vcftools software.

Selection signal analysis and functional enrichment analysis

Selection signal analysis is a common analytical method in the field of population evolution research. It is mostly used in the analysis of artificial domestication, natural selection, and other mechanisms, which are closely related to the breeding direction of species and domestication adaptation mechanism.³⁷ Calculate log2θπ ratio and Fst values using VCFtools.³⁶ (Window size 50 KB, window step 25 KB). First, the nucleotide polymorphism analysis was conducted using sliding window method. The nucleotide polymorphism ratio of Meiren yak and other yaks was calculated, and the area with significant outlier value (p < 0.01) was considered as the first condition for selecting the region. Further, the group differentiation index of Meiren yak and other yak groups was calculated using sliding window method, and the area with significant outlier value (p < 0.01) was considered as the second condition. The θπ ratio and Fst method were used for selection signal analysis to screen genomic regions.³⁸ This is a very effective method for detecting selective elimination regions, often resulting in strong selection signals.³⁹ The genome-wide variation of Meiren and Pamir yaks was analysed and the results were used for functional enrichment analysis of identified genes. The DAVID library was used to perform the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of candidate genes, and finally, the enrichment analysis map was constructed.⁴⁰

Results

Sample collection and sequencing results

In this study, 93 blood DNA samples from eight yak breeds were sequenced (the sequencing depth was 10×). Through filtering by mass value, depth, deletion rate, and minimum allele frequency, 7,666,824 high-quality SNP-labelled data were finally obtained. The data were relatively high, indicating that the sequencing data were of good quality, could reflect the genomic information of the corresponding individual, and could be used for further bioinformatic analysis. In addition, the SNP annotation results of eight yak breeds revealed that the distribution of SNP sites in each yak population was similar. They were mainly located in the intergenic region, followed by the intron region, upstream region of the gene, downstream region of the gene and exon region (Table 1).

Table 1. Annotation information of all SNPs in 8 yak populations.

Breeds		Changtai	Jiulong	Meiren	Muli	Pamier	Sibu	Tianzhu	Zhongdian
Samples		6	16	19	15	10	11	8	8
SNPs		10343725	13971741	10257729	10957636	10325426	13529944	10671202	10769073
Ts/Tv ratio		2.55	2.54	2.52	2.55	2.56	2.54	2.56	2.56
Hom/Het		0.52	0.55	0.6	0.57	0.62	0.53	0.54	0.54
SNP types
Exon	Synonymous variant	44888	61547	33982	47958	44738	59384	46066	46158
	Initiator codon variant	15	27	13	15	16	22	12	14
	Start lost	159	219	95	183	158	226	159	158
	Stop gained	745	1222	926	815	765	1252	776	795
	Stop lost	236	337	110	293	245	317	256	239
	Stop retained variant	62	80	50	65	64	75	61	63
Splice sites	Splice region variant	8305	11737	6241	9156	8288	11239	8501	8452
	Splice acceptor variant	276	486	170	338	299	497	301	286
	Splice donor variant	385	665	227	490	390	649	403	397
Intron	Intron variant	3298309	4399517	3149398	3483506	3275875	4206247	3406730	3415608
	Intragenic variant	68	103	37	90	74	94	79	78
UTR	5 prime UTR variant	8487	12597	4745	9923	8788	12362	8952	8772
	5 prime UTR premature start codon gain variant	1195	1722	961	1273	1180	1697	1262	1249
	3 prime UTR variant	20953	28415	16743	22526	20832	27821	21573	21594
Intergenic	Upstream gene variant	675057	936438	649511	716339	668501	910370	693423	695661
	Downstream gene variant	215081	297034	212150	228446	212998	291258	220750	222108
	Intergenic region	6017519	8141055	6141606	6375723	6028438	7929077	6207902	6294148

Population genomic analysis

Analysis of cluster differentiation, group structure and kinship

Phylogenetic trees can be used to describe the evolutionary relationships and evolutionary order among populations. To study the relationship and differentiation degree of the whole genome of various breeds of yaks, the phylogenetic tree was constructed. The results revealed that all yak breeds were clustered together individually, and no mixing was observed. Among them, Changtai yaks were clustered into a separate category; Kowloon, Zhongdian, and Muli yaks were grouped together; Sibu yak, Pamir yak and Meiren yak were grouped together, and these three yak groups were grouped together with Tianzhu yak (Figure 1(A), Schedule 1). Each locus in the PCA plot represents one sample, and the farther apart the two samples in the figure, the greater the difference in the genetic background of the two samples. The results indicated that a certain degree of differentiation existed among the eight yak breeds. Among them, the Muli, Kowloon and Pamir yaks exhibited a higher degree of differentiation, which was consistent with the characteristic phylogenetic tree map (Figure 1(B), Schedule 2). Further, the accurate and appropriate number of subgroups (k-value) was determined. The k-value was assumed to be between 1 and 10, and the cross-validation error for each breed belonging to the kth subgroup was calculated. In the population structure analysis, k = 2 indicated that all yaks were divided into two groups, and Pamir yaks were isolated first. k = 3 indicated that the Kowloon and the Muli yaks are separated. k ≥ 4 indicated that the analysis results of some yak populations are similar (Figure 1(C)). The results of kinship analysis indicated that individuals without three generations of kinship in the Meiren yak and Pamir yak groups met the requirements of subsequent analysis (Figure 1(D), Schedule 3). It suggests that the individual yaks used in this study can be used for subsequent selection signalling analyses in Meiren yaks and Pamir yaks.

Figure 1. Cluster differentiation and affinities of yak populations. (A) Phylogenetic tree. (B) Principal component analysis of yak populations. (C) Ancestral component analysis. (D) Analysis of genetic relationships in yak populations.

Heterozygosity, LD, continuous homozygous regions, genetic polymorphisms and population differentiation index analysis

Genetic diversity is fundamental to a population’s adaptation to its environment and evolution. The expected heterozygosity of the eight yak breeds was approximately 0.357; the observed heterozygosity of Meiren yaks was the lowest, and that of Changtai yaks was the highest. The observed heterozygosity value of eight yak groups was less than the theoretical value (Figure 2(A), Schedule 4). To evaluate the LD of yak groups, PopLDdecay was used to calculate the pairwise r2 values of each group. LD analysis indicated that Meiren and Changtai yaks had the fastest and slowest LD decay velocity, respectively (Figure 2(B)). Continuous homozygous regional analysis revealed that the quantity and length of ROH in various yak populations were highly correlated (Figure 2(D), Schedule 5). Further, the ROH length was classified into <1 Mb, 1–2 Mb and >2 Mb. Most ROH were of size <1 Mb, followed by 1–2 Mb. The number of ROH with size >2 Mb was 0, 3, 0, 6, 3, 1, 2 and 2, respectively, in Changtai, Jiulong, Meiren, Muli, Pamier, Sibu, Tianzhu and Zhongdian yaks (Figure 2(C), Schedule 6). The results of genetic polymorphisms indicated that the nucleotide diversity of Muli and Pamir yaks was low (0.0009), and that of the rest of the yak breeds was 0.001(Figure 2(F)). Therefore, to quantitatively analyse the genetic differentiation of yak populations, the genetic differentiation index of eight yak breeds was accurately calculated, which was > 0.04. It met the currently recognized genetic differentiation index standard in yak breeds, which considers that groups with a population differentiation index > 0.04 have a higher degree of population genetic differentiation (Figure 2(E)).

Figure 2. Population genomics analysis. (A) analysis of heterozygosity in yak populations. (B) Linkage disequilibrium analysis. (C,D) ROH analysis of yak populations. (E) Population differentiation index analysis. (F) Genetic polymorphism analysis.

Selection signal analysis and functional enrichment analysis

The θπ ratio and Fst method were used for selection signal analysis using defined thresholds (log2 θπ ratio > 1.901, Fst ≥ 0.272) to screen genomic regions. Values above the threshold were considered positive selections and can be used to divide groups with unique characteristics (Figure 3(A), Schedule 9). In GO enrichment analysis, these genes were mainly enriched in molecular functions (UDP-N-acetylglucosamine-lysosomal enzyme and N-acetylglucosamine phosphotransferase activity, etc.), followed by biological processes (positive regulation of amacrine cell differentiation, odontogenesis of dentin-containing tooth, regulation of amacrine cell differentiation, etc.)(Schedule 7). Simultaneously two genes (MBP2 and ACVR2B) were screened that were related to multiple physiological processes of animals adapting to hypoxia at high altitude. KEGG enrichment analysis revealed that the genes were associated with biological pathways such as folate resistance, pantothenic acid and CoA biosynthesis, α-linolenic acid metabolism, biosynthesis of unsaturated fatty acids, fatty acid degradation, fatty acid metabolism, and nitrogen metabolism (Figure 3(B), Schedule 8).

Figure 3. Enrichment analysis results. (A) Select signal analysis. Based on θπ ratio and Fst analysis to take the intersection for selection signal analysis. The horizontal coordinate is the ratio of θπ (MR/PM) and the vertical coordinate is the Fst value, which corresponds to the frequency distribution plots on the top and on the right, respectively, and the dot plots in the centre represent the corresponding Fst and θπ ratios within different windows. The top blue and green regions are the regions selected by θπ, the red regions are the regions selected by Fst, and the Middle blue and green regions are the intersection of Fst and θπ, i.e., the candidate regions. The blue dots are candidate areas for MR and the green dots are candidate areas for PM yaks. (B) KEGG enrichment analysis of genes in selected regions(the area of green dots in Figure 3(A)).

Discussion

Yak is a very important germplasm resource in the plateau. Compared with other cattle breeds, yak attains sexual maturity late in life and has low fertility.⁴¹ It has a strong adaptability to cold and hypoxic environment; therefore, it is all-year-round livestock providing meat for local herders.⁴² Among them, Meiren yak, as an extremely important genetic resource of livestock in China, has a strong physique and muscles. To better understand the evolution and genetic characteristics of Meiren yak, we performed whole gene resequencing, and the obtained data were compared with those of other yak populations. The genome characteristics of the yak populations and differences between Meiren yak and other yak breeds were studied using population genomic analysis.

The phylogenetic tree revealed that all yak groups were clearly distinguished. Changtai yak was alone clustered into one group and had the farthest distance. The results of PCA were consistent with those of phylogenetic tree, where Meiren and Tianzhu yaks were closer. This was consistent with the real-life presence of both yak breeds in Gansu Province. The optimal number of subgroups was not observed because of frequent hybridization between yak groups. Heterozygosity analysis revealed that all yak populations exhibited rich genetic diversity, and the observed heterozygosity values were less than the theoretical values. This indicated that various degrees of decline in heterozygosity were observed in yak populations. LD is affected by the inbreeding level of livestock. The higher the degree of LD, the higher the level of inbreeding.^43,44 Meiren yaks had the fastest LD decay velocity, indicating that the inbreeding level and degree of LD were low in Meiren yaks. This is consistent with the feeding mode and local breeding mode of Meiren yaks living in the Meiren grassland for a long time. Meiren yaks have a lower degree of domestication and a large LD decay distance. Strong artificial selection in livestock causes an increase in homozygosity in a particular region.⁴⁵ Analysis of ROH in eight yak breeds indicated that compared with the other seven yak breeds, the ROH of Meiren yaks were shorter and less in number. This indicated that the possibility of inbreeding in Meiren yaks was small, and no high-intensity inbreeding selection was observed. This is consistent with the results obtained from the analysis of LD in this study. The population differentiation index revealed that the average genetic differentiation index of the eight yak breeds was > 0.04. This was consistent with the results of the genomic feature map in the previous analysis. The average genetic differentiation index of Meiren yaks was 0.05; therefore, it could be used as a yak population genetic resource.

Animal groups living in various regions have experienced natural selection, artificial domestication, and migration. Because of these processes, some species have become extinct; some species have differentiated into different subgroups or subspecies, and even some new species have formed. Selection signal analysis can provide us insights into the traces that organisms have left on their genomes during evolution.^46,47 Therefore, the detection of selection signals can help to mine genes associated with economic traits of animals, understand the underlying genetic basis, and resolve the mechanisms of interpopulation differences.^48,49 GO analysis of the annotated genes in the screened genomic region in yaks revealed that they were mainly enriched in the functions such as N-acetylglucosamine phosphotransferase activity and regulation of amacrine cell differentiation. KEGG enrichment analysis revealed that the annotated genes were mainly involved in biological processes such as pantothenic acid and CoA biosynthesis, α-linolenic acid metabolism, PPAR signaling pathway, and cytokine–cytokine receptor interaction. ACCA1, as the key enzyme of fatty acid β-oxidation in peroxisome, plays an important role in regulating fatty acid metabolism.⁵⁰ Wang et al.⁵¹ in their study on sheep precursor adipocyte differentiation reported that ACAA1 deficiency promoted sheep precursor adipocyte differentiation while upregulating the expression of lipogenic marker genes; when ACAA1 was overexpressed, it repressed the expression of these markers. This suggested that ACAA1 affects lipid metabolism in sheep by regulating lipogenic marker genes. Similarly, some studies suggested that ACCA1 plays an important role in regulating lipid metabolism.^52,53

In this study, we found some regulatory genes (BMP2 and ACVR2B) that may be related to the adaptation of yaks to the low-oxygen environment on the plateau. Animals living at high altitudes, such as yaks and Tibetan antelopes, have evolved various adaptation methods to adapt to low-oxygen environments. These adaptations include increased production of red blood cells and blood vessels and alteration in functions such as deposition of muscle fat. Bone morphogenetic protein 2 (BMP2) is involved in bone growth and development. Studies have reported that it plays a role in regulating erythropoiesis and angiogenesis,^54–58 which are adaptive changes in terms of the body’s response to hypoxic environment. In response to hypoxia, animals increase the production of erythropoietin (EPO), a hormone that stimulates red blood cell production. BMP2 and EPO have been reported to act together to promote the differentiation of erythroid progenitors into mature erythrocytes.⁵⁹ Hypoxia-inducible factor 1 alpha (HIF-1α) is a transcription factor that is activated in the hypoxia reaction, which in turn stimulates the production of vascular endothelial growth factor (VEGF) to promote the production of blood vessels. BMP2 can interact with HIF-1α to promote the expression of VEGF. This promotes vascular formation.⁶⁰ Moreover, BMP signalling is thought to play a key role in the heart and nerves.⁶¹ ACVR2B, also known as activator receptor 2B, is a protein-coding gene that plays an important role in cell differentiation, growth, and apoptosis. ACVR2B can reduce cardiomyocyte death and reduce cardiac damage by adjusting cardiomyocyte metabolism in response to hypoxic environments.⁶² Furthermore, activation of ACVR2B can promote stem cell differentiation to repair damaged cells. These regulations may be the adaptations of Himalayan marmots to extreme environments.⁶³ These studies indicated that research is still limited on the genes regulating animals’ adaptability to the plateau habitat and hypoxic environment. Further research is needed to fully understand the molecular mechanism underlying the adaptation of animals to high altitude.

Conclusion

In summary, we analysed eight yak breeds through whole genome resequencing combined with population genetics to understand the evolution and domestication of these yaks. The results revealed that Meiren yaks have formed unique genetic characteristics in the long-term evolution process and is a potential yak genetic resource. Several genes that may be associated with adaptability to high altitude were obtained in Meiren yaks. This study enhanced our understanding of yak populations in the plateau and provided an important reference for the conservation of genetic resource of Meiren yaks.

Authors’ contributions

Conceptualization, T.W., X.M. and P.Y.; methodology, T.W.; software, T.W. and X.M.; formal analysis, T.W.; investigation, T.W.; resources, X.M., C.M. and C.L.; data curation, C.M., C.L. and X.W.; writing – original draft preparation, T.W.; writing – review and editing, C.M., C.L. and X.W.; visualization, P.Y.; supervision, P.Y.; project administration, P.Y.; funding acquisition, P.Y. All authors contributed to the article and approved the submitted version.

Institutional review board statement

All animal studies conducted in this study were approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (SYXK-2014–0002).

Acknowledgements

The authors would like to thank the Key Laboratory of Yak Breeding Engineering in Gansu Province for providing the experimental conditions and the teachers and brothers for their help.

Disclosure statement

The authors declare no conflict of interest.

Data availability statement

The datasets generated during the current study are available in the [FigShare] repository, [10.6084/m9.figshare.23723817].

Additional information

Funding

This research was funded by the Major Science and Technology Special Projects in Gansu Province (21ZD10NA001,GZGG-2021-1), the National Beef Cattle Industry Technology & System (CARS-37), the Agricultural Science and Technology Innovation Program (25-LZIHPS-01), the Science and Technology program of Gansu Province (20JR5RA580), the National Natural Science Foundation of China (32102524) and the Shenzhen Science and Technology Program (KCXFZ20201221173205015).