Screening of Genes Related to Milk Production Traits in Bactrian Camels in the Eastern Xinjiang Region Based on Whole Genome Resequencing Technology

1. Introduction

The ancestors of camels originated in Canada, North America [1], and later migrated to Central Asia due to changes in the natural environment, where they gradually evolved humps, giving rise to dromedary and Bactrian camels [2].

In the 19th century, camels were primarily used for transportation. Since then, their economic value has increasingly been reflected in milk production, meat production, wool production, tourism, and cosmetics [3]. Camel milk is rich in nutrients such as fat, protein, minerals, unsaturated fatty acids, and amino acids, and also contains various bioactive compounds such as lactoferrin, lysozyme, immunoglobulin, lactoperoxidase, and insulin [4]. These compounds mediate immune regulation [5], antibacterial and anti-inflammatory effects [6], anticancer properties [7,8], anti-Parkinson’s effects [9], antioxidant activity [7], neuroprotective effects [10], and therapeutic effects on diabetes. Moreover, studies have shown that camel milk does not contain β-lactoglobulin and has a high content of α-lactalbumin, which is similar to human milk. Therefore, camel milk can be considered an alternative to cow’s milk for infants and young children prone to allergies [11].

According to the Food and Agriculture Organization of the United Nations (FAO) statistics for 2023, there are approximately 1.3 million Bactrian camels worldwide, with the largest populations in China (579,700), Mongolia (473,900), Kazakhstan (264,900), and Uzbekistan (18,000) [12]. In 2024, China’s “National Catalogue of Livestock and Poultry Genetic Resources” included five local breeds of Bactrian camels: Alxa Bactrian camel and Sunite Bactrian camel from Inner Mongolia, Tarim Bactrian camel and Dzungarian Bactrian camel from Xinjiang, and Qinghai camel. According to data from the “National Statistical Yearbook of China” (ISBN 978-7-5037-8743-1), there were 579,700 camels in China in 2023, distributed across Xinjiang, Inner Mongolia, Gansu, Qinghai, and other regions. Among them, Xinjiang and Inner Mongolia had the largest populations, with 316,000 and 206,000 camels, respectively. Currently, the development of the Bactrian camel industry faces several prominent issues, including insufficient milk supply, significant variations in the nutrient and functional substance content of milk from different breeding regions, and unclear quality characteristics of the milk. Additionally, compared to other dairy animals, there has been less research on the molecular breeding of Bactrian camels. Studies on genes related to lactation traits and their regulatory mechanisms are not yet sufficiently in-depth, and a large number of genes remain functionally unverified. Therefore, selecting camel groups with high milk yield and good milk quality for breeding will help produce high-quality camel milk and promote the healthy and sustainable development of the camel milk industry.

Screening lactation-related genes and candidate gene regions in Bactrian camels is of great significance for improving their lactation traits [13]. Currently, in camel molecular breeding, studies have shown that casein genes [9] and PICALM genes (phosphatidylinositol-binding clathrin assembly protein) [14] are associated with lactation traits in dromedary camels. Gubin et al. [15] found that two SNP sites in the intergenic regions of OSBPL8, MRPL37, SSBP3, and LOC102516351 genes in Gobi Red Camels from Inner Mongolia, China, are linked to the content of conventional nutrients (fat, protein, and lactose) in camel milk. In addition, Yao et al. [16] discovered that the secretion of GnRH (gonadotropin-releasing hormone), mTOR, PI3K-Akt, and MAPK signaling pathways is related to the regulation of lactation traits in Bactrian camels. Relatively few genes have been reported to regulate lactation traits in Bactrian camels. Milk production traits are typically regulated by multiple genes. Therefore, this study employed whole-genome resequencing to analyze the population structure and identify selection signals in two Bactrian camel populations: those from Northern Xinjiang and Eastern Xinjiang. The goal was to screen for lactation-related genes and regulatory signaling pathways, providing important molecular markers for the genetic selection and breeding of Bactrian camels in the Xinjiang region.

