Genetic structure of Arabian Peninsula dromedary camels revealed three geographic groups.

Dromedary camels (Camelus dromedarius) are widespread in the desert and semi-desert areas of Africa, the Arabian Peninsula, some parts of southwest Asia and Australia. In the Arabian Peninsula, these well-adapted species have been classified based on their ecology into Desert camels, found mainly in the north and center of the Peninsula, Mountain camels, distributed along the west and south of the Peninsula, and Beach camels, populating the west to southwest of the Peninsula. Here, we aimed to investigate the genetic relationship between 386 camels corresponding to 12 dromedary populations from different geographical locations and ecology in the Arabian Peninsula with the genotyping of 17 microsatellite loci. No significant deviation was observed in heterozygosity, allelic richness, Fis (inbreeding coefficient) among the studied populations had a mean value of 0.5849, 4.808 and 0.04, respectively. A mean Fst (fixation index) value of 0.0304 was calculated for the various populations with the highest value obtained between racing Omani and Awarik camel populations (0.079). Both the neighbor-joining phylogenetic tree and the STRUCTURE analysis divided the populations into three different groups corresponding to their Arabian Peninsula geographic location (North, Central and West, South-West, and South-East of the Arabian Peninsula), rather than their ecological classification, with a high level of genetic admixture and gene flow among them. Investigating the genetic relationship of dromedary populations in the Arabian Peninsula can be considered as the first milestone to conserve this well-adapted species. The results obtained here need to be further validated using whole genome sequencing data.

Introduction

Dromedary camels (Camelus dromedarius), also known as Arabian camels, are single-humped even-toed ungulates populating the desert and semi-desert environment of the Arabian Peninsula, Africa, southeast Asia and Australia (Burger et al., 2019, Wilson, 1998). This species, in addition to the domesticated two-humped Bactrian camels (Camelus bactrianus) and the two-humped wild camels (Camelus ferus), belongs to the Old World Camel tribe (Camelini) of the Camelidea family. This tribe diverged about 16.3 million years ago from the New World camel tribe (Lamini), which inhabits the Andes plateau in South America (Burger et al., 2019).

Dromedary camels have different anatomical and physiological properties that confer their adaptability to the surrounding harsh environmental conditions. Anatomically, they have large feet with thick pads that help them to walk easily on soft sands (Dagg, 1974, Hoter et al., 2019). They also ensure air circulation around their body upon raising their sternum at a recumbent position (Ouajd and Kamel, 2009). Physiologically, dromedary camels are able to drink a large amount of water, approximately 200 L in 3 min, to withstand walking for long distances without drinking water (Ouajd and Kamel, 2009). Upon water deprivation, they can fluctuate their internal body temperature from 6 °C to 8 °C (42 °C to 34 °C), to avoid body water evaporation, a mechanism known as adaptive heterothermy (Hoter et al., 2019, Tibary and El Allali, 2020). These physiological characteristics have supported the use of this well adapted species as mean of transport over large geographic distances for trading and travel. Dromedaries were used in trading goods, such as spices and perfumes, between the southern and northern parts of the Arabian Peninsula. They have also historically connected the African continent and the Arabian Peninsula transporting ivory, wool and incenses on both sides of the Gulf of Aden (Al-Ghabban et al., 2010).

