Phenotypic differentiation of the Solidago virgaurea complex along an elevational gradient: Insights from a common garden experiment and population genetics

Abstract Plant species distributed along wide elevational or latitudinal gradients show phenotypic variation due to their heterogeneous habitats. This study investigated whether phenotypic variation in populations of the Solidago virgaurea complex along an elevational gradient is caused by genetic differentiation. A common garden experiment was based on seeds collected from nine populations of the S. virgaurea complex growing at elevations from 1,597 m to 2,779 m a.s.l. on Mt. Norikura in central Japan. Population genetic analyses with microsatellite markers were used to infer the genetic structure and levels of gene flow between populations. Leaf mass per area was lower, while leaf nitrogen and chlorophyll concentrations were greater for higher elevations at which seeds were originally collected. For reproductive traits, plants derived from higher elevations had larger flower heads on shorter stems and flowering started earlier. These elevational changes in morphology were consistent with the clines in the field, indicating that phenotypic variation along the elevational gradient would have been caused by genetic differentiation. However, population genetic analysis using 16 microsatellite loci suggested an extremely low level of genetic differentiation of neutral genes among the nine populations. Analysis of molecular variance also indicated that most genetic variation was partitioned into individuals within a population, and the genetic differentiation among the populations was not significant. This study suggests that genome regions responsible for adaptive traits may differ among the populations despite the existence of gene flow and that phenotypic variation of the S. virgaurea complex along the elevational gradient is maintained by strong selection pressure.

