Historical domestication-driven population expansion of the dung beetle Gymnopleurus mopsus (Coleoptera: Scarabaeidae) from its last refuge in Mongolia

Populations of Gymnopleurus mopsus (family Scarabaeidae), a dung beetle that displays dung-rolling behavior (i.e., a telecoprid), have recently experienced sharp declines, and many populations are now at high risk of local extinction. However, Mongolia, which constitutes a major portion of the species’ distribution, still sustains a relatively large population. Here, we used mitochondrial COI sequences to investigate the within-population genetic diversity and both the genetic and phylogeographic structures of 24 G. mopsus populations across the species’ main distribution in Mongolia. Several lines of evidence indicated that the phylogeographic structure of G. mopsus had been influenced by a recent and sudden demographic expansion. Interestingly, the expansion of Mongolia’s G. mopsus population corresponded to the advent of livestock domestication in the region, and the species’ genetic structure coincided with road networks, which presumably serve as migration routes for livestock that might mediate the beetle’s dispersal. In addition, we also found that G. mopsus possesses high levels of haplotype diversity, which is generally indicative of large effective population sizes (Ne). Overall, the present study contributes to the current understanding of G. mopsus’ demographic history and dispersal patterns and also provides valuable information for the species’ conservation and management.

imens collected in Mongolia (N = 406), as well 12 specimens from Korea, China, and France, we identified 94 polymorphic sites and a total of 230 haplotypes (Table 1, Fig. 1). Only the 14 populations that were represented by more than ten samples (N > 10) were included in downstream population genetic analyses, whereas the samples from all 27 sampling localities were used for haplotype network analysis, and those from the 24 sampling localities in Mongolia were used for AMOVA (Table 1). Sequence correspond to GenBank accession numbers MF674025-MF674381.
Estimates of haplotype diversity (h) for the 14 large (N > 10) populations ranged from 0.889 to 1.000 in the ST23 and ST24 populations, respectively, with an average of 0.974, whereas estimates of nucleotide diversity (π) ranged from 0.006 to 0.010, with an average of 0.008 (Table 1). Furthermore, the number of haplotypes in these populations ranged from ten to 29 in the ST23 and ST18 populations, respectively, which were located in the central and southernmost of Mongolia. The most common haplotype (H19) was shared by 47 individuals (11.2%) and was obtained from all but seven of the localities in Mongolia (ST02, ST05, ST08, ST09, ST11, ST12, and ST18), and the next most common haplotypes, H21, H15, H12, H30, and H02, accounted for 4.5, 3.6, 2.9, 2.9, and 2.6% of the sampled individuals, respectively. Meanwhile, 181 of the haplotypes (43.4%) were obtained from single individuals and, thus, classified as unique haplotypes (i.e., singletons). Most of the localities yielded relatively high haplotype richness (HR) values, ranging from 5.616 to 10, with the lowest HR value obtained for the ST23 population, which was located in the center of our sampling localities. COI haplotype network. A haplotype network of COI sequences was star-shaped, with the most common haplotype (H19) located in the most internal position in the network (Fig. 2), thereby indicating presumably an ancestral haplotype. Furthermore, the ancestral H19 haplotype was detected from all but seven of the sampling localities in Mongolia (ST02, ST05, ST08, ST09, ST11, ST12, and ST18), which indicates its wide distribution. Related haplotypes that radiated from H19 were found repeatedly in different parts of the network, which suggests that local population expansion was followed by radiation from the ancestral haplotype. The haplotypes of individuals from Korea (N = 7, H25, H56, H80, H81, H83, and H84) and China (N = 3, H82 and H65) were located at the edge of the network but with only one or two mutational steps separating them from the haplotypes  Table 1. Information, localities, and diversity indices of Gymnopleurus mopsus samples analyzed in the present study. Significant values (P < 0.05) are indicated in bold. C, Central; E, East; N, North; S1, South1; S2, South2; S3, South3; Superscript P in locality denotes populations that were represented by more than ten samples (N > 10); *Indicates very close localities but collection in different years. Geographic population structure. Hierarchical analyses of molecular variance (AMOVA) revealed significant spatial genetic structuring among the 24 Mongolia populations, according to three different biogeographic groupings of (1) two groups, the humid areas (North and Central populations) and desert areas (South and East populations); (2) three groups, the North and Central area (NC), South area (S), and East area (E); or 3) three groups of South populations (S1: ST06, ST07, ST08, ST09, ST10, ST11, and ST18; S2: ST12, ST13, ST19, ST20, and ST21; and S3: ST05, ST14, ST16, and ST17). These groupings were based on our a priori expectation of the spatial population differentiation of G. mopsus, based on geographic proximity to one another and to the region's main roads (Table 2). We found that 3.06 and 2.17% of the molecular variance was partitioned between and among the groups, respectively (φ CT ; P < 0.01), whereas most of the variance occurred within populations (φ ST ; 96.74 and 97.60% of the total variation, P < 0.01; Table 2). However, only 0.20 and 0.23% of the molecular variance was partitioned among populations within groups (φ SC ; P > 0.05). The greatest between-group difference (3.06%) was observed when the four populations (ST01, ST15, ST22, ST23) in the northern humid area were compared to the 20 desert populations. However, the genetic structure (φ CT = − 0.002) among the three South groups (S1, S2, and S3) was insignificant (P = 0.775), which indicated the occurrence of relatively high levels of either historical or contemporary gene flow among populations in South area.

