Phenotypic data, along with milk and blood samples, were meticulously collected from 125 lactating Alxa Desert Bactrian camels located at the Alashan League, Inner Mongolia. All camels were in a consistent lactation period, with a parity ranging from 2 to 5, and were maintained in optimal body condition within a uniform environment. They were fed the same diet: dry clover supplemented with 2 kg of grain concentrate (68% corn +12% wheat bran +20% soybean cake after oil extraction) and 30 g of table salt for each animal daily. Blood and milk samples were collected from 32 and 83 Bactrian camels at approximately 30 and 270 days postpartum, respectively. Daily milk production was estimated based on the average milk yield over three consecutive days. For compositional analysis, 25 mL of milk was sampled to determine the concentrations of key constituents such as fat, protein, and lactose using mid-infrared spectrometry (MilkoScan Minor, Foss Analytics, Hillerød, Denmark). Additionally, 10 mL of whole blood was collected from the jugular vein of each camel for genetic analysis. Blood samples were treated with TRIzol reagent (TaKaRa, United States) and stored at −80°C until RNA extraction could be performed.

Genomic DNA was collected from 398 Bactrian camels used for chip genotyping. These camels were from six key breeding regions in China, which are major areas for Bactrian camel breeding. Specifically, samples were obtained from Subei County, Gansu (GSB; 7 individuals), Alxa Right Banner, Inner Mongolia (NAY; 82 individuals), Alxa Left Banner, Inner Mongolia (NAZ; 257 individuals), Sunit Right Banner, Inner Mongolia (NSNT; 22 individuals), Siziwang Banner, Inner Mongolia (NSZW; 10 individuals), and Urad Rear Banner, Inner Mongolia (NWLT; 20 individuals). The dataset comprises three breeds of Bactrian camels: Alxa Desert Camel (n = 343), Alxa Gobi Camel (n = 22), and Sunit Bactrian Camel (n = 33). Alxa Desert Camels were sampled from GSB, NAY, and NAZ, while Alxa Gobi Camels were exclusively sampled from NWLT, and Sunit Bactrian Camels were sampled from NSNT and NSZW. The methodologies used in this study were approved by the Institutional Animal Care and Use Committee of Inner Mongolia Agricultural University, Hohhot, China.

We isolated total RNA from 125 blood samples using the TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The blood samples were homogenized in TRIzol reagent and chloroform, followed by precipitation using isopropanol. Total RNA from each sample was treated for genomic DNA contamination using the RNase-free DNase set (QIAGEN, Crawley, West Sussex, United Kingdom) and purified using the RNeasy mini kit according to the supplied guidelines (QIAGEN, Crawley, West Sussex, United Kingdom). RNA sample quality was assessed using the NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States) and the Agilent Bioanalyzer 2100 system. RNA samples exhibited 28 S/18 S ratios ranging from 1.8 to 2.0 and RNA integrity number values between 8.0 and 10.0.

For mRNA cDNA library preparation, 1.0 μg of total RNA was utilized from each sample using the TruSeq RNA Library Preparation kit v2 (Illumina, San Diego, CA, United States). Poly A-containing mRNA was enriched from the total RNA using poly-T oligo attached beads and fragmented for first-strand cDNA synthesis, followed by second-strand synthesis. The ends were repaired, and 3′ end adenylation and adapter ligation were performed for each library. Subsequently, libraries were polymerase chain reaction (PCR) amplified, validated using the Bioanalyzer (Agilent Technologies Inc., Cedar Creek, TX, United States), and finally normalized and pooled. Clustering of the index-coded samples was conducted on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. Following cluster generation, the library preparations were sequenced on the Illumina NovaSeq 6000 high-throughput sequencing platform, generating 150 bp paired-end reads.

