Seed production and dispersal limit treeline advance in the Pyrenees

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2019 The Authors. Journal of Vegetation Science published by John Wiley & Sons Ltd on behalf of International Association of Vegetation Science 1Institute of Botany and Landscape Ecology, University of Greifswald, Greifswald, Germany 2Department of Evolutionary Biology, Ecology and Environmental Sciences, Biodiversity Research Institute (IRBio), University of Barcelona, Barcelona, Spain 3Centre de Recerca Ecològica i Aplicacions Forestals (CREAF), Universitat Autònoma de Barcelona, Cerdanyola del Vallès, Spain

vancements in elevation may be explained by the decrease and even cessation of anthropogenic activities in the last century. In such cases, treelines are not expanding beyond their historical thermal limits, but recovering their past altitudinal distribution (Ameztegui et al., 2016;Cudlín et al., 2017). Regardless, the treeline advance, either through range expansion or through re-colonisation, may lead to changes in extant plant diversity and ecosystem function through displacement of alpine species and shifts in carbon and nutrient dynamics (Greenwood & Jump, 2014).
Production of viable seeds is the first step for a seed-based regeneration process at any treeline (Holtmeier & Broll, 2005). The next essential step is the dispersion of the seeds, which represents the first spatial demographic process that will determine the position of recruit establishment and potential range shifts (Nathan & Muller-Landau, 2000). Ultimately, any change in treeline position is thus related to successful seed production, viability, dispersal and germination, as well as seedling survival, at higher elevations (or northern latitudes). Some studies have reported a decrease in cone and seed production or seed germination with elevation and latitude (Brown et al., 2019;Juntunen & Neuvonen, 2006;Kruse et al., 2019;Šenfeldr & Maděra, 2011;Sirois, 2000), pointing to potential reproductive limitations at the treeline. These reproductive constraints seem exclusively related to changes in population density in some cases (Juntunen & Neuvonen, 2006), whereas in others they seem related to changes in the reproductive capacity of individual trees (Sirois, 2000), or a combination of both factors (Šenfeldr & Maděra, 2011). However, the reproductive ecology of treelines is still largely unknown (Körner, 2012), especially in relation to the dynamics of alpine treeline ecotones.
In this study, we aimed to address this gap by focusing on seed production, viability, dispersal, and germination dynamics of Pinus uncinata. Our main objective was to analyse whether seed production and dispersal can constrain the elevation advancement of the treeline in the Pyrenees, and to determine variations in seed germination dynamics between different provenances (i.e., elevation positions) along the elevation gradient under controlled environmental conditions. We hypothesised that: (a) viable seed production would decrease with elevation; (b) that seed rain and seed arrival would decrease along the elevation gradient; and (c) that seed germination rates and dynamics would mirror the provenance in the elevation gradient, with seeds from lower elevations (where temperatures are higher and growing seasons longer) showing higher germination success and a more extended germination period than seeds from trees growing at higher elevations.

| Study sites and species
Seed production was assessed at five sites in two protected areas of the Central Pyrenees, Catalonia, Spain. The study sites Gelada, Dellui, Son and Eixol are located in the Aigüestortes i Estany de Sant Maurici National Park, and the study site Capifonts is located in the Alt Pirineu Natural Park ( Figure 1, Table 1). A sixth study site, Amitges, was used for the seed rain experiment together with Gelada and Dellui ( Figure 1, Appendix S1). All the sites occupy NW-NE facing slopes. The anthropogenic alteration of forests in these areas has been minimal since the mid-twentieth century due to their protection status, as well as an ongoing reduction in livestock numbers and the abandonment of forest clearcutting (Domínguez Martín, 2001;Gracia, Meghelli, Comas, & Retana, 2011). Forest expansion has been noticeable from 1956 to the present, although this phenomenon has produced little displacement of the treeline in such near-natural treeline locations . The mean annual precipitation in climate stations located near the study sites for the period 2008-2014 is 1,340 mm, with a mean temperature of 9.8°C for the growing season months (JJA) and −3.6°C for the winter months (DJF) (Servei Català de Meteorologia). For detailed climate data, see Appendix S2.
The study species is Pinus uncinata, the dominant tree species in most of the treeline ecotones in the Pyrenees (Ninot et al., 2008). It forms dense subalpine forests between 1,700 m and 2,200 m a.s.l.

