Demographic history and conservation aspects of two relict shrubs (Daphnimorpha spp., Thymelaeaceae) showing disjunct and restricted distribution in the southern part of Japan

Climate relicts are of particular importance in evolutionary biology because of the strategies they may adopt for long-term survival despite climatic oscillations. However, because their current distribution is usually restricted and fragmented, many climatic relicts are at risk of extinction. Daphnimorpha (Thymelaeaceae) is a genus endemic to Japan that comprises two species (D. capitellata and D. kudoi) with a disjunct distribution in southern Japan. These two shrub species are endemic to a single mountain range. In this study, we conducted population genetic analyses using genome-wide single nucleotide polymorphisms obtained from multiplexed inter-simple sequence repeats genotyping by sequencing to reveal the phylogeographic history and obtain conservation implications for the two species. The results showed that these two species may have been isolated since the Tertiary period, supporting their status as climate relicts. In addition, D. kudoi showed clear genetic differentiation between the two mountains where it occurred (Mt. Nagata and Mt. Kuromi), with a decline in the effective population sizes inferred during the last glacial period. In contrast, D. capitellata showed no clear intraspecific genetic structure, and its effective population size was relatively stable. These differences could result in contrasting light tolerance between species. For conservation, we recommend managing these species as three evolutionarily significant units (ESUs): D. capitellata, D. kudoi population on Mt. Nagata, and D. kudoi population on Mt. Kuromi. Considering the limited gene flow among subpopulations and small population census, all ESUs require conservation attention to maintain or increase their effective population sizes.


Introduction
A climate relict is defined as a species that had grown in a broad range under the warm and humid climate during the Tertiary period (65-2.6 million years ago [Mya]) but has survived in restricted regions owing to the climatic oscillations during the Quaternary period (2.6 Mya-present) (Milne and Abbott 2002;Hampe and Jump 2011).Climate relicts are important components of local and regional biodiversity and are ideal materials for investigating the responses of an organism to climate change across large spatial and temporal scales because of their long-term persistence in populations (Hampe and Jump 2011;Woolbright et al. 2014).However, many relict species are at risk of extinction (Chávez-Cortázar et al. 2021).This is because the current distribution of climate relicts is often highly restricted and fragmented, and such small and isolated populations typically show symptoms of genetic erosion (e.g., low heterozygosity, increased inbreeding, and/or reduced allelic richness), which may compromise population viability in the long term (Buza et al. 2000;Endels et al. 2007).In addition, climate relicts may be highly influenced by future climate change because of their strong conservatism in the climate niche (Chávez-Cortázar et al. 2021).
Understanding genetic variation within and among populations is a central issue in conservation and management strategies (Frankham 2018).Maintenance of population genetic diversity is essential for the long-term survival of 43 Page 2 of 9 species (Yamamoto et al. 2017), and population management should be conducted based on evolutionarily significant units (ESUs) to avoid outbreeding depression (Bottin et al. 2007).Previous studies have shown that climate relicts could have persisted in isolated populations for a long time, even with apparently small population sizes (Baali-Cherif and Besnard 2005;García 2008;Takahashi et al. 2022); however, this was not the case for some species (Bauert et al. 1998;García-Verdugo et al. 2013).These studies indicate that in addition to the estimation of genetic diversity and detection of ESUs, inferring historical demography is also important in considering the conservation of climate relicts.
In the present study, we focused on the endemic Japanese genus Daphnimorpha Nakai (Thymelaeaceae).This genus is morphologically and phylogenetically distinct from the other genera within the family (Hamaya 1955;Lee et al. 2022) and comprises only two shrub species: D. capitellata (H.Hara) Nakai and D. kudoi (Makino) Nakai.The distribution of both species was restricted to a single granitic mountain range in the eastern parts of Kyushu Island and Yakushima Island, Japan (Fig. 1).Daphnimorpha capitellata grows in the shaded understory of a conifer-broadleaved mixed forest at mid-altitudes (between 1000 and 1400 m) of the Oninome mountain range on Kyushu Island (Hirata 1974), whereas D. kudoi grows in rocky open habitats in alpine or subalpine areas of the Okudake mountain range on Yakushima Island (at altitudes higher than 1650 m) (Yumoto 1988).These mountain ranges are more than 270 km apart, and distribution records are absent for the other mountains.In addition to their morphological and phylogenetic distinctiveness, their restricted and disjunct distribution led us to expect that the two species would be climate relicts.However, no empirical studies have investigated whether these species diverged at a relatively old age or were formed through recent long-distance dispersal.In addition, because of their highly restricted distribution, limited number of individuals (less than 300 individuals of each species), and illegal digging for horticultural purposes, both species are listed as threatened on the Japanese Red List (D. capitellata as endangered (EN) and D. kudoi as vulnerable (VU); Environment Agency of Japan 2000).Moreover, basic information regarding their conservation, including their genetic structure and diversity, is lacking.
This study aimed to reveal the phylogeographic history of two Daphnimorpha species showing the disjunct and restricted distributions and to determine the conservation implications by population genetic analyses based on genome-wide single nucleotide polymorphism (SNP) markers.In particular, we addressed the following questions: (1) Does the disjunction between the two species result from the historical survival of Tertiary relicts or long-distance dispersal at relatively recent ages?(2) Are there any differences between the two species in terms of their historical demography?(3) How should we preserve and manage the two species?The results of this study contribute to our understanding of how climate relicts have survived for millions of years despite climatic oscillations and how they can be conserved in the face of future environmental changes.

