Ambiguous species boundaries: Hybridization and morphological variation in two closely related Rubus species along altitudinal gradients

Abstract Although hybridization frequently occurs among plant species, hybrid zones of divergent lineages formed at species boundaries are less common and may not be apparent in later generations of hybrids with more parental‐like phenotypes, as a consequence of backcrossing. To determine the effects of dispersal and selection on species boundaries, we compared clines in leaf traits and molecular hybrid index along two hybrid zones on Yakushima Island, Japan, in which a temperate (Rubus palmatus) and subtropical (Rubus grayanus) species of wild raspberry are found. Leaf sinus depth in the two hybrid zones had narrower clines at 600 m a.s.l. than the molecular hybrid index and common garden tests confirmed that some leaf traits, including leaf sinus depth that is a major trait used in species identification, are genetically divergent between these closely related species. The sharp transition in leaf phenotypic traits compared to molecular markers indicated divergent selection pressure on the hybrid zone structure. We suggest that species boundaries based on neutral molecular data may differ from those based on observed morphological traits.

Interactions between gene dispersal (gene flow) and selection determine the formation and clines of hybrid zones, and a number of hybrid zone conceptual models have been constructed, such as the tension, bounded hybrid superiority, and mosaic hybrid zone models, to describe likely ecological processes that operate in these zones. Tension zone models assume that hybrids have low fitness relative to parental species (Barton & Hewitt, 1989). Intrinsic selection against hybrids can lead to narrow clines of neutral molecular markers (Martinsen, Whitham, Turek, & Keim, 2001). Bounded hybrid superiority zone models assume greater levels of fitness of hybrids than either parent species in intermediate habitats and extrinsic hybrid selection (Moore, 1977). Hybrid superiority can lead to the exclusive occupation by hybrids and smooth clines of the hybrid zone (Hamilton, Lexer, & Aitken, 2013). Mosaic hybrid zone models assume hybrid zones occur in heterogeneous environments in which parental species co-occur (Harrison, 1986). However, the assignment of hybrid zones to the most appropriate model is complex; for example, a narrow tension zone was found to be evident for flower color and associated gene polymorphism in Antirrhinum, but not for other single nucleotide polymorphisms (Whibley et al., 2006).
Comparison of morphology and neutral genetic marker clines may illustrate a balance between selection and dispersal (Brennan et al., 2009;Gompert, Mandeville, & Buerkle, 2017). For example, steeper clines in phenotypic traits than molecular markers that may indicate divergent selection despite gene flow (Campbell, Faidiga, & Trujillo, 2018). The limited number of hybrid zones (Abbott, 2017) may reflect problematic identification of hybrid zones and hybrids, because only early generations of hybrids may be recognized based on phenotypic traits. Alternatively, plant hybrid zones may indeed be rare, as a consequence of reproductive isolation and lower fitness of hybrids (Baack, Melo, Rieseberg, & Ortiz-Barrientos, 2015). Observation of transitions in traits and dispersal along hybrid zones can facilitate the identification of selection processes against hybrids and increase understanding of species boundary dynamics.
(raspberry), which form a closely related group (Okada, Kikuchi, Hoshino, Kunitake, & Mimura, 2020) in the subgenus Idaeobatus, diverged approximately one million years ago (Mimura et al., 2014) and are predominantly distributed in temperate and subtropical climatic regions, respectively. However, effects of climate change have driven secondary contact between these species and there are areas in which hybrid zones have formed, such as on Yakushima Island, Japan (Mimura et al., 2014). Yakushima Island is small (504.9 km 2 ), but mountainous (1,936 m a.s.l.), where R. grayanus is generally distributed in the lowlands, while R. palmatus tends to be distributed in uplands. Temporal flowering peaks of the hybrid zone shifts along the altitudes on the island, starting in February in the lowlands (R. grayanus) followed by the uplands (R. palmatus) in April. A previous TA B L E 1 Study site location and climate conditions. AHM, annual heat-moisture index calculated as (MAT + 10)/(MAP/1,000); DD, degree-days; MAT, mean annual temperature (°C); MWMT, mean warmest month temperature (°C); MCMT, mean coldest month temperature (°C); TD, temperature difference between MWMT and MCMT, or continentality (°C); MAP, mean annual precipitation (mm). Climate PC1 (90.6%) and PC2 (8.0%) are the first and second principal components, respectively, obtained based on these climate data.  (Mimura et al., 2014).