2. Materials and Methods

2.1. Selection of Experimental Population and Genomic DNA Extraction

The experimental population consisted of 106 Bactrian camels, including 51 from the eastern part of Xinjiang and 55 from the northern part of Xinjiang. The northern population was sourced from three different farms, while the eastern population came from a single farm. Both populations were raised under stall-feeding conditions. Information on the experimental population is provided in Table 1. Blood samples were collected from the carotid artery between March and October 2024. All blood samples were stored in anticoagulant tubes and immediately frozen at -80°C. Previous analyses of the nutritional components in camel milk from the eastern and northern regions revealed differences in nutrient content. Genomic DNA was extracted from the blood samples using a magnetic bead method. The concentration of DNA samples was measured using a Qubit fluorometer. The integrity and purity of DNA samples were assessed by 1% agarose gel electrophoresis. Only qualified samples were used for library preparation.

2.2. Whole-Genome Resequencing

The purified genomic DNA from Bactrian camel blood was sent to Huazhi Biotechnology Co., Ltd. for sequencing. Sequencing was performed on the MGI DNBSEQ-T7 platform by BGI Genomics. The library construction process included DNA fragmentation, end repair, adapter ligation, amplification, fragment selection, and circularization to generate high-quality DNA fragments ranging from 300 to 500 bp. The size of the library fragments was verified using a Qsep400 Bioanalyzer. After completing the library construction, sequencing was conducted on the platform.

2.3. Data Processing and Alignment

To ensure the removal of sequencing adapter sequences, low-quality bases, and undetermined bases, the raw sequencing data were filtered, checked for contamination, and assessed for quality. The quality control process included the following steps: (1) Filtering of low-quality reads, where read pairs were removed if more than 50% of the bases had a quality score of Q ≤ 20; (2) Removal of adapter sequences from the reads; (3) Filtering out reads with excessive N bases, specifically those containing more than 5 Ns; (4) Elimination of reads shorter than 100 bases.

