Deer grazing drove an assemblage‐level evolution of plant dwarfism in an insular system

Plant dwarfism, a syndrome characterised by a significant reduction in plant height and organ size, is a widely observed pattern of stress‐tolerant life‐form evolution that results from local adaptation to harsh environmental conditions. The drivers of assemblage‐level dwarfism have primarily been attributed to abiotic factors, such as low temperature, aridity, poor soil fertility or frequent fires. While biotic factors such as grazing pressure from herbivores can contribute to the establishment of plant dwarfism, these factors have rarely been tested at assemblage levels. Focusing on a dwarf plant assemblage comprising over 80 taxa on a small continental island in Japan with a high deer density, we hypothesised that historical deer grazing could also be a factor contributing to the large‐scale convergent evolution of dwarfism. To test this hypothesis, we measured the size of 1908 individual plants of 40 taxa pairs, comprising both palatable and unpalatable pairs from the island and their counterpart taxa from neighbouring regions, and sought to assess which factors (i.e. low solar radiation, estimated divergence time, low nutrient conditions and grazing pressure from deer) may have contributed to the formation of the dwarf plant assemblage on the island. We also performed genetic analysis to infer the time frames for the establishment of dwarf taxa. Statistical analyses revealed that plant size was significantly reduced mainly among the palatable taxa growing on the island, with preferential grazing by deer being identified as the most significant factor influencing plant size. Furthermore, genetic analyses revealed that dwarf ecotypes may have evolved over tens of thousands of years. Synthesis: To the best of our knowledge, this study is the first to demonstrate that interactions with herbivores can shape the assemblage‐level convergence of plant dwarfism. These findings enhance our current understanding of the formation of plant functional diversity.

The drivers of assemblage-level dwarfism have primarily been attributed to abiotic factors, such as low temperature, aridity, poor soil fertility or frequent fires.
While biotic factors such as grazing pressure from herbivores can contribute to the establishment of plant dwarfism, these factors have rarely been tested at assemblage levels.
2. Focusing on a dwarf plant assemblage comprising over 80 taxa on a small continental island in Japan with a high deer density, we hypothesised that historical deer grazing could also be a factor contributing to the large-scale convergent evolution of dwarfism.
3. To test this hypothesis, we measured the size of 1908 individual plants of 40 taxa pairs, comprising both palatable and unpalatable pairs from the island and their counterpart taxa from neighbouring regions, and sought to assess which factors (i.e.low solar radiation, estimated divergence time, low nutrient conditions and grazing pressure from deer) may have contributed to the formation of the dwarf plant assemblage on the island.We also performed genetic analysis to infer the time frames for the establishment of dwarf taxa.
4. Statistical analyses revealed that plant size was significantly reduced mainly among the palatable taxa growing on the island, with preferential grazing by deer being identified as the most significant factor influencing plant size.Furthermore, genetic analyses revealed that dwarf ecotypes may have evolved over tens of thousands of years.

Synthesis:
To the best of our knowledge, this study is the first to demonstrate that interactions with herbivores can shape the assemblage-level convergence of plant dwarfism.These findings enhance our current understanding of the formation of plant functional diversity.

