Frontiers of Biogeography Investigating elevational gradients of species richness in a Mediterranean plant hotspot using a published flora

The Apuan Alps are one of the most peculiar mountain chains in the Mediterranean, being very close to the coastline and reaching an elevation of almost 2000 m. Based on published flora, we investigated the distribution of plant species richness along the whole elevational gradient of this chain considering: (i) native species, (ii) endemic versus alien species; and (iii) functional groups of species based on Raunkiær life forms (RLF). Generalized Linear Models (GLMs) were used to analyse richness patterns along the elevational gradient, and elevational richness models versus the area of the elevational belts were fitted to test the effect of surface area. Our results showed decreasing species richness with increasing elevation. In contrast, endemic species richness increased along the elevational gradient. Alien species were mainly distributed at low elevations, but this result should be taken with caution since we used historical data. Species life forms were not equally distributed along the elevation gradient: chamaephytes and hemicryptophytes were the richest groups at high elevations, while therophytes showed highest species richness at low elevations. Our findings suggest that in the Apuan Alps there is a major elevational gradient in species composition that could reflect plant evolutionary history. Furthermore, we highlight the key role of published floras as a relevant source of biodiversity data. This paper is part of an Elevational Gradients and Mountain Biodiversity Special Issue.


Introduction
Patterns of plant species richness have always been a topic of interest in biogeography and ecology (Pianka 1966, Huston and Huston 1994, Lomolino 2001, Whittaker et al. 2001. Historically, species richness trends have been associated with available area and environmental gradients. An increase in species richness with increasing area is one of the most robust ecological generalizations (Rosenzweig 1995, Lomolino 2000, and a vast literature has been published on species−area relationships (e.g., Karger et al. 2011, Dengler et al. 2019, Moradi et al. 2020, Matthews et al. 2021. Species richness patterns have also been evaluated along several environmental gradients (Gaston 2000), such as temperature, productivity, potential evapotranspiration, etc. However, the two most investigated environmental gradients related to species richness patterns are latitude and elevation.
The elevation gradient of species richness has been commonly explained by the same factors driving the latitudinal gradient of species richness, such as climate, productivity, and other energy-related processes (Richerson and Lum 1980, Turner et al. 1987, Currie 1991, Wright et al. 1993, Panda et al. 2017, Vetaas et al. 2019. Indeed, increase in elevation leads to changes in several environmental factors, such as temperature, precipitation, evapotranspiration, soil nutrient availability, and solar radiation, potentially driving species richness. For this reason, since the pioneering studies of Alexander von Humboldt, analyses of elevational gradients have been pivotal to disentangling the causes behind broad-scale patterns of biodiversity (Fattorini et al. 2019). However, despite the large number of studies investigating elevational patterns of species richness, many questions still remain unresolved (Wang et al. 2017), and the factors determining the elevation−richness relationship and its variations across spatial scales, clades, and biomes still require clarifications. Understanding how elevational gradients affect species richness is both of theoretical and practical concern, especially in a climate change scenario, since it links biodiversity patterns to conservation (Rahbek 1995, Lomolino 2001, Thuiller et al. 2005, Loarie et al. 2009).
Within this theoretical and analytical framework, published Floras represent a highly valuable source of data for ecological and biogeographical analyses (König et al. 2019). This view is corroborated by several papers focused on island biogeography theory (Wilson 1988, Chiarucci et al. 2017), land-use changes (Aggemyr and Cousins 2012, Finderup Nielsen et al. 2019, biogeographical networks (Kougioumoutzis et al. 2014, Torre et al. 2019, large-scale drivers of species richness (Weigelt et al. 2020), and patterns of species richness along elevational gradients (Kluge et al. 2017).
In the present study, we used the detailed data on the local elevational distribution of more than 2,000 plant taxa included in a published Flora of the Apuan Alps (Central Italy) to model patterns of species richness for: (i) native species, (ii) endemic and alien species separately, and (iii) functional groups of species based on Raunkiaer life forms (RLF; Raunkiaer 1934). Specifically, we investigated how species richness varies along the elevational gradient, both for the whole flora and for endemic and alien species, as well as for different Raunkiaer's life forms. Additionally, we assessed whether or not the species richness-elevation relationship may be influenced by differences in available surface area among elevational belts. Finally, we explored the elevational ranges of endemics, aliens, and different Raunkiaer's life forms.

Study Area
The Apuan Alps are a mountain chain located in north-western Tuscany (Italy, Fig. 1), belonging to the Mediterranean biogeographical region (Cervellini et al. 2020). This mountain chain is characterized by steep slopes, especially on the south-west facing side, and broad rocky outcrops reaching up to 1946 m. a.s.l. at Monte Pisanino. The harsh morphology and the substrate variability of these mountains, made up of metamorphic and sedimentary rocks, shape local habitat heterogeneity (Tomaselli and Agostini 1994) that combined with recent geomorphologic and neotectonic events likely drive current plant biodiversity (Bedini et al. 2011).
The lowlands of the southern mountain side are dominated by evergreen sclerophyllous plants, which are replaced by Carpinus and Quercus forests at mid-elevations (Tomaselli et al. 2019). On the Fig. 1. Study area. The coloured area corresponds to the Apuan Alps according to Ferrarini and Marchetti 1994. Dark blue to light yellow represents the elevation gradient. In the bottom-left corner, we show the position of the study area within a map of Italy.  (Tomaselli and Agostini 1994). These areas host moorlands on the siliceous peaks, while the calcareous peaks are typically dominated by casmophytic (cliff dwelling) plant communities.

Dataset building
The floristic dataset used for the analyses was built by recording each taxon ("species" hereafter) listed in the three volumes of Prodromo alla Flora della Regione Apuana (Ferrarini and Marchetti 1994, Ferrarini et al. 1997, Ferrarini 2000, and its elevational distribution (minimum and maximum elevation at which each species is recorded for the area) and chorotype (native, endemic or alien, in all cases relatively to Apuan Alps). Taxonomical nomenclature was standardised according to Pignatti et al.(2017) and RLF were assigned to each species. When more than one RLF was listed for a species we retained the more durable RLF (e.g., between phanerophyte and chamaephyte we retained the former, and between therophyte and hemicryptophyte we retained the latter). Finally, we excluded all the species whose elevation distribution was uncertain (n=116), as well as all those species reported as extinct or cultivated (n=3 and n=32, respectively). Thus, 1,952 species were retained in the final dataset. Lastly, since hydrophytes and helophytes were missing from most of the elevational belts, occurring only below 500 m, we did not account for these RLFs in the models.

Data analysis
We used species richness estimated in each 50-m elevation belt to model the elevational patterns of plant richness in the study area. A total of 39 belts were obtained, ranging from 0 to 1950 m. In the text and plots, the upper elevation limit was used to indicate the corresponding belt, both in the text and in the plots. For each belt, we calculated the planar area and the overall species richness, the richness of endemic and alien species, and that of the seven RLF categories.
In order to graphically explore the elevational patterns of species richness, we plotted overall species richness and species richness for each category against the elevational belts, both as absolute values (i.e., counts of species) and as relative species richness (% of total richness).
We applied polynomial linear regression using Generalized Linear Models (GLMs) with a Poisson or Quasipoisson (in case of overdispersion) distribution of error to investigate the variation of species richness along the elevational gradient. In case of overdispersion, we used negative binomial regression models. In the same way, the relative richness of different chorotypes and RLF was modelled using GLMs with quasibinomial distribution of error. The best model among linear, quadratic and cubic was assessed based on Akaike Information Criterion (AIC) and parsimony criterion. The most parsimonious model was chosen if the decrease in AIC when adding the next order term was small. In particular, we assumed that models in which the difference in AIC is < 2 can be considered as equally supported (Burnham 1998, Burnham et al. 2011. It has been widely reported that species richness increases as a function of area Grytnes 2007, Lee et al. 2013). Therefore, to assess the species− area relationship we fitted linear models between the residuals of each GLM model as the response variable and the area of each belt as the predictor variable.
To investigate the elevational distribution of endemic and alien species as well as of RLF, we compared the lower and upper limits of elevational distribution and elevational ranges of the species. As the distributions of lower and upper limits, as well those of elevational ranges, were highly skewed, we used non-parametric tests. We compared the elevational distributions of alien and endemic species by means of Wilcoxon rank sum test. For the RLF, we first tested the differences among all groups by means of Kruskal-Wallis test. As we found significant differences, we applied pairwise Wilcox rank sum as a post-hoc test, by adjusting the p-values with a Holm correction to prevent Type I error. All the analyses were performed using in R v. 3.6.3 (R Core Team 2020) and the following R packages: sars (Matthews et al. 2019), reshape (Wickham 2007), ggplot2 (Wickham 2016), sf (Pebesma 2018), raster (Hijmans 2020) and tmap (Tennekes 2018), and MASS (Venables and Ripley 2002).

Results
Among the 1,952 species, 14 endemic or nearly endemic species and 22 alien species were included. Thus, endemics represented 0.72% and aliens 1.13% of the total flora.
Total species richness decreased along the elevational gradient (Fig. 2a), with about 58.35% of the total flora being found at the lowest elevational belt and 2.77% at the highest elevational belt (Fig. 3a). The model accounting for native species showed a decrease (D 2 =0.94) of species richness. In the lowest elevational belts more than 1,000 species were found, while at the highest belts less than 150 species occurred (Table 1).

Endemic and alien species
The endemic and alien species showed contrasting trends. The former had a unimodal positive pattern along the elevational gradient (Fig. 2b), peaking around 1200 m a.s.l. (D 2 = 0.74), whereas alien species displayed a less definite trend, with a drastic decrease of species richness starting from 200 m a.s.l. and a complete lack of recorded alien species from the 1,050 m belt up to the 1,400 m belt (D 2 = 0.75; Table 1).
The relative species richness along the elevational gradient for endemic species showed a positive relationship, with an increase of 5% along the gradient (D 2 = 0.92). In contrast, alien species did not reveal a clear trend and the relative richness never exceeded 1.2% (Table 2). Both groups fitted to a quadratic model, but alien species showed a downward convexity (D 2 = 0.37; Fig. 3b).

Raunkiaer life forms
The sharpest trends were observed for hemicryptophytes and therophytes as the majority of the species belong to these two categories (Fig. 2c). Hemicryptophyte species richness showed a bellshaped trend, with a peak at about 500 m a.s.l. (D 2 = 0.92), while therophyte species richness had a negative exponential pattern, with a rapid decrease from the sea level up to 1,200 m a.s.l., and a less steep decline at higher elevations (D 2 = 0.99).
Phanerophytes, geophytes, and chamaephytes displayed weaker trends. Phanerophyte and geophyte species richness showed an approximately linear decrease with elevation (D 2 phanerophyte = 0.99, D 2 geophyte = 0.91), while chamaephytes showed a weak increase in species richness with a plateau between the 500 m and the 1,000 m belts and decreased up to the summit (D 2 = 0.88; Table 1 and supplementary Table S1).

Effect of area
A strong negative correlation (R 2 0.7) between elevation and area of the belts was found. The linear model between the residuals of each elevational model versus belt areas showed no significant effect of the area in almost all cases. Significant effect of the area was found only for endemic species even if the linear model had an adjusted R 2 of 0.19. This result is due to the effect of the lowest elevational belt (see supplementary Figure S1), which has the largest area and, in this case, the highest residual value. All the other models had a p-value higher than 0.2 and  Table 1. Summary statistics of fitted polynomial GLMs with vascular plant species richness as a function of 50-metres elevational belts in the Apuan Alps. For all models we employed a Poisson distribution of error, except for native species where we used a negative binomial distribution of error. "Int" is the intercept, "Elev" and "Elev 2 " are first-and second-order angular coefficients. D 2 is the percentage of deviance explained by the models.  Table 2. Summary statistics of fitted polynomial quasibinomial GLMs with proportional species richness as a function of 50-metres elevational belts. "Int" is the intercept, "Elev" and "Elev 2 " are first-and second-order angular coefficients. D 2 is the percentage of deviance explained by the models. Proportions are calculated on overall species richness per belt. an adjusted R 2 lower than 0.1 (see supplementary Figure S1).

Elevational range
Endemic species showed heterogeneous elevational ranges. Indeed, some species were distributed along the entire gradient, while others had a narrow range, being centered at about 1,800 m a.s.l. On the other hand, alien species showed similar range sizes (ca. 700 m width) and almost all the species had their minimum height of occurrence near the sea level. All the comparisons between endemic and alien species showed significant differences ( Fig. 4A; p < .001). Therefore, endemic species had wider elevational ranges and higher minimum and maximum elevation of occurrence compared to alien species.
The RLF range sizes showed a strong variability among groups (Fig. 4b), which differed significantly in all comparisons (p < .001). The most common pattern was the minimum height of occurrence, indeed all the RLF included a large number of species that had their lowest limit near the sea level. As mentioned before, hydrophyte and helophyte species had very narrow elevational ranges. Within the remaining categories species were distributed very differently. Therophytes and geophytes, on average, had smaller elevational ranges, with therophytes showing the lowest elevational occurrences both for minimum and maximum limits. Chamaephytes showed the widest elevational ranges while most of the hemicryptophytes and phanerophytes had their elevational ranges varying from 500 m to 1,250 m. On the other hand, hemicryptophytes and chamaephytes showed the highest minimum limits with median maximum limits above 1,000 m a.s.l.

Discussion
Our study area hosts over 25% of the Italian flora (Bartolucci et al. 2018), supporting the view that this mountain chain can be regarded as a hotspot of plant diversity (Gentili et al. 2015, Carta et al. 2019. The large number of species is supported by the great habitat heterogeneity of the Apuan Alps (Tomaselli and Agostini 1994). Indeed, high environmental heterogeneity provides great niche diversity and facilitates species coexistence, resulting in greater species richness (López-González et al. 2015). The high number of endemics may be associated with the role of Apuan Alps as a refuge during the last glacial maximum (Médail andDiadema 2009, Gentili et al. 2015), indeed these data highlight the great biogeographical importance of the Apuan Alps for generating and preserving plant evolutionary history (Carta et al. 2019).
The observed pattern along the elevational gradient confirms the hypothesis of a linear decrease in species richness. It is currently considered that the most common pattern of plant species richness along an elevational gradient is a hump at middle heights (Colwell et al. 2004, Grytnes andMcCain 2007). However, our gradient spans less than 2,000 metres in elevation and it is not uncommon to find linear trends along relatively short gradients (Nogués-Bravo et al. 2008). Moreover, elevational gradients of species richness were usually investigated using standardised plots with equal area Vetaas 2003, Lee et al. 2013). In contrast, we used data from a published flora listing all the recorded species for each elevational belt and the effect of area could have potentially emerged. However, we note that, with the exception of the lowest one, there was only slight variation in area size across the elevational belts. Moreover, the results of linear models carried out to investigate the effect of area suggest, as in this case, that area has a secondary and not significant effect on species richness.
The decrease of species richness along the gradient suggests that environmental constraints related to mountain ecosystems may act as filters for species distributions (Körner 2003). Nevertheless, we found that the percentage of endemic species increases along the elevational gradient and their elevational ranges showed higher median values than other groups, confirming the hypothesis that endemics are more commonly found at higher elevations (Vetaas andGrytnes 2002, Steinbauer et al., 2016). These results can be explained by the biogeographical effect of the mountain region, which can be considered as an ecological archipelago for evolutionary processes (Körner 2003, AEgisdóttir et al. 2009, Sklenář et al. 2014, Rahbek et al. 2019). However, these crucial areas for biodiversity conservation have been threatened by direct and indirect effects of human activities in the last decades (IPCC 2019). Moreover, high-mountain areas have been experiencing some of the greatest increases in temperature in the last half century (Pepin et al. 2020). Thus, these habitats, extremely important for species conservation, are among the most vulnerable to current climate and land use changes (Walther et al. 2005). The endemic or rare species that are found at high elevations are more likely to be threatened by local extinctions in the near future (Thuiller et al. 2005, Pauli et al. 2007, Bellard et al. 2014. Thus, it is important to improve conservation strategies, particularly for endemic and threatened species growing in the Apuan Alps (IUCN 2020). Phylogeographic studies on endemics should be implemented to assess the genetic diversity of these species and the refuge role of the Apuan Alps in the last glacial maximum (Bedini et al. 2011). Alien species richness peaked at the lowest elevations, where human impact is stronger and human population density is higher. Currently, the greatest number of alien species in the world and the greatest knowledge about the extent of invasions are reported in economically developed countries (Pyšek et al. 2008, McGeoch et al. 2010, Early et al. 2016, Corcos et al. 2020. Due to an interplay between climate warming and direct human impact, an upward shift of the alien species is expected in the near future (Walther et al. 2009, Alexander andEdwards 2010). This shift could threaten vulnerable habitats that host endemic species, potentially leading these taxa to local extinctions (Bellard et al. 2016, Blackburn et al. 2019. The elevational patterns of species richness for different RLF show how the plant life strategies shape their elevational distribution. As demonstrated by previous studies on the dominant role of hemicryptophytes and chamaephytes in high mountain ecosystems (Körner 2003, Matteodo et al. 2013, Irl et al. 2020), we observed a positive relationship between the percentages of hemicryptophytes and of chamaephytes and elevation. This is most likely due to the climatic and topographic features of the Apuan Alps. Indeed, the highest peaks of these mountains are characterized by steep slopes and calcareous plateaus with thin soils colonised by stress tolerant communities. On the other hand, the proportion of therophytes showed that this group is dominant in lowland habitats, probably due to the high human impact occurring therein. Annuals (sometimes biannual), which become quite rare at high elevations (Körner 2003), commonly do not represent more than 2% of the total mountain flora, and become increasingly rare with increasing elevation (Bliss 1971, Billings 1974. Phanerophytes had relatively wider elevational ranges, but their distribution was limited by low seasonal mean ground temperature: a value of 6.7 °C mean temperature during the growing season is usually considered as the thermal threshold for forest growth at high elevation (Körner and Paulsen 2004). A linear decrease of phanerophyte species richness with increasing elevation has been recorded at the plot scale for the forest plant communities of Tuscany (Chiarucci and Bonini 2005), confirming the marked dependence of this RLF type on elevation.
Finally, although the use of published floras has been shown potentially to help disentangle ecological and biogeographical patterns and processes, it is not free of bias. One major limitation could be the underestimation of elevational gradient for rare species. Indeed, these species could actually show distributional gaps, both for their real absence or due to the sampling effort. Even if these limitations could affect the results, the great number of species with a lower elevational limit near the sea level and the relatively small areas of the higher belts ensure the reliability of the sampling effort. Moreover, we note that the Apuan Alps flora comes from a life's work with a great sampling effort distributed over more than 20 years (see Ferrarini and Marchetti 1994).

Conclusions
We highlighted a general decrease of plant species richness along the elevational gradient of the Apuan Alps. The observed trends are likely due to the increase of environmental constraints along the elevational gradient (Körner 2003). On the other hand, the increase in percentage of endemics at high elevations highlighted the key role of mountain ecosystems in shaping evolutionary processes (Vetaas andGrytnes 2002, Steinbauer et al. 2016). Moreover, our results confirmed that published floras provide useful information to investigate species diversity as well as biogeographical patterns and their ecological drivers (Aggemyr and Cousins 2012; Chiarucci et al. 2017;Finderup Nielsen et al. 2019).

Author Contributions
A.Ch. conceived the idea and coordinated the study; the database was designed and compiled by C.F., later expanded by L.S.; P.Z. and M.D.M. performed the analyses and produced the outputs; M.D.M. lead the writing of the manuscript with substantial contributions from P.Z.; all authors discussed the results and contributed to drafting the final manuscript.

Data Accessibility
Data are deposited on Zenodo and are available upon motivated request (doi:10.5281/zenodo.4419941).

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Table S1. Vascular plant species richness per 50-metres elevational belt in the Apuan Alps. Figure S1. Scatterplots of poisson GLMs residuals (species richness as function of elevation) vs. 50-metres elevational belts' areas.