Increase in CO 2 concentration could alter the response of Hedera helix to climate change

Abstract Increasing CO 2 concentration ([CO 2]) is likely to affect future species distributions, in interaction with other climate change drivers. However, current modeling approaches still seldom consider interactions between climatic factors and the importance of these interactions therefore remains mostly unexplored. Here, we combined dendrochronological and modeling approaches to study the interactive effects of increasing [CO 2] and temperature on the distribution of one of the main European liana species, Hedera helix. We combined a classical continent‐wide species distribution modeling approach with a case study using H. helix and Quercus cerris tree rings, where we explored the long‐term influence of a variety of climate drivers, including increasing [CO 2], and their interactions, on secondary growth. Finally, we explored how our findings could influence the model predictions. Climate‐only model predictions showed a small decrease in habitat suitability for H. helix in Europe; however, this was accompanied by a strong shift in the distribution toward the north and east. Our growth ring data suggested that H. helix can benefit from high [CO 2] under warm conditions, more than its tree hosts, which showed a weaker response to [CO 2] coupled with higher cavitation risk under high temperature. Increasing [CO 2] might therefore offset the negative effects of high temperatures on H. helix, and we illustrate how this might translate into maintenance of H. helix in warmer areas. Our results highlight the need to consider carbon fertilization and interactions between climate variables in ecological modeling. Combining dendrochronological analyses with spatial distribution modeling may provide opportunities to refine predictions of how climate change will affect species distributions.


