2.1 Ethics Statement
All of the animal procedures were performed in strict accordance with the guidelines and regulations proposed by the Animal Science Research Institute of Iran. All the animal experiments were approved by ethics committee of the Animal Science Research Institute of Iran under the number ASRI-34-64-1357-005-970180. Blood samples were collected during qualified veterinary treatment within the framework of governmental programs aimed at the animal identification, monitoring of health, and parentage confirmation of the dromedary populations included in our study. No other kind of tissue was used in this study.
2.1. Animals and Sample Collection
Yazd province with the area of 129,285 km2 (49,917 sq mi) is situated at an oasis where the Dasht-e Kavir desert and the Dasht-e Lut desert meet and located in 31° 53′ 41.28″ N, 54° 21′ 25.2″ E (31.8948, 54.357). Data on 51 herds of dromedaries were collected in 2018. The mean herd size was 89.90 ± 100.90 among 4279 camels, 16% were younger than one year, and named Hashi, while 12 % were older than ten years. The proportion of females (Arvaneh), males (Lok), male calves (Hasi), and female calves (Hashi) were 76%, 9%, 6%, and 9%, respectively. Among pregnant camels, 49% were older than ten years. The ratio of pregnant camels to all female camels was 46%. A total of 964 calving was registered between January to May 2018, distributed over 22% in January, 28% in February, 27% in March, 15% in April, and 8% in May. Generation intervals in females and males were estimated at 7.84 and 5.91 years, respectively. Among 256 male calves, we recorded 96 samples from 5 regions including: Bafgh (n = 41) Bahabad (n = 8), Khatam (n = 17), Mehriz (n = 8), and Ardakan (n = 22) (Figure 1). Characteristics of the sampled herds and the rangeland plants are presented in Table 1.
Table 1 characteristics of sampling herds and the rangeland plants
region
|
Herd size
|
Sampling site
|
n samples
|
The rangeland plants
|
Bahabad
|
400
|
31°52'29.6"N 56°01'14.7"E
|
8
|
Seidlitzia Rosmarinus, Artemisia spp, Salsola yazdiana, Tamarix tetragyna, Alhagi camelorum, Calligonum comosum, Zygophyllum spp
|
Khatam
|
112
|
30°28'42.4"N 54°12'36.7"E
|
17
|
Seidlitzia rosmarinus-Haloxylon ammodendron, Artemisia sieberi-Seidlitzia Rosmarinus, Tamarix tetragyna, Alhagi camelorum, Salsola yazdiana-Calligonum polygonoides, Zygophyllum atriplicoides
|
Ardakan
|
450
|
32°31'40.8"N 55°14'32.9"E
|
22
|
Salsola yazdiana, Haloxylon ammodendron, Artemisia sieberi , Zygophyllum spp
|
Mehriz
|
320
|
31°35'14.4"N 54°25'44.4"E
|
8
|
Haloxylon ammodendron-Zygophyllum atriplicoides, Artemisia sieberi-Salsola arbuscula,
|
Bafgh*
Herd1
Herd2
|
117
100
|
31°37'05.0"N 55°24'26.9"E
|
23
18
|
Zygophyllum spp, Seidlitzia Rosmarinus, Artemisia sieberi, Zygophyllum spp, Salsola yazdiana
|
* Two herds in Bafgh were recorded (1: local herd and 2: National Research and Development Station on Dromedary Camel herd)
2.2. Phenotypic Measurements
Data were recorded at the morning before grazing on the pasture. Camels were kept in closed area at night, which is called Garch. The animal identification was inferred via three-digit ear tags. Due to large distance in remote regions and transport difficulties, we constructed a portable weighting scale, consisting of 13 pieces of iron, a digital scale for 2000 kg, and one chain crane (Figure 2). The meta data collected for any calf included: ID number, characteristics of owner, geographical region, recording date, birth date, parental names, and body weight. We recording intervals were approximately three months, with the first record starting in the calving season, the second during the summer, and the third at weaning season at the beginning of autumn. We collected 252 body weight records from 96 calves in different times during 2018. The 18 calves belong to National Research and Development Station on Dromedary Camel (Bafgh), measured in 8 -times, the others were recorded 2 or 3-times including: Bafgh (n=164), Bahabad (n=8), Khatam (n=26), Mehriz (n=9), Ardakan (45).
For adjusting the body weight trait, the growth trend was estimated using linear regression model. Analysis of covariance (ANCOVA) of body weight was did among five sampling regions using SPSS v.22 software18. The mean of differences was compared using LSD test. Daily weight gain was calculated by two or three-times records.
2.3 DNA extraction, SNP discovery and genotyping
Blood samples were collected from the jugular vein in EDTA tubes during routine veterinary treatment on the pasture. The genomic DNA was extracted using the modified salting-out method19 and, after elution, was quantified using spectrophotometry and checked for quality on a 1% agarose gel. Finally, DNA samples were adjusted to a concentration of 50 ng/µl for subsequent steps.
The samples were genotyped-by-sequencing using two restriction enzymes combination, EcoR1 and HinF1 (New England Biolabs, Ipswich, MA, USA), and paired-end (150 bp) sequencing (10 X) on the Illumina HiSeq 2000 platform by Persian Bayangene Research and Training Center (Shiraz, Iran).
The sequence reads were mapped to the dromedary reference genome assembly on chromosome level (GCA_000803125.3;[1]) by using the BWA-MEM algorithm of Burrows–Wheeler Aligner (BWA)20;. PCR duplicates were detected by utilizing Picard tools and disregarded in downstream analyses both by GATK21 and SAMtools22. SNPs were called across the GBS data using GATK.
2.4 Population structure andd genetic diversity
A quality control (QC) steps, and genome-wide diversity (observed and expected heterozygosity), as well as admixture analyses were performed using TASSEL V5.023. Variants with a minor allele frequency (MAF) below 0.05 and call rate below 0.95 were removed. Of the 41,897 SNPs, 256 markers were monomorphic, and 27,375 markers were deleted because of MAF < 0.05. The final data set consisted of 14522 SNPs and 96 individuals. To investigate population structure, we used vcfR package23 in R for data manipulation and quality control as for producing input file objects for the other analysis, after that using ape and poppr package24 we carried out K-means clustering and discriminant analysis of principal components (DAPC), while all graphics were produced by means of RColorBrewer25.
2.5 Linkage disequilibrium analysis (LD) and SNP-based haplotype blocks estimating
TASSEL 5.026 was used to calculate the linkage disequilibrium (LD) (r2) for all pairwise loci. Haplotype blocks (HAP) were constructed in Haploview27.
2.6. Genome-wide association studies and candidate genes prediction
The association between the SNPs and the traits were tested using mixed linear models with PCA and kinship matrix in TASSEL software26. The regions of this study and age at weighting date were used as a fixed and covariate effect, respectively.
The statistical analysis model, the MLM-PCA + K analysis, was expressed as:
y = αX + βP + UZ + e
where y was phenotype value; α was the vector of SNP effects; β was vector of population structure effects based on PCA; u was vector of kinship background effects; e was vector of residual effects; X, P, Z were incidence matrix relating the individuals to fixed marker effects α, fixed principal component (PC) effects β, random group effects u, respectively. The suggestive significant Bonferroni P-value thresholds were set (-log p value > 3.9) using the GEC software tool28. The associated SNPs (-Log p value > 3) was traced in NCBI and the candidate genes were detected by blasting to the dromedary camel’s genome (GCA_000803125.3).
2.7. Bayesian genomic prediction
The estimation of Genomic Breeding values (GEBVs) was performed with the BGLR package including BRR, Bayes A, B, and C approach (nIter=100k, Burn In=5k)29. two sets of SNPs were used to predict GEBV: (1) all 14522 SNPs and (2) the 99 associated SNPs (-Log p value ≥ 2.5 from GWAS). The prediction accuracy was estimated using the average Pearson’s correlation (r) and regression (b) coefficient between the GEBVs and observed values30,31, 32. The replicated training – testing approach (10 replications) was used for evaluation of the models. We also used 3:1 size ratio of training set and validation set randomly selected from the 96 camels, which is a three-folds cross-validation, and repeated 100 times for evaluation of models by the 99 associated SNPs33,34.