| INTRODUC TI ON
Plant dwarfism, a syndrome characterised by a marked reduction in plant height and organ size, is a widespread phenomenon observed in diverse plant families and regions world-wide, particularly among alpine and Arctic plant communities (Billings & Mooney, 1968).Given that a reduction in plant size can contribute to reductions in reproductivity and competitive ability (Mitchell-Olds, 1987;Schmid & Weiner, 1993), dwarfism is considered a typical example of stress-tolerant life-form evolution that occurs as a consequence of local adaptation to harsh environments.Cold stress (Körner, 2003), strong winds (Combrinck et al., 2020), water current-associated disturbance (Itoh et al., 2005), desiccation (Fonseca et al., 2000) and frequent fires (Maurin et al., 2014;Simon & Pennington, 2012) are major abiotic factors that can potentially induce plant dwarfism, as relatively smaller plants can minimise damage provoked by environmental stressors and reduce evaporation (Wright et al., 2017).Moreover, a shorter stature could enable plants to evade the grazing activities of vertebrate herbivores and reduce the fitness costs of the associated losses (Carman & Briske, 1985).As a biotic factor, grazing pressure can promote the evolution of dwarf plants in diverse vegetation types (Fahnestock & Detling, 2000;Skaien & Arcese, 2020;Suzuki, 2008;Yin et al., 2020).Consequently, although multiple factors, both biotic and abiotic, can drive plant dwarfism (Wright et al., 2017), most previous studies have tended to focus on evaluating the influence of single factors in promoting a reduction in plant size, whereas their relative contributions have rarely been assessed (Quiroga et al., 2010).
Convergent evolution is the process whereby phylogenetically distinct lineages independently evolve similar phenotypes in response to similar environmental challenges (James et al., 2023).
The significance of convergent evolution in studies of natural selection is well-established, as parallel phenotypic evolution cannot be attributed solely to random processes such as genetic drift (Schluter et al., 2004).In particular, assemblage-level functional convergence indicates that strong selective forces have influenced assemblages and has thus served as a model for investigating relationships between selective agents and phenotypes (Rajakaruna et al., 2014).In this context, a notable plant assemblage, comprising more than 80 plant taxa that have convergently evolved as dwarf forms, is known from Yakushima Island (Hatusima, 1991), a mid-latitude continental island located 70 km south of Kyushu Island in the Japanese Archipelago (Figure 1).Although this island is relatively small (covering an area of approximately 500 km 2 ), it has mountains reaching elevations of up to 1937 m and annually receives large amounts of precipitation (>8000 mm per year).In the high-elevation area of this island (>1600 m), which is characterised by Pseudosasa owatarii grassland, wetlands and rocky granite slopes, many of the herbaceous species (within more than 60 genera of 37 families) are of a distinctively smaller size compared with normal types growing elsewhere in the archipelago (Hatusima, 1991;Masamune, 1934;Figure 2).Due to their remarkable difference in plant sizes, dwarf individuals of the island are recognised as distinct species, varieties or forms in many species (Yahara et al., 1987).Furthermore, a majority of the dwarf plants are perennial and they are renowned as ornamental plants among Japanese horticulturalists.Notably, the size of these plants remains virtually unaltered by repeated cultivation (Minehana Society, 2013), which has also been confirmed by the findings of common garden studies (Ishikawa et al., 2006;Shinohara, 2015;Tsukaya, 2005).Moreover, Fujishima et al. (1990) have reported that F 1 individuals obtained from a cross between a dwarf individuals of Ranunculus yaegatakensis and its counterpart (normal types in the mainland) were characterised by an intermediate morphology with respect to leaf size.In addition, Shinohara and Murakami (2006) have observed reductions in the sizes of cells in the leaves of dwarf species, which rarely occurs as a consequence of environmental stress (Körner et al., 1989).Collectively, these findings thus tend to indicate that the dwarf forms of these Yakushima plants might be genetically fixed, as opposed to plastic responses to environmental factors.Given the lack of any distinct karyological differences (Yamamoto et al., 2009) and the low genetic differentiation (Ishikawa et al., 2006;Shinohara & Murakami, 2006) between dwarf taxa and their counterparts, it can be speculated that the parallel evolution of dwarf forms may have occurred over a comparatively short time scale.

The assemblage-level dwarfism of plants growing on Yakushima
Island has been attracting the attention of biologists for a number of years, during which time, several ecological factors have been proposed to account for this phenomenon.Among these, a candidate is the poor nutrient status of soil in the high-elevation area associated with the parent granite rock, combined with the large volumes of precipitation falling (Hatusima, 1991;Makino, 1913).
Similarly, limited solar radiation due to the extended periods of cloudy and foggy weather and the short growing season limited by the harsh alpine-like environment have also been proposed as contributory factors (Hatusima, 1991;Yumoto et al., 2018).
Moreover, Mitsuta and Nagamasu (1984) have suggested that some of these dwarf taxa represent early divergent lineages, and that their dwarf forms are actually ancestral types of the respective groups.However, although several factors have been proposed to explain the mechanisms underlying the formation of dwarf flora, their influence is yet to be systematically investigated, and thus at present, the ecological factors contributing to parallel dwarfism remain unclear.In this study, we propose a novel alternative hypothesis, which posits that the dwarf assemblages found on Yakushima Island have evolved a short stature to evade the grazing pressure of a subspecies of the Japanese sika deer (Cervus nippon yakushimae), which is endemic to Yakushima and a small neighbouring island (Izawa et al., 1996).In the high-elevation area of Yakushima Island, this deer occurs at relatively higher densities (14.1-49.1 head/km 2 in 2008 or 2009) than deer populations on mainland Japan (0.8 head/ km 2 in 1989; Japanese Ministry of the Environment, 2022; Koda et al., 2009;Koda & Kawamura, 2012).Previous genetic studies (Goodman et al., 2001;Terada et al., 2021) have reported that C. nippon yakushimae in Yakushima Island exhibits relatively high genetic diversity and a large effective population size, with no indication of past population shrinkage.In addition, deer in Yakushima Island show strong genetic differentiation from mainland or surrounding insular populations (Terada & Saitoh, 2018).
These findings imply that the high densities of deer in Yakushima Island have likely been sustained at least since the end of Last Glacial (LG) period (about 11,500 years ago), when a land bridge connecting the island and the mainland collapsed (Kimura, 1996).
Accordingly, we speculate that plants in the high-elevation area of Yakushima Island would have been historically exposed to a high grazing pressure from these deer.Furthermore, in addition to dwarfism, we have found that several species of plants growing in the high-elevation area of the island, such as Berberis thunbergii (Berberidaceae) and Cirsium yakusimense (Asteraceae), are characterised by spines that are sharper and longer than those found in other populations or related species on the mainland.In addition, it is notable that compared with palatable species, plant species considered unpalatable to deer appear to have undergone considerably less miniaturisation on the island.On the basis of this evidence, we accordingly hypothesised that the dwarfism of taxa growing in the high-elevation area of the island represents an adaptation of this flora that has evolved to evade the grazing pressure imposed by deer.
Accordingly, in this study, we sought to determine the mechanisms associated with the formation of the endemic dwarf flora on Yakushima Island, focusing on the interaction of deer with the flora.To this end, we conducted a multi-species comparison using 40 taxa pairs, including both palatable and unpalatable taxa pairs from Yakushima and their counterpart taxa in neighbouring regions.To evaluate the influence of the candidate ecological factors hypothesised to promote dwarfism, we statistically modelled variations in plant size by considering climatic factors, a proxy of soil nutrient conditions, and deer preference.We also performed genetic analyses to gain an estimate of the time frames over which the dwarf taxa have become established and their demography.
We considered that this multi-dimensional approach would facilitate the identification of the ecological drivers of the assemblagelevel dwarf evolution that has occurred within the high-elevation area of the island, and thereby enable us to propose the evolutionary implications of the mechanisms whereby local endemics emerge.

