Tetraploids expanded beyond the mountain niche of their diploid ancestors in the mixed-ploidy grass Festuca amethystina L.

One promising area in understanding the responses of plants to ongoing global climate change is the adaptative effect of polyploidy. This work examines whether there is a coupling between the distribution of cytotypes and their biogeographical niche, and how different niches will affect their potential range. The study uses a range of techniques including flow cytometry, gradient and niche analysis, as well as distribution modelling. In addition, climatic, edaphic and habitat data was used to analyse environmental patterns and potential ranges of cytotypes in the first wide-range study of Festuca amethystina—a mixed-ploidy mountain grass. The populations were found to be ploidy homogeneous and demonstrate a parapatric pattern of cytotype distribution. Potential contact zones have been identified. The tetraploids have a geographically broader distribution than diploids; they also tend to occur at lower altitudes and grow in more diverse climates, geological units and habitats. Moreover, tetraploids have a more extensive potential range, being six-fold larger than diploids. Montane pine forests were found to be a focal environment suitable for both cytotypes, which has a central place in the environmental space of the whole species. Our findings present polyploidy as a visible driver of geographical, ecological and adaptive variation within the species.

The patterns of geographic occurrence. The pattern of distribution of diploid and tetraploid F. amethystina is parapatric: a general allopatric distribution with a few areas with closer occurrence, which could act as potential contact zones (Fig. 2). In general, the central part of the species range is occupied by tetraploids, and the marginal area by diploids. The most visible exception is the presence of tetraploid populations on the northern border of the species range, e.g., in the Polish Lowlands. Six areas were found to be potential contact zones between cytotypes of F. amethystina (marked in Fig. 2 and described in detail in Supplementary Table S5). It can also be seen that diploids of F. amethystina are associated with mountain ranges, such as the Alps, Carpathians and mountains of the Balkan Peninsula, while tetraploids occur both in the mountains and lowland areas, particularly the Polish Lowlands (Fig. 2).

