Leaf anatomy of two reciprocally non-monophyletic mountain plants (Heliosperma spp.): does heritable adaptation to divergent growing sites accompany the onset of speciation?

Evolution is driven by natural selection, favouring individuals adapted in phenotypic traits to the environmental conditions at their growing site. To shed light on ecological and (epi-) genetically based differentiation between Heliosperma pusillum and Heliosperma veselskyi, two reciprocally non-monophyletic, but morphologically and ecologically divergent species from the south-eastern Alps, we studied various leaf anatomical traits and investigated chloroplast ultrastructure in leaves of the two species grown either in their natural habitat or in a common garden. The alpine H. pusillum occurs in open, wet rock habitats, whereas its close relative H. veselskyi is restricted to dry, shady habitats below overhanging rocks in the montane belt. H. pusillum exhibited higher thickness of leaves and palisade layers as adjustments and/or adaptations to higher irradiance and a higher stomatal area index reflecting better water availability. Traits were adjusted plastically, but differed between species grown in a common garden, suggesting that the differentiation between the two species is not solely based on phenotypic plasticity but also has a genetic basis. Our study thus supports the hypothesis that differentiation between the highly interfertile species is likely driven by natural selection.


Introduction
Adaptation to environmental conditions can become manifest in various phenotypic traits, as natural selection favours individuals being able to successfully survive and reproduce in a particular environment (Coyne and Orr 2004). Habitat heterogeneity, which is particularly pronounced in mountain areas (Scherrer and Körner 2011), may trigger the formation of ecotypes within species. Ecotypes originate from differentiation of populations that are adapted in various phenotypic traits to a specific microenvironment (Hufford and Mazer 2003;Lowry 2012). Ecotypic differentiation is most likely driven by a combination of heritable and non-heritable traits (Pfennig et al. 2010;Bonduriansky et al. 2012). Despite the absence of intrinsic reproductive barriers, ecological isolation may over time lead to the formation of new species (Lowry 2012).
An example of a recent ecotypic differentiation is provided by Heliosperma pusillum and Heliosperma veselskyi (Caryophyllaceae, Fig. 1). The former has a broad distribution throughout the southern and central European mountain ranges and occurs in humid, partly sun-exposed rock crevices and screes in the upper montane to alpine zone (1700-2300 m a. s. l.; authors' personal observations). In contrast, H. veselskyi is restricted to a few scattered populations below cliff overhangs in the lower montane belt of the south-eastern Alps and northernmost Balkan Peninsula (500-1300 m a. s. l. ;Neumayer 1923;Frajman and Oxelman 2007). Its growing sites are usually characterised by low irradiance and limited water availability (Bertel et al. 2016). Due to morphological divergence, both entities have been described at the species rank and are still treated as independent species (e.g . Fischer 2008;Frajman and Oxelman 2007;Poldini 2002). H. pusillum is glabrous or has sparsely hairy leaves, whereas H. veselskyi is characterised by a dense indumentum of long multicellular glandular hairs (Janka 1858;Neumayer 1923;Frajman and Oxelman 2007;Fischer 2008). Chloroplast and nuclear low copy DNA sequence data (Frajman and Oxelman 2007;Frajman et al. 2009) as well as highly resolving restriction associated DNA markers sampled across the nuclear genome (Trucchi et al., unpublished) suggest that the two species are phylogenetically not distinct and that H. veselskyi is inextricably nested within H. pusillum. H. veselskyi thus rather represents a habitat specific ecotype, whose disjunct populations have evolved postglacially from geographically close populations of H. pusillum. For the sake of simplicity, we treat H. pusillum and H. veselskyi as species throughout the text in spite of the lack of consistent genetic divergence and their highly debatable taxonomic value.
Significant differences in environmental conditions between the habitats of H. pusillum and H. veselskyi comprise irradiation, temperature and soil water availability, and result in specific photosynthetic adjustments in both species (Bertel et al. 2016). In general, adjustments of photosynthesis are partly plastic, but also involve non-reversible traits influencing light absorption and carbon fixation at various organisational levels ranging from that of organelles within mesophyll cells over leaf anatomy up to the whole-plant level (Larcher 2003). Leaves of plants growing under low irradiance display a set of morphological, biochemical and physiological adjustments (Valladares and Niinemets 2008). They usually have a thinner mesophyll with no or only a few layers of palisade parenchyma with shorter cells and fewer chloroplasts per area with larger grana and more stromal thylakoids and fewer stomata per leaf surface area (Larcher 2003). Mesophyll conductance is often low contributing to a low overall photosynthetic capacity (Monti et al. 2009).
Low water availability often leads to the formation of more trichomes, highly waterproof cuticles and also influences stomatal density (Schulze et al. 2002). Fewer stomata on upper leaf surfaces prevent evapotranspiration, and more and shorter stomata on lower leaf surfaces allow a more precise regulation of gas exchange (Lösch 2001;Larcher 2003). As the formation of a particular leaf anatomy is both delimited by genetic constraints and the peculiar environmental growth conditions during development, investigations aiming at disentangling genetic and environmental effects of leaf structural traits ought to include plants grown at their natural growing sites as well as plants cultivated in a common garden under similar environmental conditions.
To shed light on ecological and (epi-)genetically based differentiation of leaf anatomy between H. pusillum and H. veselskyi, we address the following questions: (1) Does leaf anatomy differ between both species grown under the divergent ecological conditions of their natural habitats? A lack of differentiation would point to a weak ecological difference between habitats or minor relevance for the species, suggesting that plastic adjustment of physiological characteristics in the short term may be sufficient to compensate the divergent habitat conditions.
(2) In order to resolve the extent of phenotypic plasticity vs. (epi-)genetically based determination of leaf anatomical traits, we further asked if leaf anatomy differed between the two species grown from seeds under similar environmental conditions in a common garden. Absence of genetic differentiation between the species is expected to lead to similar leaf anatomy under identical growth conditions, whereas intrinsic genetic constraints would limit the plastic adjustment of leaf anatomy. In addition, we investigated leaf ultrastructure of both species from their natural growing sites and from common garden qualitatively by transmission electron microscopy in order to evaluate subcellular rearrangements (e.g. abundance of grana stacks), which mirror ecological and (epi-)genetic differentiation.

Study plants and experimental design
Leaf samples were both taken from plants in their natural habitat and from individuals grown from seeds in a common garden. Samples from natural growing sites originated from two geographically close and recently diverged (Trucchi et al., unpublished) populations of H. pusillum (Mt. Rudnigkofel, 46.76°N, 12.88°E, 2060 m a. s. l.) and H. veselskyi (Anetwände, 46.77°N, 12.90°E, 790 m a. s. l.) situated in the Lienzer Dolomiten, Austria. In autumn of 2013, seeds were collected from both populations and sown in pots filled with soil composed of compost, ground earth, lava, turf, quartz sand, rock flour and pumice gravel. Thirty individuals of each species were cultivated in a common garden in the Botanical Garden of the University of Innsbruck, Austria (47.27°N, 11.38°E, 600 m a. s. l.) under the same climatic conditions with adequate water supply. Environmental conditions differed between the natural growing sites of both species, as higher irradiance and temperature fluctuations were measured at the alpine site of H. pusillum. They also differed between the natural growing sites and the common garden, where high irradiance-comparable to the site of H. pusillum-and higher leaf temperatures-comparable to the site of H. veselskyi-were measured. Results of microclimatic measurements from the same experimental sites and the same vegetation period are presented in Bertel et al. (2016). Leaf samples were taken from adult plants. From here on, we use the following abbreviations for sample origin: HP G , H. pusillum grown in the common garden; HP N , H. pusillum collected in nature; HV G , H. veselskyi grown in the common garden; HV N , H. veselskyi collected in nature.

Leaf anatomical traits
For investigation of leaf anatomical traits, fully developed stem leaves from middle parts of the upright shoots were collected randomly from 12 well-developed and randomly chosen individuals of each of the four groups. These leaves had developed in the collecting season and were in a physiologically active state, permitting comparison among groups. Cross sections were obtained from the middle part of the leaf lamina using a hand-microtome and investigated with a microscope (Olympus BX50, Olympus Optical Co., Tokyo, Japan) in combination with the software cellD (Version 3.1., Olympus Optical Co., Tokyo, Japan). Several leaf anatomical traits were measured on one cross section per individual, i.e. 12 cross sections from 12 individuals per study group: total thickness of leaves (LT), thickness of palisade parenchyma (PT), thickness of spongy mesophyll (SMT), thickness of cuticle on adaxial (CT ad ) and abaxial (CT ab ) epidermis, thickness of epidermis on adaxial (ET ad ) and abaxial (ET ab ) surfaces. Trait values represent the mean of three measurements per cross section. Relative proportion of palisade parenchyma (PT rel ), spongy mesophyll (SMT rel ) and intercellular spaces (IC rel ) were estimated from five images per study group of leaf cross sections according to Kubinova (1993). Stomatal density and stomatal length on adaxial (SD ad , SL ad ) and abaxial (SD ab , SL ab ) leaf surface were measured from imprints. For imprints, nitrocellulose polish was applied on whole leaves as a thin film, and carefully removed after 5 min air drying with transparent adhesive tape, which was then placed on microscope slides. Stomatal density was determined from digital microscopic images of imprints as the mean number of stomata located within 30 randomly chosen grid cells of 100 × 100 μm in the middle section of leaves. Stomatal length was measured on digital images of imprints and is the mean of five measurements. Adaxial and abaxial stomata area indexes (SAI ad and SAI ab ) were calculated as the product of mean stomatal length and stomatal density (Ashton and Berlyn 1994;Gratani and Varone 2004;Gratani et al. 2014).

Transmission electron microscopy
Stem leaves of the middle section of vertical shoots were sampled on sunny days. Samples of both species were taken at the same day at natural sites and in the common garden with smallest possible temporal delay (HP N and HV N , 10.07.2013, between 13.00 and 15.30 h; HP G and HV G , 30.06.2014 at 11.30 h). Leaf samples were immediately fixed (2.5 % glutaraldehyde in 50 mM sodium cacodylate buffer, pH 7.0) for 2 h and rinsed in the same buffer. Samples were postfixed in buffered 1 % OsO 4 over night at 4°C and further processed as described by Holzinger et al. (2007). Ultrathin sections were investigated in a Zeiss Libra 120 TEM at 80 kV.

Statistical analyses
A principal component analysis (PCA) was applied to illustrate differentiation among the four study groups using the leaf anatomical traits LT, PT, SMT, ET ad , ET ab , CT, SAI ad and SAI ab standardised to zero mean and unit variance. As the Spearman correlation index was lower than 0.8 for any pair of characters, all characters were retained. Differences between study groups were compared by a one-factorial analysis of variance (ANOVA), conducted for each leaf anatomical trait separately. Significant differences between study groups were analysed by pairwise comparisons applying Bonferroni correction at a significance level of α = 0.05. Leaf anatomical traits grouped by species (thus pooling growing sites) were tested by a two-factorial ANOVA (p < 0.05); we only report but do not interpret the test comparing natural site and common garden pooling species as natural growing sites are strongly divergent. For all analyses, data transformations were conducted if normal distribution was not given: PT was transformed by common logarithm, IC rel by natural logarithm and SL ad to the power of four. For SMT, data were normally distributed, but errors were not homogenous; in this case Welch's test and pairwise Wilcoxon rank sum test (for comparison of study groups) were applied. All statistical analyses were computed in R 2.14.0 (R Core Team 2014). PCA was calculated using the functions dudi.pca (package ade4: Dray and Dufour 2007) and the package plotrix (Lemon 2006) was used for graphical representations.

Results
Anatomical characteristics of leaves of H. pusillum and H. veselskyi from natural growing sites as compared to the common garden are illustrated in Fig. 2. Study groups differed significantly in most leaf anatomical traits according to the one-factorial ANOVA (Table 1). Differences in leaf anatomy were most pronounced in traits related to leaf thickness (i.e. LT, PT and SMT), with highest values in HP G , intermediate values in HP N and HV G and lowest values in HV N . In congruence, a second row of palisade parenchyma (No. P) was formed in all leaves of HP G and in some leaves of HP N , but was never observed in HV N or HV G. Further differences concerned thickness of cuticles (CT), epidermal layers (ET) and stomatal characteristics (SD, SAI, SL ab ; boxplots for important traits are shown in Fig. 3, all traits are summarised in Table 1). Similarly, in the two-factorial ANOVA, the effect of species was significant for traits related to leaf thickness (LT, PT, PT rel ) and surface characteristics (ETad, CT, SD ad , SAI ad , summarised in Table 2).
The PCA showed a clear differentiation between HP N and HP G , HV N and HV G as well as between HP G and HV G (Fig. 4). The first three axes represented 76.2 % of variance (i.e. 47.9, 17.2 and 11.1 %, respectively). The traits PT, SAI ad , CT, SMT and LT correlated strongest with the first axis and SAI ab , ET ab and ET ad with the second axis.
The chloroplast ultrastructure is illustrated in Fig. 5. HP N chloroplasts were more roundish-ellipsoidal and more densely packed with starch ( Fig. 5a) than chloroplasts of HV N , which were lens-shaped containing only a few-mostly one or two-starch grains (Fig. 5b). Grana stacks contained more thylakoid membranes in HV N than in HP N . Differences in chloroplast ultrastructure between HP G and HV G were less pronounced, but in HV G plastoglobules were abundant in all examined chloroplasts, which were lacking or rarely observed in chloroplasts of all other study groups.

Discussion
Leaf anatomy of the two mountain plant species H. pusillum and H. veselskyi, which occupy strongly divergent habitats despite their inextricable genetic relationship (Frajman and Oxelman 2007;Frajman et al. 2009), is  Traits include total thickness of leaves (LT), number of palisade cell layers (No. P), thickness of palisade parenchyma (PT), thickness of spongy mesophyll (SMT), thickness of epidermis on adaxial (ET ad ) and abaxial (ET ab ) leaf surface, thickness of cuticle on adaxial (CT ad ) and abaxial (CT ab ) epidermis, stomatal density, length and stomatal area index on adaxial (SD ad , SL ad , SAI ad ) and abaxial (SD ab , SL ab , SAI ab ) leaf surfaces, as well as relative proportions of palisade parenchyma (PT rel ), spongy mesophyll (SMT rel ) and intercellular spaces (IC rel ). Differences between groups according to one-way ANOVA and multiple pairwise comparisons applying Bonferroni correction or pairwise Wilcoxon rank sum test (for SMT exhibiting unequal error variance) (p < 0.05) are indicated by different letters. Significant differences at p < 0.05 are indicated in italics

Leaf anatomical traits reflect divergent environmental conditions
Environmental conditions at the natural growing sites of H. pusillum (HP N ) and H. veselskyi (HV N ) were reflected in their leaf anatomy, which was adapted to higher irradiance (PPFD) and better water availability in HP N . Thickness of leaves (LT) and palisade parenchyma (PT) were higher in  The effect of species pooling growing sites (i.e. HP N , HP G vs. HV N , HV G , BSpecies^), of growing site pooling species (i.e. HP N , HV N vs. HP G , HV G , BSite^) and their interaction (BSpecies × Site^) were included as predictor variables for single leaf anatomical traits. As the natural growing sites of the two species are strongly divergent, the effect of growing site is not meaningful. Traits include total thickness of leaves (LT), thickness of palisade parenchyma (PT), thickness of spongy mesophyll (SMT), thickness of epidermis on adaxial (ET ad ) and abaxial (ET ab ) leaf surface, thickness of cuticle on adaxial (CT ad ) and abaxial (CT ab ) epidermis, stomatal frequency, length and stomatal area index on adaxial (SD ad , SL ad , SAI ad ) and abaxial (SD ab , SL ab , SAI ab ) leaf surfaces, as well as relative proportions of palisade parenchyma (PT rel ), spongy mesophyll (SMT rel ) and intercellular spaces (IC rel ). Significant differences at p < 0.05 are indicated in italics. To achieve normal distribution and variance homogeneity, three variables were transformed by common (log10) or natural logarithm (log) or to the power of four (^4) HP N than in HV N (Fig. 3, Table 1), reflecting the higher irradiance at the growing site of HP N . Irradiance was characterised by an almost fourfold mean daily sum of photosynthetic photon flux density during the vegetation period at the HP N site as compared with the HV N site (Bertel et al. 2016). Increase in thickness of the palisade parenchymawhich contains the majority of chloroplasts and channels direct light, enabling plants to increase net photosynthetic rates (Boardman 1977;Björkman 1981;Vogelmann et al. 1996)-was due to longer cells and in some individuals to a second layer of palisade cells (Table 1). Congruently, higher net photosynthetic rates at high irradiance were measured in HP N compared to HV N (Bertel et al. 2016). Higher values of LT in HP N could also be advantageous at its alpine growing site, as adjustment and/or adaptation to increased elevation often involves an increase in leaf thickness to compensate for the lower CO 2 pressure at increased elevation (Körner 2003;Gratani et al. 2014). LT was low in HV N as was thickness of spongy mesophyll (SMT), but its relative proportion (SMT rel ) and the relative proportion of intercellular spaces (IC rel ) were in tendency higher in HV N than in HP N (Fig. 3, Table 1). This pattern is typically observed in shade-adapted leaves, as a relatively large proportion of spongy mesophyll enhances leaf absorbance due to greater internal light scattering at gas-liquid interfaces and thus returns photons to the palisade cells (Vogelmann et al. 1996). Intercellular spaces ranged from 4 to 12 % of the leaf cross section area and appear rather underestimated, as they usually comprise 30-40 % in dorsiventral leaves of other species (Vogelmann et al. 1996). Chloroplasts ultrastructure reflected the lower PPFD at the growing site of HV N , as chloroplasts were thinner and contained fewer starch globuli (Fig. 4). Thickness of grana thylakoids seemed to be elevated in HV N as expected for shade-adapted leaves, since highly stacked grana enhance light absorption (Larcher 2003). However, the evidence remains scarce as our TEM analysis was only qualitative.
Adjustment and/or adaptation to drier growing sites is usually reflected in stomatal patterns (Larcher 2003). Under dry conditions, reduced stomatal density on upper leaf surfaces prevents evapotranspiration while more and shorter stomata on lower leaf surfaces allow a more precise regulation of gas exchange (Lösch 2001;Larcher 2003). The lower stomatal area index on adaxial leaf surfaces (SAI ad ) measured in HV N compared to HP N may be advantageous at the drier growing site of HV N as a lower SAI can limit transpiration (Lloyd and Woolhouse 1978;Dixon 1986;Gratani et al. 2014  Missing values of traits SAI ad and SAI ab were replaced by the group mean. Variables were transformed to achieve normal distribution, and scaled and centred. Confidence ellipses are defined by the centroid and the standard deviation of the cloud. Ordination axes represent 48.0 % (xaxis) and 17.2 % (y-axis) of the explained variance. Arrows in the dashed circle (r = 1) represent direction and magnitude of effects of structural traits HV N as plants growing at low irradiance develop fewer stomata (Björkman 1981). However, on abaxial leaf surfaces, stomatal area index (SAI ab ) did not differ between HV N and HP N , and stomatal density (SD ab ) was in tendency higher and stomatal length (SL ab ) lower in HV N (Fig. 3, Table 1). The higher stomatal density on adaxial leaf surfaces (SD ad ) in HP N may also be advantageous at higher elevation, as alpine plants often have more stomata and a higher proportion of stomata located on upper leaf vs. lower leaf surfaces. This most likely counteracts photosynthetic limitations resulting from lower CO 2 pressure (Körner 2003;Aryal and Neuner 2012;Kammer et al. 2015). The reduced thickness of the cuticula (CT) observed in HV N (Fig. 3, Table 1) is not necessarily directly related to an increased water permeability. The main barrier for water diffusion within the cuticula is located within a thin (1 μm) waxy band; hence, thickness may not reliably indicate overall water permeability (Kerstiens 1996;Lambers et al. 2008;Riederer and Schreiber 2001).
Different leaf anatomical traits of H. pusillum and H. veselskyi grown in a common garden indicate an (epi-)genetically based differentiation between species Leaf anatomy differed between H. pusillum and H. veselskyi when grown in the same common garden environment (HP G and HV G ; Figs. 2, 3, Tables 1, 2). Hence, the differentiation between the two species must be (epi-)genetically based to some extent. This is supported by the fact that in the common garden leaf anatomical traits (LT, PT, CT ab , SAI ad ) were differentiated in the same direction between H. pusillum and H. veselskyi as at natural sites. Both species adjusted their leaf anatomy in response to the increased irradiance in the common garden, although the increase in irradiance was minor for H. pusillum but considerable for H. veselskyi (i.e. a 1.2-vs. 4.6-fold increase in daily sum of photon flux density; Bertel et al. 2016).
Adaptation to increased irradiance was mirrored in LT and PT, which was higher in HP G than in HV G ( Table 1). The higher PT in HP G was due to a second row of palisade cells, whereas in HV G , palisade parenchyma consisted of a single cell row. The number of rows of palisade cells might be genetically fixed as not all plants have the ability to build a second layer of palisade parenchyma (Lambers et al. 2008). In contrast, a similar SMT was attained in both species (Fig. 3,  Table 1). We speculate that the different reaction norms for LT and PT reflect an adaptation to the natural growing sites, which are characterised by moderate to high irradiance at the site of H. pusillum vs. very low irradiance at the site of H. veselskyi (Bertel et al. 2016). The significant effect of species in the two-factorial ANOVA for leaf traits related to irradiance and water availability, such as LT, PT, PT rel , ET ad , CT ad , CT ab , SD ad and SAI ad , indicates an (epi-)genetically based component of trait expression (Table 2). Chloroplast ultrastructure is known to be highly plastic (Larcher 2003) and did not differ remarkably between the species. However, indications for adjustments and/or adaptations to shade in HV N are given by the larger grana stacks. The considerable amount of plastoglobules in HV G may indicate some sort of irradiation stress for H. veselskyi compared to H. pusillum, as the formation of plastoglobules usually results from the degradation of thylakoid membranes (Austin et al. 2006). Similar observations were quantitatively analysed in drought-stressed Picea plants, where the chloroplast volume occupied by thylakoid membranes increased from 27 to 37 %, and that of plastoglobules from 2 to 12 %, respectively (Zellnig et al. 2010). We are aware that such changes are reversible; however, they are a good indication of the physiological plasticity that will occur under continuously altered environmental factors.
In conclusion, our study, even if based on the comparison of a single population each of both investigated species, supports the hypothesis that differentiation between the two postglacially diverged species H. pusillum and H. veselskyi (Trucchi et al., unpublished) is likely driven by natural selection as (1) differentiation in leaf traits at natural growing sites was in congruence with environmental habitat conditions and (2) differentiation in the common garden reflected adaptation to the conditions of their growing site. Continuing studies investigating phenotypic differentiation among six population pairs and the fitness of each species in the habitat of the other will further contribute to our understanding of the functional differentiation between H. pusillum and its multiply independently evolved, rare and ecologically peculiar descendant H. veselskyi.