| Study species and morphological measurements
From 2020 to 2022, we collected morphological data of 40 plant taxa from several locations in the high-elevation area of Yakushima Island to account for local variations in plant sizes (Figure 1).We refer to these plants as the Yakushima alpine taxa.For comparison, in 2021, we also collected morphological data from populations of the 40 respective probable closest related taxa or at least comparable taxa (hereafter referred to as the counterpart taxa) from mountainous areas in the adjacent main islands of Kyushu, Shikoku and Honshu, or low-elevation areas of Yakushima Island.The counterpart taxon for each Yakushima alpine taxon was selected based on taxonomic treatments according to Yahara et al. (1987), a taxonomic database of the Japanese plants (YList: Yonekura & Kajita, 2003) or the results of recent phylogenetic studies (Fuse & Tamura, 2016;Higashi & Setoguchi, 2014).Among the 40 taxa pairs, two taxa pairs are annual (Eriocaulon and Melampyrum), one is biannual (Gentiana) and the remaining pairs are perennial species.We categorised each taxa pair as either preferable or non-preferable to deer, with reference to lists of plants considered palatable and unpalatable to Japanese sika deer (Hashimoto & Fujiki, 2014;Kyusyu Regional Forest Office, 2011).
These lists are based on previously reported field observations of deer grazing and the results of feeding experiments using deer in the Japanese Archipelago.Consequently, 30 and 10 pairs were categorised as preferable and non-preferable, respectively.Detailed information pertaining to these samples is presented in Table S1.The field works were carried out with the permissions of the Japanese Ministry of the Environment, Agency for Cultural Affairs, Forestry Agency, Kagoshima Prefecture, Miyazaki Prefecture and Ehime Prefecture.
Morphological measurements were performed for all 40 taxa pairs.Among these, for 38 taxa pairs, we collected entire specimens of mature individuals bearing flowers (or spores in the case of ferns), which were pressed and dried for subsequent measurements of plant size.For the remaining two taxa pairs (Cirsium and Chloranthus), we measured plant height from living flowering individuals in the wild, on account of either their limited population size or their large vegetative size.For the majority of taxa, we measured plant height, defined as the length from the stem base to the plant apex.For taxa lacking above-ground stems, we measured the maximum leaf length as a proxy of plant size.We also measured individual specimens deposited in three herbaria (KAG, KYO and TUS).
In total, we obtained measurement data for 1908 individuals (946 individuals of the Yakushima alpine taxa, with an average of 23.7 individuals per taxon; and 962 individuals of the counterpart taxa, with an average of 24.1 individuals per taxon).To statistically assess the morphological differences between the Yakushima alpine and counterpart taxa, we conducted a t-test for each taxa pair, using the base package in R v.4.0.5 (R Core Team, 2021).Among the 40 taxa pairs, were used for subsequent nutrient measurements and molecular analyses (Table S1).The selection of this smaller subset of taxa pairs was based on the considerations of ploidy level (avoiding non-diploid taxa pairs), sample availability (removing taxa pairs with insufficient DNA material) and DNA quality (removing taxa pairs, for which DNA extracts contained PCR inhibitors).

| Measurement of leaf nitrogen concentrations
The high-elevation area of the island is characterised by the high levels of precipitation and granitic soils, with the resultant low soil fertility plausibly leading to reduced plant sizes.To investigate whether nutrient conditions play a role in the reduction of plant size, we used leaf nitrogen concentrations of the Yakushima alpine taxa as a proxy for the contents of soil nutrients.This selection was based on the fact that within lineages, the concentration of nitrogen in leaves has been established to be positively correlated with soil fertility (contents of soil phosphorus and nitrogen) at both global and local scales (Andersen et al., 2012;Bowman et al., 1993;Ordoñez et al., 2009).
For each of the 26 assessed taxa pairs, the leaves of four to eight individuals of each taxon were randomly selected and used for nitrogen measurement.In total, we obtained measurements from 124 and 114 individuals of the Yakushima alpine and counterpart taxa, respectively.For the purposes of determining leaf nitrogen concentrations, we ground dried leaves to yield powders, the suspensions of which were then analysed using a Thermo Scientific™ FLASH 2000 analyser (Thermo Fisher Scientific, Waltham, MA, USA).Differences in leaf nitrogen concentrations between the two taxa in each pair were assessed using the Mann-Whitney U test in the base package of R.

