Whole‐genome resequencing reveals recent divergence of geographic populations of the dung beetle Phelotrupes auratus with color variation

Abstract Knowledge of population divergence history is key to understanding organism diversification mechanisms. The geotrupid dung beetle Phelotrupes auratus, which inhabits montane forests and exhibits three color forms (red, green, and indigo), diverged into five local populations (west/red, south/green, south/indigo, south/red, and east/red) in the Kinki District of Honshu, Japan, based on the combined interpretation of genetic cluster and color‐form data. Here, we estimated the demographic histories of these local populations using the newly assembled draft genome sequence of P. auratus and whole‐genome resequencing data obtained from each local population. Using coalescent simulation analysis, we estimated P. auratus population divergences at ca. 3800, 2100, 600, and 200 years ago, with no substantial gene flow between diverged populations, implying the existence of persistent barriers to gene flow. Notably, the last two divergence events led to three local populations with different color forms. The initial divergence may have been affected by climatic cooling around that time, and the last three divergence events may have been associated with the increasing impact of human activities. Both climatic cooling and increasing human activity may have caused habitat fragmentation and a reduction in the numbers of large mammals supplying food (dung) for P. auratus, thereby promoting the decline, segregation, and divergence of local populations. Our research demonstrates that geographic population divergence in an insect with conspicuous differences in traits such as body color may have occurred rapidly under the influence of human activity.