Dromedary camels are spread over different geographic regions of the Arabian Peninsula, which are associated with different ecological locations, in which dromedaries have been classified according to them. The north of the Arabian Peninsula is mainly characterized by desert areas, which also cover part of the center and east of the Peninsula. Different camel types are populating this habitat and are referred to as desert camels, e.g. Magaheem, Wodeh, Sofor and Shual. Magaheem camel type, in the Arabic language means the darker color, is widely distributed across the Arabian Peninsula. They are characterized by a dark brownish to black coat color, also known as Malah. This camel type is divided into two main subtypes based on the tribe that owns them, Ad-Dwsaria and Almarria. Ad-Dwsaria is associated the Adwasire tribe that originated from the region of the Wadi Ad-Dawasir, in the center of the Arabian Peninsula. Whereas Almarria belongs to the Almurah tribe that lives in the eastern part of Saudi Arabia and other Gulf countries. This camel population has a high potential for milk and meat production (Wardeh, 2004). Wodeh (also named Magateer) are amongst the commonest camel types in the Arabian Peninsula. The term “Wodeh” means bright color in Arabic as this camel type is characterized by white coat color. They are found in the North of the Arabian Peninsula (North of Saudi Arabia, Kuwait and Iraq). It is reared for meat and milk production, as well as for beauty contests due to its unique coat color. Sofor and Shual are common camel types in the northern part of the Arabian Peninsula. They are characterized by brownish coat color, with uniquely in Sofor a dark brownish to black color at the top of their neck, shoulder, hump and tail. Both Sofor and Shual camels are used for milk, meat and wool production.

Camels populating the west of Arabian Peninsula are mainly distributed on the mountains and along the Red Sea coast, and hence are divided into Hill or Mountain camels and Beach camels (Almathen, 2014, Almathen et al., 2018). Hadana and Awadi are the Mountain camel types populating the mountainous areas in the west and southwest of the Arabian Peninsula, respectively. Hadana camels are mainly distributed in Saudi Arabia in the mountainous area of the Al-Baha region. They are smaller in size and shorter compared to other camel types with light brownish coat color. Awadi camels are populating the Faifa mountain in the Jazan region up to the Yemen border of Saudi Arabia. They are similar to Hadana in shape with remarkable shorter hair. Mountain camels are also found in the south of the Arabian Peninsula, such as Omani camel populations in the Dhofar region of Oman.

The Beach camels, Sahlia and Awarik are found mainly near the Red Sea coast of Saudi Arabia. Awarik camels are characterized by a light brown coat color with short hair. Their names derived from a plant, they are feeding on, growing in the Jazan region called Alarak. Sahlia camels, which resemble Awarik camels but with a darker brown coat, are populating a wider geographical zone than the Awarik camels from (Almadinah) west of Saudi Arabia to the Jazan region. Besides Desert, Hill and Beach camels, there are two racing camel populations in the Arabian Peninsula, Hurra and Omani. Hurra camels are found in the northern part of Saudi Arabia in the Al-Jawf region of Tabarjal town. Morphologically, they have a fine head with a flat forehead, narrow feet and long legs. Omani racing camels, which are similar to Hurra in shape, populate mainly in the East of Oman (Al-Sharqiya) as well as in most regions of the United Arab Emirates (UAE).

All of the aforementioned dromedary classification does not follow a standard registered breeding system as in other livestock populations, such as sheep (“Jacob Sheep Breeders Association,” 2009) and cattle (“The Dexter Cattle Society,” 2019). Establishing a standard genetically informative breeding system needs to be progressed with a detailed characterization of the genetic structure and diversity of these camel types. Genetic research efforts in this field have mainly been linked with autosomal microsatellite markers (Mburu et al., 2003, Ould Ahmed et al., 2010, Almathen et al., 2016, AlAskar et al., 2020, Hossam Mahmoud et al., 2020) and partial mitochondrial DNA (mtDNA) sequences (Almathen et al., 2016), with the exception of (Bahbahani et al., 2019) and (Ming et al., 2020) studies which reported genotype-by-sequence (GBS) and whole genome sequence results, respectively. Genetic studies based on microsatellite markers indicated a level of genetic differentiation between indigenous dromedary camels. Two genetic groups present in Kenya dromedary populations, namely the Somali population and a group including the Turkana, Rendille and Gabra populations were defined by Mburu et al. (2003). Ahmed et al. (2010) also differentiated the Tunisian dromedary populations into two distinct genetic entities, namely the Nafzawa (Kebili) and the Aaradh group, the later including the Medenine and Tataouine populations. A study by Almathen et al. (2016) on dromedary populations from Africa and Asia detected a phylogeographic distinction separating the East African dromedary populations from the Arabian Peninsula one. Excluding these East African dromedary populations, high level of historical gene flow among the Asian and African populations were also revealed. A similar continental-wise genetic distinction has been observed in AlAskar et al. (2020). Likewise, a recent study by Mahmoud et al. (2020) on four dromedary populations from Saudi Arabia (Magaheem, Wodeh, Sofor and Shual) indicated a high level of genetic admixture among them.