| Molecular analyses and SNP calling
To gain an estimate of the time of divergence between the corresponding Yakushima alpine and counterpart taxa in the 26 assessed taxa pairs, we identified genome-wide single-nucleotide polymorphisms (SNPs) using the MIG-seq method (Suyama et al., 2022).
We randomly selected 10-18 individuals per taxon, and in total, 818 individuals were used for molecular analyses (414 and 404 individuals in the Yakushima alpine and counterpart taxa, respectively; Table S1).Total DNA was extracted using the CTAB method (Doyle & Doyle, 1987), and libraries for MIG-seq analysis were subsequently constructed following the protocol described by Suyama et al. (2022).For the purposes of sequencing, we used the Illumina MiSeq platform (Illumina, San Diego, CA, USA) in conjunction with a MiSeq Reagent Kit v.3 (150 cycle) (Illumina).The sequence data thus obtained have been deposited in the database of the National Center for Biotechnology Information (BioProject ID: PRJNA862386).
Forward and reverse files obtained for each sample were merged and analysed as single-end data following Suyama et al. (2022).The raw reads were filtered and trimmed using Trimmomatic v.0.32 software (Bolger et al., 2014) using the following commands: HEADCRAP:6, CROP:77, SLIDINGWINDOW:10:30 and MINLEN:51.For each taxa pair, we generated three different datasets that were separately de novo assembled: a taxa pair dataset, including individuals of both Yakushima alpine and counterpart taxa, and two datasets comprising only the respective individual taxa.The former dataset was used for estimating divergence time between the two taxa and the latter two datasets were used for estimating fluctuations in the effective population size of each taxon.These datasets were assembled using the Stacks v.2.41 pipeline (Rochette et al., 2019), with the minimum depth being set to 6 (−m 6), and default values being employed for other parameters.
For each taxa pair dataset, we extracted SNPs present in both taxa and shared by more than five individuals per taxon using the -p and -r options of the populations program in Stacks.For the singletaxon datasets, we extracted SNPs detected in at least half of the analysed individuals (−R = 0.5).Moreover, we calculated the average values of nucleotide diversity and the expected heterozygosity of each taxon using the population program in Stacks.The significance of differences between the Yakushima alpine and counterpart taxa was assessed using a permutation test implemented in the coin package of R (Hothorn et al., 2008).