| INTRODUC TI ON
Knowledge of the history of population divergence is crucial to our understanding of organisms' diversification mechanisms. The evolutionary histories of wild organisms have been shaped largely by past climatic and geographic events. For example, environmental changes that occurred during the last glacial period strongly influenced species' distribution patterns. The combined interpretation of past climatic and geographic events and molecular population genetic data has provided a powerful means to elucidate the evolutionary history of organisms over relatively long time scales such as the Quaternary (Hewitt, 2000(Hewitt, , 2004. In recent years, the impact of human activities on the evolution of organisms has gradually become clear. Human-induced environmental changes have affected the evolution of adaptive traits in many organisms (Sullivan et al., 2017) and have led to the decline and extinction of wild populations (Ceballos et al., 2015). During the early to mid-Holocene, humans began to practice agriculture and pastoralism, rapidly increasing their population size and impact on the environment (Boivin et al., 2016). It is thus important to consider the possibility that evolutionary changes in the populations of wild organisms have been influenced by not only climatic and geographic events but also human activities.
Genetic variation between local populations is determined by the balance of the homogenizing effect of gene flow and the diversifying effect of local adaptation (Endler, 1977). Geographic variation in adaptive traits often reflects genetic differences among populations due to restricted gene flow. Among many adaptive traits, coloration and pigment patterns are inferred as inter or intraspecific signals that have been used to define populations in studies of evolutionary processes (Curran et al., 2020;Hoyal Cuthill & Charleston, 2015;McLean & Stuart-Fox, 2014;Mullen et al., 2009;Stankowski et al., 2017). Recently, the development of high-throughput sequencing technology has enabled the acquisition of large amounts of fine-scale sequence variation data that provide powerful clues about the processes underlying population divergence through whole-genome resequencing (Ellegren, 2014) or reduced-representation approaches such as restriction siteassociated DNA sequencing (RAD-seq; Baird et al., 2008).
Currently, great efforts are being made to estimate past demographic changes and evolutionary events that have occurred in various wildlife populations based on large amounts of genomic data (Excoffier et al., 2013;Li & Durbin, 2011;Schiffels & Wang, 2020;Terhorst et al., 2017;Zhou et al., 2020).
Here, we focused on the geotrupid dung beetle Phelotrupes auratus (Motschulsky). It exhibits notable variation in its structural metallic body coloration, which is provided by a simple multilayer reflector structure (Akamine, Ishikawa, et al., 2011) and is categorized as red, green, or indigo based on the peak wavelength on the elytral reflectance spectrum (Akamine et al., 2008;Araki & Sota, 2021). Phelotrupes auratus is a diurnal beetle that inhabits montane forests; its food resource is the dung of large wild mammals, especially sika deer (Tsukamoto et al., 2014). In the Kinki District of Honshu, the main island of Japan, the three color forms occur separately in different geographic areas Figure 1a); however, the mechanism by which this variation formed and has been maintained remains unclear. Our previous population structure analysis using RAD markers (Araki & Sota, 2021) showed that P. auratus populations in the Kinki District diverged into five geographic groups defined by combinations of three color forms and three genetic clusters-west/red (WR), east/ red (ER), south/indigo (SI), south/green (SG), and south/red (SR)and that geographic color variation was maintained by hybrid transition zones with steep genetic clines (i.e., barriers to gene flow) or transition zones without genetic differentiation. In the latter zone type, we inferred the action of strong selection on different body colors in different regions. However, as the function of body color in P. auratus is unclear, the inference that color differentiation is maintained solely by selection despite frequent gene flow was questionable. Obscure genetic differences between sister populations may have been caused by very recent population divergence, even if they are separated by barriers to gene flow. Therefore, the lack of detection of genetic differentiation in our previous study may be attributable to the low resolution of the RAD-seq data and/ or the methodology of our population clustering analysis. Analysis using a larger genomic dataset, such as that provided by wholegenome resequencing, may provide new insights into geographic color diversification in P. auratus.
In this study, to construct robust demographic history models and estimate the dynamics of effective population sizes, gene flow, and divergence times among populations, we constructed a de novo draft genome of P. auratus and analyzed genome-wide singlenucleotide polymorphism (SNP) data obtained by whole-genome resequencing of individuals sampled from five local populations. We found that the five local populations diverged recently-within several thousand years, under the influence of human activity-and that no substantial gene flow had occurred among the diverged populations. Our findings reveal an unexpected influence of human activity on population differentiation through the divergence of a conspicuous trait in this insect species.

| Sampling and individual coloration data
Between 2019 and 2020, P. auratus adults were collected using traps baited with horse dung at six locations around Lake Biwa in the Kinki District, Honshu, Japan (Figure 1a). To obtain five or more individuals from each of the five local populations of P. auratus identified in the region, one sampling site was selected for each population based on our previous data; SG was sampled at two locations to obtain sufficient samples. A total of 27 beetles were collected ( Table 1) and fixed in 99% ethanol; one beetle collected from SI was immersed in RNAlater solution for DNA and RNA sampling and then stored at −30°C until DNA/RNA extraction. As in our previous study (Araki & Sota, 2021), we quantified the elytral colors of individuals using a spectrometer (USB2000 + UV-VIS-ES; Ocean Optics) to measure the reflection spectra of the elytral surfaces of dry specimens in a dark room and classified the color forms as indigo (<525 nm), green (525-600 nm), or red (>600 nm; Araki & Sota, 2021).

| Genomic DNA extraction and sequencing
Total genomic DNA was extracted from the flight muscle and testes using a Genomic-tip 20/G kit (QIAGEN), following the manufacturer's protocol. The quality of DNA samples was assessed with Qubit (Thermo Fischer Scientific), NanoPhotometer (Implen), and 2200 TapeStation (Agilent) instruments. For genome assembly, the genomic DNA of an SI individual (SI03) was sent to GeneBay, Inc.
(Yokohama) for library construction and long-read sequencing using the ONT PromethION system (Oxford Nanopore Technologies).
Short-read sequencing was performed with 150-bp paired-end reads in three runs on an Illumina HiSeq X Ten sequencer (Illumina) at Macrogen Japan. For genome assembly, a single sample (SI03) was sequenced using one lane of the HiSeq X Ten sequencer. For resequencing, 10 samples from the SI and ER groups and 15 samples from the WR, SG, and SR groups were separately sequenced using two lanes of the HiSeq X Ten sequencer. Library construction for SI and ER samples was performed at Macrogen Japan and that for WR, SG, and SR samples was performed at Kyoto University using the NEB Next Ultra II FS DNA Library Prep Kit for Illumina (E7805S; New England Biolabs). Adapters on raw reads were trimmed using Porechop (v0.2.4, https://github.com/rrwic k/Porechop) for ONT PromethION long reads and fastp (Chen et al., 2018) for Illumina HiSeq X short reads.

| Transcriptomic RNA extraction and sequencing
For gene prediction on the reference P. auratus genome, we obtained mRNA sequence data from an adult SI specimen (not included in Table 1) fixed in RNAlater solution. Total RNA was extracted using a Monarch Total RNA Miniprep Kit (New England Biolabs Japan), following the manufacturer's protocol. Stranded paired-end 150-bp library construction and sequencing were conducted using GENEWIZ. Adapters on raw reads were trimmed using fastp. F I G U R E 1 (a) Distribution areas and sampling sites of P. auratus in the Kinki District of Honshu, Japan. Distribution areas (encircled with different colors): red, west/red; purple, east/red; blue, south/indigo; green, south/green; and orange, south/red. For sampling sites (closed circles), see Table 1

| Read mapping and variant calling
Short reads were mapped to the draft genome with bwa-mem 0.7.17 ) and sorted with picard SortSam v2.25.7 (https:// broad insti tute.github.io/picar d/). The average read depth was 16.02, excluding SI03 (read depth, 154.11). SNP calling was conducted for each sample using bcftools 1.12 ; genotype likelihoods were computed using the mpileup command, and variants were called using the call command. Based on the average read depth of each sample, generated by samtools 1.13 , we collected SNPs with Phred quality scores >10 and depths between two-thirds (or 10 if the average depth was <15×) and twice the average depth, and removed indels. Called variants were merged into a single VCF file with bcftools merge, then filtered with maf >0.05, max missing <0.15, −m = 2, and −M = 2 to remove rare substitutions, missing abundant sites, and multiple substitution sites.

| Phylogenetic and genetic population structure analyses
A phylogenetic tree was constructed using all the obtained SNP data and the maximum-likelihood method in IQ-TREE 2.1.2 (Minh et al., 2020) with 1000 UF-bootstrap iterations (Hoang et al., 2018) and 1000 bootstrap replicates for the SH-like approximate likelihood ratio test (SH-aLRT; Guindon et al., 2010). The best-fitting substitution model was estimated using ModelFinder (Kalyaanamoorthy et al., 2017) with an ascertainment bias correction model (Lewis, 2001); as a result, PMB + F + ASC + R4 was selected. In the optimal topology, five individuals of WR were assigned to the outgroup based on the results of our previous study (Araki & Sota, 2021), although we did not assume the monophyly of WR individuals prior to the phylogenetic analysis. We also performed population genetic structural analysis using the ADMIXTURE v1.3.0 software (Alexander et al., 2009) with the number of ancestral clusters, K, ranging from 2 to 5 and principal component analyses using the pcaMethods package in the R software (Stacklies et al., 2007) based on the SNPs to assess differences in the genetic compositions of individuals among the populations.

| Demographic history
We

| Phylogenetic and population genetic structural analyses
The phylogenetic analysis revealed that individuals from WR, ER,   Figure 3). Although there was little genetic divergence among the southern populations, the genetic divergence between SR and SI + SG individuals was also shown by the PC1-PC2 plot.
These population genetic structure results were similar to our previous RAD-seq results (Araki & Sota, 2021).
We performed a coalescent simulation with fastsimcoal2 to re- and further population increases occurred in SI, SG, and SR ( Figure 5, Table S2). Although we assumed continuous gene flow between populations, the estimated effective number of migrants (N e m) was much lower than unity (10 −10 to 10 −7 ) for all cases ( Figure 5, Table S3), suggesting a lack of substantial gene flow (Slatkin, 1987).

| Demographic histories of the geographic populations
Our demographic analysis revealed the history of divergence of the five geographic populations of P. auratus. The SMC++ results   (Table S1). This result is consistent with the results of population structural analysis ( Figure 2) and our previous finding that WR diverged deeply from the other populations (Araki & Sota, 2021). However, the divergence times estimated by SMC++ were much older than those estimated by fastsimcoal2. This difference may have resulted from methodological differences between the two programs. analyses are needed to understand the differences in divergence times estimated using different approaches.

F I G U R E 4
Historical changes in effective population sizes, inferred by SMC++. Line colors represent local populations: red, west/red; purple, east/red; blue, south/indigo; green, south/green; and orange, south/red. Years ago Effective population size F I G U R E 5 Demographic history of P. auratus populations inferred by coalescent simulation using fastsimcoal2 during the past 5000 years. Horizontal bar widths indicate effective population sizes, N e (median and 95% confidence intervals are indicated by numerals); black circles and red vertical bars indicate medians and 95% confidence intervals, respectively, of divergence times.
The coalescent simulation conducted with fastsimcoal2 and the incorporation of gene flow in the model revealed that no substantial gene flow among the five populations had occurred since their divergence ( Figure 5), suggesting that the geographic color variation of P.
auratus is maintained mainly by barriers to gene flow. In our previous study, we inferred the occurrence of frequent gene flow among the three southern populations (SR, SG, and SI) and the presence of hybrid zones between WR and SG and between ER and SR (Araki & Sota, 2021), and the genetic similarity among southern populations was also supported by the population structural analysis conducted in this study. However, the fastsimcoal2 analysis revealed that the divergence of the southern populations occurred within the past 600 years, implying that an insufficient number of unique genetic variations had accumulated in each population to detect divergence between local populations through population clustering analysis, despite the establishment of geographic isolation.

| Factors restricting gene flow
Although P. auratus has dispersal potential by means of flight, they are confined to mountain forest areas and do not occur in wide open lands between lowland forests. The distribution of P. auratus among separate mountain areas in the Kinki District reflects this habitat selection pattern (Figure 1a). Forest fragmentation can also lead to a decline in the supply of food resources for P. auratus, as forests are also the main habitat of large mammals in the Japanese islands.
Therefore, forest fragmentation would have affected population divergence in this species.
The confinement of P. auratus to forest areas may also be related to the high risk of predation during flight outside of forests. The metallic coloration of beetles may be effective for predator avoidance when they fly among light and dark forest patches (Schultz, 1986).
However, when beetles fly out into open land, they are easily tracked by birds. Thus, migration between patchy forests is likely to be difficult. The cost of subsocial behavior on parents with larval nests may also restrict the time allowed for dispersal by flight. Large dung beetles that construct deep nests and lay few, large eggs are expected to have low dispersal rates (Hanski & Cambefort, 2014).
This trend is likely to occur in P. auratus, as females build nests at depths of 50-100 cm and lay only one large egg in each nest (Tateno et al., 2020); this relatively sedentary lifestyle may limit opportunities for dispersal and migration.

| Factors influencing population divergence
Our fastsimcoal2 analysis showed that P. auratus populations in the Kinki District diverged about 3800 years ago and more recently ( Figure 5). Although the Holocene climate has been relatively stable and warm (Petit et al., 1999), minor local climatic changes (± 2°C) have occurred several times (Seppä et al., 2009;Wang et al., 2011).
In Eurasia, marked climate cooling occurred 4200-4000 years ago, which may have caused the decline of several human civilizations, including the Jomon population in Japan (Yasuda et al., 2004). Japan also experienced climatic cooling 3800-3600 years ago (Kawahata et al., 2009;Yasuda et al., 2004). Declines in mammal populations due to hunting may have led to declines in the species richness and abundance of dung beetles (Andresen & Laurance, 2007). Phelotrupes auratus is thought to depend largely on sika deer dung, due to its overlapping distribution with sika deer (Toda & Akei, 2003). In addition, P. auratus likely has limited dispersal propensity due to its reproductive ecology. Thus, reduction and fragmentation of the sika deer population may have strongly affected the P. auratus population. Although P. auratus can also use dung produced by other large mammals as food resources, archeological records for the Jomon period (1200-2400 years ago) and documents mentioning local products from the mid-Edo era (1730 s) have shown that sika deer and other large mammals such as wild boars had reduced their distribution areas by the 18th century (Tsujino et al., 2010). A genetic structural analysis also suggested that the sika deer population in Japan has experienced habitat fragmentation and bottlenecks due to human activity since the mid-19th century (Goodman et al., 2001). To test our hypothesis that human activity has impacted large mammal and dung beetle populations, it is necessary to study the demographic history of large mammals that provide dung food for P. auratus with sufficient temporal and spatial resolution.

| Divergence of color morphs among local populations
Although we detected marked divergence of local populations with different coloration within short periods in this study, it remains uncertain whether different structural coloration has adaptive importance (Araki & Sota, 2021). In a previous study, a Müllerian mimicry hypothesis was proposed based on the similarity of geographic color variation patterns between P. auratus and a sympatric congener, Phelotrupes laevistriatus (Watanabe et al., 2002). However, carcasses of the these beetles are often found in the feces of mammals such as raccoon dogs (Matsuyama et al., 2006) and chickens have been observed to feed on them without any problems (Kochi, 2001), implying that the Phelotrupes beetles are not unpalatable. As mentioned previously, the metallic coloration of P. auratus may be effective for avoiding bird predation in forests, where they fly among light and dark patches (Schultz, 1986). However, color differences may not be important in this predator avoidance strategy. If differences in the structural coloration of P. auratus are not important, then color differences may have been fixed by chance. A reduction in the effective population size of P. auratus was occurring when the indigo and green forms diverged ( Figure 5), which may indicate that this color divergence was driven by genetic drift.

| CON CLUS ION
In this study, we conducted a detailed genomic analysis of genetic information from five local populations of P. auratus. Our analyses revealed that geographic color variation in P. auratus is of very recent origin and is maintained by restricted gene flow. The recent divergence of these populations may have been influenced by increased human activity in the past few thousand years, providing an interesting example of the indirect impact of human activity on trait divergence in wild insect populations. Our study implies the need for additional research conducted with a population genomics approach to reveal the historical relationships between the trophic levels, namely large mammals and dung beetles. Future studies of the demographic history of sika deer and other large mammals using genomic data will be essential to our understanding of the process of recent P. auratus population differentiation.

ACK N OWLED G M ENTS
We thank S. Yamamoto, T. Fujisawa, and S. Nomura for their helpful advice. This study was partly supported by JSPS KAKENHI (nos. 15H02637 and 18H04010).

CO N FLI C T O F I NTE R E S T
The authors declare that there is no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw sequence reads are available from the DDBJ Read Archive Teiji Sota https://orcid.org/0000-0002-0400-9147