The identification of SNPs was executed following the workflow depicted in Figure 1, with the comprehensive SNP calling pipeline script provided in Supplementary file S4. Fastp v0.23.4 (19) was employed to assess the quality of sequence reads, targeting the identification of sequencing read artifacts such as low-quality Phred scores, duplicated reads, uncalled bases (N sequences), and potential contamination. Subsequently, the filtered reads from each sample underwent individual alignment to the Bactrian camel reference genome (Ca_bactrianus_MBC_1.0) using Hisat2 v2.1.0 software (20). The aligned results were exclusively utilized for Single Nucleotide Polymorphism (SNP) calling, while insertions and deletions (indels) were excluded from the analysis due to challenges associated with accurate indel calling. Variant calling was performed for each read merging method using the “mpileup” and “call” commands from BCFtools v1.9–77-gd0cf724+ (21). Only those SNP variants were retained where the alternative allele manifested in all samples, accompanied by a Phred quality score of at least 25 and a minimum read depth of 10. We calculated the minor allele frequencies (MAFs) and missing rates using PLINK v1.90 (22), and the SNPs with MAF <0.05 were excluded. We tested the Hardy–Weinberg equilibrium (HWE) using VCFtools v0.1.11 (23) with the -hwe option and removed SNPs that severely departed from HWE (p < 0.01).

To assess the stability and reliability of the “CamelBell No. 1” SNP chip, we performed genotyping on 24 lactating camels randomly selected from the 398 DNA samples. The sample set included three pairs of duplicate samples to evaluate reproducibility. Next, the sequencing reads were first subjected to quality control using FastQC v0.11.5 (29) to assess the quality of the raw data. Low-quality reads and adapter sequences were trimmed using Trimmomatic v0.39 (30). The resulting data were aligned to the Bactrian camel reference genome (Ca_bactrianus_MBC_1.0) using BWA v0.7.17 (31). Post-alignment, the alignment rate was calculated to assess sequencing efficiency and accuracy. The alignment files were processed using SAMtools v1.17 (21) to convert, sort, and index the aligned reads. Alignment quality was further assessed by examining metrics such as the percentage of mapped reads, coverage depth, and uniformity of coverage across the genome. To ensure the chip’s precision and reliability in genotyping, we meticulously analyzed the genotype concordance rate among the duplicate samples. This analysis is vital for validating the chip’s capability to produce consistent and reproducible results. Additionally, we quantified the detection rates of all test samples across various site coverage depths, providing a comprehensive evaluation of the chip’s performance under different genomic conditions.

To evaluate the performance of the SNP array in detecting population genetic diversity and structure, we genotyped and analyzed 398 DNA samples. VCFtools v0.1.11 (23) was employed to filter the SNPs and individuals. Specifically, individuals and SNPs with a detection rate lower than 95% were excluded. SNPs with minor allele frequencies lower than 0.05 were excluded. SNPs with significant deviation from Hardy–Weinberg equilibrium (p < 0.01) in any population were excluded. The genotyping data were extracted and converted into PLINK. bed and .fam formats, and then imported into PLINK v1.90 (22) and ADMIXTURE v1.3.0 (32) software for analysis. Using PLINK, PCA is conducted to identify major sources of genetic variation across the six breeding areas. ADMIXTURE is used to infer population structure and assign individuals to genetic clusters. It provides estimates of the proportion of an individual’s genome that originates from each of K ancestral populations. The tested K was set from 2 to 9.

The GBLUP model is used to calculate GEBV as follows:

where y is the vector of corrected phenotypic values. These are the phenotypic measurements (e.g., milk yield, milk protein, milk lactose milk fat) that have been adjusted for fixed effects (e.g., breed, parity, lactation period). 1 is a vector of ones, which ensures that the overall mean ( μ ) is added to each observation, with 𝜇 being the overall mean of the phenotypic values across all individuals in the study. X β represents the fixed effects, where 𝑋 is the incidence matrix for these effects, and 𝛽 is the vector of fixed effects coefficients. 𝑍 is the incidence matrix that relates the additive genetic values (𝑢) to the phenotypic values (𝑦), with each row of 𝑍 corresponding to an individual and each column to a genetic effect. 𝑢 is the vector of additive genetic values (or genomic breeding values, GEBVs). These values represent the genetic contribution of each individual to the trait. It is assumed that 𝑢 ∼ N(0, Gσ_𝑢²), where G is the relationship matrix built with the HIBLUP v1.4.0 software (33). This matrix represents the genetic relationships between individuals based on SNP marker data. σ_𝑢² represents the additive genetic variance. 𝑒 is the vector of random residual effects, representing the variation in the phenotypic values not explained by the model. It is assumed that 𝑒∼N (0, Iσ_𝑒²), where I is an identity matrix and σ_𝑒² is the residual variance.

