Beyond geography: climatic gradients shape Reeves’s muntjac population structure in Taiwan

Beyond geography: climatic gradients shape Reeves’s muntjac population structure in Taiwan | Authorea try { document.documentElement.classList.add('js'); } catch (e) { } var _gaq = _gaq || []; _gaq.push(['_setAccount', 'G-8VDV14Y67G']); _gaq.push(['_trackPageview']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })(); Skip to main content Preprints Collections Wiley Open Research IET Open Research Ecological Society of Japan All Collections About About Authorea FAQs Contact Us Quick Search anywhere Search for preprint articles, keywords, etc. Search Search ADVANCED SEARCH SCROLL Ecology and Evolution This is a preprint and has not been peer reviewed. Data may be preliminary. 27 October 2025 V1 Latest version Share on Beyond geography: climatic gradients shape Reeves’s muntjac population structure in Taiwan Authors : YiLun Peng 0009-0003-1999-9619 , Hsuan-Wien Chen , Shih-Wei Chang , Chun-Yi Hsiao , and Hurng-Yi Wang 0000-0003-1708-8734 [email protected] Authors Info & Affiliations https://doi.org/10.22541/au.176160808.85226808/v1 Published Ecology and Evolution Version of record Peer review timeline 214 views 169 downloads Contents Abstract Information & Authors Metrics & Citations View Options References Figures Tables Media Share Abstract Studies of Taiwanese mammals reveal varied patterns of population divergence. Low-mobility species such as mole-shrews and Formosan wood mice exhibit strong north–south splits. Surprisingly, similar patterns also occur in mobile taxa like Formosan serow and sambar deer. In contrast, other mobile species, including flying squirrels, and Reeves’s muntjac, show weak or no population structure. This recurring north–south divergence across taxa with diverse ecologies suggests that shared environmental and geographic gradients, in addition to historical isolation, may underlie these patterns. If so, species occupying similar habitats may exhibit comparable genetic breaks regardless of life history traits. Prior studies based solely on mitochondrial DNA may have missed fine-scale structure in species like muntjac; high-resolution SNP data now offer improved resolution. Here, we analyzed genome-wide SNPs from 71 Taiwanese Reeves’s muntjac and comparative Chinese samples. We detected deep divergence from Chinese muntjac (~0.24 MYA), and further north–south subdivision within Taiwan (~0.06 MYA). Demographic modeling revealed a complex history involving glacial isolation and asymmetric gene flow, mainly from north to south. Within Taiwan, genetic differentiation was shaped by both geography and climate, especially temperature annual range (Bio7), with niche models showing environmental separation. Selection scans identified PLA2-associated genes, potentially linked to thermal adaptation. This is the first study to demonstrate that both geographic and environmental heterogeneity jointly contribute to mammalian divergence in Taiwan. The repeated north–south split across ecologically diverse species highlights shared climatic and topographic factors driving parallel population structure in Taiwan’s montane ecosystems. Beyond geography: climatic gradients shape Reeves’s muntjac population structure in Taiwan Yi-Lun Peng 1# , Hsuan-Wien Chen 2# , Chun-Yi Hsiao 3 , Shih-Wei Chang 4 , Hurng-Yi Wang 1,5,6,7* 1. Institute of Ecology and Evolutionary Biology, National Taiwan University, Taipei, Taiwan 2. Department of Biological Resources, National Chiayi University, Chiayi, Taiwan 3. Technology Commons, College of Life Science, National Taiwan University, Taipei, Taiwan 4.Taiwan Biodiversity Research Institute, Nantou, Taiwan 5. Graduate Institute of Clinical Medicine, College of Medicine, National Taiwan University, Taipei, Taiwan 6. Graduate Institute of Medical Genomics and Proteomics, College of Medicine, National Taiwan University, Taipei, Taiwan 7. Department of Entomology, National Taiwan University, Taipei, Taiwan # These authors contribute equally to this work * Corresponding authors. Hurng-Yi Wang, [email protected] ABSTRACT Studies of Taiwanese mammals reveal varied patterns of population divergence. Low-mobility species such as mole-shrews and Formosan wood mice exhibit strong north–south splits. Surprisingly, similar patterns also occur in mobile taxa like Formosan serow and sambar deer. In contrast, other mobile species, including flying squirrels, and Reeves’s muntjac, show weak or no population structure. This recurring north–south divergence across taxa with diverse ecologies suggests that shared environmental and geographic gradients, in addition to historical isolation, may underlie these patterns. If so, species occupying similar habitats may exhibit comparable genetic breaks regardless of life history traits. Prior studies based solely on mitochondrial DNA may have missed fine-scale structure in species like muntjac; high-resolution SNP data now offer improved resolution. Here, we analyzed genome-wide SNPs from 71 Taiwanese Reeves’s muntjac and comparative Chinese samples. We detected deep divergence from Chinese muntjac (~0.24 MYA), and further north–south subdivision within Taiwan (~0.06 MYA). Demographic modeling revealed a complex history involving glacial isolation and asymmetric gene flow, mainly from north to south. Within Taiwan, genetic differentiation was shaped by both geography and climate, especially temperature annual range (Bio7), with niche models showing environmental separation. Selection scans identified PLA2-associated genes, potentially linked to thermal adaptation. This is the first study to demonstrate that both geographic and environmental heterogeneity jointly contribute to mammalian divergence in Taiwan. The repeated north–south split across ecologically diverse species highlights shared climatic and topographic factors driving parallel population structure in Taiwan’s montane ecosystems. Keywords: Population structure, Climatic gradients, Genome-wide SNPs, Taiwan phylogeography, Adaptive divergence INTRODUCTION Understanding genetic diversity and differentiation within species is crucial for biodiversity conservation, as it underpins adaptability and ecosystem resilience. Analyzing genetic patterns helps identify hotspots of variation, informs targeted conservation strategies, and supports habitat restoration. These insights are particularly important for developing climate-resilient conservation plans (Hughes, et al. 2008; Allendorf, et al. 2010; Hoban, et al. 2021; Shaw, et al. 2025). Genetic differentiation in mammals is shaped by historical processes (e.g., glacial cycles), environmental variation (e.g., habitat and climate), and species-specific traits such as dispersal and demography (Clobert 2012; Wang and Bradburd 2014; Theodoridis, et al. 2020; Afonso Silva, et al. 2025). In this context, Taiwan provides a valuable natural laboratory, where diverse ecological settings and complex topography have given rise to a wide range of phylogeographic patterns across mammalian species. Genetic studies of Taiwanese mammals have revealed divergent patterns of population structure that appear linked to dispersal ability. Species with limited dispersal, such as the mole-shrew ( Anourosorex yamashinai ) and Formosan wood mouse ( Apodemus semotus ), show distinct northern and southern lineages (Hsu, et al. 2001; Yuan, et al. 2006), with the boundary broadly corresponding to the Wu (Dadu) and Zhoushui Rivers. Interestingly, some large mammals with no apparent dispersal restriction like Formosan serow ( Capricornis swinhoei ) and sambar deer ( Rusa unicolor swinhoii ) also exhibit similar north-south divergence (Horng, et al. 2003; Li, et al. 2023). In contrast, species with strong dispersal ability, including the white-bellied rat ( Niviventer culturatus ), Formosan Reeves’s muntjac ( Muntiacus reevesi micrurus ), and flying squirrels ( Petaurista lena and P. grandis ), exhibit weak or no population structure (Hsu, et al. 2000; Chan 2004; Oshida, et al. 2011; Tan 2017) . Distinct mitochondrial DNA (mtDNA) lineages observed in various taxa are proposed to have originated from Pleistocene glacial refugia (Oshida, et al. 2006; Yuan, et al. 2006; Wang 2012). Although glacial refugia might have promoted genetic divergence, subsequent post-glacial dispersal likely blurred their geographic patterns. For example, the Formosan serow, adapted to steep terrain, may experience fewer topographic barriers and retain multiple mitochondrial lineages within a single population. Consequently, similar north–south phylogeographic patterns across species with varying dispersal abilities suggest that historical isolation alone cannot fully account for the observed differentiation and merit further investigation. In addition to historical processes, contemporary geographic and environmental factors may also contribute to phylogeographic structure. If this is the case, species with similar distributions should exhibit comparable north-south splits, regardless of their historical background. Prior studies may have overlooked fine-scale patterns in species such as Reeves’s muntjac and flying squirrels, due to reliance on mtDNA, which reflects only maternal lineages and may miss subtle population differentiation (Hurst and Jiggins 2005). By using genome-wide SNP data, we can resolve phylogeographic patterns and identify the ecological and genetic factors shaping intraspecific divergence across Taiwan’s montane ecosystems. To further investigate the drivers of population divergence in Taiwan’s montane fauna, we analyzed the genomic variation of the Taiwanese Reeves’s muntjac using double digest Restriction-site Associated DNA sequencing (ddRAD-seq). Our results revealed a clear genetic division between northern and southern populations, with the boundary broadly corresponding to the Wu (Dadu) and Zhoushui Rivers. This north–south split is associated with environmental variables, particularly temperature annual range, implying that local climatic conditions, in addition to geographic barriers, contribute to population structure. We also noted that this genetic break mirrors those observed in other mammals, such as the Formosan serow and sambar deer, indicating that shared climatic and geographic gradients may drive parallel divergence across species in Taiwan’s mountainous landscapes. MATERIALS and METHODS Sample Collection, DNA Extraction and Sequencing A total of 71 tissue samples of Taiwanese Reeves’s muntjac ( Muntiacus reevesi micrurus ), collected across Taiwan, were provided by the Taiwan Biodiversity Research Institute (Fig. 1 and Fig S1; Table S1). Genomic DNA was extracted from muscle tissue using the FavorPrep DNA Extraction Mini Kit (Favorgen Biotech, Taiwan). ddRAD libraries were prepared using SbfI and MseI at Technology Commons, College of Life Science, National Taiwan University and sequenced on an Illumina NovaSeq 6000 (150 bp paired-end; Genomics BioSci & Tech, Taiwan). The mitochondrial cytochrome b gene was amplified using species-specific primers (F: 5′-ACCACGACTAATGATATGAAAAACC-3′; R: 5′-TGTCCTCCTTTTCTGGTTTACAA-3′) and sequenced following PCR using the KAPA Taq Kit under standard cycling conditions. Data Processing Raw reads were processed using Stacks v2.62 (Rochette, et al. 2019) with quality filtering via process_radtags (–score-limit 20). Filtered reads were aligned to the M. reevesi reference genome (GCA_008787405.2) using BWA-MEM , and loci were assembled with gstacks, retaining paired-end reads with mapping quality ≥20. SNPs were called with populations, retaining loci present in ≥80% of individuals, and exported in VCF format (complete SNP dataset). Whole-genome sequences of Chinese M. r. reevesi were downloaded (Chen, et al. 2024) and mapped to the reference genome (GCA_008787405.2) using GPU-accelerated FQ2BAM workflow in NVIDIA Parabricks (NVIDIA 2019). Variants were called with bcftools (-q 20 and -Q 20) (Danecek, et al. 2021), and sites with depth ≥5 and <mean+ two standard errors were retained. Regions homologous to the ddRADseq loci of Taiwanese muntjac, including both variant and invariant sites, were extracted to ensure consistency in downstream analyses. Mitochondrial genomes were assembled from whole-genome data using NOVOPlasty v4.3.5 (Dierckxsens, et al. 2017), aligned with MARS (Ayad and Pissis 2017), and the cytochrome b gene was extracted for phylogenetic comparison with Taiwanese muntjac. Population Structure Analysis To minimize linkage bias, SNPs with minor allele count 0.2 within 25 kb windows were removed using PLINK v1.90 (Purcell, et al. 2007), yielding a pruned dataset for population structure analyses. Population structure was inferred using the LD-pruned SNP set. ADMIXTURE v1.3.0 (Alexander, et al. 2009) was run across K = 1-10, with cross-validation to identify the optimal K. Principal component analysis (PCA) was performed using adegenet v2.1.10 (Jombart 2008), which depends on ade4 v1.7-22 (Dray and Dufour 2007), and visualization of the results using ggplot2. Pairwise F ST (Weir and Cockerham 1984) and nucleotide diversity (π) were calculated in VCFtools. Phylogenetic reconstruction and divergence time estimation Phylogenetic relationships were reconstructed using both nuclear and mitochondrial datasets. A maximum likelihood tree based on genome-wide SNPs was inferred using IQ-TREE2 v2.2.2.6 (Minh, et al. 2020), with the best-fit model selected by ModelFinder (Kalyaanamoorthy, et al. 2017) and node support assessed via 10,000 ultrafast bootstrap replicates. Mitochondrial phylogeny was inferred from 943 bp of cytochrome b sequences from Taiwanese and Chinese muntjacs, supplemented with public data. Sequences were aligned using MUSCLE (Edgar 2004) in MEGA X (Kumar, et al. 2018), and a neighbor-joining tree was constructed under the Kimura 2-parameter model. Divergence times were estimated using Bayesian inference of species trees implemented in SNAPPER (Stoltz, et al. 2021). To improve computational efficiency, 2,000 SNPs and 10 individuals per population were subsampled following the recommendations of Chen et al. (2024). Analyses were run in BEAST2 (Bouckaert, et al. 2019) for 1,000,000 MCMC iterations, discarding the first 100,000 as burn-in. Divergence times were calibrated using the median estimate for western Chinese populations (0.065 Mya; 95% CI: 0.003–0.154 Mya) and reported as medians with 95% highest posterior density (HPD) intervals. Demographic Inference of Taiwanese Reeves’s muntjac To infer the demographic history of Taiwanese muntjacs, site frequency spectrum (SFS) were generated from the complete dataset using easySFS (https://github.com/isaacovercast/easySFS), assuming a mutation rate of 2.38 × 10⁻⁹/site/generation (Yin, et al. 2021) and a generation time of 2.5 years (Pacifici, et al. 2013). Stairway Plot (Liu and Fu 2020) was used to infer independent demographic histories for northern and southern populations. To explore more complex demographic scenarios involving population size changes and gene flow, 12 models were evaluated using fastsimcoal2 v2.7.0.9 (Excoffier, et al. 2021) (Fig S3). Divergence times were constrained between 20,000 and 100,000 years, based on SNAPPER estimates. Each model was run 100 times, with 40 expectation–conditional maximization (ECM) cycles and 100,000 simulations per cycle. The best-fitting model was selected based on maximum likelihood using fsc-selectbestrun.sh, and further evaluated using AIC via calculateAIC.sh (https://github.com/speciationgenomics/scripts). Given the limited resolution of SFS for recent events, we used GONE (Santiago, et al. 2020), a linkage disequilibrium-based method, to infer changes in effective population size over the past 100 generations, assuming a recombination rate of 1 cM/Mb consistent with mammalian averages (Dumont and Payseur 2008). Landscape and Climatic Correlates of Genetic Differentiation To evaluate the drivers of genetic differentiation in Taiwanese Reeves’s muntjacs, we examined the effects of geographic distance (straight-line distance), landscape resistance (least-cost path distance), and climatic variation (environmental dissimilarity) using 46 individuals with precise coordinates (one per locality)(Table S1). Genetic distances (Euclidean) were calculated using the gl.dist.ind in dartR v2.9.7 (Gruber, et al. 2018). Geographic distances were computed as straight-line distances using distm from geosphere v1.5-18 (https://github.com/rspatial/geosphere). Least-cost path distances were generated from a resistance surface combining slope and land-use data. Slope layers were generated in QGIS (QGIS Development Team 2009) with resistance values increasing with gradient (0–10º = 1, 20–40º = 1.5, 50–60º = 4, 70–80º = 6, >80º = 10). Land-use categories (Chen, et al. 2019) were similarly assigned values reflecting presumed resistance to muntjac movement: Forest = 1, Grassland = 2, Bare land = 2, Agricultural land = 3, Inland water = 5, Built-up land = 10 (Wang 1989; McCullough, et al. 2000). Final distances were computed with gl.costdistances in dartR . Climatic data were obtained from WorldClim v2, which includes 19 bioclimatic variables representing current climate conditions (Fick and Hijmans 2017). Pearson’s correlation coefficients among the 19 bioclimatic variables were calculated using raster.cor.matrix in ENMTools v1.1.2 (Dinnage 2024). For variable pairs with |r| > 0.7, only the variable more strongly correlated with genetic distance was retained. Selected variables were used to calculate environmental dissimilarity with vegdist in vegan (Oksanen J 2024). Associations between genetic distance and geographic, resistance, and climatic distances were tested using Mantel and partial Mantel tests, with sampling year as a covariate to control for temporal effects. Environmental Drivers of Genetic Structure and Niche Differentiation Climatic variables significantly associated with genetic distance were used to quantify niche differentiation between northern and southern populations. PCA was conducted using dudi.pca in ade4 , and population niches were modeled with ecospat.grid.clim.dyn in ecospat v4.0.0 (Broennimann O 2023). Niche equivalency and similarity (ecospat.niche.equivalency.test, ecospat.niche.similarity.test) were tested with 999 permutations to assess whether the two populations occupy distinct climatic spaces. To evaluate environmental influences on genetic structure, Generalized Dissimilarity Modeling (GDM; gdm v1.5.0-9.1) (Fitzpatrick M 2022) was used to estimate the contributions of geographic distance, climate, and landscape resistance to genetic differentiation. For distribution modeling, we used Maxent v3.4.4 (Phillips and Dudík 2008) with GPS data from sampled individuals and GBIF records (2024). After removing duplicates, oceanic/erroneous points, and records with retained. We used QGIS v3.28.2 (QGIS Development Team 2009) to set spatial boundaries between 21.708°–26.5°N and 118°–123°E. We excluded highly correlated variables and retained six Bioclim layers, bio1, bio2, bio7, bio13, bio14, and bio15, for modeling (Merow et al., 2013). We also used Bioclim layers from the Last Glacial Maximum (LGM, the species’ potential distribution during glacial-period. Identification of putative loci under selection We used two complementary approaches to identify loci putatively under selection. PCadapt v4.3.5 (Privé, et al. 2020) identified SNPs strongly associated with principal components of genetic variation (FDR < 0.01), and BayeScan v2. 1 (Foll and Gaggiotti 2008) detected outliers based on allele frequency differences under a Bayesian framework, using 20 pilot runs x 5,000 iterations, burn-in = 50,000, and prior odds = 100. Genes within ±150 kb of each candidate SNP were selected as candidates. From the reference genome (GCA_008787405.2), 23,112 protein-coding genes were curated using AGAT v1.4.1 (Jacques Dainat 2024), as the background for enrichment analysis. Functional annotation was performed using EGGNOG-mapper v2.1.12 (Cantalapiedra, et al. 2021). Gene Ontology (GO) enrichment was conducted using clusterProfiler v4.10.1(Wu, et al. 2021), using an adjusted q < 0.05. Candidate gene functions were inferred by identifying Mus musculus orthologs via Ensembl (Harrison, et al. 2024), and phenotype associations explored using modPhEA(Weng and Liao 2017). Population Structure We sequenced 71 Taiwanese Reeves’s muntjac ( Muntiacus reevesi micrurus ) individuals from across Taiwan (Fig. 1, Fig. S1-S2, Table S1), generating approximately 61.9 Gb of sequencing data. The average number of bases per individual was 756 Mb, with a standard deviation of 123 Mb (Table S2). Using bwa mem and gstacks, we obtained 524,360 loci, with an average sequencing depth of 60.3x per sample (standard deviation: 9.6x, with maximum and minimum depths of 89.5x and 35.3x, respectively) (Table S3). After filtering with populations, 37,878 loci were retained, containing a total of 94,034 SNPs (complete dataset) (Table S4). To further investigate the evolutionary history of Taiwanese Reeves’s muntjac, 60 whole-genome sequences of M. reevesi reevesi from China were downloaded (Fig. S2), and genomic regions overlapping with our ddRAD-seq loci were retained and combined with the Taiwanese ddRAD-seq data for downstream analyses (see Materials and Methods). Compared to other deer species, Reeves’s muntjac exhibits lower nucleotide diversity (π) (Table S5) (Hu, et al. 2019; de Jong, et al. 2021; Combe, et al. 2022). In addition, the nucleotide diversity of Taiwanese Reeves’s muntjac (0.26%) is approximately 40% less than that of Reeves’s muntjac from China (0.42 - 0.49%). The nucleotide diversity of Reeves’s muntjac from China, estimated from whole-genome data, closely matches that of the subset overlapping with our ddRAD-seq loci (Table S5), suggesting that these loci are representative of the genome-wide diversity. As expected, both admixture and principal component analysis (PCA) revealed that Reeves’s muntjac can be divided into two major genetic lineages, Taiwan and China (Fig. 1). When K = 3, the muntjacs in China are further divided into eastern and western groups. As K increased to 4, the muntjacs collected from northern and southern Taiwan formed two distinct genetic clusters, separated by the Wu (Dadu) River and Zhoushui River, with admixed individuals observed in the intermediate region (Fig. S1). Phylogenetic analysis based on nuclear SNPs revealed patterns consistent with the Admixture results (Fig. S4). In contrast, the mitochondrial cytochrome b phylogeny showed that individuals from the same sampling locations often fell into different clades, indicating weak geographic structuring of mtDNA lineages. The F ST based on nuclear SNPs between Taiwan and China ranges from 0.069 to 0.227 (Table 1). The Taiwanese Reeves’s muntjac is genetically closest to the easternmost population of China’s muntjac. Within Taiwan, where samples were collected from locations less than 300 km apart, the F ST between the northern and southern populations is 0.036. This level of differentiation is comparable to that observed among muntjac populations in China, where F ST values range from 0.017-0.056 with an average of 0.035, despite these samples being collected from locations over 1,000 km apart. Interestingly, the F ST derived from mtDNA show similar pattern as those from nuclear SNPs (Table 1B), indicating that mitochondrial genome also retains signals of differentiation. The estimated divergence time between China and Taiwan was 0.24 million years ago (MYA) with a 95% confidence interval (CI) of 0.08 – 0.54 MYA. Within Taiwan, the divergence between northern and southern populations occurred 25,000 years ago (95% CI: 8,200 - 56,000) (Fig. 1D). Demography and divergent history According to the stairway plot analysis, both northern and southern Taiwanese Reeves’s muntjac populations experienced a population reduction beginning around 0.25 MYA (Fig 2A), which is consistent with their divergence from China’s muntjac. Subsequently, population sizes remained relatively stable until a second decline occurred approximately 100 - 150 years ago, likely associated with recent human population expansion and increased anthropogenic pressures in Taiwan (Fig. 2B). To further investigate the history of population differentiation in Taiwanese Reeves’s muntjac, we applied a coalescent-based simulation program fastsimcoal2 to estimate demographic parameters under 12 evolutionary scenarios (Fig. S3 and Table S6). The most supported model incorporates population size change as well as both early and recent gene flow between the northern and southern populations (Fig 3). According to this model, the ancestor population of Taiwanese Reeves’s muntjac was approximately 534,000 and the population divergence occurred 66,000 years ago. After divergence, the effective population size of the southern population declined to 28,407 , while that of the northern population was 57,679. Subsequently, the southern population expanded to 190,000 around 41,000 years ago, while the northern population grew to 81,000 around 32,000 years ago. The model also indicates that there was gene flow between populations after divergence. Gene flow ceased from approximately 45,860 to 31,660 years ago. Afterwards, gene flow was reestablished. However, both instances of gene flow were asymmetric, with a higher migration rate from the northern population to the southern populations. Factors that cause population differentiation Mantel test was used to assess the correlation between genetic and geographic distance, as well as the relationship between genetic distance and environmental factors among individuals. When all samples were analyzed together, a significant correlation between genetic and geographic distances was observed (r = 0.51, p < 10⁻⁴; Fig. 4). However, when analyzed separately, the northern population showed a strong correlation (r = 0.7, p < 10⁻⁴), while no such correlation was found in the southern population. Since the sampling of Taiwanese muntjac in this study spanned approximately 25 years , Mantel and partial Mantel tests were conducted to assess whether sampling time influenced genetic distance. The results indicated no significant relationship between genetic distance and sampling time ( r = 0.01, p = 0.42). Mantel test was also used to assess the correlation between genetic distance and environmental factors. The results showed that eight climatic variables, Bio2, Bio3, Bio4, Bio7, Bio14, Bio15, Bio17, and Bio19, were significantly correlated with genetic distance of Taiwanese muntjac populations (Table S7A). However, when geographic distance or least-cost path was accounted for, the effects of these environmental variables were no longer significant. This suggests substantial overlap between environmental variables and geographic or least-cost path distances. When the northern and southern populations were analyzed separately, the northern population showed a significant correlation between genetic distance and Bio3, Bio4, Bio14, Bio17, and Bio19, but these correlations were no longer significant after accounting for geographic distance or least-cost path (Table S7B). In contrast, the southern population showed significant correlations with mean diurnal range (Bio2), temperature seasonality (Bio4), temperature annual range (Bio7), and precipitation of driest month (Bio14), and these associations remained significant even after controlling for geographic distance or least-cost path (Table S7C). Generalized Dissimilarity Modelling (GDM) was used to evaluate the relative contribution of different environmental and geographic factors on genetic divergence. When the analysis included both northern and southern populations, the best model is least-cost distance which explains 44.75% of genetic differentiation (Table S9A; Fig. S5). For the northern populations, the best model considers only geographic distance which explains 59.84% of genetic differentiation (Table S9B; Fig. S6). For the southern populations, the best model contains only temperature annual range (Bio7) which explains 62% of genetic variations (Table S9C; Fig. 5). To model the distribution of the Taiwanese Reeves’s muntjac in Maxent, we first removed highly correlated variables (Table S8), retaining only those most associated with genetic distance (Table S10). The final model included four predictors, annual mean temperature (Bio1), temperature annual range (Bio7), precipitation of driest month (Bio14), and precipitation seasonality (Bio15), which contributed 3.2%, 62.6%, 13.6%, and 20.7%, respectively. The model showed strong predictive accuracy (AUC = 0.904; Fig. S8). The niche of the Taiwanese Reeves’s muntjac during the Last Glacial Maximum (LGM) appears to have contracted toward the southeastern region of the island, particularly around the southern tip (Fig. 6). To assess niche differentiation between Taiwanese muntjac populations, we conducted PCA-env using three uncorrelated bioclimatic variables (Bio2, Bio7, Bio15). The first two axes explained 96.15% of climatic variation (PC1 = 57.63%, PC2 = 38.52%) (Fig. 4A). Southern individuals occupied a narrower climatic space, while northern individuals spanned broader conditions. The niche equivalency test (D = 0.13, p = 0.05; Fig. 4B) indicated statistically distinct niches, but the niche similarity test was not significant (p = 0.86; Fig. 4C), suggesting the difference may not exceed background environmental variation. These results imply that divergence is more likely shaped by geographic separation and environmental gradients than by adaptive niche differentiation. To further test whether different Taiwanese muntjac populations exhibit niche differentiation, we used PCA to summarize variation across different environmental variables (PCA-env). Among eight bioclim variables that show correlation with genetic distances, we select mean temperature diurnal range (Bio2), temperature annual range (Bio7), and precipitation seasonality (coefficient of variation) (Bio15) for PCA-env analysis, as the pairwise correlation between these variables was ≤ explained 96.15% of the total variation along climatic gradients (PC1 = 57.63%, PC2 = 38.52%) (Fig. S7A). The southern population occupies a relatively narrow range in both dimensions of the climatic space, whereas the northern population spans a broader climatic niche along both PC1 and PC2. Based on the niche equivalency test, the D value between the northern and southern Taiwanese populations was 0.13 (p = 0.05) (Fig. S7B), indicating that the two populations occupied statistically non-equivalent climatic spaces. However, the niche similarity test was not significant (p = 0.86) (Fig. S7C), suggesting that the observed differences in niche occupancy are not greater than expected by chance. Together, these results suggest that the divergence between northern and southern populations is more likely driven by geographic separation and environmental gradients than by adaptive niche differentiation. Putative regions under selection Using PCadapt and BayeScan, 90 and 14 SNPs potentially under positive selection were identified, respectively, with 10 SNPs detected by both methods (Table S11). Within a 150 kb region surrounding these SNPs, 154 protein-coding genes were identified. Gene enrichment analysis revealed a significant Gene Ontology (GO) term, GO:0004623, which is associated with phospholipase A2 (PLA2) activity (Table S12). PLA2 is involved in various biological functions, including signal transduction, lipid mediator production, and membrane remodeling, and plays a critical role in responding to environmental stressors such as temperature fluctuations and oxidative stress (Murakami, et al. 2017). Among these, responding to temperature fluctuations is the most interesting, as the most prominent environmental factor associated with Taiwanese muntjac divergence is temperature annual range (Bio7), which is the indicator of seasonal temperature variation. The mouse orthologs of these genes are Pafah1b3, Pla2g4d, Pla2g4e, and Pla2g4f. According to phenotype enrichment analysis (modPhEA)(Weng and Liao 2017), these genes were involved in ion homeostasis (Pla2g4f) and physiological strength and integument phenotype (Pla2g4e) (Table S13). DISCUSSION We sequenced 71 individuals of Taiwanese Reeves’s muntjac ( Muntiacus reevesi micrurus ) using ddRAD-seq and incorporated 60 whole-genome sequences of M. reevesi from China to facilitate comparative analyses. Homologous genomic regions, encompassing both variant and invariant sites, were extracted from the whole-genome dataset based on ddRAD-seq loci to ensure consistency, and nucleotide diversity comparisons confirmed that the ddRAD-seq data reliably represents genome-wide diversity. Divergent history of muntjac The divergent history of Reeves’s muntjac in East Asia coincides with past climatic fluctuations and geological events. Bayesian inference indicates that the divergence between Taiwanese and Chinese muntjac populations occurred 0.24–0.25 MYA, at the onset of Marine Isotope Stage 7 (MIS 7; ~0.24–0.19 MYA), likely driven by sea-level rise during this warm interglacial period (Sato, et al. 2017). This timing is consistent with demographic analyses showing a sharp population decline in Taiwanese muntjac between 0.3 and 0.25 million years ago, likely due to a founder effect associated with the colonizing lineage. Furthermore, the nucleotide diversity of Taiwanese muntjac is only half that of their Chinese counterparts, supporting a scenario of reduced genetic diversity following colonization. The divergence between the northern and southern Taiwanese muntjac populations was estimated at approximately 66,000 years ago, during MIS 4 (∼71–57 ka), a period characterized by cold and dry climatic conditions in East Asia (Doughty, et al. 2021). This initial split was followed by gene flow between the two populations. We observed that this divergence estimate contrasts with that inferred from Bayesian inference (Fig. 1D), which does not account for post-divergence migration and instead estimates a more recent split at ~25,000 years ago. Interestingly, when migration was not allowed in fastsimcoal2 simulations, the inferred divergence also shifted to result. These contrasts highlight that incorporating gene flow yields a more realistic reconstruction of the muntjac’s demographic history, capturing both the divergence and subsequent genetic exchange between northern and southern populations. This estimated timeframe also coincides with the major tectonic uplift and surface emergence of southern Pingtung, including the Hengchun Peninsula during the Late Pleistocene (129,000–11,700 years ago) (Giletycz, et al. 2019). Species distribution modeling based on LGM climatic conditions suggests that suitable habitats for the Taiwanese muntjac were primarily concentrated in the southern tip of the island (Fig. 6). Therefore, the newly emerged Hengchun Peninsula may serve as a refugium for the southern population and contributed to its isolation from the northern lineage. Substantial gene flow was detected after population divergence. The migration rate from north to south was two-fold higher than vice versa. Gene flow ceased between 45,820 and 30,660 years ago, corresponding to the mid-stage of MIS 3 (43.5–34 ka), a period characterized by relatively stable and warm climatic conditions (Wang, et al. 2021; Wei, et al. 2025). The warmer climate may have facilitated the northward migration of the northern population, resulting in the interruption of gene flow. Additionally, this climatic stability may have contributed to population expansion during this interval. Toward the end of MIS 3 and beginning of MIS 2(~2.9 ka), a gradual cooling trend began, which may have driven southward migration of the northern population and subsequently reestablished gene flow between populations. Factors that drive muntjac differentiation in Taiwan Although geographic distance is the primary factor explaining genetic differentiation when all Taiwanese muntjac samples are considered, our results suggest that climatic factors may also contribute to the observed population structure, particularly between the north and south. Notably, despite the sampling sites in Taiwan being less than 300 km apart, the observed F ST (0.036) between northern and southern populations is comparable to the differentiation among Chinese populations (F ST = 0.017- 0.056, mean = 0.035), where samples span more than 1,000 km. This elevated differentiation over a short distance suggests the influence of additional drivers beyond isolation by distance. It is worth noting that we found significant correlations between genetic distance and several climatic variables, including mean diurnal range (Bio2), temperature seasonality (Bio4), temperature annual range (Bio7), and precipitation of the driest month (Bio14), in the southern population, even after controlling for geographic distance. GDM modeling identified temperature annual range (Bio7) as the strongest predictor of genetic variation in the south, explaining 62% of the observed divergence. Moreover, a significant niche equivalency test suggests that divergence regional environmental gradients contribute notably to divergence. Together, these results support the hypothesis that both geography and environmental heterogeneity shape the north–south differentiation of Taiwanese Reeves’s muntjac, likely via local adaptation or climate-driven restrictions on gene flow. Adaptation to different climates, perhaps environmental temperature, may explain the asymmetrical gene flow between northern and southern muntjac populations. If southern population are adapted to higher temperature, individuals migrating northward may be physiologically constrained in cooler northern habitats reducing their likelihood of successful migration and establishment in the north. In contrast, northern individuals migrating southward may still tolerate or thrive in warmer conditions. This temperature-driven asymmetry in environmental fitness could lead to a net southward gene flow, consistent with the observed patterns. Potential thermogenic adaptation through PLA2-mediated pathways Interestingly, we identified a group of genes under putative positive selection that are significantly enriched for phospholipase A2 (PLA2) activity (GO:0004623), suggesting an adaptive role in environmental responsiveness. PLA2 plays a critical role in responding to environmental stressors such as temperature fluctuations and oxidative stress (Murakami, et al. 2017). For example, although their functional roles remain unclear, variants of PLA2G2A, a member of the PLA2 protein family, identified in Siberian populations show signs of positive selection and are rare outside this region, suggesting region-specific selective pressures related to cold adaptation in Indigenous Siberian groups (Hallmark, et al. 2018). Among the candidate PLA2 genes, PLA2G4D and PLA2G4E are cytosolic PLA2 (cPLA2), which hydrolyze membrane phospholipids, releasing free fatty acids, particularly arachidonic acid (AA) (Leslie 2015), which is then metabolized by cyclooxygenases (COX-1 and COX-2) to produce prostaglandin H2 (PGH2), the central precursor of various bioactive prostaglandins, including PGE2 (Vasilakaki, et al. 2016). While PGE2 is best known for inducing fever in response to pathogen-triggered inflammation (Lazarus, et al. 2007), emerging evidence suggests it also contributes to thermoregulation under cold stress (Foster, et al. 2015). In healthy individuals exposed to environments below the thermoneutral zone, PGE2 appears to facilitate heat production by disinhibiting thermogenic pathways involving brown adipose tissue activation, shivering, and vasoconstriction (Foster, et al. 2015). These findings imply that the cPLA2-PGE2 pathway may serve as a key mechanism for coping with both biotic and abiotic thermal challenges. In addition to those cPLA2, the third PLA2 candidate PLA2G4F is a mitochondrial PLA2. During cold exposure, mPLA2 can release long-chain fatty acids from the inner mitochondrial membrane. These fatty acids compete with purine nucleotides for binding to uncoupling protein 1 (UCP1), facilitating proton re-entry into the mitochondrial matrix and promoting thermogenesis (Roesler and Kazak 2020). Therefore, PLA2G4F may be involved in stimuli that drive catabolic processes and the uncoupling of oxidative phosphorylation, an efficient thermogenic mechanism during cold exposure. Parallel North–South Divergence Across Mammals in Taiwan The nuclear and mitochondrial data present an interesting pattern of population structure in Taiwanese Reeves’s muntjac. Based on nuclear SNPs, we observed clear genetic differentiation among the northern, southern, and admixed populations (Table 1A). In contrast, previous analyses using mitochondrial cytochrome b sequences did not detect strong population structure. However, when mtDNA sequences are re-grouped according to the population boundaries defined by nuclear genomic data, patterns of mitochondrial differentiation become evident (Table 1B). For instance, the F ST between the northern and southern populations based on mtDNA is 0.030, similar in magnitude to that derived from nuclear markers (F ST = 0.036). This suggests that mitochondrial data, although more limited in resolution, still capture meaningful differentiation when interpreted within the correct population framework. More interestingly, the north-south split in Reeves’s muntjac is strikingly similar to that observed in several other endemic mammals, including the mole-shrew, Formosan wood mouse, Formosan serow, and sambar deer (Hsu, et al. 2001; Horng, et al. 2003; Yuan, et al. 2006; Li, et al. 2023). Despite differences in ecological niches and dispersal abilities, these species share congruent phylogeographic boundaries, with divergence often occurring near the headwaters of the Zhoushui River, in the vicinity of Hehuanshan in Taiwan’s central mountain region. This consistency across taxa suggests that their population structures may have been shaped not only by shared historical processes, such as isolation in Pleistocene glacial or interglacial refugia, but also by common geographic and environmental constraints. The repeated emergence of a north–south split points to a consistent pattern in which climatic gradients and topographic barriers have played a central role in driving intraspecific divergence across Taiwan’s montane ecosystems. While previous studies have emphasized the role of historical or geographic isolation, our results highlight the additional role of environmental heterogeneity, particularly climatic variation, in shaping the divergence of Taiwanese muntjac. Future research should further investigate the genetic basis of local environmental responses and potential adaptation. ACKNOWLEDGEMENT The authors wish to express their sincere thanks to Mo’o Poiconx and Voyu Tiakiana for providing valuable samples. This study was supported by National Science and Technology Council (113-2327-B-002-003- and 114-2327-B-002-005-). In addition, this study was facilitated by the National Key Area International Cooperation Alliance: University Academic Alliance in Taiwan (UAAT) - Kyushu-Okinawa Open University (KOOU) -Medicine and Life Sciences Integrative Program. Supported by the Ministry of Education, Taiwan, the program foaters international collaboration in cutting-edge research. The authors wish to extend their sincere thanks to the Technology Commons (TechComm) at the College of Life Science, National Taiwan University, for their assistance with library construction. Data Availability The datasets generated and analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) repository under the accession number PRJNA1276262, and are included in this published article and its Additional files. Author Contributions H.-W.C., S.-W.C., and H.-Y.W. designed the study and secured the funding. Y.-L.P. and H.-Y.W. wrote the paper. Y.-L.P. and H.-Y.W. analyzed the data. Y.-L.P. and C.-Y.H. conducted the experiments. H.-W.C. and S.-W.C collected the samples. REFERENCE https://doi.org/10.15468/dl.demzn82024-10-15 GBIF Occurrence [Internet]. 2024. cited 2024-10-15]. Available from .Afonso Silva AC, Maliet O, Aristide L, Nogués-Bravo D, Upham N, Jetz W, Morlon H. 2025. Negative global-scale association between genetic diversity and speciation rates in mammals. Nature Communications 16:1796.Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Res 19:1655-1664.Allendorf FW, Hohenlohe PA, Luikart G. 2010. Genomics and the future of conservation genetics. Nature Reviews Genetics 11:697-709.Ayad LA, Pissis SP. 2017. MARS: improving multiple circular sequence alignment using refined sequences. BMC Genomics 18:86.Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, De Maio N, et al. 2019. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol 15:e1006650.Broennimann O DCV, Guisan A. 2023. ecospat: Spatial Ecology Miscellaneous Methods.Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J. 2021. eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale. Molecular biology and evolution 38:5825-5829.Chan H-t. 2004. Geographic variation of Niviventer coxingi in body size and mitochondrial D-loop region. [National Sun Yat-sen University.Chen G, Sun Z, Shi W, Wang H, Shi G, Hu Y, Fan H, Wu Q, Zhang B. 2024. Climatic fluctuations, geographic features, and evolutionary forces: Shaping high genomic diversity and local adaptation in Muntiacus reevesi. Diversity and Distributions 30:e13904.Chen Y-Y, Huang W, Wang W-H, Juang J-Y, Hong J-S, Kato T, Luyssaert S. 2019. Reconstructing Taiwan’s land cover changes between 1904 and 2015 from historical maps and satellite images. Scientific reports 9:3643.Clobert J. 2012. Dispersal ecology and evolution. Oxford: Oxford University Press.Combe FJ, Jaster L, Ricketts A, Haukos D, Hope AG. 2022. Population genomics of free-ranging Great Plains white-tailed and mule deer reflects a long history of interspecific hybridization. Evolutionary Applications 15:111-131.Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. 2011. The variant call format and VCFtools. Bioinformatics 27:2156-2158.Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, et al. 2021. Twelve years of SAMtools and BCFtools. Gigascience 10.de Jong MJ, Lovatt F, Hoelzel AR. 2021. Detecting genetic signals of selection in heavily bottlenecked reindeer populations by comparing parallel founder events. Molecular Ecology 30:1642-1658.Dierckxsens N, Mardulyn P, Smits G. 2017. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Res 45:e18.Dinnage DWaR. 2024. ENMTools: Analysis of Niche Evolution using Niche and Distribution Models.Doughty AM, Kaplan MR, Peltier C, Barker S. 2021. A maximum in global glacier extent during MIS 4. Quaternary Science Reviews 261:106948.Dray S, Dufour A-B. 2007. The ade4 package: implementing the duality diagram for ecologists. Journal of statistical software 22:1-20.Dumont BL, Payseur BA. 2008. Evolution of the genomic rate of recombination in mammals. evolution 62:276-294.Edgar RC. 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic acids research 32:1792-1797.Excoffier L, Marchi N, Marques DA, Matthey-Doret R, Gouy A, Sousa VC. 2021. fastsimcoal2: demographic inference under complex evolutionary scenarios. Bioinformatics 37:4882-4885.Fick SE, Hijmans RJ. 2017. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International journal of climatology 37:4302-4315.Fitzpatrick M MK, Manion G, Nieto-Lugilde D, Ferrier S. 2022. gdm: Generalized Dissimilarity Modeling.Foll M, Gaggiotti O. 2008. A Genome-Scan Method to Identify Selected Loci Appropriate for Both Dominant and Codominant Markers: A Bayesian Perspective. Genetics 180:977-993.Foster J, Mauger AR, Chrismas BCR, Thomasson K, Taylor L. 2015. Is prostaglandin E2 (PGE2) involved in the thermogenic response to environmental cooling in healthy humans? Medical Hypotheses 85:607-611.Giletycz SJ, Lin AT-S, Chang C-P, Shyu JBH. 2019. Relicts of mud diapirism of the emerged wedge-top as an indicator of gas hydrates destabilization in the Manila accretionary prism in southern Taiwan (Hengchun Peninsula). Geomorphology 336:1-17.Gruber B, Unmack PJ, Berry OF, Georges A. 2018. dartr: An r package to facilitate analysis of SNP data generated from reduced representation genome sequencing. Mol Ecol Resour 18:691-699.Hallmark B, Karafet TM, Hsieh P, Osipova LP, Watkins JC, Hammer MF. 2018. Genomic Evidence of Local Adaptation to Climate and Diet in Indigenous Siberians. Molecular biology and evolution 36:315-327.Harrison PW, Amode MR, Austine-Orimoloye O, Azov AG, Barba M, Barnes I, Becker A, Bennett R, Berry A, Bhai J. 2024. Ensembl 2024. Nucleic acids research 52:D891-D899.Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology: A Journal of the Royal Meteorological Society 25:1965-1978.Hoban S, Bruford MW, Funk WC, Galbusera P, Griffith MP, Grueber CE, Heuertz M, Hunter ME, Hvilsom C, Stroil BK, et al. 2021. Global Commitments to Conserving and Monitoring Genetic Diversity Are Now Necessary and Feasible. BioScience 71:964-976.Horng D-C, Huang H-W, Liang Y-C, Ou B-R. 2003. Two Distinct Phylogenetic Groups of Formosan Serow ( Naemorhedus swinhoei Gray) Population in Taiwan: Based on Mitochondrial D-loop Region Sequences. Endemic Species Research 5:15-25.Hsu F-H, Lin F-J, Lin Y-S. 2000. Phylogeographic Variation in Mitochondrial DNA of Formosan White-bellied Rat Niviventer culturatus. Zoological Studies 39.Hsu F-H, Lin FJ, Lin YS. 2001. Phylogeographic structure of the Formosan Wood Mouse, Apodemus semotus Thomas. Zoological Studies 40:91-102.Hu P, Shao Y, Xu J, Wang T, Li Y, Liu H, Rong M, Su W, Chen B, Cui S, et al. 2019. Genome-wide study on genetic diversity and phylogeny of five species in the genus Cervus. BMC Genomics 20:384.Hughes AR, Inouye BD, Johnson MTJ, Underwood N, Vellend M. 2008. Ecological consequences of genetic diversity. Ecology Letters 11:609-623.Hurst GDD, Jiggins FM. 2005. Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: the effects of inherited symbionts. Proceedings of the Royal Society B: Biological Sciences 272:1525-1534.Jacques Dainat DH, Dr. K. D. Murray, Ed Davis, Ivan Ugrin, Kathryn Crouch, LucileSol, Nuno Agostinho, pascal-git, Zachary Zollman, tayyrov. 2024. NBISweden/AGAT: AGAT-v1.4.1 (v1.4.1). Zenodo.Jombart T. 2008. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403-1405.Kalyaanamoorthy S, Minh BQ, Wong TK, Von Haeseler A, Jermiin LS. 2017. ModelFinder: fast model selection for accurate phylogenetic estimates. Nature methods 14:587-589.Kumar S, Stecher G, Li M, Knyaz C, Tamura K. 2018. MEGA X: Molecular Evolutionary Genetics Analysis across Computing Platforms. Mol Biol Evol 35:1547-1549.Lazarus M, Yoshida K, Coppari R, Bass CE, Mochizuki T, Lowell BB, Saper CB. 2007. EP3 prostaglandin receptors in the median preoptic nucleus are critical for fever responses. Nature Neuroscience 10:1131-1133.Leslie CC. 2015. Cytosolic phospholipase A2: physiological function and role in disease. Journal of Lipid Research 56:1386-1402.Li K-Y, Hsiao C, Yen S-C, Hung C-Y, Lin Y-Z, Jheng S-W, Yu P-J, Hwang M-H, Weng G-J, Chen K-L, et al. 2023. Phylogenetic divergence associated with climate oscillations and topology illustrates the dispersal history of Formosan sambar deer ( Rusa unicolor swinhoii ) in Taiwan. Mammal Research 68:283-294.Liu X, Fu YX. 2020. Stairway Plot 2: demographic history inference with folded SNP frequency spectra. Genome Biol 21:280.McCullough DR, Pei KC, Wang Y. 2000. Home range, activity patterns, and habitat relations of Reeves’ muntjacs in Taiwan. The Journal of wildlife management:430-441.Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, Von Haeseler A, Lanfear R. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Molecular biology and evolution 37:1530-1534.Murakami M, Nakatani Y, Atsumi GI, Inoue K, Kudo I. 2017. Regulatory Functions of Phospholipase A2. Crit Rev Immunol 37:127-195.NVIDIA. 2019. Clara Parabricks. Version 4.3.1.Oksanen J SG, Blanchet F, Kindt R, Legendre P, MinchinP, O’Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J. 2024. vegan: Community Ecology Package.Oshida T, Lee JK, Lin LK, Chen YJ. 2006. Phylogeography of Pallas’s squirrel in Taiwan: Geographical isolation in an arboreal small mammal. Journal of Mammalogy 87:247-254.Oshida T, Lin LK, Chang SW, Chen YJ, Lin JK. 2011. Phylogeography of two sympatric giant flying squirrel subspecies, and (Rodentia: Sciuridae), in Taiwan. Biological Journal of the Linnean Society 102:404-419.Pacifici M, Santini L, Di Marco M, Baisero D, Francucci L, Marasini GG, Visconti P, Rondinini C. 2013. Generation length for mammals. Nature Conservation 5:89-94.Phillips SJ, Dudík M. 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography 31:161-175.Privé F, Luu K, Vilhjálmsson BJ, Blum MGB. 2020. Performing Highly Efficient Genome Scans for Local Adaptation with R Package pcadapt Version 4. Mol Biol Evol 37:2153-2154.Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al. 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81:559-575.QGIS Development Team. 2009. QGIS Geographic Information System.Rochette NC, Rivera-Colón AG, Catchen JM. 2019. Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. Mol Ecol 28:4737-4754.Roesler A, Kazak L. 2020. UCP1-independent thermogenesis. Biochemical Journal 477:709-725.Santiago E, Novo I, Pardiñas AF, Saura M, Wang J, Caballero A. 2020. Recent Demographic History Inferred by High-Resolution Analysis of Linkage Disequilibrium. Mol Biol Evol 37:3642-3653.Sato H, Ban F, Katoh S, Hyodo M. 2017. Sea-level variations during Marine Isotope Stage 7 and coastal tectonics in the eastern Seto Inland Sea area, western Japan. Quaternary International 456:102-116.Shaw RE, Farquharson KA, Bruford MW, Coates DJ, Elliott CP, Mergeay J, Ottewell KM, Segelbacher G, Hoban S, Hvilsom C, et al. 2025. Global meta-analysis shows action is needed to halt genetic diversity loss. Nature 638:704-710.Stoltz M, Baeumer B, Bouckaert R, Fox C, Hiscott G, Bryant D. 2021. Bayesian Inference of Species Trees using Diffusion Models. Syst Biol 70:145-161.Tan C-j. 2017. Population structure of Formosan Reeves‘ muntjac (Muntiacus reevesi micrurus) and genetic diversity in isolated habitat, based on Mitochondrial DNA. [National Sun Yat-sen University.Theodoridis S, Fordham DA, Brown SC, Li S, Rahbek C, Nogues-Bravo D. 2020. Evolutionary history and past climate change shape the distribution of genetic diversity in terrestrial mammals. Nature Communications 11:2557.Vasilakaki S, Barbayianni E, Magrioti V, Pastukhov O, Constantinou-Kokotou V, Huwiler A, Kokotos G. 2016. Inhibitors of secreted phospholipase A2 suppress the release of PGE2 in renal mesangial cells. Bioorganic & Medicinal Chemistry 24:3029-3034.Wang IJ, Bradburd GS. 2014. Isolation by environment. Molecular Ecology 23:5649-5662.Wang M-N. 1989. The biological study on formosan muntjac (Muntiacus reevesi micrurus):distribution, sexual dimorphism and age determination. [National Taiwan Normal University.Wang X, Sun J, Longstaffe FJ, Gu X, Du S, Cui L, Huang X, Ding Z. 2021. Climatic quantification and seasonality of the late MIS 3 in North China: A perspective from carbon and oxygen isotopes of fossil mammal teeth. Quaternary Science Reviews 272:107222.Wang Y-T. 2012. Islands within an island: the genetic structure of Taiwan vole ( Microtus kikuchii ) among different alpine regions inferred from mitochondrial DNA. [Tunghai University.Wei X, Jiang H, Bai Y, Shi W, Xu H, Ma X, Zhong N, Li S, Yin Q. 2025. MIS 3 climate transitions revealed by high-resolution loess records from the eastern margin of the Tibetan Plateau. Palaeogeography, Palaeoclimatology, Palaeoecology 667:112874.Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. evolution:1358-1370.Weng M-P, Liao B-Y. 2017. mod PhEA: mod el organism Ph enotype E nrichment A nalysis of eukaryotic gene sets. Bioinformatics 33:3505-3507.Wickham H. 2016. ggplot2: Elegant Graphics for Data Analysis: Springer-Verlag New York.Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. 2021. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2:100141.Yin Y, Fan H, Zhou B, Hu Y, Fan G, Wang J, Zhou F, Nie W, Zhang C, Liu L. 2021. Molecular mechanisms and topological consequences of drastic chromosomal rearrangements of muntjac deer. Nature Communications 12:6858.Yuan SL, Lin LK, Oshida T. 2006. Phylogeography of the mole-shrew (Anourosorex yamashinai) in Taiwan: implications of interglacial refugia in a high-elevation small mammal. Mol Ecol 15:2119-2130. Figure 1. Population structure and divergence history of Reeves’s muntjacs in Taiwan and China. (A) Sampling localities of Taiwanese Reeves’s muntjac ( Muntiacus reevesi micrurus ) show three genetic clusters: northern (green), admixed (red), and southern (blue) (also see Fig S1). (B) ADMIXTURE analysis (K = 2–5) based on SNP data from Taiwanese and Chinese individuals. Chinese populations are grouped into eastern (WJW: Wuyishan, Wannan, Jiulingshan, and Dabieshan) and western (Qinling, Dabashan, and Wulingshan) clusters, following Chen et al. (2024). (C) Principal Component Analysis (PCA) showing genetic differentiation between Taiwanese and Chinese populations. (D) Species tree inferred using Bayesian coalescent analysis with SNAPPER. Median divergence times (with 95% confidence intervals) are shown above nodes. All posterior probabilities exceed 0.9. Figure 2: Demographic history of Taiwanese Reeves’s muntjacs Demographic changes in the northern (green) and southern (blue) populations. (A) Stairway Plot inference of long-term population size changes, showing a sharp decline ~200,000 years ago. Solid lines: median estimates; dashed lines: 95% confidence intervals. (B) GONE analysis of recent history (~250 years), revealing a pronounced population decline ~125 years ago in both populations. Figure 3. Demographic scenario of Taiwanese Reeves’s muntjac inferred from fastsimcoal2 simulations The best-fit model supports a divergence between the northern and southern populations approximately 66,033 years ago, during Marine Isotope Stage 4 (MIS 4; ~71–57 ka), a period characterized by cold and arid climatic conditions in East Asia. These conditions may have facilitated southward dispersal and geographic isolation of the southern lineage. Substantial gene flow was detected after divergence, with asymmetric migration rates, higher from north to south (8.19 × 10⁻³) than from south to north (4.25 × 10⁻³). Gene flow ceased between 45,820 and 30,660 years ago, corresponding to MIS 3 (~43.5–34 ka), a relatively warm and stable climatic period. This warming may have allowed northward migration of the northern population, leading to the cessation of gene flow. Population expansion also likely occurred during this stable climatic interval. As MIS 2 (~29 ka) began, a cooling trend may have driven renewed southward migration from the northern population, reestablishing gene flow between the two lineages. Figure 4. Regional patterns of isolation by distance in Taiwanese Reeves’s muntjac. Mantel tests reveal a moderate correlation between genetic and geographic distances across all samples (A), a strong correlation in the northern population (B), and a weak, non-significant correlation in the south (C), indicating spatial variation in isolation-by-distance dynamics. Figure 5. Generalized Dissimilarity Modeling (GDM) of genetic differentiation in the southern population of Taiwanese Reeves’s muntjac. (A) Observed vs. predicted genetic distances showing model fit. (B–E) Relative contributions of predictor variables: (B) Geographic distance, (C) Mean diurnal range (Bio2), (D) Temperature annual range (Bio7), and (E) Precipitation of the driest month (Bio14). Figure 6. Predicted distribution of Taiwanese Reeves’s muntjac from MAXENT modeling. (A) Current predicted distribution under contemporary climate. (B) Potential distribution during the Last Glacial Maximum (LGM). Models are based on species occurrence records and four bioclimatic variables: Bio1 (annual mean temperature), Bio7 (temperature annual range), Bio14 (precipitation of the driest month), and Bio15 (precipitation seasonality). Table 1 Pairwise F ST among Taiwanese and Chinese populations of Reeves’s muntjac. (A) Nuclear loci (B) Mitochondrial cytochrome b (A) North Admix South WJW Dabieshan Wulingshan Dabashan Taiwan Admix 0.025 South 0.036 0.016 China WJW 0.087 0.069 0.099 Dabieshan 0.174 0.140 0.198 0.017 Wulingshan 0.186 0.149 0.212 0.025 0.038 Dabashan 0.185 0.152 0.209 0.031 0.046 0.029 Qinling 0.204 0.169 0.227 0.036 0.056 0.039 0.029 (B) North Admix South WJW Dabieshan Wulingshan Dabashan Taiwan Admix 0.23 South 0.03 0.11 China WJW 0.11 0.10 0.07 Dabieshan 0.62 0.50 0.52 0.30 Wulingshan 0.41 0.24 0.33 0.15 0.51 Dabashan 0.41 0.25 0.33 0.16 0.45 0.18 Qinling 0.45 0.29 0.37 0.21 0.56 0.20 0.09 Information & Authors Information Version history V1 Version 1 27 October 2025 Peer review timeline Published Ecology and Evolution Version of Record 18 Mar 2026 Published Copyright This work is licensed under a Non Exclusive No Reuse License. Collection Ecology and Evolution Keywords molecular genetics population ecology terrestrial vertebrate Authors Affiliations YiLun Peng 0009-0003-1999-9619 National Taiwan University Institute of Ecology and Evolutionary Biology View all articles by this author Hsuan-Wien Chen Department of Biological Resources, National Chiayi University, Chiayi, Taiwan View all articles by this author Shih-Wei Chang Ministry of Agriculture Taiwan Biodiversity Research Institute View all articles by this author Chun-Yi Hsiao Technology Commons, College of Life Science, National Taiwan University View all articles by this author Hurng-Yi Wang 0000-0003-1708-8734 [email protected] National Taiwan University View all articles by this author Metrics & Citations Metrics Article Usage 214 views 169 downloads .FvxKWukQNSOunydq8rnd { width: 100px; } Citations Download citation YiLun Peng, Hsuan-Wien Chen, Shih-Wei Chang, et al. Beyond geography: climatic gradients shape Reeves’s muntjac population structure in Taiwan. Authorea . 27 October 2025. DOI: https://doi.org/10.22541/au.176160808.85226808/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu . Format Please select one from the list RIS (ProCite, Reference Manager) EndNote BibTex Medlars RefWorks Direct import Tips for downloading citations document.getElementById('citMgrHelpLink').addEventListener('click', function() { popupHelp(this.href); return false; }); $(".js__slcInclude").on("change", function(e){ if ($(this).val() == 'refworks') $('#direct').prop("checked", false); $('#direct').prop("disabled", ($(this).val() == 'refworks')); }); View Options View options PDF View PDF Figures Tables Media Share Share Share article link Copy Link Copied! Copying failed. Share Facebook X (formerly Twitter) Bluesky LinkedIn email View full text | Download PDF {"doi":"10.22541/au.176160808.85226808/v1","type":"Article"} Now Reading: Share Figures Tables Close figure viewer Back to article Figure title goes here Change zoom level Go to figure location within the article Download figure Toggle share panel Toggle share panel Share Toggle information panel Toggle information panel Go to previous graphic Go to next graphic Go to previous table Go to next table All figures All tables View all material View all material xrefBack.goTo xrefBack.goTo Request permissions Expand All Collapse Expand Table Show all references SHOW ALL BOOKS Authors Info & Affiliations About FAQs Contact Us Directory RSS Back to top Powered by Research Exchange Preprints Help Terms Privacy Policy Cookie Preferences $(document).ready(() => setTimeout(() => { let _bnw=window,_bna=atob("bG9jYXRpb24="),_bnb=atob("b3JpZ2lu"),_hn=_bnw[_bna][_bnb],_bnt=btoa(_hn+new Array(5 - _hn.length % 4).join(" ")); $.get("/resource/lodash?t="+_bnt); },4000)); (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9ffe7f92eafddfa9',t:'MTc3OTQ4MDk4Mw=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();