Transect
The smooth clinal changes in population genetic structure found in the previous study (Mimura et al., 2014)

km
In this study, we aimed to determine effects of dispersal and selection on phenotypic and molecular trait variations in R. grayanus and R. palmatus along altitudinal gradients within their hybrid zone on Yakushima Island to elucidate species boundary processes. Specifically, we compared leaf phenotypic trait clines to the population hybrid index, which was measured using neutral genetic markers (Mimura et al., 2014), to test for selection effects on leaf traits that are used in species determination. We also evaluated trait differences among parental species and natural hybrids subjected to the same, controlled environmental pressure, using a common garden experiment.

| Study sites and Rubus populations
We recorded leaf traits from nine Rubus populations (gYK, ab2, ab4, ab5, pYK, sr1, sr2, sr3, sr4) distributed along two transects located on trails at Anbo (ab, gYK, and pYK) and Shiratani (sr) on Yakushima Island, from where we had previously studied (Mimura et al., 2014; Table 1, Figure 1a). We also used the molecular hybrid index from the ten populations, including the population ab3, that we had previously estimated in Mimura et al., (2014). Distance between populations at the lowest (R. grayanus) and highest altitudes (R. palmatus) along the transects was estimated using Google Earth (Table 1; Anbo: 43-1,100 m a.s.l., 8.72 km; Shiratani: 100-760 m a.s.l., 3.72 km). Data for ten climate variables related to temperature and precipitation at the study sites were obtained from ClimateAP (Wang, Wang, Innes, Seely, & Chen, 2019; Table 1); the variables were highly correlated, so we estimated principal components (PC) in R (R Core Team, 2016) to generate independent climate PC variables that were analyzed using linear regression, with 95% confidence intervals for significant coefficients.

| Leaf traits
New leaves tend to be produced in Rubus species continuously through the growing season, so at the beginning of August 2017, we collected three fully developed leaves from current-year branches of the nine study populations of R. grayanus and R. palmatus on the island (n = 151; average of 16.8 individuals per population). All leaves were scanned on the sampling day for measurement of leaf area, length and width of leaf blade, and sinus width (shortest width between sinuses of a blade) using ImageJ (Schneider, Rasband, & Eliceiri, 2012). The depth of sinus was calculated as width of sinus divided by width of blade, and leaf length:width ratio was calculated as leaf length divided by width ( Figure 1b). Leaf dry mass per area (LMA) was estimated from six leaf disks (10 mm in diameter) from the leaves. The mean values of replicated leaves per plant were used for the following data analyses. In addition to morphological leaf traits, we also measured chlorophyll and carotenoid contents in fresh leaves, where three leaf disks were collected from the leaf samples and placed in 1.5 ml of 80% aqueous acetone at 4°C overnight to extract chloroplasts. Then, absorbance of the colorless leaf disks was measured at 470, 646.8, and 663.2 nm using a spectrophotometer (UV-2450, Shimadzu Corp.), and concentrations of chlorophyll a, chlorophyll b, and carotenoids were estimated according to Wellburn (1994).