For the reference genome, the study selected the genome sequence (GCA_048773025.1) from the NCBI database (https://www.ncbi.nlm.nih.gov/), with a reference genome size of 2.5 Gb. The clean reads obtained after quality control were aligned to the reference genome. Sentieon was used to detect variant sites in each sample, followed by joint calling to analyze the gVCF (genomic Variant Call Format) files of all samples collectively, yielding variant results for each individual in the population.

2.4. Variant Detection

The SNPs obtained from the joint analysis were subjected to initial hard filtering based on the following criteria: QD < 2.0, FS > 60.0, MQ < 40.0, SOR > 3.0, MQRankSum < -12.5, and ReadPosRankSum < -8.0. Subsequently, Vcftools software was used for further filtering of the SNPs with the following criteria: missing rate > 0.1 and MAF < 0.05.

2.5. Population Structure Analysis

Principal component analysis (PCA) was performed using PLINK based on the SNPs detected through whole-genome sequencing. Genetic distances were calculated, and a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method.

2.6. Detection of Selection Signals

The composite likelihood ratio (CLR) was employed to detect selection signals within populations. The detection principle was based on the Site Frequency Spectrum (SFS). In addition to detecting selection signals within individual populations, signals were also explored across different populations to identify molecular regions associated with adaptation to diverse ecological environments and the formation of high-yield phenotypes. Three methods were used for this purpose: fixation index (Fst), nucleotide diversity (π), and cross-population extended haplotype homozygosity (XP-EHH). Furthermore, KEGG pathway analysis and Gene Ontology (GO) enrichment analysis were conducted.

3. Results

3.1. Data Quality Control

3.1.1. Raw Data Preprocessing

To ensure data quality, the raw data were filtered and assessed for quality before analysis. After filtering low-quality data using FASTP, the total amount of clean data generated from sequencing for all samples was 2770.3 Gb, with an average of 26.13 Gb per sample. The average Q30 value across all samples was 96.64%, and the average GC content was 40.77%. These results indicate that the sequencing data are authentic and reliable. The data quality for some samples are shown in Table 2.

3.1.2. Reference Genome Alignment

The alignment rate of sequencing data (Clean reads) to the reference genome (MappedReadsRate), sequencing depth (MapDepth), and coverage (Cov1X) are key indicators reflecting the quality of samples, library construction, sequencing, and the reference sequence. The results show that the average MappedReadsRate reached 99.01%, indicating a high alignment success rate, with a coverage depth of 96.14%, which meets the requirements for alignment analysis. The alignment results of some samples’ sequencing data to the reference genome are shown in Table 3.

3.1.3. Genetic Variant Analysis

A total of 6,451,453 SNPs were detected across the entire genome. The distribution of these SNPs is shown in Figure 1. The majority of mutations are located in intergenic regions (3,003,030, 46.55%), followed by intronic regions (2,422,399, 37.55%). A smaller proportion of SNPs are found in the 1 kb upstream region of genes (111,624, 1.73%) and the 1 kb downstream region of genes (124,037, 1.92%), as well as in exonic regions (63,231, 0.98%). A small number of SNPs are located at splice sites (the 2 bp near the exon/intron boundary within introns) (277, 0.004%), and in regions that are both the 1 kb upstream region of one gene and the 1 kb downstream region of another gene (9,838, 0.15%). The remaining SNPs are distributed in other regions (717,017, 11.11%).

Functional annotations were performed on the SNPs located in exonic regions. The results revealed that the majority of SNPs were synonymous mutations (32,739, 54.16%) and non-synonymous mutations (27,289, 45.14%). A smaller proportion of SNPs were mutations that introduced stop codons (366, 0.61%) or removed stop codons (60, 0.10%).

3.2. Population Genetic Structure Analysis

3.2.1. Principal Component Analysis

To investigate the genetic relationships and clustering patterns among the two major Bactrian camel populations in Eastern and Northern Xinjiang, as well as among different camel populations from various farms within the same region, we performed Principal Component Analysis (PCA) using PLINK. The first three principal components, which explained 10.12%, 7.11%, and 6.71% of the total variance, were selected for analysis. These three components effectively distinguished the Bactrian camel population from Group 3 in Northern Xinjiang and the Bactrian camel population from Group 4 in Eastern Xinjiang from the other populations (Groups 1 and 2). However, the three populations in Northern Xinjiang showed significant admixture and could not be clearly separated by clustering, indicating a close genetic relationship. The PCA results showed that the overall genetic structure could differentiate Bactrian camels from Northern and Eastern Xinjiang, although there was some admixture between the two regions. The PCA results are shown in Figure 2a.

3.2.2. Phylogenetic Tree Construction

Genetic distances were calculated using PLINK, and a phylogenetic tree was constructed using the Neighbor-Joining (NJ) method (Figure 2b). In conjunction with the PCA results, the phylogenetic tree shows that the Bactrian camels from Eastern Xinjiang and Northern Xinjiang are distinctly separated into different clusters, indicating a relatively distant genetic relationship but with some gene flow between them. In contrast, the Bactrian camels from the three camel farms in Northern Xinjiang are closely related.

3.3. Detection of Genomic Selection Signals

3.3.1. Single-Population Selection Signals

We treated the 55 camels from Northern Xinjiang and the 51 camels from Eastern Xinjiang as two distinct populations. The SWEED software was used to detect selection signals within each of these populations. This method relies on the composite likelihood ratio (CLR) test to identify selection signals within populations, based on the site frequency spectrum (SFS). The genome was divided into 10-kb windows. For single-population selection signal detection in the camels from Northern and Eastern Xinjiang, we selected the top 5% of regions. This process identified 5,017 genes in the Northern Xinjiang population and 4,405 genes in the Eastern Xinjiang population, accounting for approximately 12.93% and 11.35% of the genome, respectively. The Manhattan plots of selection signals within the Bactrian camel populations from Eastern and Northern Xinjiang are shown in Figure 3.

For the 4,405 candidate genes from the 51 Bactrian camels in Eastern Xinjiang and the 5,017 candidate genes from the 55 Bactrian camels in Northern Xinjiang, we performed GO and KEGG enrichment analyses using the DAVID software. Selected GO and KEGG terms are presented, with the visualization results shown in Figure 4. The GO analysis results indicate that the 51 Bactrian camels from Eastern Xinjiang are enriched in 20 significant GO terms (Figure 4a), covering three major categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Among them, 447 genes are enriched in protein binding. After PDR correction of the P-values, the terms signal transduction, axon, and protein binding show highly significant enrichment. To understand the metabolic pathways in which the candidate genes of the Eastern Xinjiang Bactrian camels are involved, we performed KEGG pathway enrichment analysis on the 4,405 candidate genes. A total of 20 significantly enriched pathways were identified in the Eastern Xinjiang Bactrian camel population (Figure 4b). Among them, 67 genes were enriched in the Calcium signaling pathway. The three pathways of Glutamatergic synapse, Long-term depression, and Calcium signaling pathway reached highly significant P-values. The GO analysis results show that the 55 Bactrian camels from Northern Xinjiang are enriched in 20 significant GO terms (Figure 4c), covering three major categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Among them, 634 genes are enriched in the plasma membrane. After PDR correction of the P-values, the terms axon guidance, cell adhesion, and intracellular signal transduction show highly significant enrichment. To explore the metabolic pathways in which the candidate genes of the 55 Bactrian camels from Northern Xinjiang are involved, we conducted KEGG pathway enrichment analysis on the 5,017 candidate genes. A total of 20 significantly enriched pathways were identified in the Northern Xinjiang Bactrian camel population (Figure 4d). Among them, 73 genes were enriched in the Calcium signaling pathway. The P-values for four pathways—Calcium signaling pathway, signaling pathway, Pancreatic secretion, and Axon guidance—reached highly significant levels.

3.3.2. Inter-Population Selection Signals

To identify the molecular regions underlying adaptation to different geographical environments and the formation of superior milk quality traits in Bactrian camels from different populations in Northern and Eastern Xinjiang, China, we employed three methods—fixation index (Fst), genomic heterozygosity (π), and cross-population extended haplotype homozygosity (XP-EHH)—to detect inter-population selection signals. Candidate genes were identified and subjected to pathway enrichment analysis.

3.3.2.1. Fst-Based Method

Using the fixation index (Fst) method, we divided the genome into 50-kb windows with a step size of 25 kb and selected the top 5% of differentiated regions as candidate regions. A total of 2,391 candidate genes were identified between the Bactrian camels from Eastern Xinjiang and those from Northern Xinjiang (Figure S1a). We performed GO and KEGG enrichment analyses on these differentiated genes. We selected the top 20 GO terms for display (Figure S2a). The GO enrichment results showed significant enrichment in three major categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Three highly significant GO terms were identified: adherens junction, phosphatidylinositol phospholipase C activity, and cysteine-type endopeptidase inhibitor activity, mainly concentrated in CC and MF. A large number of candidate genes were enriched in protein binding. The KEGG enrichment analysis results (Figure S2b) showed that the majority candidate genes were enriched in Pathways in cancer. Three significant pathways (P < 0.05) were identified: Axon guidance, Pathways in cancer, and Ras signaling pathway.

3.3.2.2. π-Based Method

Using the genomic heterozygosity (π) method, we divided the genome into 50-kb windows with a step size of 25 kb and selected the top 5% of differentiated regions as candidate regions. Between the Bactrian camel populations from Eastern Xinjiang and Northern Xinjiang, a total of 2,334 differentiated candidate genes were identified (Figure S1b). We performed GO and KEGG enrichment analyses on these differentiated genes. The GO enrichment results (Figure S2c) showed that 135 candidate genes were enriched in the ATP binding term. The candidate genes were highly significantly enriched (P < 0.01) in three terms: one-carbon metabolic process, plasma membrane tubulation, and regulation of axonogenesis. The KEGG enrichment analysis results (Figure S2d) indicated that the majority of candidate genes were enriched in the Axon guidance pathway, which was the only highly significant pathway (P < 0.01).

3.3.2.3. XP-EHH-Based Method

Using the cross-population extended haplotype homozygosity (XP-EHH) method, the top 5% of differentiated regions were selected as candidate regions. In the 51 camel samples, a total of 3,982 genes were identified in the selected regions (Figure S1c). GO and KEGG enrichment analyses were performed on these differentiated genes.The top 20 GO terms were selected for display. The analysis results showed that the plasma membrane term had the highest number of enriched genes, reaching 522. A total of 17 GO terms were highly significantly enriched (P < 0.01), as shown in Figure S2e. The KEGG enrichment analysis results (Figure S2f) indicated that the Calcium signaling pathway had the highest number of candidate genes. Axon guidance was the only highly significant pathway (P < 0.01).

3.3.2.4. Joint Detection by Fst, π, and XP-EHH

Given Xinjiang’s vast land area and diverse geographical environments, the long-term differences in natural environments and selective breeding practices between Eastern and Northern Xinjiang have led to significant changes and differences in the physical characteristics and economic traits of Bactrian camels in these regions, which enhances the production efficiency of Bactrian camels. To identify candidate regions in the Bactrian camel population of Eastern Xinjiang compared to those of Northern Xinjiang, we used three different population selection signal methods—Fst, π, and XP-EHH (Figure 5a)—to detect inter-population selection signals. These three methods collectively identified 459 overlapping genes. By intersecting these genes with those related to Bactrian camel lactation, we annotated eight genes associated with lactation traits in Bactrian camels (see Table 4), namely DUOX1, DUOX2, CPQ, CGA, PLCG1, FYN, GNRHR, and TRHR.

To understand the signaling pathways in which the candidate genes are involved in camels, we performed GO and KEGG analyses on the 459 candidate genes. The results of the GO enrichment analysis (Figure 5b) showed that 11 pathways were included in BP (Biological Process), 4 pathways were included in CC (Cellular Component), and 5 pathways were included in MF (Molecular Function). Eight pathways were highly significant, with the thyroid hormone generation term showing a particularly significant difference. Additionally, a large number of genes were enriched in the plasma membrane term. The results of the KEGG enrichment analysis (Figure 5c) showed that the Neuroactive Ligand-Receptor Interaction pathway had the highest number of enriched genes (15). The pathways related to Axon guidance, Neuroactive Ligand-Receptor Interaction, and Thyroid Hormone Synthesis showed high significance levels. Among these, the Thyroid Hormone Synthesis pathway is associated with lactation traits.

4. Discussion

4.1. Differences in Nutritional Composition of Camel Milk from East and North Xinjiang Regions

Research suggests that the nutrient and functional substance content in camel milk is influenced by various factors, such as breed [17], lactation period [18], feeding methods, sampling season [19], and genetics. Additionally, the tools and methods used for milk testing and analysis [20] may also affect the measured nutrient content in the milk. In this study, samples were collected from camel farms in Hami City, East Xinjiang, as well as Yining County and Huocheng County in North Xinjiang. The camels in both regions were of the Xinjiang Junggar Bactrian breed. However, the two regions are approximately 1,500 kilometers apart, with significant differences in geography, climate, and feeding management practices. At the start of the study, the conventional nutritional components of camel milk from both regions were analyzed. The results revealed that the average values for conventional nutritional components in North Xinjiang camel milk were 4.44% milk fat, 4.04% milk protein, 5.83% lactose, and 15.00% total solids. In contrast, the average values for East Xinjiang camel milk were 5.55% milk fat, 4.11% milk protein, 5.76% lactose, and 16.86% total solids. Based on these results, it can be preliminarily concluded that camel milk from East Xinjiang has slightly higher nutritional content than that from North Xinjiang.

4.2. Molecular Marker Progress

This study identified eight genes: DUOX1, DUOX2, CPQ, CGA, PLCG1, FYN, GNRHR, and TRHR. Existing research suggests that DUOX1 and DUOX2 are primarily involved in the thyroid hormone synthesis signaling pathway. There is cross-regulation between thyroid hormones and prolactin (PRL), and hypothyroidism can reduce sensitivity to PRL. In dairy cows, thyroid hormone levels (T3/T4) are positively correlated with the content of conventional nutrients in milk, within the normal range [21]. This suggests that the DUOX1 and DUOX2 genes may play a role in the lactation process of Bactrian camels.

Moreover, the expression level of the DUOX1 gene is significant in normal human keratinocytes [22]. In Iranian Markhoz goats, both the DUOX1 and DUOX2 genes are located within the quantitative trait locus (QTL) region associated with shearing weight at one year of age [23]. Through transcriptome sequencing and protein interaction network analysis of Bactrian camel populations with and without mastitis, the PLCG1 gene was identified as a major immune core gene in Bactrian camel mastitis [24]. Additionally, the PLC gene family has been shown to play a role in targeted therapies for breast cancer [25]. In this study, the PLCG1 gene was enriched in multiple pathways, with the axon guidance signaling pathway reaching a highly significant level (P < 0.01). In the study of lactation mechanisms in Holstein dairy cows, it was found that the axon guidance signaling pathway influences lactation performance [26], further indirectly validating that the PLCG1 gene could serve as a candidate gene for lactation traits in Bactrian camels. Hussain Bahbahani et al. [27] conducted whole-genome sequencing of dromedary camels and found that the FYN gene plays a role in fat metabolism and energy expenditure, consistent with the results of this study. The GNRHR gene, which encodes the gonadotropin-releasing hormone receptor, plays a significant role in mammalian reproduction [28] and regulates hormone levels, such as estrogen and progesterone. These hormones significantly impact the development of mammary glands and lactation regulation in mammals. Therefore, the GNRHR gene is likely associated with lactation traits in Bactrian camels, suggesting that the findings of this study are highly reliable. The roles of the CGA and TRHR genes in regulating camel phenotypic traits and their mechanisms of action have not yet been reported.

4.3. Pathway Selection and Analysis in This Study

In this study, the enriched pathways of the differentially selected candidate genes between the East Xinjiang Bactrian camel and the North Xinjiang Bactrian camel primarily focus on axon guidance, neuroactive ligand-receptor interactions, and thyroid hormone synthesis. In mammals, signaling pathways associated with lactation traits include RAP1, Wnt, MTOR, PI3K-AKT, and MAPK pathways, which are mainly involved in mammogenesis, mammary gland development, and lactation [29,30,31]. However, compared to dairy cows and other mammals, Bactrian camels exhibit significantly different mammary gland structures and lactation patterns. There are also notable differences in the signaling pathways involved. The findings of this study provide support for this perspective. Additionally, based on previously measured differences in the nutrient content of milk from East and North Xinjiang Bactrian camels, the camels were categorized into two groups. The relatively small sample size of the Bactrian camels selected is a limiting factor. Nonetheless, this study offers valuable biological insights into the unique lactation mechanisms of Bactrian camels and lays the groundwork for further research on the nutritional characteristics of Bactrian camel milk.

5. Conclusions

This study utilized whole-genome resequencing technology to analyze the lactation-related genetic differences between Bactrian camel populations in East Xinjiang (from a single camel farm) and North Xinjiang (from three camel farms across two counties and cities), identifying a total of 6,451,453 SNPs. Principal component analysis revealed a clear distinction between the two populations. GO and KEGG enrichment analyses indicated that, compared to North Xinjiang, the candidate genes in East Xinjiang Bactrian camels were primarily enriched in the plasma membrane and thyroid hormone synthesis pathways. Whole-genome selection signals identified eight candidate genes potentially associated with lactation traits: DUOX1, DUOX2, CPQ, CGA, PLCG1, FYN, GNRHR, and TRHR. These genes may serve as key markers for future molecular breeding.