Altitudinal-climatic gradients.
Our results indicate that diploids tend to prefer higher altitudes than tetraploids, and they are not found in altitudes below 500 m. In contrast, tetraploids tend to occur at lower www.nature.com/scientificreports/  All populations were found to exhibit homogeneous ploidy. The grey areas indicate an altitude above 500 m. The distribution range of F. amethystina is an updated version of that identified by Kiedrzyński et al. 30 . The numbered potential contact zones between cytotypes are described in Supplementary Table S5.  Table S3). Both cytotypes occur in climates that have a relatively higher level of precipitation, both annually and in the warmest quarter (Fig. 3B), and therefore favour mountain regions and temperate subcontinental climate (Fig. 4). The marginal climates for diploids (high-altitude localities in Pindus Mts., Rila Mts.-see Fig. 4) are characterised by higher temperatures of the driest quarter; this is characteristic of sub-Mediterranean regions with a dry period during the warmer part of the year. Lowland localities of tetraploids are characterised by relatively   According to the analysis including Maxent niche modelling, the most influenced predictors in models performed for F. amethystina diploids were annual precipitation (Bio12), annual mean temperature (Bio1) and precipitation in the warmest quarter (Bio18) ( Table 1), while the most influenced in models performed for F. amethystina tetraploids was the mean temperature of the driest quarter (Bio9). In addition, Bio12, for diploids, and Bio9, for tetraploids, achieved above 50% contribution when only climatic predictors were used. When the geology layer was added to models, it was essential as a predictor and even the most influenced variable in the case of tetraploids (47.1%). For diploids, the influence was more evenly distributed between parameters in all tested models (Table 1).
Response curves for the most influenced climatic variables show that the diploids occur most likely in areas with annual precipitation higher than 1000 mm, precipitation of the warmest quarter of around 400 mm and mean annual temperature not higher than 3-5 °C (Supplementary Figs. S2, S3). Tetraploids have a higher optimal temperature range (Bio1 and Bio7), but the most influenced predictor, i.e. the mean temperature of the driest quarter, should be not higher than 0 °C ( Habitat and parent material of soils along with altitude. The altitudinal gradient of cytotype localities is also reflected in their habitat conditions. Both cytotypes inhabit medium-montane, calcareous pine forests and associated grasslands (  Table S3, not shown in Fig. 5A).
The frequency of cytotypes on different types of soil parent material (geological units) also varies according to altitude (Fig. 6A). Tetraploids of F. amethystina demonstrate broader occurrence with regard to geology (units no. 1-5 in Fig. 6A, Supplementary Fig. S5), while diploids have more specific preferences (units 1-3). Flysch and clastic rocks and glacial deposits in lower altitudes are unique to tetraploids. In addition, on metamorphic rocks and fluvial deposits, diploids tend to occur at higher altitudes than tetraploids. A similar pattern is also visible in calcareous rocks, which is the most common parent material of soils for both cytotypes (Fig. 6A, Supplementary Fig. S7).
The analysis of the habitat, soil parent material and altitudinal gradients (Figs. 5A, 6A, Supplementary Figs. S6, S7), indicates that the focal environment for both F. amethystina cytotypes consists of medium-montane pine forest and grasslands on calcareous soils, located in the "central" part of bioclimatic niche of the species (Figs. 5B, 6B). Other habitats and parent materials of soils occupy rather marginal climatic conditions and altitudes.
Potentially suitable areas. The presence-absence prediction maps of suitable areas for the F. amethystina cytotypes ( Fig. 7) were created based on MTSS thresholds, which for diploids was 0.601 and for tetraploids 0.484.
Areas with potentially suitable conditions for diploids appear to be smaller (MTSS predicted area 81 × 10 3 km 2 ) and concentrated in mountain ranges (Fig. 7). Tetraploids have a more extensive potential range, being six-fold larger (517 × 10 3 km 2 ).
For both cytotypes, models showed areas where they can potentially occur, but they were not found there (Fig. 7). For example, diploids have potentially suitable areas in the Dinaric Mountains or are more widely distributed in the Eastern and the Southern Alps. Tetraploids tend to be more extensively distributed in the Eastern Table 1. Percentage contribution of predictors used in Maxent niche modelling in diploids and tetraploids of Festuca amethystina. Models with two sets of predictors were used and compared in the analysis: only climatic and climatic + geology (parent material of soil) variables. www.nature.com/scientificreports/ and Western Carpathians and their northern forelands, as well as in the southern and eastern part of the Alps and the Bohemian Massif.

Discussion
Cytotype distribution. Our results indicate that the tetraploids of Festuca amethystina have a geographically broader distribution range than the diploids. The differences between the distribution are represented by the presence of tetraploids in lowland localities, especially north of the Carpathians in the Polish territory. However, our modelling of the potential distribution indicated some regions where further examples of both cytotypes may be found by more extensive studies. Polyploids can repeatedly originate in different parts of the species range 8,40 . Hence, the diffuse distribution of tetraploids and adjacent diploid localities indicated in our analysis suggest the relict character of localities or that multiple polyploidization events may have occurred in the past. Moreover, it is also possible that tetraploids from various locations could have arisen from different parents. These scenarios should be examined in future cytogenetical studies. F. amethystina is perennial, and can persist vegetatively over extended periods of time in discrete habitats; therefore, as demographic processes may be ineffective due to the effect of minority cytotype exclusion, mixedploidy populations are possible 41 . However, no such results were obtained in the present study: none of the observed regions demonstrated any sympatric occurrence of F. amethystina cytotypes. This suggests that gene flow between cytotypes, if any exists, is currently restricted. This situation could be analogous to the closelyrelated F. norica, where diploids, tetraploids and hexaploids tend to occur in different parts of the Eastern Alps, and probably in most parts of the cytotype ranges populations are ploidy homogeneous 42 .
Altitude, climatic and habitat preferences of cytotypes. Based on the current state of knowledge, it can be concluded that, in general, polyploids tend to populate higher altitudes than diploids 13 , in terms of the diploid-polyploid ratio and the percentage of polyploids in the mountain flora 43 . However, our present findings indicate a somewhat reversed pattern, where diploids and tetraploids both occur in common medium-montane altitudes, and diploids have rather high-altitude and tetraploids low-altitude localities.
Our example represents an interesting pattern of cytogeography, where the parental diploids occupy current mountain habitats, and the polyploids reduce competition to diploids through living in nowadays lowland environments, which can be stressful for mountain plants. Our present findings support the altitudinal segregation of cytotypes with regard to environmental factors: the tetraploids extended their occurrence to general warmer climates, more diverse parent material of soils and habitats. However, in this case, climate seasonality has a visible impact. Tetraploids tend to occur in areas with relatively cold and dry winters; in the studied area characterised by a subcontinental temperate climate. www.nature.com/scientificreports/ Similarly, some studies have identified other polyploids that have expanded their niches and ranges beyond those of their progenitors [51][52][53] . The possibility of occurrence of F. amethystina tetraploids in current lowland habitats is connnected with such of cases, when polyploids occupy more temperate and more mesic environments than their diploid ancestors. Such polyploids often have favourable morphological features such as larger leaves, larger stomata and slower growth 54 . Similar ploidy-dependent patterns in plant habit and size have been described for F. amethystina in our previous study 37 . Tetraploids in lowlands as glacial relicts. The occurrence of mountain species in lowland areas is often considered a relict of previous cold periods. As such, the occurrence of F. amethystina in lowland areas has been interpreted as a glacial relict since early biogeographical studies 55 . Our findings indicate that nowadays speciesrich lowland oak forests, occasionally also beech forests, or associated meadow and grassland communities harbour only tetraploid populations. Therefore, it should be considered in what sense these lowland populations can be considered as a relict. The 'glacial relicts' should then represent the remnants of tetraploid larger ancestral range that occurred during glaciations. The current disrupted range of tetraploids, with many local islands of occurrence, supports this interpretation. However, it appears that this distribution pattern considered on the background of the Pleistocene glacial ranges, is the result of migrations that took place during earlier glaciations, not just the last ice age. These issues may become clearer after examining the phylogeography of the species as a whole.
Considering lowland populations as a relict has further implications for understanding adaptations in such polyploid complexes. Assuming that the habitats in which tetraploids occur play a refugial role, their microclimate should be decoupled from regional trends and demonstrate lower species competition 56 . Tetraploids occurring in lowland forests under the tree canopy can be less exposed to climatic extremes, such as drought and high temperatures 57 . Further experimental studies will aim to discover whether tetraploids have a higher tolerance to lowland conditions than their diploid progenitors, as suggested by our analysis. Alternatively, the occurrence of tetraploids in such regions may be made possible by the specific properties of the refugial habitats rather than the adaptation of plants per se.
Mentioned habitats play a refugial role for many other plants, but are shrinking as a result of the changing anthropogenic regime and, probably, anthropogenic climate change 29,58-60 . As a result, many lowland populations of F. amethystina are now considered extinct 61 . If this trend continues, lowland tetraploid populations, which probably harbour extensions of the species ecological niche, are under threat. The extension of species ecological niche demonstrated by lowland populations of tetraploids is in line with the assumption that diploids are older than tetraploids and have a mountain distribution and niche; as such, this extension appears to be the ecological constitution of tetraploids, which have apart mountain locations also lowland populations. www.nature.com/scientificreports/

Conclusions
Our data indicate that the cytotype distribution of pure-ploidy populations of F. amethystina follows a parapatric pattern. Even populations of cytotypes located in the same region are separated by several dozen kilometres, occur at different altitudes, or exist in different habitats. Tetraploids tend to occur at lower altitudes, on more varied geological units and in more varied habitats than diploids. Our results also identify a focal environment suitable for both cytotypes, which has a central place in the environmental space. The range expansion of tetraploids of F. amethystina probably occurred as a result of survival in lowland refugia during the Pleistocene climatic oscillations. The studied complex represents an example of polyploid expansion or adaptation to lower elevations than their montane diploid ancestors. In those localities, polyploids currently occur in refugial habitats, where regional climate extremes can be mitigated by topography or vegetation structure. Any direct link between the expansion of polyploids and their adaptations to drier and warmer lowland conditions needs to be determined in further experimental studies.  Table S1). In most cases, to confirm the ploidy of the sampled populations, at least five accessions were collected from plants in different parts of each particular locality. Each accession consisted of several leaves from an individual plant; these were either placed in plastic bags and dried in silica gel, or in paper envelopes and allowed to air dry at room temperature. The geographical coordinates of the centre point of the sampled localities were recorded in the field using a GPS Garmin 60CSX.