Cross-validation (CV) is usually used to obtain a reliable and stable model, and to evaluate the quality of the model (34). To evaluate the accuracy of genomic prediction, we utilized five-fold cross-validation. The dataset was divided into five approximately equal-sized groups. In each iteration of the cross-validation, four groups were used as the training set to estimate model parameters, while the remaining group was used as the validation set to test the model’s predictive performance. Prediction accuracy was calculated as the correlation between the predicted estimated breeding values (EBVs) and the actual phenotypes in the validation set, divided by the square root of the heritability estimated in the validation population. The standard error (SE) was calculated as the standard deviation of the five calculated reliability values from the five-fold cross-validation, divided by the square root of five.

A total of 19,177 confident SNPs were detected from the RNA-seq data of lactating Bactrian camels. Initially, transcriptomes of peripheral blood from 125 selectively chosen individuals, exhibiting consistent lactation period out of 1,243 Bactrian camels, were sequenced. This sequencing yielded approximately 3,575 million paired-end reads, averaging 29 million paired-end reads per individual sample. After quality control, the base effectiveness rate stood at 97.88%. The Q30 ratio surpassed 93.14%, and the GC content exceeded 52.08%. Remarkably, about 91.72% of the reads were accurately mapped to the Bactrian camel genome (Ca_bactrianus_MBC_1.0), with nearly 80.54% of the reads from each individual uniquely aligning with the camel genome. Detailed alignment information for each sample is tabulated in Supplementary Table S1. Subsequently, a comprehensive SNP calling was conducted in the reference genome from 125 transcriptome datasets, unveiling 9,790,170 SNPs. Only those SNP variants where the alternative allele appeared in all samples and had a Phred quality score of at least 25 and a minimum read depth of 10 were retained. Minor allele frequencies (MAFs) and missing rates were calculated, with SNPs having MAF <0.05 excluded. Additionally, we tested for Hardy–Weinberg equilibrium (HWE) using the -hwe option and removed SNPs that significantly deviated from HWE (p < 0.01). Among these, 19,177 SNPs were confidently ascertained across the entire cohort.

For the lactation performance traits of milk yield, milk protein, milk lactose, and milk fat, both Exon Wide Selection Signature (EWSS) and Genome Wide Association Study (GWAS) methodologies were employed to generate a comprehensive SNP set. Initially, the calculation of F_ST was applied to 19,177 SNPs to discern associations with extreme phenotypic milk traits. Groupings according to extreme values of milk traits are shown in Table 1. The statistical analysis of other milk components in each extreme value group is shown in Tables 2–5. Using a threshold of F_ST > 0.2, we identified 178 loci as outliers associated with various milk traits. Specifically, 138 loci were related to milk yield, 28 to milk protein, 30 to milk lactose, and 6 to milk fat. These loci are depicted in Figure 2. These identified loci serve as a partial reference for the selection of SNPs associated with lactation traits. Simultaneously, leveraging the 19,177 SNPs identified within the transcribed regions, we also employed GWAS to associate them with four continuous lactation performance traits across 125 subjects, thereby providing an additional reference set for the screening of pertinent loci. Therefore, significant associations were found with 923 SNPs for milk yield, 980 SNPs for milk protein, 955 SNPs for lactose, and 756 SNPs for milk fat traits (FDR <0.05), shown in Figure 3. During the analysis, we calculated the genomic inflation factor (λ) to assess the potential inflation of test statistics due to population structure or cryptic relatedness. The results showed a genomic inflation factor of 0.945, indicating minimal inflation and suggesting that the test statistics were not significantly affected by population structure or cryptic relatedness (Figure 4A). The overlap of SNPs associated with different traits is shown in Figures 4B,C. Overall, the amalgamation of both methods furnished us with a reference set comprising 2,960 SNP loci, which are associated with four lactation traits and dispersed within the transcribed regions of 1,395 genes, shown in Supplementary Table S2. A total of 81 overlapping sites were identified between F_ST and GWAS analyses. The functional enrichment analysis conducted on the identified genes indicated a significant overrepresentation in several biological processes (BP), particularly in areas related to the immune system, organonitrogen compound metabolism, and regulation of metabolic processes. Additional processes such as vesicle-mediated transport, positive regulation of biological processes, protein metabolism, and cellular localization establishment were also notably enriched. Furthermore, genes involved in catabolic processes were highlighted as part of the biological process category. To provide a comprehensive overview, the top 10 items from each category—Biological Processes (BP), Cellular Components (CC), Molecular Functions (MF), and Kyoto Encyclopedia of Genes and Genomes pathways (KEGG) were visually summarized in Figure 4D.