| Seed production along the elevation gradient
We carried out the sampling from mid-October to mid-November 2010. At each of the five study sites, we established a transect along the elevation gradient where we set four plots corresponding to different elevation positions and plant communities: dense forest at mid-elevation within the subalpine belt (DFM), dense forest at maximum elevation within the subalpine belt (DFH), scattered trees (ST), and krummholz zone (KZ) (Figure 1, Table 1). The three highest sampling plots were close to each other and encompass the treeline ecotone. At each plot we sampled 5-6 mature trees and collected 5-13 cones per tree. In total, we collected 635 cones across plots and study sites. The minimum distance between sampled trees was 50 m, and they were selected to include the local forest heterogeneity (e.g., tree size). To include variability linked to position within the canopy, cones were collected at different orientations across the canopy.
We applied a cold treatment to the collected cones. We placed them in growth chambers for one week at 4°C, one month at −18°C, and 2.5 months at temperatures ranging between −10°C and 9°C, to simulate natural conditions. We triggered cone opening with a moderate heat treatment that consisted of an initial temperature of 30°C that was gradually increased over 8 hr up to 55°C, which was maintained for 10 min. This temperature does not affect the viability of Pinus uncinata seeds (Escudero, Sanz, Pita, & Pérez-García, 1999). We mechanically opened those cones that were not completely open after the heat treatment to extract all their seeds. For each cone, we measured the maximum length (from base to apex) and assessed the number of empty seeds (without embryo) and the number and weight of fully developed seeds (hereafter full seeds).
Additional information on the assessment of cone production at two of the study sites can be found in Appendix S3.

| Germination experiment
We carried out a germination experiment to test the potential viability and germination dynamics of all full seeds under comparable environmental conditions. We used a growth chamber under con-  ANADON-ROSELL Et AL. et al., 2002) and are close to those found in the subalpine belt in the most favourable period for germination in early summer (Ninot et al., 2008). We sowed all full seeds from the same cone in one pot (7.5 cm × 7.5 cm × 8 cm), or in two pots when cones had more than 50 full seeds, always at 1-cm depth. We used a sterile substrate of 1:1 perlite and vermiculite, watered weekly to saturation with deionised water and let the pots drain freely. We recorded seedling emergence and removed seedlings twice a week to avoid possible competition or inhibition. The germination trial lasted 121 days.

| Seed rain and seed dispersal along the altitudinal gradient
To explore in more detail seed production and dispersal of Pinus uncinata, we performed a seed rain experiment at three of the study sites, Gelada, Dellui and Amitges (Figure 1, Appendix S1). To account for the interannual variability in seed production, we carried out this experiment during three consecutive years. At each site we established one transect with four plots (elevation positions) along the elevation gradient. Plots corresponded to the subalpine dense forest (SDF), subalpine fragmented forest (SFF), krummholz zone (KZ), and alpine grassland (AG). Note that in Gelada and Dellui, seed rain plots do not correspond to the plots used in the seed production experiment (Appendix S1). In Amitges, due to its characteristics, the KZ and AG plots could not be established. Between 28 September and 5 October 2011 we installed six seed traps per plot at each site.
Seed traps were 1 m 2 and were made of artificial tuft grass. In each plot, the six traps were installed forming a circumference of 7-8 m in diameter. From 2012 to 2014 we collected all seeds in each trap during summer and counted them at the laboratory. Since dispersal distances of Pinus uncinata are frequently within 10 m from the parental tree (Camarero et al., 2005;Vitali et al., 2019), we established a circular area adding a supplementary radius of 10 m from where the circumference of seed traps was located. Within this circular area (seed trap area + external circumference), we counted all the trees present, we measured their diameter at breast height (DBH), and classified them as reproductive trees (DBH > 5 cm), potentially reproductive trees (non-reproducing sexually mature trees with DBH > 5 cm), non-reproductive trees (DBH < 5 cm), and dead trees.
We used these data to calculate reproductive tree density and the mean DBH of reproductive trees in the area covered by the seed traps at each study site and elevation position (see Appendix S4).