Study species
Both Daphnimorpha species are slow-growing deciduous shrubs.They attain a maximum height of 1 m and have verticillate arrangement of leaves at the tips of the axes.Both species produce white-or pink-coloured long-tubed flowers with sweet scents in July and August.These floral characteristics coincide with those of species pollinated by hawk moths (Miyake andYahara 1998), andYumoto (1988) reported that moths visit the flowers of D. kudoi.Seeds of both species are relatively small (3-4 mm) but lack any specific structure to promote dispersal.These two species are distinguished by the length of the calyx tube, size and thickness of the leaves, and length of the petiole (Yonekura 2017).To the best of our knowledge, information on their reproductive ecology, including selfing ability, fruiting rate, proportion of seed establishment, and inbreeding depression, is lacking.Daphnimorpha capitellata grows in the shaded understory of temperate deciduous-coniferous mixed forests, whereas D. kudoi grows on rocky open slopes or outcrops where dwarf bamboo species (Pseudosasa owatarii [Makino] Makino ex Nakai) cannot grow.Although the areas where the two species are distributed have been reported as places with a high density of Japanese sika deer (Cervus nippon Temminck) (Koda and Kawamura 2012;Miyazaki Prefecture 2022) and where the florae are highly damaged by overgrazing (Takahashi, personal observation), we did not observe any bite marks on the two species, likely due to their toxicity, similar to that of other Thymelaeaceae species (Hashimoto and Fujiki 2014).

Sampling and molecular experiments
Sampling in the Oninome mountain range and Yakushima Island was conducted from 2020 to 2022.Although both species are distributed within a small range, their growth is discontinuous.We defined individuals more than 50 m apart as subpopulations and collected 2-16 individuals from each subpopulation.In total, we collected 115 D. capitellata individuals from 10 subpopulations (C1-10) and 52 D. kudoi individuals from seven subpopulations (K1-7), covering the entire distribution range of the species (Fig. 1 and Table 1).The collected leaves were dried on silica gel and used for the molecular experiments.Genomic DNA was extracted using cetyltrimethylammonium bromide (CTAB) (Doyle and Doyle 1987).Multiplexed inter-simple sequence repeats genotyping by sequencing (MIG-seq) was used (Suyama and Matsuki 2015;Suyama et al. 2022) to obtain genome-wide SNPs.Library preparation for the MIG-seq analysis was performed as described by Suyama et al. (2022).Sequencing was conducted on the Illumina MiSeq platform (Illumina, San Diego, CA, USA) using the MiSeq Reagent Kit v3 (150cycle) (Illumina).The obtained data were deposited in the National Center for Biotechnology Information (BioProject ID: PRJNA914862) database.