The phenotypic plasticity determines short-term morphological and physiological responses to environmental fluctuations (Bradshaw, 1965;Castillo et al., 2014;Sultan, 2000). However, different environmental selection pressures may lead to long-term evolutionary changes in phenotypes adapting to local environments (Linhart & Grant, 1996). As a result, adaptive traits are genetically fixed over generations by selection. Phenotypic variation due to genetic differentiation along environmental gradients often underlies formation of ecotypes or intraspecific taxa, which are at initial stages of speciation (Fukuda, 1974;Tateoka, 1983). Thus, study on geographic variation of phenotypes is important to understand mechanisms of evolution and intraspecific diversity (Endler, 1977). A powerful way to detect genetic-based phenotypic differentiation is common garden experiments, which have been applied to many plant species (Bertel, Buchner, Schonswetter, Frajman, & Neuner, 2016;Clausen, Keck, & Heisey, 1948;Kawano, 1974;Kruckeberg, 1967;Scheepens & Stöcklin, 2013;Vitasse et al., 2014).
Mountain ecosystems are ideal to study adaptive differentiation of plant species, because phenotypic variation is often found along elevational gradients, accompanied by drastic changes of environmental conditions at short geographic distances. Temperature and air pressure are lower at higher elevations, while wind velocity and ultraviolet radiation are greater at higher elevations (Friend & Woodward, 1990). Plant species respond to these elevational changes based on their morphological and physiological plasticity and through adaptation. For example, deciduous plant species maintain a positive carbon balance at the level of individual leaves at high elevations with a short growing season by decreasing the leaf construction cost and increasing the assimilative capacity (Kudo, 1996(Kudo, , 1999Oleksyn et al., 1998). Plant height also generally decreases with increasing elevation because of resource limitation due to severe environmental conditions, that is, short growing season, prolonged snow cover, strong wind, and shallow soil (Clausen et al., 1948;Körner, 2003;Mizuno, 1991;Natori, 1964;Takahashi & Yoshida, 2009). Furthermore, seed traits may change depending on elevation. The germination rate is high in high elevations for certain species (Vera, 1997), and germination and subsequent seedling survival rates are often greater for larger seeds (Jakobsson & Eriksson, 2000;Moles & Westoby, 2004). Therefore, sharp environmental changes along elevational gradients promote differentiation of plant traits to adapt to local environmental conditions. However, genetic differentiation along elevational gradients may be prevented by gene flow because gene flow homogenizes the genetic structure among populations (Lenormand, 2002;Slatkin, 1985). Gene flow is expected to occur along elevational gradients within a mountain because geographic distance is short between populations within a mountain (Matter, Kettle, Ghazoul, & Pluess, 2013).
However, even if the elevational range is narrow, differentiation of flowering phenology is caused by the difference in such parameters as timing of snow melting (Kudo, 2000). Nonsynchronous flowering timing among populations along an elevational gradient acts as a barrier against gene flow (Hirao, Kameyama, Ohara, Isagi, & Kudo, 2006). Furthermore, the gene flow in plants with entomophilous flowers along elevational gradients is affected by the activity of pollinating insects (Byars, Parsons, & Hoffmann, 2009). Thus, evaluation of gene flow and genetic structure behind geographic variation of phenotypes is important to improve our understanding of mechanisms of adaptation and evolution of plants.
The Solidago virgaurea L. complex (Asteraceae) is a perennial herb species and is widespread in northern Eurasia from the temperate to subarctic zones (Hayashi, 1978a;Makino, 1989). This species grows in various vegetation types from lowlands to the alpine zone (i.e., riverside, forest floor, alpine meadow) and shows large morphological variation. Previous studies reported phenotypic variation of the S. virgaurea complex along elevational gradients. For example, plant height, leaf mass per area (LMA), and number of flower heads per plant decrease with increasing elevation, while leaf nitrogen and chlorophyll concentrations, diameter of involucres and number of florets per flower head increase with elevation (Nishizawa, Kinoshita, Yakura, & Shimizu, 2001;Sakurai & Takahashi, 2017;Takahashi & Matsuki, 2017).
The S. virgaurea complex shows large morphological variation along elevational gradients and so the S. virgaurea complex in central Japan is classified into two subspecies, a lowland type of S. virgaurea L. subsp.
Although many previous studies focused on morphological variation of reproductive organs (e.g., shape of involucral scales) to identify the two subspecies in Japan, they did not show whether morphological variation of the S. virgaurea complex along elevational gradients was caused by genetic differentiation and adaptation or phenotypic plasticity (Hayashi, 1976(Hayashi, , 1978a(Hayashi, , 1978bKitamura, 1957;Takasu, 1975). The two subspecies of the S. virgaurea complex collected at nine elevations between 1,260 m and 2,670 m a.s.l. on Mt. Hakusan in central Japan show no differentiation, based on genetic analyses of fluorescence in situ hybridization (FISH) and random amplified polymorphic DNA (RAPD) (Nakamura, Miyamoto, Murata, Yamagishi, & Fukui, 1997). However, the RAPD method may not be suitable to analyze population genetic structure because most of RAPD bands show dominant inheritance and the reproducibility is low (Kagaya, 2005). The shape of involucral scales used to classify the two subspecies continuously changes along elevational gradients and is often intermediate between the two subspecies at the boundary between their distribution areas (Hayashi, 1976(Hayashi, , 1977(Hayashi, , 1978bNakamura et al., 1997). Therefore, the taxonomic treatment of the intraspecific taxa is still controversial. Sakaguchi and Ito (2014) advocated the necessity to examine the population genetic structure and gene flow in populations of the S. virgaurea complex.
In this study, we hypothesized that phenotypic variation of the S. virgaurea complex along an elevational gradient is caused by genetic differentiation. We made a seed germination experiment and a common garden experiment using seeds collected from nine populations of the S. virgaurea complex from 1,597 m to 2,779 m a.s.l. in central Japan to confirm the genetic basis of phenotypic variation, and we also analyzed the population genetic structure by using codominant microsatellite markers with high reproducibility and polymorphism. Specifically, we answer the following questions to clarify the hypothesis:  (Table 1). At each elevation, seeds were collected from 10 individuals, except for locality at 1,713 m a.s.l. where they were collected from only five individuals (pop_1700, Table 1). Seeds of each individual were put in an envelope and stored at 3°C until the germination experiment. Twenty seeds per individual were weighed, and then, the mean seed weight was calculated for each maternal plant.
Although the lowest sampling site was 1,597 m a.s.l. in this study, individuals at the lowest site were not hybrids between the two subspecies and were typical S. virgaurea subsp. asiatica (Nishizawa et al., 2001;Takahashi & Matsuki 2017). In addition, as a practical matter, it was difficult to find S. virgaurea subsp. asiatica in elevations lower than 1,597 m a.s.l. because of human disturbances. Therefore, we could not sample at elevations lower than 1,597 m a.s.l., but our samples included the two subspecies.

| Germination experiment
A germination experiment was done in April 2013. We used 10 seeds per maternal plant for the experiment (i.e., total of 100 seeds per elevation, except pop_1700 for which 50 seeds were used).
Seeds were placed on two filter papers in a petri dish, which were kept wet with de-ionized water. The petri dishes were placed in an incubator (Hitachi, CRB-41LA, Tokyo). Temperature and light conditions in the germination experiment were according to Kondo (1990): Temperature was set to 22.5°C throughout the experiment; the daylight duration was 14 hr. Germination was hardly observed after 8th day of the experiment, so the experiment was stopped at 17 day.
Germinated individuals were counted every day during the experiment. The germination rate per population was calculated from the number of germinated individuals until the end of the experiment.
Germinated individuals were used for the common garden experiment. The mean temperatures of August and January were 24.7 and −0.4°C, respectively. The annual mean precipitation was 1,031 mm, with most precipitation in summer. Although the elevation of our common garden was 650 m, it was ideal to set up common gardens at the same elevational range where samples were taken (1,597-2,779 m a.s.l.).

| Common garden experiment
However, as a practical problem, it was impossible to make such ideal common gardens, and we had to make a common garden in the university campus (650 m a.s.l.). However, the elevational differences in the phenological patterns of the S. virgaurea complex at the common garden reflected those in natural populations (see Results). Therefore, the common garden experiment at 650 m a.s.l. invalidates the results of this study.
The plants were watered once a day or once every 2 days, were fertilized once a week (N-P-K = 6-10-5, HYPONeX Japan Corp. Ltd., Osaka), and were relocated every week to reduce possible positioning effects during the first and second growing seasons. No bolting individuals appeared by the end of the first growing season (i.e., all plants remained as rosettes without flowers). In September 2013 (147th day after sowing), the rosette diameter in two perpendicular directions (one was the maximum width) was measured for each individual plant, and the above-ground part was harvested. The rosette area was estimated as an ellipse. Each rosette leaf was separated into lamina and petiole, and the laminae of each individual plant were scanned.
We measured the total lamina area of each plant using free software  137°33′00.14 dry mass divided by the lamina area. The above-ground dry mass of each individual plant in the first growing season was calculated as the sum of the total lamina dry mass and petiole dry mass.
All laminae of each individual plant were ground into powder to measure nitrogen and chlorophyll concentrations. Leaf nitrogen concentration (%) was measured using FLASH2000 (Thermo Fisher Scientific Inc., Waltham). Leaf chlorophyll was extracted using dimethylformamide (4 ml). The absorbance of samples extracted from leaf samples at 663.8 and 646.8 nm was measured by using a spectrophotometer (UVmini-1240, Shimadzu, Kyoto) and was substituted into Porra's equations (Porra, Thompson, & Kriedemann, 1989)  Stem height was measured for all bolting individuals every week.
The above-ground part of each individual was harvested after the finish day of flowering. The above-ground part of each individual was separated into stem, leaf, and reproductive organs, and each organ was stored in a separate envelope. Each organ was oven-dried at 80°C for 48 hr and weighed. Seeds of the S. virgaurea complex tend to disperse from the plant by wind soon after seed maturation, and involucres also tend to drop from the plant after the release of seeds. The dispersed seeds and dropped involucres were not included in measuring plant biomass. Therefore, the sum of stem mass and leaf mass, except for seeds and involucres, was regarded as the above-ground biomass in this study. The below-ground parts (i.e., roots) were also dug up after the harvest of the above-ground part. Soil on roots was washed out.
The below-ground biomass was weighed after oven-drying at 80°C for 48 hr.

| Molecular analyses
In June 2014 (second growing season of the common garden experiment), several leaves were sampled from each plant in the common garden for molecular analysis. In September 2014, 11-18 leaves of the S. virgaurea complex were additionally sampled from each of the nine study populations in the field. The leaf sample of each population was stored in an envelope and was kept at 3°C until DNA extraction.

| Common garden experiment
A generalized linear mixed model (GLMM) was used to analyze morphological and physiological traits and flowering phenology in the common garden among the nine populations with different elevations of provenance sites at which seeds were originally collected. The germination rate per population and mean seed mass per maternal plant were regressed against the elevation of provenance sites by GLMM.
Maternal plants were treated as a random effect because the germination rate and seed mass may differ between maternal plants due to genetic differences, even at the same elevation. The GLMM was also used to analyze relationships of the elevation of provenance sites with morphological and physiological traits, that is, stem height, aboveground biomass, below-ground biomass, rosette area, ratio of belowground biomass to total biomass (above-and below-ground biomass),

| Population genetic analyses
To assess the population genetic structure and the gene flow among the nine populations of the S. virgaurea complex along an elevational gradient, we analyzed the genotyping data obtained with the 18 microsatellite markers. Deviation from the Hardy-Weinberg equilibrium was estimated within a population by using the probability test for Genepop 4.4.3 (Raymond & Rousset, 1995 (Peakall & Smouse, 2006). Allelic richness (A R : El Mousadik & Petit, 1996) per locus in each population was calculated using FSTAT 2.9.3.2 (Goudet, 1995), and then, the mean value and standard error were calculated. One-way ANOVA evaluated whether A R differs between the nine populations.
Genetic differentiation between populations was estimated, based on the coefficient of F ST . Negative F ST values were converted into zero in this study. Pairwise F ST values between each pair of populations were computed using Arlequin 3.5.2.2 (Excoffier, Laval, & Schneider, 2005). Analysis of molecular variance (AMOVA) was used to partition the genetic variance among populations using Arlequin 3.5.2.2 (Excoffier et al., 2005). Principal coordinate analysis (PCoA) was also conducted to detect genetic differentiation among the nine populations using GenAlEx 6.502 (Peakall & Smouse, 2006).

| Germination experiment
The germination rate showed a sigmoid curve against elevations of provenance sites at which the seeds were originally collected ( Figure 1a, Table 2, p < .001). Seed mass was also greater for seeds of higher elevations of provenance sites (Figure 1b, p < .01).

| Common garden experiment
In the first growing season, rosette area did not correlate with the elevation of provenance sites (Figure 2a, Table 2). However, F I G U R E 1 Relationships of elevation of provenance sites with (a) germination rate and (b) seed mass of the Solidago virgaurea complex. Table 2  In the second growing season, stem height and above-ground biomass significantly decreased with elevation of provenance sites (p < .001, Figure 3a,b, Table 2). Below-ground biomass did not correlate with elevation of provenance sites (Figure 3c), and so the proportion of below-ground biomass in the total biomass significantly increased with elevation of provenance sites (p < .001, Figure 3d). The   Therefore, these results indicate that the significant correlation between the elevational distance and pairwise F ST values may be caused by high correlation between geographic distance and elevational distance.

| Seed size and maternal effects
Seeds of studied plants of the S. virgaurea complex from higher elevations showed greater seed mass and germination rates than those from lower elevations. A trade-off relationship exists between the number of seeds and seed mass (Primack, 1978). Species with smaller seeds produce a larger number of seeds and disperse small seeds in a wide area, which contributes to avoidance of sibling competition (Augspurger, 1984;Howe & Richer, 1982). The germination rate  and seedling establishment are greater for species with larger seeds (Westoby, Leishman, & Lord, 1996). Some studies showed the increase of seed size at high elevations with severe environmental conditions for seedling establishment because of a large storage of carbohydrate (Baker, 1972;Mariko, Koizumi, Suzuki, & Furukawa, 1993;Moles et al., 2007;Pluess, Schütz, & Stöcklin, 2005). Thus, the strategy to increase seed germination and seedling establishment is thought to be advantageous at high elevations with a short growing season.
Possibly, the phenotypic variation measured in the common garden experiment partly reflected the environmental maternal effect because all seeds used were obtained from natural populations (Monty, Lebeau, Meerts, & Mahy, 2009). Growth environments of maternal plants influence phenotypes of their offsprings through morphological and physiological plasticities of maternal plants (Roach & Wulff, 1987).
For example, environmental conditions of zygotes often affect phenotypes of the sporophyte (Kirkpatrick & Lande, 1989;Schmid & Dolt, 1994). These environmental maternal effects are often transmitted to offsprings through seed mass or size. Growth environments of maternal plants also affect the seed size; the seed size positively correlates with the germination rate and the following plant growth (Schmid & Dolt, 1994;Weis, 1982). However, the environmental maternal effect is marked in the early stages of the life history only and generally decreases with progress of individual growth (Schmid & Dolt, 1994;Wulff & Bazzaz, 1992). Therefore, phenotypic variations measured in the common garden experiment over 2 years are suggested to be caused by genetic differentiation rather than by the environmental maternal effect.

| Ecological interpretation of phenotypic variation in the common garden experiment
In the first growing season of the common garden experiment, rosette area did not correlate with elevation of provenance sites. However, above-ground biomass and LMA decreased, and leaf nitrogen and F I G U R E 3 Relationships of elevation of provenance sites with (a) stem height, (b) above-ground biomass, (c) below-ground biomass, (d) below-ground allocation (%) (proportion of below-ground biomass to the total biomass), (e) number of flower heads per individual of the Solidago virgaurea complex in the second growing season of the common garden experiment. Table 2 shows the results of generalized linear mixed model. A regression line is not shown in the graph (c) because of no statistical significance (p > .05) chlorophyll concentrations increased with increasing elevation of provenance sites. These observed patterns corresponded to those related to populations of the S. virgaurea complex along an elevational gradient in the field (Takahashi & Matsuki, 2017;our unpublished data). Leaf nitrogen concentration positively correlates with the lightsaturated maximum photosynthetic rate (Reich, Walters, Ellsworth, & Uhl, 1994). The increase in chlorophyll a/b ratio indicates that there is more chlorophyll a having a higher light condensing capability than chlorophyll b that is the reaction center chlorophyll and is regarded as an adaptation to increase the photosynthesis rate in high light conditions (Peng, Wu, Xu, & Yang, 2012). The S. virgaurea complex is a deciduous perennial herb species, and the stem and leaves wither before winter. Therefore, the positive carbon balance at the level of individual leaves can be maintained by decreasing the leaf structural cost and increasing the maximum photosynthetic rate at high elevations with a shorter growing season (Kudo, 1996(Kudo, , 1999Oleksyn et al., 1998).
In the second growing season of the common garden experiment, the stem height and the above-ground biomass decreased with an increase in elevation of provenance sites. These variation patterns corresponded to those observed in the field (Takahashi & Matsuki, 2017) because of resource limitation due to severe environmental conditions, that is, short growing season, prolonged snow cover, strong wind, and shallow soil (Clausen et al., 1948;Körner, 2003;Mizuno, 1991;Natori, 1964;Takahashi & Yoshida, 2009). The proportion of below-ground biomass in the total biomass increased with elevation of provenance sites. Natori (1964) and Shibata, Kinoshita, and Arai (1975) also reported higher ratios of below-ground biomass to above-ground biomass in the S. virgaurea complex at higher elevations. Probably, the S. virgaurea complex reserves photosynthate in the below-ground part to grow soon after the start of the next growing season at high elevations with a short growing period.
Although some morphological traits, such as LMA and stem height, continuously changed with elevation of provenance sites in the common garden experiment, the morphology of involucre discontinuously changed at the boundary between the distribution areas of the two subspecies (1,950 m a.s.l.). Flowers of the S. virgaurea complex are entomophilous (Sakurai & Takahashi, 2017). However, the number of pollinators decreases at higher elevations with cooler conditions (Blionis, Halley, & Vokou, 2001;Maad, Armbruster, & Fenster, 2013).
Size of flower heads of the S. virgaurea complex is larger in higher elevations (Nishizawa et al., 2001;Takahashi & Matsuki, 2017). Large flower heads may increase the reproductive success by increasing the chance to attract pollinators (Brody & Mitchell, 1997).

| Elevational differences of flowering phenology
Each stage of flowering phenology started earlier for plants at higher elevations of provenance sites, indicating that flowering phenology differs along the elevational gradient. This result corresponds to previous studies (Sakurai & Takahashi, 2017;Shibata & Terauchi, 1990) in which also reported earlier flowering of the S. virgaurea complex at higher elevations in the field. Early flowering is adaptive to high elevations because plants must produce mature seeds by the end of the short growing period (Sandring, Riihimäki, Savolainen, & Agren, 2007).   of the S. virgaurea complex in the subalpine zone (1,597 to 2,779 m a.s.l.), suggesting that plants grown in the common garden attained the effective accumulated temperature necessary to flower earlier than natural populations in the subalpine zone because the common garden location was warmer than the subalpine zone. Therefore, the elevational differences in phenological patterns in the common garden reflected those in natural populations; nevertheless, there were elevational differences between the common garden and natural populations of the S. virgaurea complex.

| Local adaptation under existence of gene flow
Although neutral loci had only very low genetic differentiation among the nine populations, many morphological and physiological traits genetically differed among the nine populations at the common garden. Each population showed a high genetic diversity, suggesting that stochastic fluctuation of population size has not been dominant.
Populations would rather have been stable along the elevational gradient and connected genetically under strong gene flow, as promoted by pollen dispersal by various insects, such as hoverflies and butterflies, and seed dispersal by anemochory (Kawano, 1988;Sakurai & Takahashi, 2017 (Kubota et al., 2015).
Although flowering of the S. virgaurea complex starts approximately 2 weeks later at elevations below 400 m, flowering lasts more than 1 month at each elevation (Sakurai & Takahashi, 2017). Thus, gene flow is presumed to occur between neighboring elevations, following the stepping stone model of population structure (Kimura & Weiss, 1964

| CONCLUSION
In this study, S. virgaurea complex populations along an elevational gradient were shown to be linked by substantial gene flow between neighboring populations by molecular analysis using neutral microsatellite markers. However, many morphological and physiological traits and flowering phenology showed genetic differentiation along elevations of provenance sites in the common garden experiment. These results suggest that only the genome regions of adaptive traits may display differentiation due to strong selection pressures despite the existence of gene flow. Our findings provide an example of plant micro evolution that genetically maintains adaptive traits to their local environments, even in narrow geographical ranges and under gene flow that could homogenize local adaptation.
In our study area (Mt. Norikura), S. virgaurea subsp. asiatica and S. virgaurea subsp. leiocarpa are distributed below and above ca.
1,950 m a.s.l., respectively, based on morphological characteristics of flower heads (Nishizawa et al., 2001). However, the molecular analysis provided by study, using 16 microsatellite markers, did not support the differentiation of the S. virgaurea complex into the two subspecies. Likewise, Nakamura et al. (1997) reported no genetic differentiation between the two subspecies on Mt. Hakusan in central Japan, using the RAPD and FISH methods. Considering their result together with our findings, it is plausible that the two subspecies are not genetically differentiated from each other, at least in terms of neutral loci. Therefore, the widespread alpine taxon of S. virgaurea subsp.
leiocarpa would be an ecotypic entity that arose from ancestral species via ecological adaptation to alpine environments. However, more comprehensive genetic analysis of the two subspecies (or ecotypes) in multiple mountain ranges is needed to understand the origins of alpine subspecies of subsp. leiocarpa in Japan.