To explore the variation in SNP density across different chromosomes, a box plot was generated (Figure 5A). The median distances and the ranges vary significantly, with some scaffold showing a higher concentration of closely spaced SNPs. To mitigate the incidence of false positives in genotyping, SNPs were judiciously selected for inclusion on the chip, guided by criteria such as SNP quality, polymorphism status, SNP density, and the extent of association with specific traits. In brief, we submitted 2,960 previously identified candidate SNPs associated with lactation traits to the Targeted Capture Sequencing Probe Design System, and 2,681 SNPs passed evaluation. Next, considering the location and function of these SNPs, we eliminated one of the two sites within less than 100 bp from each other, focusing on retaining sites supported by existing literature. Thus, 1,002 SNPs were finally manually selected for the “CamelBell No. 1” array. For those SNP site, there were 460 SNPs located in exon region, 54 located in intron region and 10 located in intergenic region, shown in Figure 5B, and the details are in Supplementary Table S3. The association of these SNPs with lactation traits is visually represented in Figure 5C, with 201 SNPs being correlated with two or more lactation traits.

The efficacy of the SNP array was assessed through the genotyping of 24 DNA samples from Alashan League, Inner Mongolia. Inclusion of three pairs of replicate samples facilitated the evaluation of the chip’s detection performance. A genome alignment rate exceeding 99% was achieved for all samples, culminating in an average alignment rate of 99.77%, shown in Figure 5D. Separate testing of three duplicate samples yielded a typing consistency of 100% for each pair, underscoring the chip’s detection stability (Figure 5E). Comprehensive statistics pertaining to the detection rates across all test samples revealed that with a target site coverage depth exceeding 5×, an average site detection rate of 99.72% was maintained. Remarkably, even at a coverage depth surpassing 20×, the detection rate remained robust at 99.07%, shown in Figure 5F.

To assess the genetic diversity of the core breeding populations in of dairy Bactrian camels in northwest China, we collected 398 camels from six representative breeding areas, and conducted SNP genotyping using the SNP array. Firstly, we calculated the minor allele frequency (MAF) distribution of the population and plotted the histogram, as shown in Figure 6A. The MAF data from SNP array genotyping revealed that 95% of SNPs had a minor allele frequency greater than 0.1. Subsequently, we estimated the genetic structure and distance, observing minimal genetic divergence among samples from different breeding areas (Figures 6B,C). These results indicate that camels from six breeding areas share a similar genetic background regarding lactation functionality.

Genomic prediction (GBLUP) accuracy for the four milk traits was assessed using genotype information from 914 markers that passed the quality control filters (MAF >0.01, Call Rate >0.9). A total of 398 camels were randomly divided into training (80%) and validation (20%) sets for cross-validation, and this process was repeated five times to ensure robust evaluation. The box plot of prediction accuracy values is shown in Figure 7. The genomic prediction accuracies for the four milk traits vary, with milk lactose showing the highest and most consistent accuracy, followed by milk yield and milk fat, while milk protein exhibits the lowest and most variable accuracy. The mean Pearson correlation coefficients of the four milk traits were 0.27, 0.30, 0.34, and 0.44, respectively, indicating modest linear relationships between the predicted Genomic Estimated Breeding Values (GEBVs) and the actual phenotypes. Among these traits, milk fat exhibited the strongest correlation. The standard error (SE) values of the estimated correlation coefficients for all traits were less than 0.10. The mean regression coefficients for the four milk traits were 0.88, 0.76, 1.14, and 0.69, respectively. These values suggest that while predictions for milk yield are relatively unbiased, predictions for milk protein and milkfat tend to underestimate the actual genetic values, and predictions for milk lactose tend to overestimate them. The SE values of the estimated regression coefficients for all traits were less than 0.06.