Genetic differentiation (Fst) and isolation by distance (IBD). The genetic differentiation between
populations was estimated using pairwise Fst statistics (φ ST ), which ranged from −0.019 to 0.079 ( Table 3). Most of the pairwise differences (φ ST values) were relatively small, even between distant populations (e.g., the southernmost and northernmost sampling localities, which were separated by 670 km), which generally indicated low levels of inter-population genetic differentiation. As a result, after sequential Bonferroni correction applied, only 18 of the 91 comparisons among the 14 large (N > 10) populations indicated significant differentiation (P < 0.05; Table 3). The highest φ ST value was observed between the ST23 and ST18 populations (Fst = 0.079, P < 0.01). However, neither the ST07 nor ST10 population was significantly differentiated from any of the other populations, which suggests the occurrence of either historical or ongoing gene flow between these and the other populations. Isolation-by-distance (IBD) analysis detected a significant positive relationship between genetic distance (Fst) and geographic distance (decimal degrees; R 2 = 0.087, P = 0.001) between the 14 large Mongolia populations (Fig. 3).  Table 1, Fig. 4). For example, both the Tajima's D and Fu's Fs statistics for neutrality were negative for all but the ST23 population, which was located in the central area of all the main roads meet in Mongolia (Fig. 1). Moreover, the neutrality test was significantly violated when the all of the 14 populations were pooled as a single data set (Table 1), with Tajima's D and Fu's Fs values of −1.76936 (P = 0.005) and −24.64263 (P = 0.001), respectively, which suggests that a sudden historical population expansion occurred following a population bottleneck. These significant negative values also indicate that the hypotheses of population stasis, selective neutrality, and population equilibrium can be rejected. A mismatch distribution was clearly smooth and unimodal for all of the 14 large populations, further supporting the model of sudden population expansion 18 (Table 1, Fig. 4). The SSD (sum of squared deviations) statistic indicated that the observed frequency distribution fit the distribution expected under the model of demographic  (15)  7(2) 32 (2)   8(3)   9(3) 33 (3) 41 (3) 85 (3) 86 (3) (10) 1 (5)  4(4) 48 (4) 26 (2) 28 (2) 71 (2) 56 (2) 42 (2) 36 (2) 95 (    Historical demographic reconstructions, based on Bayesian skyline plot (BSP) analysis, further supported a model of sudden population expansion (Fig. 5). Indeed, the Mongolia populations appear to have undergone a population expansion ~2500 years ago, followed by continuous increases in effective population size for a long period afterward (Fig. 5). This estimate might indicate that the timing of G. mopsus expansion corresponded to both the expansion of human populations and advent of large mammal domestication in Mongolia.

Discussion
To the best of our knowledge, the present study is the first to examine the current population genetic structure and historical population demography of G. mopsus spanning major distributional ranges across Mongolia using mitochondrial COI sequences. The high frequency of unique haplotypes (singletons) and limited nucleotide differences between haplotypes could reflect rapid and recent population expansion, following a population bottleneck. Otherwise, such rare mtDNA variants would be eliminated by genetic drift. The significant negative Tajima's D and Fu's Fs statistics of the pooled G. mopsus populations also reject the hypothesis of the constant population size assumed under mutation-drift equilibrium 20 . In addition, the star-shaped topology of the haplotype network is another line of evidence that suggests the sudden demographic expansion of dung beetles, as in previous studies on flies 21,22 , moths 23 , leaf beetles 24 , and marine gastropods 25 , as well as domestic cattle and goats 26,27 . The most common and presumably ancestral G. mopsus haplotype (H19) was located in the center of the haplotype network and was connected to a large number of infrequent haplotypes by only a few mutational steps. This pattern is generally interpreted to indicate a population that has recently expanded in size from a small number of founders following a population bottleneck 28 . The single bell-shaped mismatch distribution for the COI haplotypes also supported the hypothesis of population range expansion 29 , whereas the hypothesis of sudden population expansion was further supported by the abrupt increase in effective population sizes, as suggested by BSP analysis (Fig. 5). Thus, all the analyses of COI sequences from G. mopsus in Mongolia indicated sudden population expansion.
The historical demographic expansion of domesticated animals, which has resulted from archeological or anthropogenic events (e.g., domestication), has been well documented in a number of taxa 21,27,30,31 . We speculate that the recent population expansion observed in G. mopsus is related to the population dynamics (e.g., expansion) of large mammals, such as cattle, horses, and camels, that provide food and breeding ground for G. mopsus and, therefore, limit its abundance and distribution. The expansion of G. mopsus populations is estimated to have occurred ~2500 years ago, based on BSP analysis with estimated substitution rate. Such estimates depend on which molecular clock rates are applied, since different lineages are often subject to different mitochondrial substitution rates 32 . However, the hypothesis of a recent history of expansion was supported by many other analyses, as well, and the population expansion of G. mopsus was estimated to have occurred a few thousand years ago, based on the BSP analysis that accounts for the variation in our time estimation. Interestingly, our estimated timing of G. mopsus expansion seems to correspond to the advent of livestock domestication in Mongolia. Thus, the domestication of cattle and livestock was a turning point in the histories of both human and dung beetle populations 33 . The earliest domestication of wild cattle in northeastern Asia, including Mongolia, North China, Korea, and Japan, has been suggested to have occurred between 4000 and 5000 years ago 34 . The similar population structure of livestock and G. mopsus in Mongolia might further support the hypothesis of livestock-mediated expansion of G. mopsus populations. Recent population expansion in various domesticated animals, including cattle 26 and goats 27,35 , has been supported by bell-shaped mismatch distributions or star-shaped haplotype networks, much like the ones presented for G. mopsus populations in the present study. Furthermore, other population genetics studies of cattle populations have also reported unusually high genetic diversity within domesticated populations in Mongolia 35 , as well as domestication-driven population expansion 34 . Moreover, such studies have also observed that the most common and presumably most ancestral haplotypes are located at the centers of the haplotype networks and that a number of rare haplotypes originate from that haplotype with only a few mutational steps.
In a previous study of cattle dung-using Helictopleurus species (H. neoamplicollis and H. marsyas) in Madagascar, researchers used COI sequences to document a rapid population expansion that followed a resource shift 36 . The resource shift from primate droppings to cattle dung resulted in extended geographic ranges and population expansion with low genetic diversity and was associated with the historic introduction of cattle ~1500 years ago 36 . Indeed, the two Helictopleurus species exhibited considerably lower COI haplotype diversity than more wild, forest-dwelling dung beetles with more varied diets 36 . Another common cattle dung-using species, H. quadripunctatus, possesses a more intermediate diet, in that it is less specific, and, compared with forest-dwelling species, exhibits more similarities to the G. mopsus populations in Mongolia, particularly in regards to its high haplotype diversity, low nucleotide diversity, and significant IBD. In fact, like H. quadripunctatus, G. mopsus feeds on a broad range of food sources (e.g., sheep and camel dungs). Thus, a resource shift by G. mopsus from wild animal dung to domesticated cattle droppings during the domestication of cattle in Mongolia could account for the species' demographic expansion. The contemporary cattle population in Mongolia is known to have originated from hybridization with wild yak populations ~7300 years ago and has probably undergone rapid population expansion for the last 3600 years 37 .
The dispersal potential (dispersal capacity) of a species is a critical factor in characterizing its present-day population genetic structure. There is evidence that the dispersal and movement patterns of dung beetles are strongly associated with those of mammals 13,14 . Overall, the weak genetic structuring observed among the G. mopsus populations of Mongolia, as indicated by low Fst values, suggests the occurrence of either contemporary or historical gene flow. In particular, the significant relationship between the genetic and geographic distances of population pairs (IBD) indicates that the geographic proximity of populations significantly influences the spatial genetic structure of G. mopsus, probably resulting from the fact that its dispersal ability is associated with the geographic closeness.
The movement processes of dung beetles might be stochastic since individuals could move short distances between dung pats, as well as long distances between pastures, during foraging. Maximum flight distance of 1.5-2 kilometers per day has been reported in a previous study 38 . Thus, the dispersal and migration distances of this small beetle are likely limited, especially when considering the long distances between certain populations (e.g., ~670 km). However, if the migration of G. mopsus is promoted by other factors, such as livestock dispersal (i.e., herd movement during pasturing), the high level of connectivity between the G. mopsus populations becomes much more conceivable. Pasturing is a popular type of livestock farming (i.e., nomadic livestock husbandry) in Mongolia, and farmers often move their cattle herds among grasslands on a seasonal basis, following road networks. Such long-distance dispersal, either by long-distance flight or long-distance animal-mediated movement, could explain the low levels of inter-population genetic differentiation observed in the present study.
The AMOVA indicated weak, but detectable genetic structuring among certain G. mopsus population groups (i.e., North and Central, East, South area) that were clustered by the major roads that are seasonally used for moving cattle (Table 2). However, the species' complicated dispersal system (i.e., short-distance foraging movement and long-distance dispersal), which is the result of dung availability and possibly other biological factors, could lead to variability in the species' dispersal. Therefore, further investigation of G. mopsus foraging behavior and dispersal potential is needed, especially elucidate the influence of dung beetle dispersal on population connectivity in Mongolia.
In the present study, we also identified an extraordinarily high level of genetic diversity, which could indicate extremely large effective population sizes (N e ) and, possibly, high levels of historical and ongoing gene flow. Large N e values typically indicate that populations have undergone historical demographic expansion but are now stable, thereby suggesting that G. mopsus originated in Mongolia. Our haplotype network analysis, which included individuals from multiple regions, including South Korea, China, and France, also supports the hypothesis that G. mopsus originated in Mongolia. Indeed, the haplotypes from South Korea and China were connected to those from Mongolia by only one or two mutational steps, and the haplotypes from France were not much more removed. It is possible that Mongolia is the geographic origin of G. mopsus expansion and that the species radiated (or migrated) to both the East and West during the demographic processes of population expansion.
The high genetic diversity of the remnant G. mopsus populations in Mongolia might allow the species to adapt to changes in local environments, including changes in climate, vegetation, and dung availability. Mongolia includes both arid and semi-arid regions, where climate change is generally more drastic than in other regions 17 . In particular, the largest genetic differentiation was detected between populations from the northern humid areas (grassland or forest-steppe areas) and southern dessert area, which were characterized by different vegetation types, i.e., grassland or forest-steppe areas and desert, respectively. A previous study 39 reported that there was substantial year-to-year variation in the insect community of Mongolia. For example, in one year, the number of coleopteran species at a region near our ST23 locality was reduced by 60% 39 . The ST23 population, which was located at the center of a main road and was the most urbanized of our sampling localities, possessed the lowest genetic diversity of the Mongolia populations and yielded only negative Tajima's D and Fu's Fs values. This could indicate urbanization effects, which would indicate that the species is at a high risk of population size reduction in Mongolia. However, the widespread COI haplotypes and the pattern of star-shaped haplotype network may still indicate that the species is highly adaptive. Therefore, in order to understand the relationship between the genetic diversity and ecological adaptive potential of G. mopsus, future studies should investigate the response of G. mopsus to diverse environmental dynamics.
To design effective conservation and restoration management plans for G. mopsus in South Korea, understanding the genetic characteristics (e.g., genetic diversity, genetic structure, and phylogenetic or phylogeographic relationships) of potential source populations, like those of Mongolia, will be crucial 40  efforts of endangered and regionally extinct species often require the introduction of individuals from source populations, whose genetic and evolutionary lineages are comparable to those of the recipient populations. Therefore, both the population genetic structure and phylogenetic relationships of source populations should be considered carefully when developing reintroduction programs 41 . The findings of the present study highlight the value of the Mongolia populations for G. mopsus conservation and suggest that the populations could be used as source populations for restoration efforts in South Korea. Additionally, this species has been suggested to have a very large distributional range in Palearctic region 42 . Therefore, the observed phylogeographic trends by molecular data in this study need to be corroborated in a wider spectrum of its geographic distribution in future research, and it will also allow for exploring the conservation values of other remnant and potentially important populations of this species.
The genetic data presented by the present study provide valuable information about the phylogeographic structure, evolutionary history, and the genetic structure of G. mopsus populations across Mongolia. Additional ecological studies could further contribute to the effectiveness of dung beetle restoration in South Korea. The results of the present establish a basis for future research regarding the effective conservation and restoration of this highly cherished dung beetle species.  (Table 1) because the populations in those regions are now regionally extinct. Especially, the species in the western part of Europe (e.g., France) is almost extinct or in great danger of extinction compared to one in the eastern part of Europe (e.g., Greece) 43 .

Methods
Genomic DNA was isolated from a femur of each specimen using the DNeasy Blood & Tissue Kit (Qiagen, USA), and polymerase chain reaction (PCR) was used to amplify a 658 bp fragment of the mitochondrial cytochrome oxidase I gene (COI), using 20-μl reactions, universal forward (LCO1490: 5-GGT CAA CAA ATC ATA AAG ATA TTG G-3) 44  Haplotype network and genetic diversity analyses. The COI sequences were aligned using Clustal W in BioEdit v.7.0.1 46 and checked manually, and the haplotypes of the 418 individuals were determined using Neighbor-Joining algorithms in MEGA 7.0 47 . A haplotype network was then inferred using HAPSTAR v0.7 48 .
The total and population genetic diversity indices, including the number of haplotypes (HN), haplotype diversity (h), and nucleotide diversity (π), were estimated using ARLEQUIN v3.5 49 , and Tajima's D and Fu's Fs statistics were also analyzed to test for neutrality. Ten populations that had insufficient sample sizes (N < 10) were excluded from all subsequent analyses, except the AMOVA. Haplotype richness was calculated using a rarefaction method in CONTRIB v1.02 50 that corrects for unequal sample sizes among populations.
Population genetic differentiation and isolation by distance (IBD). The fixation index (Fst) 51 was calculated for population pairs using ARLEQUIN v3.5, and the geographical distance between sampling localities was measured in decimal degrees using ArcGIS Desktop 10.5 (ESRI, USA) 52 . The correlation between genetic distance (Fst) and geographic distance (decimal degrees), i.e., the Mantel test (spatial autocorrelation multivariate approach) 53 , was used to test for isolation by distance (IBD) in GenAlEx v6.502 54 .