In this study, we aim to investigate the genetic relationships and differentiation using 17 microsatellite loci in 12 dromedary populations representing different geographic regions and ecological classifications at the Arabian Peninsula. Our main objective is to assess whether or not these populations living in different ecological and geographic areas may be distinguished at the genetic level.

Materials and methods

Sample collection

A total of 386 blood samples were collected representing dromedary camel populations from different geographical locations in the Arabian Peninsula (Table 1 and Fig. 1A). These unrelated samples, based on questioning the camel owners, were collected during routine veterinary treatments from different camel owners in three countries (Saudi Arabia, Oman and UAE). The blood was sampled from the jugular vein and stored in 5 ml Magic® storage buffer tubes (Biogen Diagnóstica, Spain).

Table 1

The sample size, geographical locations and ecological classifications of the Arabian Peninsula dromedary camel populations included in the study.

Population	Sample size	Geographical location	Ecological classification
Hurra	33	North of Arabian Peninsula	Desert
Sofor	31	North of Arabian Peninsula	Desert
Shual	35	North of Arabian Peninsula	Desert
Magaheem	48	Center of Arabian Peninsula	Desert
Wodeh	26	North of Arabian Peninsula	Desert
Sahlia	27	West of Arabian Peninsula	Beach
Hadana	50	West of Arabian Peninsula	Hill/ Mountain
Awarik	33	South-West of Arabian Peninsula	Beach
Awadi	26	South-West of Arabian Peninsula	Hill/ Mountain
United State of Emirates (UAE) (Racing camels)	31	South-East of Arabian Peninsula	Desert
Omani P (Production camels)	21	South-East of Arabian Peninsula	Hill/ Mountain
Omani R (Racing camels)	25	South-East of Arabian Peninsula	Desert

Fig. 1

A) Geographical locations of sampled dromedary camels. NCW (green color): samples from the north (Hurra, Wodeh, Sofor and Shual), center (Magaheem) and west (Sahlia and Hadana) of the Arabian Peninsula. SE (red color): samples from the southeast of Arabian Peninsula (Omani and UAE dromedary populations). SW (blue color): samples from the southwest of Arabian Peninsula (Awarik and Awadi). B) Unrooted neighour-joining phylogenetic tree of the dromedary samples translated from the pairwise fixation index (Fst) values. C) STRUCTURE plots of K = 2 to K = 5 for the Arabian Peninsula dromedary populations. Hur: Hurra, Sfr: Sofor, Shl: Shual, Mghm: Magaheem, Wdh: Wodeh, Sah: Sahlia, Had: Hadana, Awrk: Awarik, Awd: Awadi, UAE: Camels from UAE, Omani P: Production Omani camels and Omani R: Racing Omani camels.

DNA extraction and microsatellite genotyping

Genomic DNA was extracted from blood using the DNeasy® Blood and Tissue Kit (Qiagen, Hilden, Denmark). Samples were genotyped with the 17 microsatellite markers (STR) used from the previous study (Almathen et al., 2016). All PCRs were carried out in a total volume of 2 μl containing 5 – 10 ng genomic DNA, 0.2 µM of each primer and the QIAGEN® Multiplex PCR master mix (Qiagen, UK). PCR cycles for all multiplexes entailed the initial denaturation step at 95 °C for 10 min, followed by 35 cycles, each consisting of denaturation at 94 °C for 1 min, annealing at 55 °C for 90 s, extension at 72 °C for 90 s and a final step at 72 °C for 10 min. A volume of 1 µl of 1:10 water diluted PCR product was mixed with a loading mix, containing formamide (Applied Biosystems) and Rox 500 size standard (Applied Biosystems), and then denatured for 3 min at 95C° and analyzed in a Genetic Analyzer ABI 3730 DNA sequence (Applied Biosystems). Electropherograms were evaluated using GeneMapper v. 3.7 (Applied Biosystems).