| Statistical analyses
We used linear and generalised linear mixed-effects models to analyse our data. For all the study variables, we performed model selection following the Akaike information criterion (AIC) to select the most suitable model (Zuur, Ieno, Walker, Saveliev, & Smith, 2009).
When two models did not differ in more than two AIC units, we applied a model average function (Bartoń, 2019). TA B L E 1 Altitude (m a.s.l.) and UTM coordinates (ED50, zone 31T) of each sampling plot at each study site for the seed production and germination analysis Note: For each study site, the complete toponym is given in parentheses.

Journal of Vegetation Science
ANADON-ROSELL Et AL.
We used generalised linear mixed-effects models to assess the effects of the elevation gradient on the number of seeds and number of full seeds per cone (fitted with a Poisson distribution) and on the proportion of full seeds per cone (fitted with a binomial distribution).
We used linear mixed-effects models fitted with the restricted maximum-likelihood method (REML) to test the effects of the elevation gradient on cone length and weight of full seeds. For all these variables, we used the elevation position as a fixed-effect factor and site and tree identity as random factors as a starting point for model selection. Additionally, we explored the relationship between cone length and number of full seeds per cone with Pearson correlation tests and exponential regression.
To analyse differences in the germination success (i.e., seed viability) at the end of the germination experiment (seeds germinated out of seeds sown), we used generalised linear mixed-effects models fitted with a binomial distribution. We included elevation position and site as fixed effects and tree identity as a random factor prior to model selection.
For the analysis of germination dynamics, we used meta-analytic random-effects models, which allow incorporating random factors in log-logistic models of germination data (Keshtkar, Mathiassen, Beffa, & Kudsk, 2017;Ritz, Pipper, & Streibig, 2013). For this, we carried a two-step approach. First, we fitted a three-parameter log-logistic model to the cumulative germination data according to Ritz et al. (2013): where F(t) is the cumulative seed germination at time t, d is the upper limit parameter denoting the proportion of seeds that germinated during the experiment out of the total number of seeds sown, t 50 is the time when 50% of the seeds that germinated during the experiment (d) have germinated, and b is the slope of F at time t = t 50 .
In the second step, we separately fitted the meta-analytic random-effects model to our parameters of interest. Therefore, we analysed the estimates of t 50 , b and d obtained with the event-time model (first step) with a linear mixed-effect model where tree was defined as a random factor and elevation position, site and the interaction elevation position × site were defined as fixed effects. Site was included here as a fixed effect to specifically investigate potential effects of the different provenances (i.e., elevation positions) of the seeds on the dynamics of germination. The estimated standard errors for the parameters of interest were also included in the models. We made pairwise comparisons to find significant differences between the factor levels.
We used a zero-inflated negative binomial (ZINB) model to assess the effects of the elevation gradient on seed rain per area.
We used these models to account for data overdispersion and zero inflation (their AIC values were lower than in generalised linear models fitted with Poisson and negative binomial distributions).
We used site, elevation position, study year and the interactions position × year and site × year as fixed-effect factors in the count part of the model. In the binomial part, we used the three factors without interactions following the model selection procedure.
We excluded the alpine grassland plots from all the analyses on seed rain because no seeds were found in these plots at any of the study sites and years.
We tested for normality and homoscedasticity of the residuals and applied data (sqrt-and log-) transformations when necessary to reach these assumptions in the linear regression models. We also used varIdent structure when the homogeneity of residuals was not reached (Zuur et al., 2009). When fixed terms were selected in the final model, we used post-hoc Tukey HSD tests to identify significant differences between fixed term levels.
For all statistical analyses we used r software v. 3.6 (R Core

| Seed production
The number of total seeds per cone decreased 65% along the elevation gradient from DFM to KZ. Values were significantly lower at KZ (16.4 ± 0.9) than in the other three elevation positions (DFM = 46.8 ± 1.8, DFH = 36.6 ± 1.7, ST = 33.3 ± 1.7, Tukey posthoc test p < 0.05). The number of full seeds per cone also decreased along the gradient, from 37.5 ± 1.8 full seeds per cone at DFM to 14.0 ± 0.8 at KZ (Table 3, Figure 2a), which represents a 63% reduction. Post-hoc tests showed significant differences between KZ and the forest sites (DFM and DFH, p < 0.05) but not between KZ and ST (p = 0.171). The elevation position showed no significant effect on the proportion of full seeds to empty seeds per cone and on the weight of the full seeds (mean across sites and elevation positions of 9.15 ± 0.2 mg/seed) ( Table 3).
The elevation gradient also had an effect on cone length, which decreased with increasing elevation (Table 3). Post-hoc tests showed that at the two higher elevations, KZ (3.92 ± 0.05 cm) and ST (4.35 ± 0.06 cm), cone length was significantly lower than in the two forest plots, and within these two, the lower subalpine forest DFM (4.79 ± 0.05 cm) showed longer cones than the higher subalpine forest DFH (4.43 ± 0.06 cm). Correlation analyses showed a significant exponential correlation between cone length and total full seeds per cone (no. viable seeds = exp (0.67 + 0.57 × cone length), p < 0.001 for both parameters;

| Seed germination
A mean of 65% of the sown seeds across the four assessed elevation positions germinated after 121 days (Figure 3, Appendix S6).
For the analysis of the germination success (proportion of germinated seeds at the end of the experiment) we selected two generalised linear mixed-effects models with a difference in AIC lower than two units. Elevation, study site, and their interaction were included in one model whereas only site was included in the other.
Model average only showed a significant effect of study site on the germination success (Table 4, Appendix S6), with the largest differences between Dellui and Eixol (post-hoc Tukey test run for each model, p < 0.01).
The meta-analytic random-effects models did not show consistent patterns in the germination dynamics between elevation plots across study sites for any of the estimated parameters ( Figure 3). However, when looking at the study sites independently we found significant had a lower t 50 than DFH (estimated difference of 26.6 ± 11.9 days, p = 0.025). At the study site Capifonts, the scattered trees (ST) had a higher t 50 than (KZ; estimated difference of 25.3 ± 9.7 days, p = 0.009) and higher than DFM (estimated difference 25.2 ± 11.6 days, p = 0.030). Finally, Dellui KZ had a lower t 50 than DFM (estimated difference 21.9 ± 11.2 days, p = 0.05, respectively).

| Seed rain and dispersal
Seed rain significantly decreased along the elevation gradient. The number of seeds recovered at each elevation position was significantly lower than at the elevation directly below within the elevation gradient, and it was zero at alpine grasslands (Table 5, Figure 4, p < 0.001 after Tukey post-hoc tests between each elevation plot).
The study site and the study year also had a significant influence on seed rain values.

Journal of Vegetation Science
ANADON-ROSELL Et AL.

| D ISCUSS I ON
Our study investigated the reproductive and dispersal patterns of the treeline-forming species Pinus uncinata along the treeline ecotone in the Pyrenees. Our findings on seed production, viability, and dispersal provide evidence of reproduction as a limiting factor for treeline advance.

| Seed production and viability at the treeline
Our results show a clear decrease in seed production of Pinus uncinata along the elevation gradient from the dense subalpine forest to the krummholz zone. Our analysis suggests that such decrease is directly related to the size of the cones.
Heat sum fluctuations, thermal limitations, and shorter growing seasons at high elevations may affect both cone production and development, and may have negative effects on fertilisation and embryo growth processes (Almqvist, Bergsten, Bondesson, & Eriksson, 1998;Sirois, 2000). Climate may act in synergy with other factors such as forest structure, resulting in lower pollen availability and an increased self-fertilisation (selfing) with the decrease in forest density along the elevation gradient (Iwasaki et al., 2013;Smith, Hamrick, & Kramer, 1988). We did not find a reduction in the proportion of full seeds with elevation, which would be the first negative outcome of selfing in other Pinus species (Iwasaki et al., 2013). Contrarily, our results suggest that seed development on the Pinus uncinata alpine treeline could be more influenced by low pollen supply at higher elevations (Iwaizumi & Takahashi, 2012). The decrease in the production of full seeds found along the gradient in our study sites is not accompanied by a reduction in the number of forest treeline transect in northern Québec, Canada. This author did not find differences in cone production along the transect (similar to our results), and suggested that the decrease in seed production and viability was related to a low viability of the pollen at the treeline.

TA B L E 3
Results of the statistical models for the parameters related to seed production Note: The estimated coefficient and standard error, z-or t-value, and p-value are shown for the fixed terms, including the interaction between two fixed terms when selected in the final model. The estimated variance and standard deviation are shown for the random effects in the mixed-effects models.

ANADON-ROSELL Et AL.
Šenfeldr and Maděra (2011) also found a clear decrease in the reproductive output of Norway spruce (Picea abies) in a former pastoral timberline ecotone in the Czech Republic, namely fewer seeds per cone, fewer cones per ha, and fewer cones per fertile specimen at the upper part of the ecotone than at the lower part. Another study in northern Finland found slight differences in seed anatomical maturity between timberline and treeline populations for Scots pine (Pinus sylvestris) but no differences in the number of cones per tree (Juntunen & Neuvonen, 2006). Overall, it seems a common pattern to find a decrease in reproductive output related to seed set of the treeline-forming species across elevation or latitudinal gradients, whereas patterns related to cone production and the proportion of viable seeds produced seem more variable among species.
Our findings strongly suggest that the production of viable seeds is the first limiting factor for forest expansion at the Pyrenean treeline. Even under the ongoing climate warming (IPCC, 2018), and unlike the situation in other treeline-forming species (Kullman, 2002), this suggests that the reproduction of Pinus uncinata at the treeline is still strongly limited by the present climate conditions. We acknowledge that fluctuations in the reproductive output due to inter-annual climate variations may occur (see section 4.2, Seed rain and dispersal, below), but we believe that the reported patterns and the relationship between the parameters studied (e.g., proportion Contrary to our expectations, we did not find a significant effect of elevation provenance on seed viability (i.e., germination success) and germination dynamics across sites, although the influence of elevation was marginally observed when assessing the sites separately. The lack of a consistent effect on the germination success by the elevation provenance reported here contrasts with other studies showing that germination is an adaptive trait that can vary largely with elevation (Castanha, Torn, Germino, Weibel, & Kueppers, 2013;Rehfeldt, 1989). Sirois (2000) showed that, 30 days after sowing, the percentage of total germination from seeds originating in trees closer to the treeline was lower than in those from the dense boreal forest. Results of Šenfeldr and Maděra (2011) reflected a similar pattern after 21 days, but they did not find significant differences. The different study species and the longer length of our germination trial (121 days) make it difficult to compare between these studies and our experiment. We did find, however, large differences in the germination success of seeds from

| Seed rain and dispersal
The seed rain experiment showed large differences between study sites and study years, emphasising the importance of local and inter-annual conditions for the reproduction dynamics of Pinus uncinata. However, despite the seed rain being much larger in 2012 than in the following years (2013 and 2014), we found a clear decrease in the number of seeds recovered in the seed traps along the elevation gradient that was constant through time. No seeds were found in the alpine grassland beyond the treeline, indicating that seed arrival, which is the first step for an eventual forest colonisation of alpine grasslands, is a limiting factor for treeline advance. Although the area covered by our seed traps is relatively small for concluding that no seeds arrive at alpine grasslands, our findings strongly suggest that important limitations for seed arrival beyond the krummholz zone exist. We acknowledge the possibility that seed predation influences our results, since it has been found to largely affect regeneration beyond the upper elevation limit in some species (Brown & Vellend, 2014;Castro, 1999  of dispersal constraints for Pinus uncinata seed (Camarero et al., 2005;Vitali et al., 2019).
The decrease in seed recovery with elevation was not significantly correlated with reproductive tree density, although the number of reproductive trees per area also decreased along the elevation gradient (Appendix S4). We believe these results need to be taken cautiously because there was a large variation in tree density among sites in the two subalpine forest positions (SFF and SDF), and the forest structure of the study sites was complex. There is, however, a clearly lower density of reproductive trees in the krummholz zone, which suggests that the lower number of seeds per area recovered there may be driven both by a lower reproductive capacity of individual trees and a decrease in tree density, as reported in other studies (Juntunen & Neuvonen, 2006;Šenfeldr & Maděra, 2011). Contrarily, the decrease in seed recovery from the subalpine dense forest (SDD) to the subalpine fragmented forest (SFF) may be mostly driven by a decrease in the reproductive capacity of individual trees, as shown by Sirois (2000). We found a positive correlation between the mean DBH of reproductive trees and the number of seeds recovered per area (Appendix S4), similarly to results found by Vitali, Camarero, Garbarino, Piermattei, and Urbinati (2017) on cone production. This reinforces the notion that the performance and characteristics of individual trees play an important role in explaining the reproductive output at treeline regardless of the forest structure.
We suggest that seed dispersal is of utmost importance in limiting treeline advance in the Pyrenees. Forest movement upslope will mainly depend on the establishment of seedlings germinating from nearby, scattered trees located at high elevation and rarely from long-distance dispersed seeds (Camarero et al., 2005;Kruse et al., 2019;Nathan, 2006). This may explain why, at the regional scale, ecotone densification (i.e., between the forest and the tree limit) is a more common response to global change in the Pyrenean treelines (Batllori et al., 2010; than treeline advance. Johnson, Gaddis, Cairns, and Krutovsky (2017) argued that seed dispersal was extensive enough in mountain hemlock (Tsuga mertensiana) in Alaska to ensure the advancement to higher elevations. However, this species has an average dispersal distance of 73 m, and has a long-distance dispersal capacity of up to 450 m, values much greater than those reported for Pinus uncinata (Camarero et al., 2005;Vitali et al., 2019). Therefore, the biology of the treeline-forming species may be of paramount importance in their reproductive dynamics, precluding any generalisation in terms of the reproduction-related limiting factors for treeline advance at the global scale. Additionally, seedling recruitment and the eventual treeline advance will not only depend on the arrival of the seeds, but also on the subsequent biotic and abiotic processes affecting seedling establishment and survival Nathan & Muller-Landau, 2000), such as the environmental harshness and scarcity of favourable sites for establishment (Smith, Germino,

TA B L E 5
Results of the zero-inflated negative binomial model for seed rain F I G U R E 4 Seeds collected in the seed traps along the elevation gradient (SDF = subalpine dense forest; SFF = subalpine fragmented forest; KZ = krummholz; AG = alpine grassland).
Data are plotted across the three study sites (Dellui, Amitges and Gelada) and the three study years (2012)(2013)(2014). Capital letters indicate significant differences between elevation positions. The asterisk for AG indicates that no seeds were found in any of the seed traps, and thus it was excluded from the analysis. See Appendix S9 for the number of seeds recovered in the seed traps at each elevation position per study year and site [Colour figure can be viewed at wileyonlinelibrary.com]

| CON CLUS IONS
Our study provides the first evidence of reproductive output as a limiting factor for the advance of the treeline in the Pyrenees. The decrease in seed production along the elevation gradient together with a lack of seed arrival at higher elevations (i.e., alpine grassland) strongly suggests that seed production and dispersal are constraining the ongoing rates of treeline advance in the Pyrenees. By contrast, the lack of clear differences in germination success and germination dynamics from the subalpine forests to the krummholz zone seems to indicate a limited role of genetic constraints at the site level, although genetic variability may have a more important role in treeline dynamics at the regional scale (differential site dynamics). Overall, despite potential positive effects of warming on seed production, our findings suggest that the poor dispersal capacity of Pinus uncinata seeds to the alpine grasslands beyond the treeline will slow down the upslope advancement of the forest. This, together with other climate change-related events such as increased droughts, which may imply higher seedling mortality rates, may reduce the expected climatic sensitivity of treeline position to warming climates.

ACK N OWLED G EM ENTS
We would like to thank Christian Ritz, from the University of Copenhagen, for his statistical advice on meta-analytic randomeffects models. We are grateful to Albert Ferré for providing the map and the picture in Figure 1. We also thank Oriol Grau, Francesc Talavera, Cecília Roma, Xavier Jou, Gerard Fuentes, Andreu Passola, Enric Cuadras and Xavier Trullàs for their help in the field. Finally, we are grateful to Matteo Garbarino and an anonymous reviewer for their helpful comments on the manuscript.

AUTH O R CO NTR I B UTI O N S
MT, JMN, EC and EB conceived the research idea. MT and JMN collected cones, seeds and data and performed the laboratory work.
AAR performed the statistical analyses. AAR wrote the manuscript with contributions from MT and EB. All the authors discussed the results and commented on the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data supporting this publication is stored at https ://github.com/ anadon-rosel l/Anadon-Rosell_et_al_2019_JVS.