| INTRODUC TI ON
Ecosystems are changing rapidly in response to increasing temperature, altered precipitation patterns, and increasing CO 2 concentrations ([CO 2 ]), all of which are influencing the abundance and distribution of many plant and animal species across the globe (Lindner et al., 2014). Positive carbon fertilization effects have been shown in experiments (Bader et al., 2013), but it is unclear whether increased carbon availability could affect the growth of natural populations. Additionally, climate change drivers may show complex long-term interactions that can modify their effect on plant performance and species interactions (e.g. Lindner et al., 2014). However, the main tool to model species response to climate change at large scales, spatial distribution models (SDMs), seldom account for increased [CO 2 ], or variable interactions, likely because this is methodologically challenging (Norby & Luo, 2004). As a consequence, it is unclear whether carbon fertilization might interact with other climate change drivers to influence the outcome of climate models.
Lianas are common in forests across the world and influence forest ecosystem functioning (Schnitzer & Bongers, 2002;Tymen et al., 2016) and carbon storage capacity (van der Heijden, Powers, & Schnitzer, 2015). They affect key ecological processes such as tree mortality, susceptibility to wind damage, compositional turnover, and species diversity (Allen, Sharitz, & Goebel, 2007;Körner, 2004;Ladwig & Meiners, 2010;Schnitzer & Bongers, 2002. Lianas may also be one group that could strongly benefit from climate change and they significantly increased in abundance and productivity in the last decades (Phillips et al., 2002;Schnitzer, 2015;Schnitzer & Bongers, 2002). The drivers behind this trend remain contentious. Liana increase was first reported in the neotropics (Phillips et al., 2002) but since then, contrasting results have been shown in temperate (Heinrichs & Schmidt, 2015) and African forests (reviewed in Schnitzer, 2015), where there have been no long-term changes in liana abundance. European lianas, on the other hand, have attracted limited attention, other than as invasive species in North American forests (Heinrichs & Schmidt, 2015). This is likely due to the low abundance and economic impact of European liana species in forestry, to the point that they have even been dismissed as "mainly decorative" (Silvertown, 2008). However, lianas are still an important part of European forest biodiversity and their loss from forests could have cascading effects on the biodiversity of other groups. Therefore, understanding the response of European lianas to climate change is important to predict future trends in forest biodiversity.
The most important and widespread of the liana woody species in Europe is English ivy, Hedera helix. There have been a handful of works studying H. helix's growth, effects on host performance, and host preference (Nola 1997, Schnitzer & Bongers, 2002Garfi & Ficarrotta, 2003;Heuzé, Dupouey, & Schnitzler, 2008;Castagneri, Garbarino, & Nola, 2013); however, it yet remains unknown how H. helix might respond to the new environmental conditions brought by recent and future changes in climate. Three main mechanisms have been suggested for the observed expansion of liana species in recent years: higher resource availability (which includes both increased atmospheric carbon concentration and increased nutrient deposition), changes in temperature and precipitation, and increased levels of disturbance (see Figure 1; further discussed in Schnitzer, 2005Schnitzer, , 2015. Most of these effects, however, are still poorly explored, and only the positive correlation of liana density with disturbance is well documented in both tropical (e.g. Ledo & Schnitzer, 2014) and temperate environments (e.g. Allen et al., 2007). Some studies have explored the effects of carbon fertilization on liana F I G U R E 1 Main hypotheses for liana expansion under changing climate. "+": suggested positive effects, "−": suggested negative effects. Dashed lines are relationships untested in the literature. The present work focuses on the relationships: Climate-LIANA, ↑CO 2 -LIANA, and the interactive effect Climate-↑CO 2 growth (Hättenschwiler & Körner, 2003;Marvin, Winter, Burnham, & Schnitzer, 2015;Schnitzer, 2015;Zotz, Cueni, & Körner, 2006), but these results remain highly controversial as carbon fertilization experiments have shown results ranging from no response (Marvin et al., 2015) to over 60%-100% biomass stimulation (Hättenschwiler & Körner, 2003;Zotz et al., 2006). The effect of changing temperature is also uncertain, although it has been presumed that an increase in temperatures will likely facilitate lianas to survive winter conditions (Schnitzer, 2005). There are no comparative assessments of the importance of each mechanism and their interactions, which would be necessary to accurately forecast how these species will respond to future climatic change.
The response of H. helix to changing environment will likely depend on the action of several climate change drivers, which can also be affected by the interactions between them. In the case of H. helix, we expect both temperature and precipitation to be critical factors to define the species distribution. We also expect H. helix to profit from increased carbon availability. However, given the high drought stress that Mediterranean plants experience in summer, we also expect summer water availability to influence whether the plant can profit from the increased carbon availability. Specifically, we aim here to explore the following questions: (a) How might H. helix respond to future climate change? To address this point, we used presence records to model the species climatic niche using SDMs, to obtain predictions of suitable habitat under current and future climatic conditions. (b) Does H. helix benefit from increasing atmospheric CO 2 concentrations ([CO 2 ])? To explore this, we used data from a study site in Italy to look at the long-term interactive effects of [CO 2 ] and temperature on the performance of H. helix, and its host tree, Quercus cerris using tree rings and assessing cavitation risk. (c) Could carbon fertilization effects substantially modify the climateonly model habitat predictions? As higher carbon availability cannot currently be accounted for in the SDM, we explored the potential effects by creating two future carbon fertilization scenarios for H. helix and comparing them with the previously developed, climateonly model.

| Species distribution modeling
We constructed SDMs from publicly available occurrence data for one of the main liana species in European forests, H. helix. We used occurrence data from the global biodiversity information facility (GBIF, http://www.gbif.org), where we downloaded 18,265 presence records of the species. In our search, we used the species name and limited the results to "observation" data after 1900, to avoid including herbarium specimens, specimens preserved in botanical gardens, or machine observations. Repeated presence points (i.e. those with exactly same coordinates) were eliminated to avoid overfitting problems. Climate information was obtained from the Bioclim dataset (http://www.worldclim.org/bioclim), which includes 19 biologically meaningful variables derived from monthly temperature and rainfall, and is commonly used in SDM studies (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005 (Martin et al., 2011).
We used Maxent software (Phillips & Dudík, 2008) to calculate the SDMs. Maxent is a presence-only modeling software that calculates habitat distributions with a maximum entropy algorithm.
We selected the three most relevant climatic variables, based on a preliminary Jack-knife and variable importance analyses, ecological knowledge of the species, and after excluding highly correlated variables (r > 0.7 or VIF >3), to reduce overfitting or collinearity problems. For each model, we randomly split the dataset into 75% (training) and 25% (test) subdatasets (Fielding & Bell, 1997

| Dendrochronological methods
In our second analysis, we explored the long-term growth trend of H. helix individuals and their response to increasing carbon availability and changing climate in an Italian case study. We sampled tree rings in H. helix and Q. cerris in temperate submontane forests in Abruzzo region, central Italy (Supporting Information Figure S1), close to the Riserva Naturale Guidata Abetina di Rosello. For a detailed vegetation and fauna description, see Pirone et al. (2005). We targeted H. helix rather than other liana species because it is known to form distinctive and regular tree rings that have been cross-dated successfully in the literature, and it is long-lived (in our study area, the oldest individual had 79 dated years). It is also important that in our study area, H. helix allocates all its growth to one principal stem rather than dividing it between several small ones, thus providing a better proxy for the total biomass of the individual.
In August 2014, we sampled two increment cores per individual from 30 randomly spaced trees and 32 lianas dispersed over the whole study area (Supporting Information Figure S1a).
We selected dominant trees to maximize the climatic signal (Schweingruber, 1966) and for similar reasons, avoided sampling the trees on which our sampled lianas were growing to prevent confounding effects of lianas on their hosts. We targeted Q. cerris in patches where it dominates the canopy in order to avoid effects of interspecific competition. Hedera helix individuals growing as a continuous grass-like layer were not sampled to avoid pseudoreplication problems. Sampling, cross-dating, and measurement followed standard dendrochronological methods (Bräker, 2002;Schweingruber, 1966;Speer, 2010). Despite challenging ring recognition in some H. helix (Supporting Information Figure   S1b,c), tree-to-tree agreement in growth patterns allowed accurate cross-dating of both species (Supporting Information Figure   S2a, Table S2). A small percentage of cores with growth anomalies or undetectable rings were discarded, but at least one core was maintained for all but one individual. The tree-ring data are available online in the International Tree Ring Data Bank depository (accessible at: https://www.ncdc.noaa.gov/data-access/paleoclimatology-data/datasets/tree-ring with access number ITAL045 and ITAL046 for Q. cerris and H. helix, respectively).
We used long-term climate records from the Climate Research Unit (CRU), University of East Anglia (available at: http://climexp. knmi.nl/) and atmospheric [CO 2 ] records from the longest recording measurement, Mauna Loa station (available at: http://www. esrl.noaa.gov/). Each individual tree ring series was transformed to basal area increment (BAI) before analyzing. BAI provides a better quantification of biomass production in ever-increasing individuals (Schweingruber, 1966): where R i represents the radius of the tree in a given year, i. We used linear mixed models to model the response of BAI to annual changes in precipitation, temperature, age, and [CO 2 ]. Tree identity and series temporal autocorrelation were included as random terms in the model. We included all two-way interactions between the variables, and BAI was square rooted to meet normality and heteroscedasticity assumptions. We did model selection by backward deletion, based on likelihood ratio tests and comparisons between models using ANOVA and Akaike's information criterion (AIC), to find the most plausible, minimal, model (Crawley, 2013). It is important to notice that in the case of lianas, BAI may underestimate the annual increase in biomass, as lianas are characterized by a stronger increase in length, rather than in diameter, meaning that our results may underestimate liana response to increased carbon availability.
To visualize the interaction between [CO 2 ] and maximum temperature on secondary growth for both species, we calculated the change in slope in the growth-[CO 2 ] correlation with regularly increasing maximum temperature (Figure 3). We also showed the overall influence of this interaction on the growth time series for each species by plotting predicted growth over time underrecorded [CO 2 ] and constant [CO 2 ]. When assumed to be constant, [CO 2 ] was assigned the average value of the complete available record. Tree age was also considered constant in the model predictions ( Figure 4).
To reinforce our findings and better understand the observed growth trends, we complemented the tree ring modeling with an assessment of the xylem cavitation risk (Supporting Information Figure S3). In vascular plants, cavitation refers to the breakage of the water transport column by the entrance of air bubbles into the xylem due to strong water tension inside the conduits (Tyree & Zimmermann, 2002). We estimated the vulnerability to cavitation for each tree ring using anatomical measurements of stems from a subsample of individuals of both species (Lens et al., 2011). We collected 20μm slides of seven H. helix and 11 Q. cerris, stained them with 1% safranin-astra blue (Arbellay, Fonti, & Stoffel, 2012), and analyzed them using WinCell Pro V 2004a to obtain a time series of vulnerability per tree. We chose the oldest H. helix individuals and the youngest Q. cerris individuals to minimize age differences between species. The xylem vulnerability index (VI) was calculated is calculated following White's equation, and V f is the total number of earlywood vessels counted per unit (Arbellay et al., 2012;Carlquist, 2001;Lens et al., 2011).

| Effects of carbon fertilization on H. helix model predictions
To illustrate how a positive effect of increasing [CO 2 ] could affect model predictions of future habitat suitability, we defined three potential CO 2 -effect scenarios for future H. helix habitat distributions ( Figure 4). These were as follows: "no fertilization," "slight carbon fertilization," (simply assuming a 10% increase in habitat suitability across Europe following carbon fertilization), and "strong carbon fertilization" (a 20% increase in habitat suitability) scenarios. These numbers are necessarily arbitrary and used for explorative purposes only, as there is no available method to translate growth responses into changes in habitat suitability. They are also different to the previously described modeling approach that reconstructs growth with and without CO 2 effects. A 10%-20% increase in habitat suitability seems a realistic assumption given the strong response to [CO 2 ] observed in our long time series, as well as the strong, unsaturated response of H. helix to carbon fertilization in previous experimental studies. Hättenschwiler and Körner (2003) recorded up to 100% increase in biomass and 137% in stem length under 660 ppm of CO 2 , without signs of carbon saturation, and Zotz et al. (2006) similarly observed a 30%-60% increase in length and biomass increment in shoots of H. helix under elevated CO2. All these analyses were carried out in the R software (R Core Team 2015).

| Climate-based spatial distribution models
The SDMs showed that changes in temperature and precipitation alone will likely alter the amount of habitat suitable for H. helix in Europe, requiring large range shifts for the species to track its climatic niche (Figure 2). The H. helix climatic niche was mainly defined by precipitation in the coldest quarter of the year (51% permutated variable importance), annual mean temperature (42.6% variable importance), and precipitation in the wettest quarter (6.4% variable importance). The positive correlation with mean temperature is in line with the biogeographical history of H. helix, which evolved from tropical families under warmer climates (Metcalfe, 2005). Our models project a slight increase in H. helix habitat by 2050 and a slight decrease by 2070 (Figure 2), but there is a large shift in the species range toward the north and east in the future (Figure 2). A future increase in the mean temperature in eastern and northern parts of Europe might make these regions suitable for frost-sensitive liana species in the near future (Gianoli, 2015;Schnitzer, 2005). Also, a positive correlation with winter precipitation supports the hypothesis that evergreen lianas may profit more from warmer conditions in winter and early spring than their broadleaved counterparts, as they can benefit from high light conditions before tree bud break occurs (Schnitzer, 2005 (Fitzpatrick & Keller, 2015). A particularly relevant factor to define the realized distribution in the case of lianas is host preference. However, it is unlikely that this factor plays an important role in the distribution of H. helix, as it seems to have low host specificity and preferences are driven more by size and bark characteristics (Castagneri et al., 2013) or forest successional stage (Ladwig & Meiners, 2010) than responses to particular tree species.

| Long-term climate interactions using tree rings
While our models predict that changes in temperature and precipi-  (Körner, Morgan, & Richard Norby, 2007;Zotz et al., 2006).  Figure   S3a,b). Q. cerris growth, on the other hand, has not increased over time. The trends in vulnerability to cavitation support these findings: Increasing cavitation risk with increasing temperature suggests that the inability of Q. cerris to respond to increasing [CO 2 ] concentrations is likely related to its greater sensitivity to drought (Supporting Information Figure S3c,d). The increased growth of H. helix, on the other hand, is not accompanied by an increase in vulnerability to cavitation, further suggesting that H. helix seems to take greater advantage of increasing [CO 2 ] than its tree host, as tree growth seemed to have been limited by the simultaneous rise in temperature and water stress.
In recent years, there has been a large controversy about whether tree ring measurements can be used for evaluating historical growth F I G U R E 2 Climate-only models for the suitable habitat of Hedera helix predict slight changes in the total habitat area but large range shifts toward north and east in future projections. (a-c) Species distribution models show the change in area suitable for H. helix from the present to 2050 and 2070 patterns (Bowman, et al. 2012;Brienen, Gloor, & Ziv, 2016;Brienen, Gloor, & Zuidema, 2012) because of the so-called big-tree selection bias and slow-grower survivorship bias that may confound climate growth trends with population dynamics or sampling strategies . A central aspect of both biases is the unbalanced relationship between age and annual growth in the sampled population, resulting in a spurious positive trend in growth produced by the different composition of slow-and fast-growing trees in each sampled age class (Brienen et al., 2016; more info in Supporting Information). We found no significant correlation between age and mean annual growth in our samples (Supporting Information Figure S4). Therefore, we do not expect these biases to influence our results. We also do not expect substantial disturbance or logging to have occurred in our study area, which has been protected since 1997. Although we cannot completely rule out the possibility that our forests are recovering from some old disturbance, the well-distributed samples across the landscape and their consistency make this unlikely.

| Potential consequences for the spatial distribution models
We illustrate how the positive effects of [CO 2 ] on H. helix might affect its future distribution, using three carbon fertilization scenarios.
Even small carbon fertilization effects could drastically change the predictions of how H. helix will cope with drier conditions in central and southern Europe, as well the extent to which it will expand northwards ( Figure 5). It is noteworthy that, due to the nature of the interaction between the temperature and carbon fertilization effects, it is likely that our model underestimates suitability in the colder areas of Europe while overestimating it for the hotter southern European localities, as [CO 2 ] effects were found to be lower at high temperature (Supporting Information Figure S5). can become limiting to plant growth, and therefore, no long-term fertilization effect would be expected (Bader et al., 2013). We found no saturating effect in the carbon fertilization curve; however, as our approach is not manipulative, carbon fertilization effects are limited to those increases which have occurred in the last decades.
The growth response that we observed is therefore in response to much lower [CO 2 ] than those used in fertilization experiments but it is in line with reports from other fertilization experiments with temperate lianas that did not find a reduction in carbon fertilization effect after a 80% increase in [CO 2 ]. This is a much larger increase than what has been observed in the last decades (Hättenschwiler & Körner, 2003

| CON CLUS IONS
We found that climate-only models predicted slight changes in the total area suitable for Hedera helix in the future, accompanied by large shifts in the species range toward the north and east of Europe. These models, however, could not account for the effects related to an in-

ACK N OWLED G M ENTS
We thank Michael Scherer-Lorenzen and Filippo Bussotti for assisting in the tree core sampling procedure and the Ispettorato Generale Zavala, Michael Scherer-Lorenzen, and three anonymous reviewers for helpful comments on an earlier version of this manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
RDM conceived and led the project. RDM, JBC, and FS designed the experimental sampling and performed the dendrochronological and ecophysiological analyses. RDM and EA performed the statistical analyses. All authors contributed to data interpretation and writing in consecutive versions of the manuscript.

DATA ACCE SS I B I LIT Y
Presence data: publicly available data from Global Biodiversity