Materials and methods
Material for four localities was obtained from herbarium sheets deposited in Munich (M) and Lodz (LOD) herbaria (Supplementary Table S1). The herbarium materials were located based on coordinates contained in the source or those determined based on the locality description.
The altitude of the centre point of the localities (Supplementary Table S3), indicated by the GPS coordinates, was determined from Google Earth using our own script.
Ploidy level determination. The ploidy level of the sampled plants was estimated using flow cytometry (FCM). The nuclear DNA content was measured in the dry leaves of accessions. The samples were prepared for flow cytometric analysis either by using PI according to Rewicz et al. 37 (Supplementary Table S1) or DAPI, according to Šmarda et al. 27 (Supplementary Table S1). The present analysis included also FCM measurements of 108 accessions from our previous study mentioned above, and 328 taken for the present one.
In addition, further information on the ploidy level of F. amethystina was obtained from previously published studies; however, only data with indicated geographic locations were used: coordinates, locality names or forest complex names which can be precisely localized by topography or distribution of vegetation patches on local maps. Our search revealed twelve additional localities with specified ploidy levels (Supplementary Table S2). Therefore, our analysis included 71 localities with specified ploidy levels.

Current bioclimatic parameters.
Thirty arc second (~ 1 km) resolution raster data was used, incorporating 19 bioclimatic variables from the WorldClim database 62,63 . Firstly, all 19 bioclimatic variables were assigned to each locality from the raster cells. The values were extracted according to the coordinates of localities in the ArcGIS DesktopTM 9.2: Spatial Analyst tools, Extract values to point tool [ESRI Inc. 1999, Redlands, CA, USA].
The procedure of variable selection for analysis assumed that the main climatic variables of precipitation and temperature (Bio1, Bio 12) should be included, as well as some additional values. Variables should not be closely and positively correlated to each other were then identified, as such correlations may impede the interpretation of the principal component and range modelling analyses. We used the correlation matrix to choose the additional climatic values which show seasonality of climate. Correlation matrix was constructed using 'corrplot' package in R 64 (Supplementary Fig. S1). No single threshold of correlation coefficient was used in the procedure. For temperature variables, Pearson's correlation coefficient can be accepted lower than 0.60 to the main Bio1 variable; for even highly correlated precipitation variables, the correlation coefficient should not exceed 0.85 of the main Bio12 variable. Finally, the distribution of vectors of climatic variables was checked on the Principal Components Analysis (PCA). The selected values should show as far as possible the different directions of vectors in the PCA diagram. Finally, we chose six variables: three climate parameters for temperature (Bio1-annual mean temperature, Bio7-annual temperature range, Bio9-mean temperature of the driest quarter) and three for precipitation (Bio12-annual precipitation, Bio18-precipitation of the warmest quarter, Bio19-precipitation of the coldest quarter).
General habitats in localities. Each locality was classified as one of the following habitat types: (i) subalpine grasslands, (ii) montane pine forests and grasslands, and (iii) lowland oak forests and meadows, and (iv) lowland beech forests and grasslands. The assignment is based on our field observation/descriptions and supported by published materials (Supplementary Table S3). In the case of materials from herbaria, the habitat was assigned based on the information/descriptions on the herbarium sheets.
Raster with the parent material of soils. The  www.nature.com/scientificreports/ converted into groups of major geological units (Supplementary Table S3) to standardise any differences in geological unit interpretation between different countries. To account for any possible inaccuracy in the raster on the local scale, the specified geological data was analysed according to our localities, and the pixel assignment was changed in the raster if necessary. Finally, the geology of the localities was assigned based on the converted raster using ArcGIS DesktopTM 9.2: Spatial Analyst tools Extract values to point tool [ESRI Inc. 1999, Redlands, CA, USA] (Supplementary Table S3).

