Genome-wide diversity and global migration patterns in dromedaries follow ancient caravan routes

We performed double-digest restriction site associated DNA (ddRAD) sequencing on 122 dromedary DNA samples from 18 countries (Supplementary Data 1) representative of the species distribution range. We included one Bactrian camel to test for potential interspecific hybridization, as this continues to be a widespread practice in Central Asia that might have started as early as pre-Roman times¹¹. Higher numbers of reads mapping to the Bactrian camel were detected in three individuals from Iran and in six from Kazakhstan (see “Methods”), and we decided to remove these samples from downstream analysis due to potential introgression from Bactrian camel (Supplementary Data 2). After stringent filtering for genotype and individual missingness, minor allele frequency and relatedness, the final dataset consisted of 95 dromedaries and 22,721 SNPs present in at least 75% of the individuals.

Moderate genome-wide diversity and low population structure

With 22,721 SNPs, we estimated expected (H_E = 0.27 ± 0.17; mean ± SD) and observed (H_O = 0.25 ± 0.17) heterozygosities in the global dromedary population (n_pop = 17; n_ind = 95). Separating the samples according to their continental origins, both Asian (n_ind = 49, H_E = 0.27 ± 0.17/H_O = 0.25 ± 0.17) and African dromedaries (n_ind = 46, H_E = 0.26 ± 0.17/H_O = 0.25 ± 0.18) showed similar genomic diversity. The mean H_E (t = −2.2641, df = 45,398, P = 0.02) and inbreeding coefficients (t = −2.5159, df = 43,024, P = 0.01) were higher in Asian than African dromedaries, but mean H_O (t = −1.2791, df = 45,385, P = 0.2) was not different between continents, according to Welch′s t test. Complete diversity and inbreeding values are given in Supplementary Table 1. In comparison with other domestic species, i.e., sheep (H_E = 0.22–0.32)²⁶ or cattle (H_E = 0.24–0.30)²⁷, we consider the genome-wide diversity in dromedaries as moderate at the best. Several bottlenecks during the last glacial period (see demographic analysis below, and Fitak et al.²⁴) and during domestication left modern dromedaries with a minimum of only six maternal lineages¹ and limited genome-wide diversity. This will have implications on future intensification of breeding and genomic selection in dromedaries from regions with increasing desertification.

In general, the genome-wide differentiation within the global dromedary population was very low. Analysis of Molecular Variance (AMOVA) revealed that most of the variation, ~94.3%, is explained within individuals (Supplementary Table 2). Allelic richness (AR) was similar between countries (AR = 0.25–0.27) with exception of Kenya which was lower (AR = 0.21). The pairwise fixation index between African and Asian individuals was very low (F_ST = 0.006; P < 0.001), and indices between dromedaries from different countries (if significant at all) were lowest in geographically close populations (e.g., Libya/Algeria: F_ST = 0.0002) and increased with geographic distance (Pakistan/Tunisia: F_ST = 0.0328) (Supplementary Table 3).

We screened for loci deviating from neutraliy using BayeScan 2.1²⁸ and identified sixteen F_ST outliers to be putatively under selection (false discovery rate (FDR) <0.05) between African and Asian dromedaries (Supplementary Fig. 1). We found it reasonable to investigate the biological functions of those genes harboring the SNPs as they might relevant for the adaptation of dromedaries to their respective environments. We found SNPs in two genes, CALN1 and TREM1, which are responsible for calcium ion binding and amplifying inflammatory responses triggered by bacterial and fungal infections, respectively (scaffold:SNP-location; JWIN01030783.1:128274 and JWIN01033764.1:729703). In addition, we examined (potentially linked) regions 200 kbp upstream and downstream of the F_ST-outlier loci and detected fifty-three genes related to a number of biological functions (Supplementary Data 3). Interestingly, around one fifth of the detected genes had functions related to the immune system hinting to an adaptive process in response to different pathogens in the respective environments. Other protein coding gene functions were related to pathways such as circadian rhythm, (ga)lactose, metabolism, reproductive or various cellular and developmental processes. A full list of genes is presented in Supplementary Data 3. Signatures of selection related to photoperiod, metabolism, immunity and growth have also been observed in chicken²⁹, sheep³⁰, and cattle (TGFB3³¹).