Microsatellite loci and dromedary populations summary statistics

Summary statistics were calculated at both microsatellite loci and dromedary populations levels using the summary function of adegenet package (Jombart, 2008) in R software version 3.6.2 (R Development Core Team, 2013). At loci level, the number of alleles per loci and the observed heterozygosity (Ho) were calculated. Deviation from Hardy-Weinberg equilibrium was also assessed per loci using the hw.test function of the pegas package (Paradis, 2010) in the R software. After correcting for multiple testing using Bonferroni correction, a p-value of (0.05/17 = 0.0029) was set as a threshold to define loci not in Hardy-Weinberg equilibrium. The total number of alleles, mean Ho, allelic richness and Fis values were calculated per population using the heirfstat package in the R software. The non-parametric Kruskal-Wallis one-way analysis of variance was used to test for differences among the dromedary populations for the mean heterozygosity, mean allelic richness and Fis values. All these values are not normally distributed as indicated by Shapiro-Wilk test (p-value < 0.05).

Genetic differentiation and molecular phylogeny

Pairwise genetic differentiation, expressed by Weir and Cockerham Fst value (Weir and Cockerham, 1984), was assessed among the different dromedary populations using the pairwise.WCfst function implemented in the hierfstat package of R software. These pairwise Fst values were used to build a neighbor-joining phylogenetic tree for the dromedary populations using the nj function in the ape package of R software.

STRUCTURE analysis

The genetic structure and the extent of admixture on the genotyped dromedary samples were assessed by the Bayesian clustering method implemented in STRUCTURE software version 2.3.4 (Pritchard et al., 2000) using the mixed ancestry admixture model. Ten independent simulations with independent allele frequencies were run for K = 1 to K = 12 using 50,000 iterations and a burn-in period of 10,000 Markov Chain Monte Carlo (MCMC) steps. CLUMPP software was used to concatenate the outputs of the multiple runs to be displayed in R software. The Delta K (ΔK) approach (Evanno et al., 2005) implemented in STRUCTURE HARVESTER (Earl and vonHoldt, 2012) was used to define the optimal number of clusters.

Results

The mean number of alleles per microsatellite locus is 7.882 alleles, with a minimum of two alleles observed in three loci and a maximum of 22 alleles found at a single locus. The mean Ho among microsatellite loci is 0.5859, with a minimum value of 0.2883 and a maximum value of 0.7775 (Table 2). After Bonferroni correction, a total of four loci were found to be deviated from Hardy-Weinberg equilibrium (Table 2).

Table 2

A description of the microsatellite loci used in the study including: the number of alleles in each locus, observed heterozygosity (Ho), and χ2 test P-value for the deviation from Hardy-Weinberg equilibrium.

Locus	Number of alleles	Ho	P-value
CMS09	8	0.704960836	0.787221882
CMS121	12	0.767624021	0*
CMS13	9	0.646753247	0.756965568
CMS15	10	0.7265625	0.31555063
CMS18	4	0.380829016	0.003556515
CMS25	7	0.603133159	0.362461874
CMS50	13	0.777479893	3.35e-09*
CVRL01R	22	0.772845953	0.109515848
CVRL04R	4	0.596354167	0.070702244
CVRL05R	12	0.663185379	5.62e-05*
CVRL06R	5	0.307291667	0.926751348
CVRL08	2	0.308510638	0.010110995
LAC66	7	0.688311688	0.049290151
VOLP10	7	0.746073298	0.950538007
VOLP32	2	0.288311688	0.134416688
YWLL44	8	0.545931759	3.40e-07*
YWLL59	2	0.435695538	0.068049145

Loci significantly deviated from Hardy-Weinberg equilibrium after Bonferroni correction.

