This is the submitted version of the article: Bagaria Morató, Guillem; Rodà, F.; Pino i Vilalta, Joan. «Extinction and coloni- sation of habitat specialists drive plant species replacement along a Mediterranean grassland[U+2010]forest succession». Journal of Vegetation Science, Vol. 30, issue 2 (March 2019), p. 331-340. DOI 10.1111/jvs.12722 This version is avaible at https://ddd.uab.cat/record/216994 under the terms of the license 1 Extinction and colonization of habitat specialists drive plant species replacement along a 2 Mediterranean grassland-forest succession. 3 4 Running title: β-diversity along succession 5 6 Guillem Bagaria, Ferran Rodà & Joan Pino 7 8 Bagaria, G. (corresponding author, g.bagaria@creaf.uab.cat)1,2, 9 Rodà, F. (ferran.roda@uab.cat)1,2, 10 Pino, J. (joan.pino@uab.cat)1,2 11 1 12 CREAF, Cerdanyola del Vallès 08193, Catalonia 2 13 Univ Autònoma Barcelona, Cerdanyola del Vallès 08193, Catalonia 14 15 Funding information: 16 Ministerio de Ciencia e Innovación, Spain (CGL2012-33398) 17 Spanish Consolider-Ingenio 2010 programme (CSD2008-00040) 18 Ministerio de Educación, Cultura y Deporte, Spain (AP2009-4599); GB 19 20 1 21 Abstract 22 Questions: Land-use change causes shifts in species richness, which can be delayed. However, 23 beta-diversity patterns and especially the relative role of species replacement and nestedness in 24 these situations with time-lagged extinctions and colonizations remain unknown. We aim to (1) 25 quantify beta-diversity change, species replacement and nestedness for vascular plants along a 26 grassland-forest succession with time-lagged biodiversity change for more than 50 years; (2) check 27 its consistency between all species, grassland specialists and forest specialists, and (3) identify the 28 role of forest encroachment relative to other drivers. 29 Location: Prades Mountains, Catalonia (NE Iberian Peninsula). 30 Methods: We sampled 18 sites representing a gradient in past and current grassland area and 31 connectivity, and in forest encroachment intensity, to obtain plant composition of all species, 32 grassland specialists and forest specialists. We quantified overall beta-diversity and its components 33 at each species classification group along the forest encroachment gradient and other drivers. Then, 34 we used general linear models to study (1) the change rate of beta diversity along the forest 35 encroachment gradient and (2) the relative importance of the drivers in explaining beta diversity. 36 Results: Following forest encroachment gradient, we found an overall noticeable species 37 replacement, while nestedness was the main component for habitat specialists. Landscape 38 differences contributed to explain most compositional differences (both nestedness and 39 replacement), while soil characteristics and geographic distance had a more restricted contribution. 40 Conclusions: Species replacement due to environmental sorting occurred along succession, 41 triggered by selective grassland specialists' extinctions and selective forest specialists' colonizations. 42 Nonetheless, historical landscape characteristics, current landscape characteristics and geographic 43 distance modulate plant extinctions and colonizations, suggesting biological inertia, mass effects 44 and habitat isolation, respectively. Partitioning beta-diversity into nestedness and replacement 45 components and exploring the extinction and colonization patterns of habitat specialists groups 2 46 might provide relevant insights on the drivers and processes of community shift after land-use 47 change. 48 49 Keywords 50 Beta-diversity; Habitat change; Habitat specialists; Nestedness; Species replacement; Vascular 51 plants 52 53 Nomenclature Bolòs et al. (2005) 54 55 β-diversity after forestation 56 57 Introduction 58 Land-use change is one of the main causes of biodiversity changes, driving extinctions and 59 colonizations worldwide, which result in richness gains or losses (Fahrig 2003; Jackson & Sax 60 2010; Baan et al. 2012), species replacement (Fahrig 2003; Jackson & Sax 2010) and biotic 61 homogenization or differentiation (Olden 2006). Recently, there has been a wide interest in 62 investigating both immediate and time-lagged effects of habitat loss and transformation with a focus 63 on species richness, while compositional changes remain much more overlooked (but see Larrea & 64 Werner 2010; Arroyo-Rodríguez et al. 2013; Kopecký et al. 2013; Trentanovi et al. 2013; Lindborg 65 et al. 2014). Biodiversity change and its time-lags, namely extinction debt (Tilman et al. 1994; 66 Kuussaari et al. 2009) and colonization credit (Jackson & Sax 2010; Cristofoli et al. 2010) have 67 been frequently identified following habitat change (for counterexamples see Adriaens et al. 2006; 68 Cristofoli et al. 2010). However, little is known about β-diversity patterns, largely ignoring whether 69 the payment of these debts and credits follows any pattern along the habitat change gradients. 70 Investigating β-diversity patterns is particularly interesting across large environmental changes 3 71 (Mori et al. 2018), e.g. in ecological communities along succession. In this situation, a species 72 replacement is expected, promoted by extinctions of species of the former habitat type and 73 colonizations of species specific to the new habitat type, and resulting in a directional change in 74 composition. Composition differences can only be due to species replacement and nestedness 75 between sites, or a combination of both. Species replacement means a substitution of species by 76 others, while nestedness implies net species gain or loss. Both changes usually occur together but 77 are due to different ecological processes (Legendre 2014) and, then, they must be investigated 78 separately in order to unravel the processes behind composition changes. Baselga (2010) proposed 79 additively partitioning total β-diversity into species replacement and species nestedness between 80 sites, in order to account for the different processes leading to them: (i) environmental sorting or 81 spatial and historical constraints are the processes resulting in a pattern of species replacement, and 82 (ii) selective extinction/colonization, habitat nestedness or interspecific variation in the breadth of 83 environmental tolerance are the processes which result in a pattern of species loss or gain (Leprieur 84 et al. 2011). Species replacement due to environmental sorting has shown as a major process behind 85 total community change along a gradient of environmental conditions (Leibold et al. 2004). 86 However, little is known about the relative contribution and drivers of both species replacement and 87 nestedness to species compositional change following habitat change, especially along succession. 88 During the last century and especially over the past 70 years, there has been an important decline 89 in European semi-natural grasslands, which experienced encroachment by forest (Eriksson et al. 90 2002; Adriaens et al. 2006). In the Mediterranean Basin, semi-natural grasslands also experienced a 91 widespread loss and fragmentation which was predicted to result in a decline of grassland species 92 and a spread of forest species (Debussche et al. 1999), and evidence for a partly paid extinction debt 93 and colonization credit, and a tendency to total richness decrease was already found (Bagaria et al. 94 2015). 95 Regarding habitat specialists, Bagaria et al. (2015) pointed out strong richness gradients which 4 96 suggest that nested species gain or loss might be more important than replacement for explaining 97 changes in specialist species composition. Selective extinctions or colonizations are likely to occur 98 on the basis of the breadth of environmental tolerance of species (Kopecký et al. 2013), with few or 99 no replacement (substitution of some species by others) within each species group along a habitat 100 change. However, Guardiola et al. (2013) suggested that grassland specialists' extinctions in the 101 same habitats were idiosyncratic following habitat loss, thus determining a poorly nested pattern for 102 grassland specialists' composition along the forest encroachment gradient. Forest specialists' 103 colonizations, on the other hand, are expected to occur in a more ordered pattern, as they tend to be 104 faster for very mobile and long-dispersed plants, while plants with large seeds, low fecundity and 105 unassisted dispersal tend to be poor colonizers (Verheyen et al. 2003; Ozinga et al. 2005; Hermy & 106 Verheyen 2007; Svenning et al. 2009). Moreover, both historical and current landscape might 107 modulate species extinctions and colonizations (Bagaria et al. 2015). Forest connectivity, among 108 other factors, enhances forest specialists' colonization (Honnay et al. 2005), but for some species 109 showing limited dispersal, the signal of historical landscape composition might be found after 110 several decades or even more than a century (Vellend et al. 2006; Başnou et al. 2016). If habitat 111 change strongly influences species composition, resulting in a highly nested pattern, a biotic 112 homogenization occurs (Kopecký et al. 2013), but habitat change can also result in extinction and 113 colonization-driven biotic differentiation, through idiosyncratic species gains and losses, when 114 fragmentation is high and seed dispersal is limited (Arroyo-Rodríguez et al. 2013; Lindborg et al. 115 2014). 116 In the present work, we investigated β-diversity patterns regarding species replacement and 117 nested species loss or gain in a forest encroachment process that drives strong habitat change from 118 Mediterranean semi-natural grasslands to forests, which is a representative successional gradient of 119 a common process (Debussche et al. 1999). We selected 18 sites across a grassland-forest gradient 120 resulting from a forest encroachment process that started more than 50 years ago (Bagaria et al. 5 121 2012). Our main aims were: i) to quantify β-diversity and its components (species replacement and 122 nestedness) along the forest encroachment gradient; ii) to assess its consistency between species 123 groups (all species, grassland specialists and forest specialists), and iii) to quantify the relative 124 importance of forest encroachment, soil characteristics, landscape and geographic distance on 125 species replacement and nestedness. 126 127 Methods 128 Study area 129 The study system comprises Mediterranean semi-natural calcareous grasslands and forests resulting 130 from grassland encroachment in Prades Mountains in southern Catalonia (NE Iberian Peninsula; 131 41º14'-41º23'N, 0º56'-1º10'E). These grasslands are dominated by hemicryptophytes and 132 chamaephytes up to 0.5 m high, and vegetation cover ranges from 50% to 80% (see Guardiola et al. 133 2013 for details). A forest-grassland mosaic persisted in the study area until the decade of 1940s, 134 which was grazed by sheep and goats. During the second half of the 20th century, the semi-natural 135 grasslands in the area experienced an important decrease due to livestock grazing reduction that led 136 to forest encroachment (Guardiola et al. 2013), resulting in relatively open, mixed forests of Scots 137 pine (Pinus sylvestris L.) and holm oak (Quercus ilex L.). The study sites are located between 800 138 and 1200 m a.s.l. and have mean annual temperature around 11 °C, mean annual precipitation 139 around 750 mm, and marked summer drought (Digital Climatic Atlas of Catalonia; 140 www.opengis.uab.cat/acdc). 141 Site selection 142 We selected 18 sites, each with a diameter of 200 m, in a set of calcareous plateaus. All sites were 143 located in a small area within the mountains (11 x 7 km) to prevent compositional variation due to a 144 change in climatic conditions. Sampling within sites was conducted without distinguishing habitat 145 patches, because habitat transitions in these Mediterranean systems show a high fuzziness (very 6 146 gradual and irregular transitions), and the delimitation of habitat boundaries is not reliable (Bagaria 147 et al. 2012). Site selection was conducted using historical (1956) and recent (2009) orthophotomaps 148 and other topographic and thematic maps (e.g. wildfires, soil types), according to the following 149 criteria: (1) it was required that site current grassland cover and grassland cover change ranged 150 considerably (Fig. 1); (2) forest encroachment had to be spontaneous (i.e. not due to forest 151 plantation); (3) no crop fields could exist in the past nor nowadays within them; (4) sites could not 152 have been burned in the previous 25 years; (5) slope was required to be lower than 20º, and (6) the 153 distance between them had to be at least 1 km in order to avoid landscape overlap and increase data 154 independence. After preliminary selection over orthophotomaps, a preliminary field survey was 155 used to discard unsuitable areas (see Bagaria et al. 2015 for details). In the 18 selected sites, current 156 grassland cover ranged from 16 to 70%, while grassland cover change varied from 61% loss to 8% 157 increase. 158 Sampling species composition 159 In spring and early summer 2011, all vascular plant species were recorded in 13 circular plots of 3 160 m in diameter at each of the 18 selected sites (see Appendix S1 for details). Plots were arranged in a 161 regular grid with a distance about 33 m between each other. Using regional floras (Rivas-Martínez 162 et al. 2001; Bolòs et al. 2005) and expert advice, all species were classified as grassland specialists, 163 i.e. plants growing mainly in calcareous grasslands and open shrublands; forest specialists, i.e. 164 mainly growing in forests; and other (either generalist species or specialists from other habitats, 165 e.g., ruderal or rocky habitats). From 246 recorded taxa, there were 138 grassland specialists and 54 166 forest specialists (see Appendix S2). In order to study β-diversity changes and colonization and 167 extinction patterns along the forest encroachment gradient, tables of species composition (presence- 168 absence) at site level were obtained for each of the different plant groups (all species, grassland 169 specialists and forest specialists). 170 β-diversity measures 7 171 We calculated β-diversity metrics for species presence-absence between all pairs of sites for the 172 three different species groups (i.e. all species, grassland specialists, forest specialists) in order to 173 investigate changes in community composition along environmental and spatial gradients. β- 174 diversity and its two additive components (nestedness and species replacement) were calculated as 175 proposed by Baselga (2010), using the function 'designdist' in the vegan R package (Oksanen et al. 176 2013). In this framework, for any pair of sites, total β-diversity (βsor) is calculated by the Sørensen 177 dissimilarity index: b+c β sor= 178 2 a+b+c (1) 179 The species replacement component of β-diversity (βsim) is calculated by the Simpson dissimilarity 180 index: min (b , c) β sim= 181 a+min(b , c) (2) 182 Finally, the nestedness-resultant dissimilarity (βsne; a measure of richness differences between 183 nested sites) is formulated as: max (b , c)− min (b , c) a β sne= × 184 2 a+b+c a+min (b , c ) (3) 185 where a, b and c are the shared species between both sites, species unique to the first site, and 186 species unique to the second site, respectively (Krebs 1999),.Total β-diversity is the sum of the two 187 components, replacement and nestedness-resultant dissimilarity; βsor = βsim + βsne. 188 Habitat and landscape characteristics 189 We obtained local-habitat variables that describe site characteristics and the process of forest 190 encroachment, which is expected to be the main driver of plant β-diversity changes through species 191 replacement and nestedness. Soil depth was calculated as the mean of three measures made at each 192 plot and a soil sample up to a depth of 10 cm was taken to analyze soil pH, organic matter content 193 and clay proportion. Solar radiation (MJm-2day-1) was calculated as the product between site spring 8 194 solar radiation (obtained from the Digital Climatic Atlas of Catalonia; www.opengis.uab.cat/acdc) 195 and the proportion of Photosynthetic Active Radiation (PAR) reaching understory plants (measured 196 using a ceptometer at plot centre; AccuPAR LP-80; Decagon Devices, Inc.). Grazing intensity from 197 livestock and/or wild ungulates (mostly roe deer) was measured at each plot, based on excrement 198 groups and eaten plants, in an ordinal way (0, no grazing signs; 1, one grazing sign; 2, more than 199 one grazing sign). Tree cover in a 10-m radius was obtained for each plot as the proportion of forest, 200 from reclassified orthophotomaps of 2009 (pixel size of 1 m). In order to reduce redundancy among 201 local variables related to habitat, a principal components analysis (PCA) on the Pearson correlation 202 matrix was conducted for tree cover, soil depth, pH, organic matter proportion, clay proportion, 203 solar radiation and grazing (taken as quantitative) at plot scale. The first axis of the PCA explained 204 38.6% of the total variance of the table (Fig. 2), and it shows a gradient which was closely related 205 (negatively) with forest encroachment, from forest habitat in the negative side (high tree cover, soil 206 organic matter, clay texture and soil depth) to grassland habitat in the positive side (high solar 207 radiation, soil pH and grazing). The second axis explained 16.6% of the variance and corresponded 208 to a gradient in soil characteristics, from high soil depth and pH on the positive side to high organic 209 matter in the negative side. Then, site centroids (from all plots within a site) for these two axes were 210 used in the statistical analysis as forest encroachment and soil characteristics proxies. Scores from 211 environmental variables on the first and second PCA axes are shown in Appendix S3. 212 In addition, geographic coordinates and connectivity in both the historical and current landscapes 213 were obtained for each site. Landscape connectivity measures were calculated at 500 m radius from 214 each site centre, using maps produced by combining land-cover layers and reclassified 215 orthophotomaps, for both historical (1956) and current (2009) contexts. The combined maps 216 presented fine-scale transitions from grassland to forest (from reclassification of orthophotomaps) 217 and well-defined crop fields, bare soil, and urban areas (from digitized land-cover maps; see 218 Bagaria et al. 2015 for details). These landscape connectivity measures were species-group specific, 9 219 i.e., grassland/forest ratio for all species, grassland ratio (grassland to non-grassland ratio) as 220 connectivity for grassland specialists, and forest ratio (forest to non-forest ratio) as connectivity for 221 forest specialists. 222 Statistical analyses 223 We performed a set of general linear models relating β-diversity to Euclidean distances between 224 sites in (a) forest encroachment (site centroids along axis 1 of the PCA); (b) soil characteristics (site 225 centroids along axis 2 of the PCA); (c) current landscape connectivity (log-transformed to improve 226 normality); (d) historical landscape connectivity (log-transformed), and (c) geographic coordinates. 227 For each of the plant groups (i.e. all species, grassland specialists, forest specialists), and β- 228 diversity measures (total β-diversity [βsor], nestedness-resultant dissimilarity [βsne] and species 229 replacement [βsim]), the composition change rate along forest encroachment (slope of the regression 230 of a β-diversity measure against the forest encroachment gradient) was obtained in order to compare 231 the slope (standardized coefficients) and strength (contribution to total R2) of the relationship 232 between species groups and β-diversity components. 233 To investigate the relative contribution of forest encroachment and the other drivers affecting 234 nestedness or species replacement (i.e. soil characteristics, historical landscape, current landscape 235 and geographic distance), we decomposed total R2 of each model into the contribution of each 236 predictor using relaimpo R package (Grömping 2006) and 'pmvd' metric (Feldman 2005) in the 237 models for the two different β-diversity components. Confidence intervals were calculated using 238 1000 bootstrap repetitions and the bias-corrected and accelerated method. The significance level for 239 each predictor was corrected using error degrees of freedom equal to the number of sites (18), 240 because of the non-independence of pair-wise site comparisons (Qian & Ricklefs 2012). No model 241 selection was performed for these models because the aim was to study the relative contribution of 242 all the included factors to beta-diversity, and to compare it between models. All statistical analyses 243 were conducted using R (R Core Team 2016). 10 244 245 Results 246 β-diversity along the forest encroachment gradient 247 For the whole species composition, a noticeable rate of change along the forest encroachment 248 gradient was found for total β-diversity (standardized coefficient ± SE = 0.48 ± 0.06; R2 249 contribution=0.36) and replacement (0.44 ± 0.06; R2 contribution=0.29), but not for nestedness 250 (0.18 ± 0.08; R2 contribution=0.04). In the case of grassland specialists, the total β-diversity change 251 rate along the forest encroachment gradient was similar to that of all species (0.47 ± 0.06; R2 252 contribution=0.33), but it was mainly due to nestedness (0.35 ± 0.07; R2 contribution=0.19), rather 253 than to species replacement (0.27 ± 0.07; R2 contribution=0.08). For forest specialists, we found 254 rather weak total β-diversity (0.20 ± 0.08; R2 contribution=0.06), but noticeable nestedness (0.25 ± 255 0.08; R2 contribution=0.06) change rates and no offect of species replacement (-0.08 ± 0.08; R2 256 contribution=0.01) along the forest encroachment gradient (Fig. 3). The effect of forest 257 encroachment on β-diversity is significant in all nine models, except for forest specialists' 258 replacement. 259 Drivers of nestedness and species replacement 260 For all species composition, the model for the β-diversity component accounting for species 261 nestedness between sites (βsne) explained 16% of the variance, while historical landscape (R2 262 contribution=0.07), forest encroachment (R2 contribution=0.04) and soil characteristics (R2 263 contribution=0.03) were positively correlated and had a similar contribution to total R2. Both 264 current landscape and geographic distance were not significantly associated with βsne (Table 1; 265 Appendix S4a). As expected, the model for all species replacement (βsim) accounted for higher 266 amount of total variance, 49%. Forest encroachment (R2 contribution=0.29) and current landscape 267 (R2 contribution=0.16) explained the majority of this variance, while geographic distance (R2 268 contribution=0.03), historical landscape and soil characteristics had almost no effect on species 11 269 replacement, although geographic distance was significant (Appendix S4b). All the predictors were 270 positively correlated with βsim. 271 For grassland specialists, the model for the nestedness component of β-diversity (βsne) between 272 sites explained 38% of the variance, and showed that forest encroachment (R2 contribution=0.19), 273 historical landscape (R2 contribution=0.12) and current landscape (R2 contribution=0.04) were the 274 variables significantly and positively associated with species nestedness for this group. Soil 275 characteristics and geographic distance effects were not significant (Appendix S4c). Regarding 276 species replacement (βsim), the model for grassland specialists explained 32% of the variance, and 277 included current landscape as the most important driver (R2 contribution=0.18), followed by forest 278 encroachment (R2 contribution=0.08) and historical landscape (R2 contribution=0.05) (Appendix 279 S4d). All these predictors but historical landscape dissimilarity were positively correlated with βsim. 280 Again, geographic distance and soil characteristics effects were not significant. 281 For forest specialists, the model for the nestedness component of β-diversity (βsne) explained 282 25% of the variance. Current landscape was the best predictor of this component (R2 283 contribution=0.17), followed by forest encroachment (R2 contribution=0.06), and being both 284 positively correlated with βsne. In contrast, soil characteristics, historical landscape and geographic 285 distance were non-significantly associated with βsne (Appendix S4e). Finally, the model for the 286 species replacement component of β-diversity (βsim) only explained 9% of the variance, and only 287 geographic distance (R2 contribution=0.05) was significant, being positively correlated with βsim 288 (Appendix S4f). 289 290 Discussion 291 Our study has explored changes in β-diversity and its components (i.e. species nestedness and 292 replacement) across forest encroachment in Mediterranean grasslands. Species replacement for the 293 whole community is driven by ordered extinctions and colonizations of grassland and forest 12 294 specialists, respectively, along the forest encroachment gradient. This is consistent with the richness 295 loss for grassland specialists and richness gain for forest specialists without total richness change 296 found by Bagaria et al. (2015) at plot scale in the same sites, after 50 years of forest encroachment. 297 Also Kopecký et al. (2013) found a similar pattern of nested extinctions, after forest management 298 cessation, corresponding to light-demanding species. On the other hand, Guardiola et al. (2013) 299 found an idiosyncratic (non-nested) extinction of grassland specialists at local scale in the same 300 region of our study in grassland patches that decreased in size. However, strong habitat 301 transformation similar to the present study had not occurred in their study sites. 302 Forest encroachment, but also landscape connectivity, are the overall primary drivers of β- 303 diversity for both species replacement and nestedness, while soil characteristics only played a 304 significant role in all species nestedness, and geographic distance was significantly related with 305 species replacement of both all species and forest specialists. A small effect of the historical 306 landscape on nestedness of all species was found, indicating that part of the nowadays weak nested 307 composition pattern still relies on grassland/forest ratio 50 years ago. Current landscape markedly 308 contributes, altogether with forest encroachment, to the replacement of total species, suggesting that 309 the replacement of grassland specialists by forest specialists following forest encroachment can be 310 modulated by the surrounding landscape. However, to elucidate the mechanisms under these 311 patterns, the study of nestedness and replacement of grassland and forest specialists is needed. 312 Historical and, to a lesser extent, current landscape distances affect nestedness for grassland 313 specialists, and they complement nestedness for this group along forest encroachment. Hence, there 314 also exists an ordered loss of grassland specialists (nestedness) as grassland prevalence in the 315 landscape decreases. In the same vein, many studies stressed the role of landscape connectivity in 316 promoting richness of grassland specialists in a grassland fragmentation context (e.g. Adriaens et al. 317 2006; Lindborg et al. 2014). On the other hand, the effect of historical landscape suggests that 318 richness is not yet in equilibrium with the current environmental characteristics, indicating the 13 319 existence of an extinction debt (Helm et al. 2006; Kuussaari et al. 2009; Krauss et al. 2010; Jackson 320 & Sax 2010; Bagaria et al. 2015) caused by biological inertia of grassland specialist plants, as found 321 for a grassland specialist species in the study area (Bagaria et al. 2018). Although habitat conditions 322 strongly drive nestedness of grassland specialists between sites, grassland specialists' replacement 323 was more dependent on current landscape than on the forest encroachment gradient, indicating that 324 different grassland specialists appear for different current landscape compositions regardless of 325 habitat characteristics. Plant populations in the surrounding landscapes likely reinforce site 326 populations by mass effects (Mouquet & Loreau 2003; Leibold et al. 2004), providing new 327 individuals to otherwise non-viable populations. Moreover, since species substitutions occur along 328 this gradient, species with different biological characteristics and/or competitive and dispersal 329 abilities might be favoured by different landscape compositions, as shown in previous works 330 (Lindborg 2007; Saar et al. 2012). Historical landscape distance also affects, to a lesser extent, 331 grassland specialists' replacement but with a negative relationship, indicating that sites similar in 332 historical grassland availability show higher replacement. This suggests that sites with very low 333 connectivity in the past are species-poor analogues of sites with high past connectivity, as supported 334 by the nestedness component. In contrast, idiosyncratic extinctions might be occurring depending 335 on other site characteristics (Marini et al. 2013), until the extinction debt will ultimately be paid. 336 Current landscape is the main predictor of nestedness for forest specialists, although forest 337 encroachment is also significant. Thus, forested landscapes might act as propagule sources able to 338 increase forest species number at site scale in a nested way, likely starting with high dispersability 339 plants and ending with those having large seeds, low fecundity or unassisted dispersal (Verheyen et 340 al. 2003). However, this process is expected to last at least some more decades, since ancient forest 341 species are known to rely on the amount of ancient forest in the landscape and not in current forest 342 connectivity (Vellend 2003). Contrary to grassland specialists, no historical effects were found, 343 probably because relevant historical effects act at longer time scales. Forest specialists' replacement 14 344 remains unexplained by forest encroachment or landscape, and only geographic distance has a weak 345 effect on it, indicating a distance decay of similarity (sensu Nekola & White 1999). That might 346 occur because colonization is still in progress even in the most forested sites and landscapes 347 (Jackson & Sax 2010; Bagaria et al. 2015) and dispersal limitation plays an important role 348 (Jamoneau et al. 2012) due to a long history of forest management and exploitation of 349 Mediterranean forests (e.g. Grove & Rackham 2001; Guirado et al. 2007, 2008; Başnou e⁠t al. 2016). 350 In our study system, some biotic homogenization is likely to occur for grassland specialists as 351 forest encroachment increases, increasing similarity in species composition for this group between 352 sites, which is supported by the nestedness of grassland specialists' composition along the forest 353 encroachment gradient. However, replacement, although weaker, is also significant, which 354 advocates for being cautious about this possibility. For forest specialists, the nestedness rate along 355 forest encroachment is weaker, and replacement occurs as distance increases, likely due to habitat 356 isolation, which has been proposed as a cause of biotic differentiation by seed dispersal limitation 357 between populations (Jamoneau et al. 2012; Arroyo-Rodríguez et al. 2013). A specific study would 358 be needed in order to fully investigate homogenization/differentiation patterns and their causes in 359 this system. 360 In conclusion, the patterns of β-diversity analysis for specialist species groups give insight into 361 the processes shaping communities after habitat change. As a result of the forest encroachment 362 process that affected the Mediterranean Basin during the last century (Debussche et al. 1999), a 363 noticeable species replacement occurred with few richness change (few overall nestedness along the 364 forest encroachment gradient), and it was driven by the mainly selective extinction of grassland 365 specialists and the weaker selective colonization of forest specialists, which partly compensated 366 each other. Therefore, an environmental sorting exists for the complete plant community, driven by 367 two species groups with opposite habitat requirements. However, extinction of grassland specialists 368 and colonization of forest specialists are also modulated by landscape characteristics, suggesting 15 369 that they might be delayed or accelerated depending on current habitat connectivity, and even a past 370 connectivity signal remains for grassland specialists. 371 The additive decomposition of β-diversity has proven useful to disentangle community assembly 372 processes (e.g. Leprieur et al. 2011; Marini et al. 2013), and here we depicted the patterns of species 373 replacement and nestedness, and their underlying processes, that drive total compositional change 374 after forest encroachment. Likewise, investigating only the complete community without 375 differentiating specialist species groups would not have provided insight into the processes leading 376 to total species replacement. Although previous studies showed the usefulness of investigating 377 richness gradients for habitat specialists compared to all species after grassland loss and 378 fragmentation (e.g. Adriaens et al. 2006; Öster, Cousins & Eriksson 2007; Öckinger et al. 2012), 379 they were limited to former habitat specialists, and no replacement patterns were investigated. Also, 380 studying the β-diversity patterns for different plant groups has already proven useful to elucidate 381 community assembly mechanisms in a plant invasion context (Marini et al. 2013). Hence, 382 partitioning total β-diversity into its replacement and nestedness components for different species 383 groups allowed clarifying how the extinction and colonization processes occur after land-use 384 change, being selective extinction and colonization the causes under plant replacement for the 385 whole community which, in turn, is driven by environmental sorting. 386 The processes under community change in this system could well hold for other communities 387 experiencing strong habitat transformation due to forestation. The study of a well delimited system 388 and habitat transition allowed us to deal with reduced species pools and environmental variability, 389 since larger geographic gradients might mask diversity patterns (e.g. Guardiola et al., 2013). 390 However, β-diversity patterns might be different if different processes occur, like random 391 extinctions or colonizations resulting in high replacement within habitat specialists’ groups, or 392 strong total richness gain or loss resulting in nestedness for the whole community. Then, more 393 studies are needed to improve the knowledge on the complete process of compositional shift 16 394 (nestedness and replacement) after land-use change in different systems, focusing on the 395 environmental and species characteristics that might accelerate, delay or modify extinctions and 396 colonizations. 397 398 Acknowledgements 399 We acknowledge contributions by M. Guardiola for fieldwork and plant identification support, L. 400 Sáez for aid in plant identification, G. Esparza for fieldwork assistance, J.M. Ninot for advice on 401 plant specialist groups, P. Vicente for digitizing the orthophotomap of 1956 and performing the 402 training points for orthophotomaps' reclassification, R.J. Lewis for advice in performing beta 403 diversity decomposition. 404 405 References 406 Adriaens, D., Honnay, O., & Hermy, M. 2006. No evidence of a plant extinction debt in highly 407 fragmented calcareous grasslands in Belgium. Biological Conservation 133: 212–224. 408 Arroyo-Rodríguez, V., Rös, M., Escobar, F., Melo, F.P.L., Santos, B.A., Tabarelli, M., & Chazdon, 409 R. 2013. Plant β-diversity in fragmented rain forests: testing floristic homogenization and 410 differentiation hypotheses. Journal of Ecology 101: 1449–1458. 411 Baan, L., Alkemade, R., & Koellner, T. 2012. Land use impacts on biodiversity in LCA: a global 412 approach. The International Journal of Life Cycle Assessment 18: 1216–1230. 413 Bagaria, G., Helm, A., Rodà, F., & Pino, J. 2015. Assessing coexisting plant extinction debt and 414 colonization credit in a grassland-forest change gradient. Oecologia 179: 823–834. 415 Bagaria, G., Pino, J., Rodà, F., & Guardiola, M. 2012. Species traits weakly involved in plant 416 responses to landscape properties in Mediterranean grasslands. Journal of Vegetation Science 417 23: 432–442. 418 Bagaria, G., Rodà, F., Clotet, M., Míguez, S., & Pino, J. 2018. Contrasting habitat and landscape 419 effects on the fitness of a long-lived grassland plant under forest encroachment: do they 420 provide evidence for extinction debt? Journal of Ecology 106: 278–288. 421 Baselga, A. 2010. Partitioning the turnover and nestedness components of beta diversity. Global 422 Ecology and Biogeography 19: 134–143. 423 Başnou, C., Vicente, P., Espelta, J.M., & Pino, J. 2016. Of niche differentiation, dispersal ability and 424 historical legacies: what drives woody community assembly in recent Mediterranean forests? 17 425 Oikos 125: 107–116. 426 Bolòs, O., Vigo, J., Masalles, R.M., & Ninot, J.M. 2005. Flora manual dels Països Catalans. Pòrtic, 427 Barcelona, ES. 428 Cristofoli, S., Piqueray, J., Dufrene, M., Bizoux, J., & Mahy, G. 2010. Colonization credit in 429 restored wet heathlands. Restoration Ecology 18: 645–655. 430 Debussche, M., Lepart, J., & Dervieux, A. 1999. Mediterranean landscape changes: evidence from 431 old postcards. Global Ecology and Biogeography 8: 3–15. 432 Eriksson, O., Cousins, S.A.O., & Bruun, H.H. 2002. Land-use history and fragmentation of 433 traditionally managed grasslands in Scandinavia. Journal of Vegetation Science 13: 743–748. 434 Fahrig, L. 2003. Effects of habitat fragmentation on biodiversity. Annual Review of Ecology, 435 Evolution, and Systematics 34: 487–515. 436 Feldman, B. 2005. Relative Importance and Value. 437 Grömping, U. 2006. Relative Importance for Linear Regression in R : The Package relaimpo. 438 Journal of Statistical Software 17: 1--27. 439 Grove, A.T., & Rackham, O. 2001. The nature of Mediterranean Europe: An ecological history. 440 Yale University Press, New Haven, CT, US. 441 Guardiola, M., Pino, J., & Rodà, F. 2013. Patch history and spatial scale modulate local plant 442 extinction and extinction debt in habitat patches. Diversity and Distributions 19: 825–833. 443 Guirado, M., Pino, J., & Rodà, F. 2007. Comparing the role of site disturbance and landscape 444 properties on understory species richness in fragmented periurban Mediterranean forests. 445 Landscape Ecology 22: 117–129. 446 Guirado, M., Pino, J., Rodà, F., & Başnou, C. 2008. Quercus and Pinus cover are determined by 447 landscape structure and dynamics in peri-urban Mediterranean forest patches. Plant Ecology 448 194: 109–119. 449 Helm, A., Hanski, I., & Pärtel, M. 2006. Slow response of plant species richness to habitat loss and 450 fragmentation. Ecology Letters 9: 72–77. 451 Hermy, M., & Verheyen, K. 2007. Legacies of the past in the present-day forest biodiversity: a 452 review of past land-use effects on forest plant species composition and diversity. Ecological 453 Research 22: 361–371. 454 Honnay, O., Jacquemyn, H., Bossuyt, B., & Hermy, M. 2005. Forest fragmentation effects on patch 455 occupancy and population viability of herbaceous plant species. New Phytologist 166: 723– 456 736. 457 Jackson, S.T., & Sax, D.F. 2010. Balancing biodiversity in a changing environment: extinction debt, 458 immigration credit and species turnover. Trends in Ecology & Evolution 25: 153–160. 459 Jamoneau, A., Chabrerie, O., Closset-Kopp, D., & Decocq, G. 2012. Fragmentation alters beta- 18 460 diversity patterns of habitat specialists within forest metacommunities. Ecography 35: 124– 461 133. 462 Kopecký, M., Hédl, R., & Szabó, P. 2013. Non-random extinctions dominate plant community 463 changes in abandoned coppices. Journal of Applied Ecology 50: 79–87. 464 Krauss, J., Bommarco, R., Guardiola, M., Heikkinen, R.K., Helm, A., Kuussaari, M., Lindborg, R., 465 Öckinger, E., Pärtel, M., Pino, J., Pöyry, J., Raatikainen, K.M., Sang, A., Stefanescu, C., Teder, 466 T., Zobel, M., & Steffan-Dewenter, I. 2010. Habitat fragmentation causes immediate and time- 467 delayed biodiversity loss at different trophic levels. Ecology Letters 13: 597–605. 468 Krebs, C. 1999. Ecological methodology. Addison-Wesley, California, US. 469 Kuussaari, M., Bommarco, R., Heikkinen, R.K., Helm, A., Krauss, J., Lindborg, R., Öckinger, E., 470 Pärtel, M., Pino, J., Rodà, F., Stefanescu, C., Teder, T., Zobel, M., & Steffan-Dewenter, I. 2009. 471 Extinction debt: a challenge for biodiversity conservation. Trends in Ecology & Evolution 24: 472 564–571. 473 Larrea, M.L., & Werner, F.A. 2010. Response of vascular epiphyte diversity to different land-use 474 intensities in a neotropical montane wet forest. Forest Ecology and Management 260: 1950– 475 1955. 476 Legendre, P. 2014. Interpreting the replacement and richness difference components of beta 477 diversity. Global Ecology and Biogeography 23: 1324–1334. 478 Leibold, M.A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J.M., Hoopes, M.F., Holt, R.D., 479 Shurin, J.B., Law, R., Tilman, D., Loreau, M., & Gonzalez, A. 2004. The metacommunity 480 concept: a framework for multi-scale community ecology. Ecology Letters 7: 601–613. 481 Leprieur, F., Tedesco, P.A., Hugueny, B., Beauchard, O., Dürr, H.H., Brosse, S., & Oberdorff, T. 482 2011. Partitioning global patterns of freshwater fish beta diversity reveals contrasting 483 signatures of past climate changes. Ecology Letters 14: 325–334. 484 Lindborg, R. 2007. Evaluating the distribution of plant life-history traits in relation to current and 485 historical landscape configurations. Journal of Ecology 95: 555–564. 486 Lindborg, R., Plue, J., Andersson, K., & Cousins, S.A.O. 2014. Function of small habitat elements 487 for enhancing plant diversity in different agricultural landscapes. Biological Conservation 169: 488 206–213. 489 Marini, L., Bertolli, A., Bona, E., Federici, G., Martini, F., Prosser, F., & Bommarco, R. 2013. Beta- 490 diversity patterns elucidate mechanisms of alien plant invasion in mountains. Global Ecology 491 and Biogeography 22: 450–460. 492 Mori, A.S., Isbell, F., & Seidl, R. 2018. β-Diversity, Community Assembly, and Ecosystem 493 Functioning. Trends in Ecology & Evolution. doi: 10.1016/J.TREE.2018.04.012 494 Mouquet, N., & Loreau, M. 2003. Community patterns in source-sink metacommunities. American 495 Naturalist 162: 544–557. 19 496 Nekola, J.C., & White, P.S. 1999. The distance decay of similarity in biogeography and ecology. 497 Journal of Biogeography 26: 867–878. 498 Öckinger, E., Lindborg, R., Sjödin, N.E., & Bommarco, R. 2012. Landscape matrix modifies 499 richness of plants and insects in grassland fragments. Ecography 35: 259–267. 500 Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R., O’Hara, R.B., Simpson, G.L., 501 Solymos, P., Stevens, M.H.H., & Wagner, H. 2013. vegan: Community ecology package. R 502 package version 2.0-10. 503 Olden, J.D. 2006. Biotic homogenization: A new research agenda for conservation biogeography. 504 Journal of Biogeography 33: 2027–2039. 505 Öster, M., Cousins, S.A.O., & Eriksson, O. 2007. Size and heterogeneity rather than landscape 506 context determine plant species richness in semi-natural grasslands. Journal of Vegetation 507 Science. doi: 10.1658/1100-9233(2007)18[859:SAHRTL]2.0.CO;2 508 Ozinga, W.A., Schaminée, J.H.J., Bekker, R.M., Bonn, S., Poschlod, P., Tackenberg, O., Bakker, J., 509 & Van Groenendael, J.M. 2005. Predictability of plant species composition from 510 environmental conditions is constrained by dispersal limitation. Oikos. doi: 10.1111/j.0030- 511 1299.2005.13632.x 512 Qian, H., & Ricklefs, R.E. 2012. Disentangling the effects of geographic distance and 513 environmental dissimilarity on global patterns of species turnover. Global Ecology and 514 Biogeography 21: 341–351. 515 R Core Team. 2016. R: A language and environment for statistical computing. 516 Rivas-Martínez, S., Fernández-González, F., Loidi, J., Lousa, M., & Penas, A. 2001. 517 Syntaxonomical checklist of vascular plant communities of Spain and Portugal to association 518 level. Itinera Geobotanica 14: 5–341. 519 Saar, L., Takkis, K., Pärtel, M., & Helm, A. 2012. Which plant traits predict species loss in 520 calcareous grasslands with extinction debt? Diversity and Distributions 18: 808–817. 521 Svenning, J.-C., Baktoft, K.H., & Balslev, H. 2009. Land-use history affects understorey plant 522 species distributions in a large temperate-forest complex, Denmark. Plant Ecology 201: 221– 523 234. 524 Tilman, D., May, R., Lehman, C.L., & Nowak, M.A. 1994. Habitat destruction and the extinction 525 debt. Nature 371: 65–66. 526 Trentanovi, G., Lippe, M., Sitzia, T., Ziechmann, U., Kowarik, I., & Cierjacks, A. 2013. Biotic 527 homogenization at the community scale: disentangling the roles of urbanization and plant 528 invasion. Diversity and Distributions 19: 738–748. 529 Vellend, M. 2003. Habitat loss inhibits recovery of plant diversity as forests regrow. Ecology 84: 530 1158–1164. 531 Vellend, M., Verheyen, K., Jacquemyn, H., Kolb, A., Van Calster, H., Peterken, G., & Hermy, M. 20 532 2006. Extinction debt of forest plants persists for more than a century following habitat 533 fragmentation. Ecology 87: 542–548. 534 Verheyen, K., Honnay, O., Motzkin, G., Hermy, M., & Foster, D.R. 2003. Response of forest plant 535 species to land-use change: a life-history trait-based approach. Journal of Ecology 91: 563– 536 577. 537 538 539 Supporting Information 540 Additional Supporting Information may be found in the online version of this article: 541 542 Appendix S1. Details on sampling design procedure. 543 Appendix S2. List of plant species found in the survey. 544 Appendix S3. Environmental variables scores on PCA axes. 545 Appendix S4. Figure showing variance decomposition on models for each species group and β- 546 diversity component. 21 547 Table 1. Standardized coefficients and their standard error for each predictor in the linear model for 548 each combination of species group (all species, grassland specialists, forest specialists) and β- 549 diversity component (total, βsor; richness differences, βsne; species replacement, βsim) between pairs 550 of sites. Total R2 of each model is given and significant predictors at P<0.05 are marked in bold 551 (*** P<0.001, ** P<0.01, * P<0.05). 552 Species β-diversity Predictor Total group index ΔForest ΔSoil ΔCurrent ΔHistorical Geographic R2 encroachment characteristics landscape landscape distance All species βsor 0.48±0.06*** 0.09±0.05 0.35±0.06*** 0.12±0.05* 0.19±0.05** 0.58 βsne 0.18±0.08* 0.17±0.08* -0.04±0.08 0.23±0.08* 0.12±0.08 0.16 βsim 0.44±0.06*** 0.05±0.06 0.37±0.06*** 0.06±0.06 0.16±0.06* 0.49 Grassland βsor 0.47±0.06*** 0.10±0.05 0.42±0.05*** 0.06±0.05 0.14±0.05* 0.61 specialists βsne 0.35±0.07*** 0.12±0.07 0.18±0.07* 0.31±0.07*** 0.10±0.07 0.38 βsim 0.27±0.07** 0.02±0.07 0.39±0.07*** -0.26±0.07** 0.08±0.07 0.32 Forest βsor 0.20±0.08* -0.01±0.08 0.30±0.08** 0.00±0.08 0.22±0.08* 0.20 specialists βsne 0.25±0.08** -0.12±0.07 0.37±0.07*** -0.08±0.07 0.01±0.07 0.25 βsim -0.08±0.08 0.12±0.08 -0.12±0.08 0.09±0.08 0.19±0.08* 0.09 553 22 554 Figure 1. Relation between current site grassland cover and site grassland cover change (current 555 minus historical grassland cover). Both measures are only weakly and non-significantly correlated 556 (R2 = 0.18; P = 0.08). 557 23 558 Figure 2. Plot corresponding to the first and second axes from the PCA of habitat characteristics in 559 the 234 sampled plots performed on the correlation matrix. Axis 1 explained 38.6% of the total 560 variance and corresponds to a habitat gradient negatively related with forest encroachment, from 561 forest in the negative side to grassland in the positive side. Axis 2 accounted for 16.6% of the 562 variance and is consistent with a gradient in soil characteristics, from high soil depth and pH on the 563 positive side to high organic matter in the negative side. Circles show sampled plot positions, while 564 arrows correspond to variables. 565 24 566 Figure 3. β-diversity gradients along the forest encroachment distance gradient for the different 567 components: (a) total β-diversity, βsor; (b) nestedness, βsne, and (c) species replacement, βsim. The 568 three species groups are depicted with different symbols: all species (dots, solid line), grassland 569 specialists (triangles, dotted line) and forest specialists (crosses, dashed line). β-diversity values are 570 partial residuals, and each symbol is a pair of sites. 571 25 572 Figure 1 573 26 574 Figure 2 575 27 576 Figure 3 577 578 28