To understand the low genome-wide differentiation in dromedaries across their global range, we investigated population structure and admixture between populations. We projected the genetic variation of each dromedary on the first three axes inferred from a principal component analysis (PCA) and incorporated continental information (Africa/Asia) for each sample (Fig. 1). Principal component 1 (PC1) clearly separated African from Asian dromedaries, while PC2 and PC3 split Kenyan individuals from the rest of Africa and identified a single population from Saudi Arabia grouping closer to African than Asian individuals, although showing some cross-continental admixture (Fig. 1 and Supplementary Fig. 2). This separated population belongs to a specific breed, Hadhana, and is one of the twelve recognized dromedary ecotypes in Saudi Arabia, limited to mountain regions in the South of the Arabian Peninsula, Al-Baha³². In this case, the geographic accessibility might have an important role in the observed genetic distinctiveness. A possible explanation for the close relationship of Hadhana and African dromedaries might be the historic sea route from Jiddah in Saudi Arabia to Aydhab and Port Sudan. On the western coast of the Red Sea existed a trading route connecting the Horn of Africa to Petra and Damascus via Port Sudan, Aydhab and Myos Hormos, near today’s Kosseir (Fig. 2)^6,7. In general, the Asian dromedary population showed higher genetic variability, although the genetic variation explained by the three first axes was rather low with only 5.3% (Fig. 1). While this could be a sign for ancestral variation (the Arabian Peninsula was a center of domestication¹, we cannot discard the hyphothesis of post-domestication movements of camels or multiple origins of the founder populations as this would have left similar signals in the genomes.

We next inferred potential ancestry and admixture among Asian and African dromedaries using unsupervised genetic clustering in ADMIXTURE³³ (Fig. 3). Based on the lowest cross-validation error, the best clustering solution was 1 (Supplementary Fig. 3), which suggests a panmictic dromedary population and reflects the low genetic differentiation of 0.6% among individuals from different continents. Increasing the numbers of potential ancestral populations (K) from two to seven confirmed the already observed differentiation between African and Asian dromedaries (K = 2), the clustering of the Saudi Arabian Hadhana breed with Africa (K = 4), the separation of Kenyan and Hadhana individuals, and the higher number of distinct clusters on the Asian continent (K = 7). We find a more homogenous gene pool in African animals with the exception of the East African group¹ (Figs. 1, 2 and Supplementary Fig. 2), represented in our dataset by the two Kenyan dromedaries. This can be a consequence of a random founder effect followed by lack of gene flow due to geographical, physiological (e.g., Trypanosome infestation) and/or cultural barrier, i.e., dromedaries in East Africa were dominantly used for milk production rather than transport or riding¹. There is a need to proceed with comprehensive analyses about the potential nature of natural and/or anthropogenic obstacles for gene flow between East African and other dromedaries.

Effective migration rates along ancient caravan routes

To formally test our qualitative observations of weak population structure among African and Asian dromedaries (Figs. 1, 3), we visualized the global spatial population structure using the Estimated Effective Migration Surfaces (EEMS) method³⁴. Based on a stepping-stone model, EEMS detected a corridor of significantly higher effective migration rates than the overall mean along the Mediterranean coast, connecting the Northern parts of Africa and the Arabian Peninsula until the border of the Arabian Desert (Fig. 2). This pattern shows a continuous gene flow throughout the coastal dromedary populations, and a lower than average migration in the inland desert populations. A known trading route which fits this observed effective migration pattern bordered the Mediterranean coast connecting Northwestern Africa to the North of the Arabian Peninsula from where caravans traveled toward Southern Asia along the Silk Road (Fig. 2)^1,7. The introduction of the dromedary into Northern Africa via the Sinai from Roman Egypt started in the early first millennium BCE and intensified in the Ptolemaic period^6,15. From there, dromedaries migrated along the Mediterranean coast, as archeological evidence dates their presence in Northwest Africa to the fourth to seventh century CE (Late Antiquity/Early Middle Ages)^1,6. Even earlier dispersal of taurine cattle along the Northern coast of Africa and through the Mediterranean sea to Europe was described during the Bronze age³⁵.

It is clear that camels, unlike other domesticated species, were able to penetrate deep into the Saharan desert and to connect trans-Saharan cultures. West Sahara belonged to an Islamic trading network classified as one of the major gold suppliers in the ninth to tenth centuries CE³⁶ (Fig. 2). Tadmekka, a territory located in the Southwestern Saharan desert and governed by the Tuareg, was operational by the eighth century CE and was one of the earliest towns established in the region where cross-Saharan camel caravans traded³⁷. These trades prolonged at least until the fourteenth century when Timbuktu, which similar to Tadmekka hosted large groups of Islamic traders, engaged in coin-based exchange economies across the Sahara³⁶.

While modern dromedaries along the western part of the Silk Road are still well connected today, the Iranian and Afghanistan deserts seem to present obstacles of effective migration. As EEMS assumes uniform migration rates the observed “barriers”, could however be assessed as areas of lower population density with fewer migrants exchanged per generation, producing an effective “barrier” to gene flow³⁴. The three main parallel itineraries of the Incense Routes through the Arabian Desert connecting the South of the Arabian Peninsula with the Levant⁷ also showed lower than average migration rates (Fig. 2), which could be interpreted as lower population density or potential sampling gaps. These trading routes were essential during historical periods, not only for exchanging luxury products (e.g., incense or gems), but for trading everyday local products³⁸.

The strongest barrier detected in our dataset concerned dromedaries from the Horn of Africa, which had the lowest genetic effective migration rates (Fig. 2). Geographical isolation due to the Ethiopian highlands, which might disrupt gene flow with northern populations, and the Golf of Aden would be the most likely explanation for the observed pattern. Genetic differentiation of livestock populations in East Africa has been described previously^1,22,39.

Late Pleistocene population decline and medieval expansions

To complete our understanding of the moderate genome-diversity observed in the global dromedaries, we investigated the demographic history and inferred effective population size (N_E) over time with an Extended Bayesian Skyline approach (EBS)⁴⁰ (Fig. 4). As genetic differentiation between Asian and African dromedary populations was low (0.6%) but highly significant (P < 0.001), we tested whether the dromedary populations from the two continents experienced a similar demographic history. First, we investigated the global population and second, African and Asian individuals separately. With the latter approach, we accounted for a potential confounding effect of population structure for the inference of N_E⁴⁰. Due to the observed substructure in Kenyan and Hadhana individuals, we excluded these two populations from the continental groups.

Irrespective of the continental origin, all inferences showed similar patterns of an initial population expansion from one million ybp until ~700,000 ybp (Fig. 4). Our genome-wide population approach confirmed previous N_E estimates based on a single dromedary whole genome sequence²³ using pairwise sequentially Markovian coalescent⁴¹. This population expansion coincides with two remarkable periods: the middle Pleistocene transition (1.25–0.70 million ybp) characterized by climatic cycles, and the Galerian Mammal Age (1.2–0.60 million ybp), which influenced the distribution and evolution of biota and resulted in some species being adapted to arid, cold climates^42,43. Moreover, this timeframe also overlaps with the maximal diversity of the family Camelidae (early Galerian), supporting the adaptation of the dromedary ancestor to environmental changes with an expansion of its population during the middle Pleistocene transition^23,44. Population expansion was followed by a drastic decline in N_E beginning 700,000 ybp until the dromedary population collapsed during the last glacial period (LGP; 100,000–20,000 ybp)⁴⁵. This is a finding shared by previous Old World⁹ and New World camelid²³ studies and those focusing on Late Quaternary Megafauna⁴⁶.

Conversely, no bottleneck was picked up by any of the EBSPs during the time scale when dromedaries are predicted to have been first domesticated ~3000–4000 ybp^15,16. Previous BSP analysis using mtDNA likewise did not show a population decline during the time of domestication¹. It is possible that the detection of a bottleneck with the EBSP analysis related to domestication has been superimposed by the drastic decrease in N_E ending ~30,000 ybp (Fig. 4). Similar demographic changes were observed in alpacas²³, where three population bottlenecks were detected throughout the cold conditions of the LGM in South America, yet no bottleneck was visible during the domestication period.

After the Pleistocene bottleneck, the dromedary population slowly increased until reaching a stable N_E around 300 ybp, with a higher N_E present within Asia than in Africa. Demographic inferences based on mtDNA sequences described slightly earlier expansion of the maternal lineages around 600 ypb¹ associated with the rise of the Ottoman empire and the conquest of Constantinople (1453 CE), followed by the extension to Southern Asia and the Red Sea coasts⁴⁷.

Source: Ecology - nature.com