Among the dromedary populations analyzed, the Awadi population has the lowest number of alleles over the microsatellite loci (73 alleles), while the Magaheem camels have the highest number of alleles with a total of 95 alleles. The mean Ho value of the dromedary populations is 0.5849 ± 0.041. The population-specific Ho values are not significantly different from each other (Kruskal-Wallis test P-value = 0.058), with the smallest value in the Awarik population (0.4982) and the largest one in the Hurra population (0.6329). The mean allelic richness is calculated to be 4.808 ± 0.36 with the minimum value in Awarik population (4.039) and the maximum value in Wodeh population (5.22). No significant difference in allelic richness is observed among the different dromedary populations analyzed (Kruskal-Wallis test P-value = 0.86). Wodeh camels have the lowest Fis value (- 0.0012), while camels from UAE show the highest value (0.106), with the mean Fis value of 0.04 ± 0.036 calculated among the different dromedary populations (Table 3). No significant differences are observed in the Fis values among the dromedary populations (Kruskal-Wallis test P-value = 0.806).

Table 3

Allelic richness, observed heterozygosity (Ho) and inbreeding coefficient (Fis) values of the dromedary populations including in the study.

Population	Allelic richness	Ho	Fis
Hurra	4.996	0.633	0.008
Sofor	5.086	0.621	0.015
Shual	4.801	0.609	0.018
Magaheem	5.027	0.615	0.044
Wodeh	5.221	0.622	−0.001
Sahlia	4.805	0.577	0.071
Hadana	4.761	0.571	0.084
Awarik	4.039	0.498	0.068
Awadi	4.168	0.528	0.032
UAE	4.922	0.558	0.106
Omani (production)	5.118	0.599	0.006
Omani (racing)	4.749	0.588	0.032

Molecular phylogeny and genetic differentiation

The unrooted neighbor-joining phylogenetic tree divides the dromedary populations into three main clades; 1) Awarik and Awadi camels, 2) Omani production camels and camels from UAE, and 3) Hadana, Magaheem, Hurra, Sofor and Shual camels. The Omani racing camels are genetically closer to the Omani camels (clade 2), whereas both of the Wodeh and Sahlia camels are closer to Hadana and Magaheem (clade 3) than the other clades (Fig. 1B).

The mean pairwise Fst value among the dromedary populations is 0.0304 ± 0.021. The highest genetic differentiation is observed between racing Omani camels and Awarik camels (Fst = 0.079), whereas the lowest differentiation is between Wodeh and Shoel camels (Fst = 0.001) (Fig. 2 & Supplementary Table S1).

Fig. 2

Heatmap of the pairwise (fixation index) Fst values between the dromedary populations. OmnR: Racing Omani camels, OmnP: Production Omani camels, UAE: Camels from UAE, Awd: Awadi, Awrk: Awarik, Had: Hadana, Sah: Sahlia, Wdh: Wodeh, Mghm: Magaheem, Shl: Shual, Sfr: Sofor and Hur: Hurra.

A substantial degree of genetic admixture is observed as indicated by the different K values of the STRUCTURE analysis on the dromedary populations. Based on the ΔK approach of Evanno et al. (2005), the optimal number of clusters is at K = 2 (Supplementary Table S2). At this K value, two genetic ancestries can be observed corresponding on one side to the dromedary populations from the North, Center and West, and on the other side the dromedary population from the South-East and South-West of the Arabian Peninsula. At K = 3, a third genetic ancestry background is observed in relation to the Awarik and Awadi camel populations in the southwestern part of the Arabian Peninsula. At K = 4, Hadana and Hurra camels emerge with separate genetic background ancestries. Separate genetic ancestries were also revealed for racing camels from Oman and UAE at K ≥ 5 (Fig. 1C).

Discussion

The genetic structure and relationship of dromedary camels from different geographical locations of the Arabian Peninsula have been investigated here by genotyping with microsatellite markers. The analyzed dromedary populations revealed a genetic distinction based on their geographical distribution, rather than their ecological classification, with a high level of genetic admixture observed among them.