Gradient analysis.
The distribution of the studied cytotypes was then analysed according to climatic and altitudinal gradients using principal component analysis (PCA). In the PCA, the localities were treated as ' cases' and six of the bioclimatic parameters of localities as 'variables' . Unconstrained, linear PCA analysis was performed using Canoco software; the analysis was run together with default centring and standardizing of the general environmental data 65 . Altitude was related to the multidimensional PCA space of the climatic parameters using the Generalized Additive Model. The model parameters were as follows: response variable-altitude, predictors-PCA axes, response distribution-Gaussian, link function-identity. Stepwise selection of models was performed using the AIC (Akaike information criterion) and different Term Smoothness values. The lowest AIC was demonstrated in the model with Term Smoothness = 5.0 for both predictor 1 and predictor 2. Ordination, modelling, and graphic plots were performed using Canoco 5 software 65 . The analysis of the cytotype occurrences along the altitude gradient, and according to habitat and geological units, was performed using ridgeline plots in the ggridges package in R 66 . These are partially overlapping line plots that can be useful for visualizing changes in distributions along environmental gradients.
Modelling of potentially suitable areas. The predictive modelling of the distribution and niche analysis was performed using Maxent (version: 3.4.0), a machine-learning method which uses environmental variables presented in raster layers to predict the suitability of distributions or habitats [67][68][69] . Based on the maximumentropy principle, Maxent finds the most probable distribution for a species that maximizes the entropy 70 .
The raster layers used in the analysis were restricted to the background area; these included the whole range of the studied localities and covered a broad area of Europe. The layer with the parent material of soils was treated as categorical, and the six bioclimatic layers as numerical. Geological data were used to increase the realism and spatial precision of the results, as shown in our previous studies on F. amethystina 29,30 . The raster with the parent material of soils was transformed to match cell sizes with the bioclimatic data. After transformation into an ESRI ASCII format, the obtained data set was used in Maxent modelling. All the above transformations and point generations and extractions were done in ArcGIS DesktopTM 9.2: Spatial Analyst tools and Conversion tools [ESRI Inc. 1999, Redlands, CA, USA].
Our models were fitted using hinge features with default regularization parameters. All of the models were fitted on the basis of the full data sets, and tenfold cross-validation was used to estimate errors around the fitted functions and predictive performance of the held-out data 70 . Model selection. Because the studied localities are not equally distributed across the study area, the models based on all localities were compared with a model based on spatially filtered data. Spatial filtering is an approved solution when small sample sizes are used with respect to the study area, as demonstrated in a previous study of F. amethystina 30 . To select the best model, we used specific evaluation measures, such as AIC (Akaike information criterion), BIC (Bayesian information criterion), and AUC (Area Under the Receiver Operating Characteristics Curve) 69 . AIC and BIC indexes were calculated in the ENMTools software 71 , and AUC was calculated in Maxent. Detailed results of model selection are included in Supplementary Table S4. Generally, models with a filtered sample-set resulted in the best model quality measurements and were used in further analyses.
Modelling of potential ranges. To determine the influence of differences between niches on potential distribution of cytotypes, the prediction maps were developed. The Maxent prediction map provides a geographical visualization of the prediction model, which can be interpreted as the predicted distribution of the most suitable areas for particular cytotypes 70 .
The potential distribution of suitable areas for both cytotypes was modelled in the current climate, using the geology layer as a predictor. The prediction maps were generated in ASCII file format and were visualized in ArcMap 9.2 software (ESRI Inc. 1999-2008, Redlands, CA, USA).
We generated binary (presence, absence) maps for each cytotype. To obtain binary predictions of the climate + geology suitability of ENMs, MaxEnt's logistic probability of occurrence output was converted to a binary mode (presence-absence output) using the maximum training sensitivity plus specificity logistic (MTSS) threshold. By maximizing the proportions of actual positives and negatives that were correctly identified, MTSS-based predictions are believed to be the most accurate forecasts of the potential presence or absence of species 72  www.nature.com/scientificreports/