| Demographic inference
For demographic analyses, we calculated two-and one-dimensional folded site frequency spectra (SFS).To eliminate paralogous loci using populations software, prior to estimating SFS, we excluded loci with an observed heterozygosity of >0.6 for all datasets.In addition, SNPs characterised by high lineage disequilibrium values (R 2 > 0.4) and those deviating from a Hardy-Weinberg equilibrium (p < 0.01) were excluded using Plink v.1.9(Chang et al., 2015).SFS was estimated using easySFS (https:// github.com/ isaac overc ast/ easySFS).To extract as many SNPs as possible for estimating SFS, we employed a down-projection method.For taxa pair datasets, a two-dimensional folded SFS was projected down to 10 gene copies in each taxon.One exception in this regard was Lycopodium clavatum var.nipponicum in the counterpart population, for which we obtained data from only 10 individuals, and thus we set a downprojection value of 8 for this population.For taxon datasets, the projection value of each taxon was set to half of the analysed number of individuals.
For each taxa pair, a simple branching model for two populations (Figure S1) was fitted to the two-dimensional folded SFS.Simulations and parameter estimations were conducted using fastsimcoal2.7 (Excoffier et al., 2021).Given that mutation rates appear to increase with increasing genome size (Lynch, 2010), we employed a taxonspecific mutation rate calculated from the respective genome sizes, which were inferred from the C-values of study species or related species derived from the Plant DNA C-values Database (https:// cvalu es.scien ce.kew.org/ , accessed on 2022/1/30).We calculated mutation rates using the following regression equation described by Lynch (2010): log 10 (m) = −0.81+ 0.68log 10 (G), where m denotes the mutation rate (/10 9 bp/generation) and G denotes the genome size in megabases.Furthermore, the generation times of each taxa pair were determined with reference to published articles or cultivation information available from the internet.We adopted the same prior intervals of parameters for all taxa pairs (Table S2).For each model, we initially ran 20 independent iterations with 100,000 coalescent simulations and 20 optimisation cycles and obtained point estimates of demographic parameters of the model based on the highest maximum composite likelihood.Subsequently, the 95% confidence intervals (CIs) of parameters were obtained using a parametric bootstrapping procedure, creating 100 pseudo-observed SFS using the point estimates of parameters and repeating the estimation of these datasets using the same simulation settings.We inferred the relatively recent demographic history of each taxon in all analysed taxa using Stairway plot v.2.1.1 (Liu & Fu, 2015).These inferences were based on 200 spectral iterations, with the mutation rates and generation times being the same as those used for the fastsimcoal2 analysis.All analytical settings are shown in Table S3.

| Statistical analyses of the morphological data
To gain an indication of which factor(s) (i.e.low solar radiation, high annual precipitation, low nutrient conditions and grazing pressure from deer) may have influenced dwarfism in the high-elevation area of the island, we adopted a model selection approach using a linear model of the variation in plant size between Yakushima alpine and counterpart taxa.The relative plant size of each taxa pair was calculated by dividing the mean plant height or leaf size of the Yakushima alpine taxon by that of its counterpart taxon.Furthermore, we logtransformed these values to ensure a normal distribution of the data and used these values as the response variable.As explanatory variables, we calculated the differences of annual precipitation, mean annual temperature and sunshine duration for the collection localities of the Yakushima alpine and the counterpart taxa as climatic factors.Values for these parameters were obtained from 1 km-mesh data of climatic variables of Japan, which were downloaded from the website of the Ministry of Land, Infrastructure, Transport and Tourism (https:// nlftp.mlit.go.jp/ ksj/ gml/ datal ist/ KsjTm plt-G02.html, accessed on 2022/10/20).To determine the effects of deer grazing pressure, we considered the preference of deer for the plants in each taxa pair (preferable or non-preferable) as the explanatory variable.Furthermore, using data obtained for the estimated divergence time and leaf nitrogen concentrations, we attempted to control for the effects of lineage-specific demography and soil nutrient conditions on plant dwarfism.The ratios of nutrient contents for each taxa pair were obtained by calculating the mean leaf nitrogen concentration in the Yakushima alpine taxon divided by that in the counterpart taxon.Furthermore, as it has been suggested that plant defence strategy can be also influenced by resource availability (Endara & Coley, 2011;Tomlinson et al., 2016), we included an interaction term between deer preference and the ratios of nutrient contents to the candidate models.In total, we analysed 80 candidate models based on the combination of explanatory variables, including a null model that included no explanatory variables.Given limitations regarding the availability of genetic and leaf nitrogen concentration data for certain taxa pairs, models were constructed using the data of a subset of 26 taxa pairs.For the purposes of model construction, we used the R package lme4 (Bates et al., 2015), with the best model being selected based on the Akaike information criterion (AIC) value of the model.The effects of explanatory variables in the best model were inferred based on an analysis of variance (ANOVA) implemented in the base package of R and the estimated coefficient value.
To determine the taxa pair-specific changes in plant size in the high-elevation area of the island, we fitted a linear mixed model using the plant size of all individuals (whole individual dataset).The leaf size or plant height of each individual was set as the response variable, and deer preference of taxa pairs (preferable or nonpreferable), locality (counterpart or Yakushima alpine taxa) and their interaction were used as the explanatory variables.In these models, we also included the random intercept and random slope of taxa pairs.The taxa pair-specific changes in plant size estimated using the models were visualised using the plot_model function in the R package sjPlot (Ludecke, 2021).
Finally, we performed the phylogenetic generalised least squares (PGLS) analysis, which can assess the influences of the candidate factors taking into account the phylogenetic relationships between taxa (Symonds & Blomberg, 2014).Initially, we determined correlations between the relative plant size and deer preference in the original dataset of 40 taxa pairs.Subsequently, we performed model selection using the dataset comprising 26 taxa pairs, consistent with the LM analysis, evaluating a total of 80 candidate models.A phylogenetic tree including 40 taxa was generated based on megatree, comprising 31,749 angiosperm species (Zanne et al., 2014), using phylomatic software (Webb & Donoghue, 2005).PGLS analysis was performed using the pgls function in the R package caper (Orme et al., 2013).

| Morphological variation and nutrient status of the Yakushima alpine taxa
Our morphological measurements indicated that among the 40 taxa pairs, 33 pairs where characterised by significant difference between the Yakushima alpine and counterpart taxa with respect to plant size (Figure 3 and Table S1).A majority (28 pairs) of these 33 pairs were preferable taxa pairs, for each of which, Yakushima alpine taxa were characterised by a smaller plant size.Regardless of the mean size of the corresponding counterpart taxa, the mean size of preferable plants among the Yakushima alpine taxa was approximately 20-50 mm.For the preferable taxa, the ratio of the size of the Yakushima alpine taxa to the counterpart taxa was 0.39 on average, and the corresponding ratio for the non-preferred taxa was 0.97.The measured leaf nitrogen concentration within four of the 26 taxa pairs indicated that the Yakushima alpine taxa were characterised by higher mean values than their counterpart taxa (Figure S2).
Moreover, we obtained mean relative leaf nitrogen concentration of 1.02 and 0.91 for preferable and non-preferable taxa pairs, respectively (Table S4).

| Estimated divergence times and population size fluctuations
In 10 of the 26 analysed taxa pairs, the Yakushima alpine taxa showed significantly higher genetic diversity, in terms of both nucleotide diversity and expected heterozygosity, with only six pairs showing an opposite tendency, whereas we detected no significant differences with respect to the diversities of the remaining 10 pairs (Table S5).
Our demographic analysis revealed that the 95% CI of divergence times for 19 taxa pairs overlapped with that of the LG period (11,500-116,000 years ago; Figure 4a and Table S6).Contrastingly, the 95% CI for one of the taxa pairs (Melampyrum) was included in the Holocene period (present-11,500 years ago), whereas those for the remaining six taxa pairs were traced to a period prior to the LG period.Especially, Scutellaria and Thalictrum were diverged at an older period (95% CI: 283,034 and 368,088 years ago,respectively).Furthermore, Stairway plot2 analysis revealed that both the Yakushima alpine and counterpart taxa have undergone similar demographic histories, with population sizes increasing or remaining stable during the LG period, with subsequent gradual or rapid declines (Figure 4b and Figure S3).

| Linear mixed model analyses
Our linear model analyses based on the 26 taxa pair datasets indicated that models including deer preference had lower AIC values (Table S7).The best model included deer preference and divergence time (AIC = 40.3),with the second-best model including deer preference and differences in annual precipitation, mean annual temperature and annual duration of sunshine (AIC = 40.9).Applying the best model, ANOVA revealed that deer preference provided a significant explanation of the measured variation in plant size (p < 0.001), whereas the effect of divergence time was not significant (p = 0.076) (Table S8).The estimated coefficient of deer preference (nonpreferable) was −0.88 (p < 0.001), whereas that of divergence time was −1.4E-06 (p = 0.076; Table 1).These results would thus tend to indicate that among the assessed potential contributory factors, deer preference was closely associated with variation in the size of plants growing in the high-elevation area.Linear mixed model analysis of the individual-based datasets revealed that all assessed terms (deer preference, locality and their interaction) could explain the observed variation (Table S8).The estimated coefficient of the interaction term was significant (p = 0.005), whereas those obtained for deer preference (p = 0.722) and locality (p = 0.247) were nonsignificant (Table 1).In addition, the taxa pair-specific random effect revealed that for the preferable taxa, there was a positive correlation F I G U R E 3 Boxplots of the measured plant sizes for each taxa pair.Red and blue boxplots indicate the Yakushima alpine and counterpart taxa, respectively.A dashed vertical line indicates a boundary size of deer grazing as implied by Suzuki (2008).An asterisk adjacent to the taxa pair name denotes a non-preferable taxa pair.
between the random intercept and the slope of the taxa pair term (Figure S4).These findings thus indicate that Yakushima alpine taxa paired with the preferable large-sized counterpart taxa tend to be characterised by an appreciable degree of dwarfism.Moreover, the PGLS assessment revealed that the detected levels of dwarfism were significantly correlated with the preference of deer (r = −0.677,p < 0.001; Figure 5).Among the 80 models, the model including deer preference and divergence time showed the lowest AIC value with both variables significantly explaining the variation in plant size (p = 0.001 for deer preference, and p = 0.034 for divergence time; Tables S9 and S10).The estimated coefficients of deer preference (non-preferable) and of divergence time were −0.94 (p = 0.001) and −1.6E-06 (p = 0.034), respectively (Table 1).

| DISCUSS ION
The replicated independent evolution of multiple similar phenotypes provides evidence to indicate that evolutionary change is driven by a deterministic process, and thereby contributes to TA B L E 1 Estimated values of coefficients within the best models of LM and phylogenetic generalised least squares (PGLS) analyses.determining how phenotypes interact with and respond to the surrounding environments (James et al., 2023).In this context, previous studies have reported that herbivore grazing pressure can promote the evolution of plant dwarfism, although most of the relevant studies in this regard have tended to focus on graminoids with similar life forms in grassland ecosystems, in which livestock has been grazing for hundreds or even thousands of years (Carman & Briske, 1985;Fahnestock & Detling, 2000;Yin et al., 2020).
Contrastingly, in the present study, we sought to assess assemblage-level dwarfism in a community comprising plants from approximately 40 families growing in a natural insular ecosystem (Hatusima, 1991), the finding of which indicated that significant reduction in plant size has occurred in taxa deemed preferable to sika deer (Figure 3).Moreover, we found that regardless of taking into consideration abiotic factors (Table 1, Tables S7 and S8) and phylogenetic relationships (Figure 5, Table 1, Tables S9, and S10), deer preference consistently explained the observed variation in plant size.On the basis of the statistical evidence obtained, it is reasonable to assume that as a consequence of intense grazing pressure, plants growing in the high-elevation area of Yakushima Island have undergone an evolutionary transition to dwarf forms, primarily to enable them to evade the detrimental effects of deer grazing.A small plant size or creeping pattern of growth, as seen in most of the dwarf taxa colonising the high-elevation area in the island (Figure 2), would presumably make these plants less susceptible to the grazing habits of deer and/or contribute to reducing biomass loss if grazed upon (Kotanen & Bergelson, 2000;Martin et al., 2015).Interestingly in this regard, our measurements of plant sizes (Figure 3) and estimated taxa pair-specific random effects (Figure S4) tended to indicate that a majority of the preferable taxa have undergone reductions in size to similar mean sizes of less than 50 mm, which corresponds to a threshold level below which deer are unable to browse effectively, owing to anatomical limitations of the mouth structure and ingestive behaviour (Hamilton et al., 1995;Ozaki et al., 2007).The measured plant sizes in the high-elevation area of the island were almost consistent with the reported sizes (20-50 mm) of herbaceous plants adapted to high grazing pressure of C. nippon in other regions (Ishikawa et al., 2006;Suzuki, 2008).It is notable, however, that in dwarf plants, the reduction in leaf numbers and leaf areas contributes to a decrease in total photosynthetic production, consequently leading to reduced reproductivity (Körner, 2021).
Moreover, smaller plant heights can reduce the competitive ability for light (Gaudet & Keddy, 1988) and decrease the ability to F I G U R E 5 Phylogenetic relationships of 40 taxa pairs and the relative plant size of the Yakushima alpine taxa.An asterisk adjacent to a taxa pair name and dark grey bar denote a non-preferable taxa pair.
disperse propagules (Sakaguchi et al., 2020).Therefore, the minimum sizes of Yakushima alpine taxa may reflect a trade-off between the avoidance of grazing and minimising the associated adverse effects of dwarfism.
Given our findings in this study, the question arises as to why such dwarf plant assemblages appear to have developed only on the small island of Yakushima.Intuitively, it would seem reasonable to assume that intense deer grazing within the closed environment of an insular system is a plausible contributory factor.In this regard, it is speculated that the endemic deer population inhabiting the island may have become geographically isolated at least 11,500 years ago when the glacial land bridge connecting the island to the mainland was submerged as a consequence of rising sea levels concomitant with the period of post-glacial climatic warming (Harunari, 2001).
Deer on the island could have maintained stable populations with a large effective population size (Goodman et al., 2001).The small body size of this deer on Yakushima Island is conceivably an outcome of prolonged food limitation, owing to high population densities that presumably developed in the absence of any controlling predatory species (Terada, 2012).Although deer generally forage selectively on preferable plant species, at high population densities when sources of preferential forage may become limited, they would of necessity switch to feeding on less preferable palatable plants (Takatsuki & Kajitani, 2019).Consequently, on this particular island, the deer population may have imposed a stronger selective pressure on local plant assemblages than they would perhaps normally do on mainland Japan.A further factor that may have facilitated the large-scale convergent evolution seen on this island is the isolation of plant populations from congeneric taxa.In general, sea straits not only act as strong barriers to the migration of plants with low dispersal ability but also limit the opportunities of migration for taxa with a capacity of long dispersal ability, particularly to small islands (MacArthur & Wilson, 2001).The findings of our genetic analysis indicated that a majority of the Yakushima alpine taxa may have diverged from the corresponding mainland lineage during the LG period (Figure 4a), implying that the evolution of the dwarf forms seen today has occurred over a period of tens of thousands of years.
Following the LG period, the submergence of the land bridge would have eliminated or reduced the opportunities for secondary contact between the Yakushima alpine taxa and their mainland counterparts.Moreover, similar to other alpine plants in the Japanese Archipelago (Ikeda, 2022), it is probable that the Yakushima alpine taxa have retreated to the high-elevation area in response to rising temperatures, whereas most counterpart taxa may have retreated to similar areas on the mainland, which are located at distances of at least 180 km from the island.Consequently, this geographical isolation may have effectively prevented an introgression of the alleles that are non-adaptive in a deer-abundant environment and promoted the evolution of dwarfism within the island, especially for taxa with low dispersal ability.
The long life cycles of Yakushima alpine taxa could have also contributed to the formation of dwarfism.Plant defence strategy against herbivores is closely linked to resource-based trade-offs between defence and growth or reproductive investment (Herms & Mattson, 1992;Skarpe & Hester, 2008).In annual species, as plants need to allocate all resources to reproduction (Cates & Orians, 1975), reduced plant sizes may directly impact their reproductive success.On the other hand, as perennial plants can allocate resources to growth and maintenance for survival by storing them to underground structures (Friedman, 2020), they may be more tolerant to resource limitations caused by plant size reduction compared with annual species.
We also speculated that demographic history may have contributed to the large-scale convergent evolution seen on Yakushima Island.In this regard, insular populations tend to experience founder effects and genetic drift that are typically associated with small isolated populations (Laurance, 2008).Notably, however, on the basis of our Stairway Plot2 analysis, we failed to detect any evident decline in the population size of the Yakushima alpine taxa, at least during the LG period (Figure 4b), and their current population genetic diversity was found to be comparable to that of the counterpart taxa (Table S5).The relatively high or moderate genetic diversity among Yakushima populations has similarly been reported in previous phylogeographic studies (Aoki et al., 2014;Kimura et al., 2014) and palaeontological evidence indicates that temperate forests may have become established on the island during the LG period (Tsukada, 1984).This evidence accordingly implies that Yakushima Island may have served as a glacial refugium within the southern region of mainland Japan, and that the Yakushima alpine taxa may have retained a stable population size within the island even during the glacial period.In addition, as compared with annual species, perennial species tend to exhibit less temporal variation in population size over generations (García et al., 2008).Given that, with the exception of strongly selected alleles, genetic drift can overcome the deterministic effects of natural selection (Lande, 1976), we speculate that the stability of the populations of Yakushima alpine taxa during the LG period and their long life cycles may have contributed to mitigating the adverse effects of genetic drift.Thus, overall, we believe that the large-scale convergent evolution of dwarf taxa in the high-elevation area of Yakushima Island may be a consequence of the combined effects of the strong selective pressure of herbivore grazing to a smaller plant size, isolation from their counterparts and population stability associated with their biogeographical history.
In the face of herbivores, plants have evolved various defencing strategies such as spinescence, low nutrient content, physical toughness and chemical contents to resist or avoid herbivory (Hanley et al., 2007;Skarpe & Hester, 2008;Tomlinson et al., 2016).In the high-elevation area of the island, although many palatable plants may avoid deer grazing by reducing their plant sizes, we should exercise caution as this study did not consider other defence strategies in our statistical modelling.Among

| CON CLUS IONS
Plant dwarfism represents one of the most common examples of convergent evolution, and assemblage-level dwarfism has been attributed to strong abiotic stressors such as low temperature, aridity, strong winds and frequent fires (Billings & Mooney, 1968;Maurin et al., 2014).While grazing pressure from herbivores is also recognised as a potential driver of evolution of dwarf forms (Skarpe & Hester, 2008), its contribution to assemblage-level evolution had been untested.To the best of our knowledge, this study is the first to indicate that dwarf plant assemblages may have evolved in response to herbivore grazing pressure.This large-scale convergence towards dwarf forms could be linked to the island's unique characteristics: remarkably high density of deer resulting from the predator-free environment and oceanic isolation from the mainland.Given the comparable convergence of plant size among species in approximately 40 plant families and their endemism in Yakushima Island, we believe that our findings will make a valuable contribution to enhance our evolutionary ecology, insular plant evolution, island syndrome, plant defence strategy, plant dwarfism, plant-herbivore interaction F I G U R E 1 Maps of southern Japan (a) and Yakushima Island (b) with sampling locations of the Yakushima alpine and counterpart taxa (in red and blue, respectively).The light blue area indicates the coastal line during the last glacial maximum (LGM; 21,000 years ago).

F
I G U R E 2 (a) Photographs of a wetland grazed by deer in the high-elevation area of Yakushima Island (a-1), dwarf plants (a-2: Solidago minutissima and a-3: Angelica longiradiata var.yakushimensis).(b) A photograph of representative taxa pairs: Yakushima alpine taxa (left) and its counterpart taxa (right).Dennstaedtia and Melampyrum are non-preferable taxa pairs, where the other taxa pairs are preferentially grazed by deer.

F
Results of genetic analyses based on genome-wide single-nucleotide polymorphisms detected using MIG-seq analysis.(a) Estimated divergence times between the Yakushima alpine and counterpart taxa (mean and 95% CI of the parameter).An asterisk adjacent to the taxa pair name denotes a non-preferable taxa pair.Grey areas indicate glacial periods according to Akiyama (2009).(b) Historical fluctuations in the effective population size of each taxon of the Yakushima alpine taxa (upper panel) and the counterpart taxa (lower panel).Each line represents a single taxon and the grey area indicates the last glacial period (11,500-116,000 years ago).
the 40 taxa pairs, Rubus and Cirsium pairs have spines.In Cirsium, the counterpart taxon (C.suffultum) not only has sharp spines (>10 mm) but is also considered to be non-preferable taxa to deer (Kyusyu Regional Forest Office, 2011).In Rubus, although the counterpart taxon (R. illecebrosus var.illecebrosus) has been eaten by deer regardless its spines (<5 mm; Hashimoto & Fujiki, 2014), our morphological measurement showed that Yakushima alpine taxon (var.yakusimensis) is not as miniaturised as other preferable taxa pairs (Figure3), implying the possibility that var.yakusimensis could avoid deer grazing by other strategy, especially by developing density and length of spines.Therefore, further investigation of plant defence strategies is required for more comprehensive understanding of plant evolution driven by herbivory on Yakushima Island.
current understanding of convergence evolution driven by plantherbivore interactions and formation mechanisms of plant diversity in insular systems.AUTH O R CO NTR I B UTI O N S Daiki Takahashi, Yoshihisa Suyama, Keitaro Fukushima, Hiroaki Setoguchi and Shota Sakaguchi designed the research.Daiki Takahashi and Shota Sakaguchi conducted the field survey and analysed the data.Daiki Takahashi, Yoshihisa Suyama and Keitaro Fukushima conducted the experiments.Manuscript writing was led by Daiki Takahashi with contributions from all authors.