In summer 2001, four willow ramets of each genet were tagged with basal plastic slips. The shoot growth of tagged ramets was monitored every year. The shoot lengths of each tagged willow ramet were measured in early August (at the end of the active growing season) to an accuracy of 0.5 cm. In cases where an individual ramets died, a new ramet of the same genet were tagged. Initially, 320 ramets were tagged in 2001, while the total number of individual ramets tagged was 626. In addition, the level of insect herbivory was estimated using the following defoliation classes: 0 no observable insect herbivory, 1 (0–1%), 2 (1–5 %), 3 (6–25 %), 4 (25–100%) or 5 (leaves totally defoliated) (den Herderet al . 2004). Herbivory shows remarkable among-year variation in 2001–2018, but no trend (Supporting Information).

Climate data were obtained from the nearest Finnish Meteorological Institute weather station (Kilpisjärvi Kyläkeskus) located 10 km from the experimental site (480 m a.s.l.). The weather station provides daily records on mean temperature. These climatic station data were used to calculate thermal sum for the season, when shoot growth takes place. Therefore, daily mean temperatures exceeding + 5 °C in May–July were summed up to produce a “growing degree days - GDD” index, hereafter denoted as Climate . Thermal sum shows remarkable among-year variation in 2001–2018, but no clear trend (Supporting Information). Note that although there are many other potential climate-related variables that we could have considered (e.g. length of growing season, growing degree days with different threshold, mean July temperature), these all correlated with the thermal sum as described above – we therefore chose to only consider this commonly used variable in our analyses.

Empirical dynamic modelling

We used EDM to model the relationships between willow shoot growth, climate and herbivory. EDM is a recently developed technique that combines two tools (Deyle et al . 2016) – convergent cross mapping (CCM) and sequential locally weighted global linear maps (S-mapping). CCM is a nonparametric method for tracking information flow among variables, and is a powerful tool for identifying causal links in systems with complex dynamics (Sugihara et al . 2012; Clarket al . 2015). Because the direction of information flow is directional in causal interactions (i.e. “forcing” processes impart information onto “response” processes, whereas the reverse is not true), CCM is also often successful at identifying the direction of causal interactions based on observational data. S-mapping is a related technique that can be used to extract information about the strength and direction of these relationships (Sugihara 1994). These methods are particularly useful for our purposes, as they can be applied withouta priori knowledge of the functional form of the relationship governing interactions and apply even for highly complex or chaotic dynamics. They also provide two major advantages over more classic methods, such as mixed-effects linear models or GLM’s. First, in complex ecological systems, results from CCM tend to be more closely aligned with actual causal forcing than is true for regression (Sugiharaet al . 2012; Clark et al . 2015). This feature allows us to generate putative interaction networks among species based on observational data. Second, S-mapping captures effects of dynamic, nonlinear relationships among variables. Thus, S-mapping allows us to analyse the strength and direction of interactions among species even when they are complex, time-varying, and context dependent. The analyses were employed using rEDM package and EDMhelper packages in R (Yeet al . 2018; Karakoç et al . 2020).

Based on a time-delayed application of Convergent Cross Mapping (CCM) (Fig. 1), we chose to quantify effects of climate (Climate ) and insect herbivory (Herbivory ) on mean shoot length of topmost part of the ramets (Shoot ). In CCM analysis, causal forcing is indicated when predictive skill increases significantly with the total number of historical observations used to make predictions (i.e. “library length”). We identified significant forcing of ramet dynamics by both climate and herbivory (Fig 1). Because multi-trophic interactions often involve time-delayed effects, we tested for both immediate and lagged effects. Analyses suggested that effects of insect herbivores on shoot growth were immediate (i.e. time-lag = 0), whereas climate was best represented by a time-lagged effect (time-lag = -1). All variables were standardised to zero mean and unit standard deviation before analysis, to better allow for comparison of effect sizes among variables. Note that although “best practices” dictate that 30–40 sequential observations are usually necessary for CCM and, to a lesser extent, S-mapping, this limit can be met either by single long time series, or using composites that have been stitched together from multiple shorter time series (Hsieh et al . 2008; Clark et al . 2015). For all analyses presented here, within each treatment, time-series from all ramets and replicates were stitched together, resulting in time series of > 640 measurements for each treatment.