The observed genetic admixture between the different Arabian Peninsula dromedary populations analyzed is reflected by the non-significant differences in allelic richness, observed heterozygosity, and inbreeding coefficient among them. This level of genetic admixture can be attributed to the historical use of dromedary camels in cross-continental transportation and trading, which was associated with high level of gene flow. Dromedary camels were considered as the main livestock species in connecting the south with the north of the Arabian Peninsula to transfer goods, such as spices and perfumes, through a route known as “incense road” (Epstein, 1984). Such connection has been widespread linking the different parts of the Arabian Peninsula, and even Africa, with the west of the Arabian Peninsula especially during the annual pilgrimage “hajj” when Muslims move to Makkah, which is located on the west of the Arabian Peninsula (Wilson, 1998). Genetic admixture on dromedary camel populations has been previously observed in dromedary camels from Saudi Arabia (Mahmoud et al., 2020), the Arabian Peninsula (AlAskar et al., 2020, Almathen et al., 2016), and in the African dromedary populations, such as in Sudan (Bahbahani et al., 2019) and Morocco (Piro et al., 2020).

In parallel to the genetic admixture noticed among the dromedary populations, a signal of the genetic distinction was also observed between them through both the phylogenetic tree and the STRUCTURE analyses that is related to their geographical distribution. Although the optimal cluster value of the STRUCTURE analyses, K = 2, differentiated northern, central and western dromedary populations from the southern populations, both the phylogenetic tree and K = 3 of the STRUCTURE analyses divided the dromedary populations in the Arabian Peninsula into three main geographical groups (Fig. 1A). The first group was composed of populations widespread in the north, center and west of the Peninsula; Wodeh, Sofor, Shual and Hurra (North), Magaheem (Center), Sahlia and Hadana (West). The second group comprised populations from the southwest of the Peninsula; Awarik and Awadi. Whilst, the third group is the Omani camels populating the southeast of the Arabian Peninsula. This level of grouping has also been observed previously in Almathen et al., 2016, AlAskar et al., 2020.

The second group of dromedary camels; Awadi and Awarik, is populating the Jazan region at the southwest of the Arabian Peninsula either on the top of mountains, as in Awadi, or along the coast of the red sea, as in Awarik, with the possibility to see some Awadi camels along the coast leading to a level of gene flow among them. These camel populations are isolated from the Hadana camels, which are populating the west of the Arabian Peninsula. This is due to the specific distribution of the Hadana camels at the mountains of Al-Baha region, which is about 500 km away from Jazan, and hence a separate genetic ancestry background appears for this group at K = 4. This type of population isolation, reflected by the overall lower level of Ho and allelic richness, and higher level of inbreeding compared to the dromedary populations from the north and center of the Arabian Peninsula, which have a greater level of movement increasing the chance of cross-breeding. Similarly, the Sahlia camels did not show shared genetic ancestry with other southwest Arabian Peninsula camel populations due to the widespread of Sahlia camels from Almadinah region on the west of the Arabian Peninsula to Jazan region, which covers more than 1000 km.

Separate genetic ancestry backgrounds corresponding to racing camel populations have appeared at K = 4 (Hurra camels) and K = 5 (Racing camels from Oman). These distinct genetic backgrounds might be attributed to the continuous selection for the racing ability trait from these populations for the past 200 years, which may lead to the formation of distinct camel breeds.

Conclusion

In conclusion we have investigated the phylogeography of dromedary populations from different geographical locations and ecological classifications of the Arabian Peninsula. The main outcome of this study revealed a genetic differentiation among the dromedary camels from the North, Central and West, South-West, and South-East of the Arabian Peninsula. Moreover, a signal of breeds formation with distinct genetic backgrounds has also been observed on the racing dromedary camel populations from the north and south of the Arabian Peninsula. These investigations need to be further validated using more powerful genomic tools, such as genotyping SNP arrays and whole genome re-sequencing data.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.