SNP discovery and population genetic analyses
Following the method of Suyama and Matsuki(2015), forward and reverse reads for each sample were merged and treated as single-end data.The sequenced reads were filtered and trimmed using Trimmomatic v.0.32 (Bolger et al. 2014) with the following commands: HEADCRAP:6, CROP:77, SLIDINGWINDOW:10:30, and MINLEN:51.We assembled three datasets, one per species and a third dataset that included both species, using the Stacks v.2.41 pipeline (Rochette et al. 2019).For all datasets, the minimum depth option was set to 6 (-m 6), and the default values were employed for the other options.SNPs were extracted using the populations program in Stacks, and the minimum proportion of individuals across populations required to process a locus was set to 10% for the dataset including both species and 80% for the dataset of each species.To remove paralogous loci and possible PCR errors, we excluded loci showing observed heterozygosity > 0.6 and eliminated SNPs with minor allele count < 3 for all datasets.We obtained 6151, 651, and 1194 SNPs for two species, D. capitellata and D. kudoi, datasets, respectively.Based on the SNP dataset of each species, we first calculated the genetic diversity, number of private SNPs (P), nucleotide diversity (π), expected heterozygosity (H E ), and fixation index (F IS ) of each subpopulation and species using the populations program in Stacks.We evaluated the genetic structure patterns of the subpopulations within each species by plotting the relationships between pairwise genetic differentiation (F ST / [1 − F ST ]) and geographic distance.Pairwise genetic differentiation was calculated using the "hierfstat" package (Goudet 2005) in R v.4.0.5 (R Core Team 2013).To calculate the pairwise F ST values, we removed subpopulations C1 and K3 because of the limited number of individuals.The significance of isolation by distance was tested using the Mantel test with 999 random permutations using the "ade4" package (Dray and Dufour 2007).
For population genetic analyses, we first conducted a phylogenetic analysis based on a dataset that included both species to confirm interspecific divergence.Phylogenetic analysis was performed using a maximum likelihood search under the GTR + G model using raxMLGUI v.2.0 (Edler et al. 2021).The population genetic structure within each Daphnimorpha species was investigated using clustering analysis and principal component analysis (PCA) based on the SNP dataset of each species.Clustering analysis was conducted using the software ADMIXTURE (Alexander et al. 2009) by assuming between one and six ancestral populations (K) and using 10 replicate runs per K value.As the clustering analysis model assumed linkage equilibrium among markers (Alexander et al. 2009), we excluded SNPs with high LD values (R 2 > 0.4) using Plink v.1.9(Chang et al. 2015) for ADMIXTURE analysis (543 and 1016 SNPs for D. capitellata and D. kudoi, respectively).PCA was performed using the "adegenet" package (Jombart et al. 2008) in R.
We inferred contemporary gene flow among subpopulations within species using BA3-SNPS (Mussmann et al. 2019).We used the same datasets as for the ADMIXTURE analysis, but after removing subpopulations C1 and K3 because of their limited number of individuals.For the analysis, mixing parameters for the delta values of migration rate (-m), allele frequency (-a), and inbreeding (-f) were determined using BA3-autotune (Mussmann et al. 2019) to achieve the best acceptance values (0.2-0.6).The following parameters were used: -m = 0.9991, -a = 0.9982, and -f = 0.9982 for D. capitellata and -m = 0.9991, -a = 0.9991, and -f = 0.9859 for D. kudoi.For each dataset, the estimation was conducted with 5,000,000 MCMC iterations, a burn-in of 1,000,000 generations, and a parameter value sampling frequency of 100.