| Cline analysis and migration estimation
The leaf trait data and hybrid index were fitted to a series of one-dimensional geographic cline models (Gay, Crochet, Bell, & Lenormand, 2008;Szymura & Barton, 1986, 1991 using the Metropolis-Hastings Markov Chain Monte Carlo algorithm, implemented in the R package "hzar" (Derryberry, Derryberry, Maley, & Brumfield, 2014). For leaf trait data, individual measurements were transformed to normal distribution, when needed, and scaled to values between zero and one for each transect, before clines were fitted to the population means (Derryberry et al., 2014) in five models with contrasting exponential tails (none fitted, left or right tail only, mirror tails, both tails) for best model selection. A null model estimation of no cline was calculated against distance and altitude.
For molecular clines, we used molecular hybrid index that had been previously estimated in Mimura et al. (2014). They calculated population means of maximum likelihood estimates of hybrid index (Buerkle, 2005) for the ten populations (16 samples each population) on Yakushima Island based on the partial sequences of seventeen nuclear functional genes, which were all statistically neutral.
Calculation of the molecular hybrid index was fitted to the fifteen cline models and null model that represented possible combinations of three trait intervals (fixed to 0 and 1, observed values, estimated values) and the five different tails described above (Derryberry et al., 2014). Each model was run with chain length of 100,000 steps with 10,000 burn-in steps, and the best models were selected based on corrected Akaike information criteria (AICc).
Rates and directions of gene dispersal among populations along the hybrid zones were estimated using the sequence data (Mimura et al., 2014), on which the hybrid index was based, using Bayesian inference in MIGRATE-N (version 3.6.11; Beerli & Felsenstein, 2001).
Although flowering time in the hybrid zones may overlap between neighboring populations, it shifts with altitude, so movement of pollinators across the hybrid zones was assumed to be limited to among in-flower patches; thus, we assumed a stepping-stone migration modeling approach along each transect to allow migration rates to vary between neighboring populations. The Bayesian inference was performed with one single long run of static heated chains (temperatures: 1.0, 1.5, 3.0, and 10,000) and, after discarding the first 25,000 trees, we sampled and recorded every 100th of 1,000,000 steps in each of the ten replicates. Then, posterior probabilities of theta and mutation scaled immigration rates (m/μ) along the transects were estimated.  (Mimura et al., 2014). In 2016, the plants were transplanted into larger plastic pots (30 cm in diameter, 12.8 L), and in spring 2020, three fully developed fresh leaves per individual were collected and scanned on the sampling day for measurement of the same leaf traits using ImageJ (Schneider et al., 2012) and analysis of chlorophyll and carotenoid content, as described above (Section 2.1.1). A Wilcoxon rank sum test, with Bonferroni p-value adjustment, was performed to compare the leaf traits among the groups (two parental species and hybrids). Analysis of Variance (ANOVA) was also conducted for each trait with the three groups as a fixed effect and the hybrid zones as a random effect. To test correlation among all traits, we used pairwise Spearman rank correlation analysis.

| Transect climate and environmental characteristics
Transect distance increased linearly with altitude ( Figure 1c, Table 1), where the transect at Shiratani was steeper than at Anbo. Climate variables tended to be correlated (average of pairwise correlations, |r| = 0.88) and contributed to climate PC1 that explained 90.6% of total variation in climate conditions (Table S1). Greater PC1 values indicate cooler and wetter conditions and were linearly related to altitude. Climate PC2 explained 8.0% of the total variance in climate conditions (lower mean annual precipitation and greater difference in temperature between the mean warmest and coldest months).

| Leaf trait differentiation under controlled environmental conditions
There were differences in four phenotypic leaf traits between the parental species under controlled environmental conditions in the common garden experiment (Figure 2), where leaf width was greater

| Traits and hybrid index cline models
We selected traits for the cline analysis from the common garden experiment that differed between the parental species (LMA, length:width ratio, sinus depth, and carotenoid content) and were noncorrelated (r < 0.70 between trait pairs); the population hybrid index and field transect data for these leaf traits were fitted to the clinal models by distance and altitude. Estimates of the best models for distance (Figure 3a, b) and altitude (Figure 3c, d) were similar for the two transects; the null models were not selected for the tested traits or the hybrid index. There was a gradual and smooth transition in the hybrid index along the transects (red lines in Figure 3), whereas there were clear transitions in leaf morphology (sinus depth, leaf Length:width, and LMA). Cline centers of the hybrid index were located at 5.38 and 2.83 km along the Anbo and Shiratani transects, similar to those of the sinus depth, dry weight, and Length:width ratio leaf traits (Table 2). However, cline widths for these leaf traits were narrower than the molecular hybrid index in both hybrid zones and cline centers of carotenoid content shifted toward R. grayanus (Table 2). Although transect length differed between Anbo and Shiratani, the centers of the transition of the three leaf traits were similar and occurred around 600 m a.s.l. (Figure 3c, d); at this altitude, the models estimated greater likelihood of phenotypes with R. grayanus-like sinus depth.

| Estimation of migration rates between neighboring populations
Based on the assumption of stepping-stone migration, estimated migration rates along the hybrid zones were asymmetric at the altitudinal extremes, where migration rates between neighboring populations were greater from lower to higher altitudes at low altitudes and the reverse at higher altitudes ( Figure 4). Thus, the populations at c. 600 m a.s.l. on the two transects were dominated by immigrants.

| Hybrid zone structure and gene dispersal
The smooth clines in neutral molecular hybrid indices in the two hybrid zones in this study indicated a lack of strong dispersal barriers between the parental species. Survival and fitness of hybrids characterize hybrid zone structure; for example, lower fitness of post-F 1 generations was found to create a tension zone in Senecio (Brennan et al., 2009). In some cases of bounded hybrid superiority zones, there may be an absence of later generations of hybrids, but the presence of F 1 hybrids, as has been reported for Rhododendron (Milne, Terzioglu, & Abbott, 2003) and Utricularia (Kameyama, Toyama, & Ohara, 2005), or higher levels of fitness of later generations of hybrids in habitats intermediate to parental species and a smooth molecular cline, such as in Picea (Hamilton et al., 2013).
Previously, we found that the Rubus hybrid zones largely comprised backcrossed and later generations of hybrids, based on interspecific heterozygosity, indicating hybrid survival and reproduction (Mimura et al., 2014), and an artificial crossing experiment between the parental species produced viable seeds (Mimura et al. unpublished data). Although hybrid fitness was not directly evaluated, the smooth molecular cline indicates that strong intrinsic selection against hybrids is unlikely in the hybrid zones.  TA B L E 2 Parameters estimated by best-fitting cline models for molecular hybrid index and leaf trait clines with distance from populations at the lowest altitudes (R. grayanus populations). The estimated cline center, measured as the distance from the population at the lowest altitude (R. grayanus), and cline width (1/ maximum slope) were estimated in R using "hzar" (Derryberry et al., 2014); two log-likelihood units of the parameter estimates are shown in parentheses.
study, which was based on statistically neutral molecular sequences (Mimura et al., 2014), should reflect neutral dispersal via seeds or pollen among neighboring populations. We observed no apparent declines in current population density at the centers of the hybrid zones, although overall, the population density of R. palmatus was lower than R. grayanus, perhaps due to increased patchiness in the distribution of R. palmatus at higher altitudes. Asymmetric dispersal indicates that hybrid population expansion and colonization from the upper and lower extremes of the hybrid zones, in which parental populations initially acted as source populations, were followed by frequent backcrossing toward spatially closer parental species.

| Drastic changes in leaf traits at species boundaries
Phenotypic variation is driven by genetic and environmental pressure. Although the environmental conditions at the common garden experiment contrasted with those at the island study sites, they were standardized and controlled, to allow for analysis of genetic differences in Rubus phenotype and separation of phenotypic plasticity responses. There was genetic divergence between parental species leaf traits ( Figure 2) and the cline analysis estimated that leaf sinus depth of the hybrids at 600 m a.s.l. on the two transects was more similar to R. grayanus than R. palmatus (Figure 3). These findings were supported by the common garden experiment that showed similar estimated clinal patterns in phenotypic differences in the leaf traits. Thus, modeling of our field data indicated that leaf traits, including those used for species identification, were mostly the product of genetic differences, rather than phenotypic plasticity. It is possible, however, that although the common garden experiment indicated contrasting photosynthetic strategies between parental species, differences among populations in the field may have been overestimated, because phenotypic plastic responses of chlorophyll content to different light conditions have been reported (Gratani, Covone, & Larcher, 2006).
The narrower leaf trait clines than molecular clines may be due to neutral diffusion or selection. Sharp phenotypic clines under neutral evolution may occur at relatively early stages of secondary contact, with limited migration rates between species (Gompert & Buerkle, 2016). The coalescent analysis estimated that hybrid zones in this study were formed following secondary contact more than c. 8,000 years ago, and later generations of hybrids have been observed (Mimura et al., 2014), indicating that the narrower phenotype than molecular clines did not emerge purely from neutral diffusion.
The narrow tension zone of some leaf traits and smooth molecular clines in this study support the hypothesis that these traits were under extrinsic divergent or disruptive selection, despite the presence of gene flow. A sharp change in phenotypic traits compared to neutral markers has also been interpreted as a result of divergent or disruptive selection (Campbell et al., 2018;Scordato et al., 2017). It is possible that selection may act only on parts of genomic regions responsible for traits, while other regions may be unaffected in hybrid zones (Gompert et al., 2017;Whibley et al., 2006).
Altitudinal gradients may elicit strong environmental selection pressure in hybrid zones (Abbott & Brennan, 2014;Brennan et al., 2009), including from effects of variation in temperature, precipitation, and light availability (Körner, 2007

Anbo
Shiratani are yet to be clarified, but may be related to development of veins for support and water supply (Givnish & Kriebel, 2017), efficient photosynthesis in cool temperatures at the beginning of growing seasons (Baker-Brosh & Peet, 1997;Royer & Wilf, 2006), and side-effects as a consequence of bud-packing (Edwards, Spriggs, Chatelet, & Donoghue, 2016). A recent study suggested the direct advantage of leaf teeth to carbon gain in early spring (Zohner, Ramm, & Renner, 2019

| CON CLUS IONS
We observed genetic and morphological differences in leaf traits between two closely related, but ecologically distinct Rubus spp.
in two hybrid zones on Yakushima Island, Japan. We repeatedly observed that the leaf traits, some of which indicated genetic differences between the two species, had much narrower cline widths than the molecular hybrid index. The results indicate extrinsic divergent selection pressure on Rubus leaf morphology at the hybrid zones, although it was difficult to identify an environmental factor as a selective pressure. Our results show that hybrid traits may quickly resemble those of their parents in the absence of strong reproductive barriers, and species boundaries may be less obvious than expected from the use of morphological traits in plant identification. This could explain why a low number of hybrid zones have been observed to date (Abbott, 2017), despite frequent hybridization in plants (Mallet, 2005). Hybridization and formation of hybrid zones increase with ongoing climate change and human activities, and our results indicate that genes penetrate closely related species at the margins of their ranges, where typical morphological differences can be retained. Hybridization may lead to increased risk of extinction (Todesco et al., 2016), but also increased chances of adaptive introgression, such as reported in a recent genomic study of evolutionary rescue in the context of water pollution (Oziolor et al., 2019). Thus, careful observation at species boundaries is necessary to identify and predict the outcomes of secondary contacts of species, especially under ongoing climate change conditions.

ACK N OWLED G EM ENTS
We are grateful for field assistance by T. Inoue and H. Ikeda. We also thank anonymous three reviewers and the Associate Editor for the insightful comments that improved the manuscript. This work was partially supported by JSPS KAKENHI Grant Number JP17K07571 and 25840161 to MM.

CO N FLI C T O F I NTE R E S T S
None declared.