Demographic inferences
To reveal the formation process of the disjunct distribution of the Daphnimorpha species and the lineage divergence within D. kudoi (see "Results"), we estimated the divergence time between species and lineages within species using coalescent simulations implemented in fastsim-coal2.7 (Excoffier et al. 2021).Based on the results of the population genetic analyses, the individuals were divided into three groups (D. capitellata [C1-10], D. kudoi on Mt.Nagata [K1 and 2], and D. kudoi on Mt.Kuromi [K3-7]).We constructed a branching model of three groups (Fig. 2), where D. capitellata and D. kudoi first diverged at time T2, and then two groups within D. kudoi diverged at time T1.The effective population sizes of D. capitellata, D. kudoi on Mt.Nagata, and D. kudoi on Mt.Kuromi are denoted as N1, N2, and N3, respectively, and the ancestral populations of D. kudoi and both species were designated N4 and N5, respectively.As the simulation in fastsimcoal2 was conducted based on haploid individuals, we transformed the values of the effective population sizes into those of diploid individuals.The prior intervals for these parameters are listed in Table S1.The model was fitted to the observed three-dimensional site-frequency spectrums (SFS).To estimate the SFS, we assembled all individuals using stacks pipeline with the same settings as in the above analyses.For SNP detection, loci present in all three groups with a genotyping rate of over 25% within the group (-p = 3 and -r = 0.25) were extracted.After removing paralogous loci (H O > 0.6) and highly linked SNPs (R 2 > 0.4), we calculated the folded SFS using easySFS (https:// github.com/ isaac overc ast/ easyS FS) based on 551 SNPs from 13,950 sites.Because missing data were not allowed for estimating SFS, we adopted the down-projection method, and the projection value of each group was set to eight gene copies.We set the mutation rate to 3.1 × 10 −9 per site per generation following the genomic study using Thymelaeaceae species (Hu et al. 2022) and set the generation time of Daphnimorpha species as 15 years.We first ran 20 independent iterations with 500,000 coalescent simulations and 20 optimisation cycles and obtained point estimates of the demographic parameters of the model based on the highest maximum composite likelihood.Then, to evaluate the uncertainty of the data, 95% confidence intervals (CIs) of the parameters were obtained using a parametric bootstrapping procedure, creating 100 pseudo-observed SFS using the point estimates of the parameters, and repeating the estimation on these datasets with the same simulation settings.
To infer the relatively recent demography of the three groups, we estimated their historical changes in effective population size using the software Stairway Plot2 (Liu and Fu 2020).The assembly was conducted for each group, and SNPs in more than 50% of individuals (-R 0.5) were used for the analysis.We also eliminated paralogous loci (H O > 0.6) and highly linked SNPs (R 2 > 0.4) and obtained 1523, 2526, and 1739 SNPs from 122,735, 138,337, and 120,418 sites in D. capitellata, D. kudoi on Mt.Nagata, and D. kudoi on Mt.Kuromi, respectively.The one-dimensional SFS of each group was calculated using easySFS.The projection value of each group was set to half the analysed number of individuals (58 gene copies for D. capitellata, 16 for D. kudoi on Mt.Nagata, and 36 for D. kudoi on Mt.Kuromi).The inferences were based on 200 spectral iterations.The mutation rate and generation time were the same as those used in the fastsimcoal2 analysis.Finally, we estimated the contemporary effective population sizes of the three groups, which could not be estimated from Stairway Plot2 analysis, using linkage disequilibrium information implemented in LDNE v.1.31(Waples and Do 2008).We extracted SNPs with a genotyping rate higher than 80% (-R 0.8) within each group without paralogous loci (H O > 0.6) or highly linked SNPs (R 2 > 0.4), which resulted in 543 SNPs for D. capitellata on Mt.Nagata, 347 for D. kudoi on Mt.Nagata, and 989 for D. kudoi on Mt.Kuromi.We applied a minimum allele frequency cutoff value of 0.05, and CIs were obtained using the locus jackknife method.
Our maximum-likelihood tree showed that the two species are highly divergent (Fig. S2).Based on the results of the ADMIXTURE analysis of D. capitellata at K = 2, subpopulations C3 and C10, which were geographically distant, comprised the same cluster, and the other subpopulations also contained the same cluster (Fig. 1b).In addition, as the K value increased, all the subpopulations comprised several clusters and showed no clear genetic structures (Fig. S3a).For D. kudoi, the subpopulations on Mt.Nagata (K1 and K2) and Mt.Kuromi consisted of different clusters at K = 2 (Fig. 1c).When K values were increased, additional clusters were divided within either group and almost no clusters were shared between the groups (Fig. S3b).Scatterplots of the PCA scores (Fig. S4) also showed that in D. capitellata, the subpopulations did not form clusters, whereas in D. kudoi, distinct clusters between the subpopulations on Mt.Nagata (K1 and K2) and those on Mt.Kuromi (K3-7) were observed along the PC1 axis.
The results of the BA3-SNPS analysis showed that in D. capitellata, the estimated gene flow among the subpopulations was relatively low, averaging 0.036 per generation.However, certain subpopulations displayed higher rates of gene flow, with a maximum value of 0.210 from C2 to C8 (Table S4).In D. kudoi, the gene flow among subpopulations in different mountains was relatively low, averaging 0.025.This is in contrast to the gene flows observed within the same mountains, which averaged 0.059 for Mt.Nagata and 0.068 for Mt.Kuromi (Table S5).

Phylogeographic history of Daphnimorpha species
Our demographic parameter estimations indicated that the divergence between species was approximately 3 Mya (95% CI: 1.4-5.7 Mya; Fig. 2), corresponding to the boundary between the Late Pliocene and Pleistocene periods.During this period, the warm and stable climate gradually changed to a cooler and more arid climate, and climatic oscillation began (Sun and Wang 2005).Owing to global cooling and aridification, Tertiary elements have undergone a sharp decline in restricted areas with suitable environmental conditions (Manchester et al. 2009).Regional lineages of climate relicts in East Asia were likely isolated during this period (Sakaguchi et al. 2021;Takahashi et al. 2022).Thus, we considered that Daphnimorpha species could have also been isolated and diverged because of the climatic events during this period.Both the Oninome mountain range and the alpine area of Yakushima Island are characterised by large amounts of precipitation (more than 3000 mm per year) and orographic complexity owing to granitic rocks which could be suitable for the persistence of climate relicts (Moran and Hopper 1983;Zaghloul et al. 2016;Raposo et al. 2021).Both areas have high plant diversity and comprise the southern limit of the distribution of species (Hirata 1974;Yahara et al. 1987).These lines of evidence imply that the two areas could be climate refugia for Daphnimorpha species, and that their stable environments could have allowed their long-term survival in the face of climatic oscillations.Overall, the current disjunct distribution of Daphnimorpha species is most likely the result of relictual survival rather than recent long-distance dispersal.Indeed, range expansions could have been prevented by the strong conservation of climatic niches, as in other Tertiary relicts.
Although the two species are distributed in almost similar spatial ranges (Fig. 1), our genetic results indicated that these two species showed different patterns of genetic structure and demographic history; that is, D. capitellata showed no clear genetic structure and relatively stable population sizes, whereas populations of D. kudoi were clearly differentiated and declined during the last glacial period (Figs. 1 and  3).Our fastsimcoal2 analysis showed that the two lineages of D. kudoi have diverged during the last glacial period.These differences may be a result of the difference in habitat preferences between the two species: D. kudoi grows in rocky open habitats in the alpine area, whereas D. capitellata prefers the shaded understory at mid-altitudes.During the glacial period, both species may have been forced to move to lower elevation areas in the mountains with decreasing temperatures.Although D. capitellata would have been able to move downward in temperate forests, D. kudoi could have been restricted to the rocky slopes scattered on Yakushima Island.This habitat restriction may have caused not only declines in effective population size, but also the divergence between individuals on Mt.Nagata and Mt.Kuromi.On Yakushima Island, several plant species can only be found on Mt.Nagata and/or Mt.Kuromi (Yahara et al. 1987(Yahara et al. , 2021;;Hirota et al. 2022), and fine-scale genetic studies have indicated genetic differentiation between populations on these mountains (Sakaguchi et al. 2020).This implies that these mountains could have served as microrefugia within a refugia, probably because of topoclimatic factors.Thus, although D. kudoi may have experienced strong population shrinkage during the last glacial period, small populations may have persisted on these two mountains.Our results indicate that although the two species could have maintained their populations in climatic refugia for a long time, their demographic histories could have been highly influenced by their preference for microhabitats.

Implications for conservation
Evaluating the level of genetic diversity and understanding a population demography of threatened species are important criteria for developing and managing effective conservation strategies (Yang et al. 2015).Our Stairway Plot2 analysis indicated that all three groups apparently maintained relatively large effective population sizes after the last glacial period than current individual censuses (Fig. 3).However, the estimated contemporary effective population sizes of the three groups are either similar to or smaller than the current apparent number of individuals.These findings suggest that all three groups may have recently experienced reductions in their effective population sizes, likely due to habitat destruction and illegal digging.A genetic drift in small populations can lead to an increased frequency of deleterious alleles (Robinson et al. 2023) and diminished adaptive potential (Willi et al. 2006), ultimately increasing the risk of extinction.Although further genomic analyses and common garden experiments are required to fully assess the presence of deleterious variants and adaptive potential of the two species, considering the limited number of seedlings and saplings observed in the wild (Takahashi, personal observation), it is plausible that Daphnimorpha populations face the risk of further decline in the future.The estimated genetic diversity and contemporary effective population size were higher in D. kudoi than in D. capitellata (Table 1).However, considering the high extinction risk of alpine plants in the face of climate change (Grabherr et al. 2010) and their restricted distribution, both D. capitellata and D. kudoi require conservation attention.Our results showed a clear genetic structure within D. kudoi (Mt.Kuromi vs. Mt.Nagata) but not within D. capitellata.Although the exact levels of genetic distinctiveness required for the definition of ESUs are debatable (Moritz 1984), considering their genetic differentiation, geographic isolation, and limited contemporary gene flow, we considered that in D. kudoi, the populations on Mt.Nagata and Mt.Kuromi should be regarded as different ESUs.Therefore, these three ESUs should be managed separately to conserve Daphnimorpha species.
In threatened species with fragmented distributions, a low level of connectivity among sub-populations can increase the risk of extinction (Ellstrand 1992).The results of the BA3-SNPS analysis showed that in all three ESUs, the gene flow among the subpopulations was relatively low (Tables S4  and S5).Although we should be cautious that we could not analyse all wild individuals and our data consisted of a limited number of SNPs, the result implies that contemporary gene flow would be limited among most subpopulations.In both species, the subpopulations are discontinuously distributed and several subpopulations consisted of fewer than 10 individuals.Such small and isolated subpopulations often result in a loss of genetic diversity and a subsequent reduction in offspring fitness (Oostermeijer et al. 1994).As the loss of small subpopulations may trigger a further decrease in connectivity among subpopulations and all three ESUs showed small contemporary effective population sizes, we emphasise the conservation of both large and small Daphnimorpha subpopulations to preserve the overall genetic diversity of the species.For this purpose, monitoring of fruiting rates and seedling establishment of subpopulations is needed.If necessary, hand cross-pollination and/or establishment of artificial populations, which can allow connecting subpopulations, could further guarantee the long-term persistence of these endemic climate relicts.

Fig. 1 a
Fig. 1 a Location of two mountain ranges in the southern part of the Japanese Archipelago where the two Daphnimorpha species are distributed.Orange areas show granitic soil areas.b Sampling locations of D. capitellata subpopulations (C1-10) and results of clustering

Fig. 2
Fig. 2 Population branching model used for fastsimcoal2 analysis and estimated parameters with 95% CIs.N1-5 indicate effective population sizes in diploid individuals and T1-2 indicate divergence times (years)

Fig. 3
Fig. 3 Historical changes in effective population sizes of D. capitellata (a) and D. kudoi on Mt.Kuromi (b) and Mt.Nagata (c) estimated from Stairway Plot2 analysis.The x-axis shows the time before the

Table 1
Locality information and the genetic diversity of two Daphnimorpha species estimated from MIG-seq data N, number of samples used for MIG-seq analysis; P, the number of private SNPs; π, nucleotide diversity; H E , expected heterozygosity; F IS , fixation index.Nucleotide diversity was calculated from whole sites including monomorphic sites.Expected heterozygosity and fixation index were calculated from only variable sites.Voucher specimens were deposited in TUS(D.capitellata; DT171, and D. kudoi; DT3762)