The role of environmental conditions in regulating long-term dynamics of an invasive seaweed

The mechanisms underpinning long-term dynamics and viability of invader populations in the receiving environment remain largely unknown. We tested the hypothesis that temporal variations in the abundance of a well-established invasive seaweed, Caulerpa cylindracea, in the NW Mediterranean, could be regulated by inter-annual fluctuations in environmental conditions. Abundance data of C. cylindracea, sampled repeatedly between 2005 and 2020 at the peak of its growing season (late summer/early fall), were related to interannual variations in seasonal seawater temperature, wind speed and rainfall recorded during different growth phases of the alga, in both subtidal and intertidal habitats. In both habitats, higher peak of C. cylindracea cover was associated with lower seawater temperature in spring and summer, when the seaweed exits the winter resting phase and starts a period of active growth. In addition, the peak abundance of subtidal C. cylindracea was positively associated with higher autumn wind speed intensity and spring daily total precipitation. Our study reveals the importance of seasonal and interannual variation of abiotic factors in shaping temporal patterns of abundance of C. cylindracea, in both subtidal and intertidal habitats. Identifying the factors underpinning invasive population temporal dynamics and viability is essential to predict the time and conditions under which an invader can thrive, and thus guide management strategies aimed to containing invasions under current and future climates.


Introduction
Biological invasions are among the major threats to biodiversity and ecosystem functioning (Strayer 2012;Simberloff et al. 2013). The potential of non-native species to establish and spread outside their native range is regulated by several factors, including biotic and abiotic characteristics of recipient systems, disturbance regimes, invader traits and propagule attributes (Stachowicz et al. 2002;Levine et al. 2003;Lockwood et al. 2013;Bulleri et al. 2020;Monnet et al. 2020). A large research effort has investigated the drivers of invasion success, but less attention has been devoted to identify the mechanisms underpinning temporal dynamics of invasive species once they have successfully gone through the establishment phase (Strayer et al. 2006;Strayer 2012).
Abiotic factors play a pivotal role in regulating community structure and species range limits. Longterm information on population dynamics can help to identify climatic and environmental factors determining the abundance and distribution of species, including invasive ones (Bell et al. 2015;Miloslavich et al. 2018;Brown et al. 2020b;Burrows et al. 2020;Rilov et al. 2020). This may be particular useful under future climates, where seawater warming and acidification, as well as the increase in the intensity and frequency of extreme events (e.g., heat waves, storms) are predicted to cause shifts in species distribution (Molinos et al. 2016;Chefaoui et al. 2019;Beas-Luna et al. 2020).
The establishment and spread of non-native species generally occur in areas that are climatically similar to those of their native range (Bomford et al. 2009), however fluctuations of abiotic conditions (e.g., temperature, precipitation, resource availability) among years can alter invader performance (e.g., reproduction, growth, foraging activities), thus influencing their population dynamics (Hobbs and Mooney 1991;Krushelnycky et al. 2005;Murphy et al. 2017;Forsström et al. 2018;Menke et al. 2018;Chefaoui et al. 2019;Brown et al. 2020a). For instance, the abundance and distribution of the invasive Argentine ant, a common urban pest of Mediterranean-type ecosystems, is primarily driven by interannual variations in winter precipitation, which may enhance their survival, reproduction and spread due to higher soil moisture (DiGirolamo and Fox 2006;Heller et al. 2008;Menke et al. 2018). Likewise, a multi-year monitoring data of the introduced crab, Rhithropanopeus harrisii, in the Baltic sea, indicated that the recruitment and growth of juveniles increased during warmer summers (Forsström et al. 2018). Indeed, in Mediterranean-type ecosystems, most of the nonnative plants causing drastic alterations to native communities are annual and exhibit marked interannual fluctuations (Mooney et al. 1986;Valliere et al. 2019). Within this context, identifying the environmental drivers regulating temporal patterns of abundance and distribution of non-native species is crucial for pinpointing areas most vulnerable to invasion and to predict how invader populations will respond to future climates.
The macroalga, Caulerpa cylindracea, is one of the most widespread invaders in the Mediterranean Sea, being able to colonize both rocky and sandy bottoms, from the intertidal habitat to about 40 m depth (Piazzi et al. 2016). C. cylindracea is characterized by strong seasonal fluctuations in abundance, with a fast-growing phase in summer and a winter resting phase, primarily driven by temperature (Ruitton et al. 2005). In particular, in shallow waters, it starts developing erect fronds in early spring, peaks in abundance in September and dies back from November to March, with the loss of fronds and ramuli. This seaweed is often more abundant in wave-sheltered sites (Vaselli et al. 2008;Iveša et al. 2015) and preferentially colonizes complex biogenic substrata, such as those formed by algal turf and dead Posidonia oceanica rhizomes (hereafter also referred to as dead matte), that may favour the trapping and anchoring of its fragments under wave stress (Ruitton et al. 2005;Bulleri et al. 2011). In addition, positive effects of nutrient supply, during both the active growth and the resting phase, have been documented on C. cylindracea cover at the peak of the growing season (Gennaro et al. 2015;Uyà et al. 2017;Bulleri et al. 2018). To the best of our knowledge, however, no study has investigated how environmental variables, which vary greatly among seasons and years, influence the long-term invasion dynamics of well-established C. cylindracea populations.
Here, we used data on C. cylindracea abundance, sampled repeatedly between 2005 and 2020, to assess whether interannual variations in environmental conditions affect the alga at the peak of the growing season, in both intertidal and subtidal habitats in the NW Mediterranean. The set of environmental conditions included seawater temperature, wave intensity and rainfall regime, which are likely to play a role in determining patterns of abundance and distribution of the alga (Piazzi et al. 2016).
Caulerpa cylindracea displays a broad distribution in its native range, spanning from temperate water of the western Australia to warmer water in the northern Australia (Verlaque et al. 2003;Sauvage et al. 2013). In shallow Mediterranean water, this species alternates a fast-growing phase in summer to a resting phase in winter, characterized by a major biomass regression. Thus, we expected that warmer seawater outside the active growth period could promote an earlier recovery during the spring recovery phase, subsequently resulting in a greater summer peak abundance. An increase in summer seawater temperature could instead negatively affect C. cylindracea cover at the peak of the growing season, as the performance and survival of this species have been shown to decrease at high seawater temperature (Verbruggen et al. 2013;Santamaría et al. 2021). In addition, an increase in wave intensity during the decay phase (autumn) could enhance the spread of C. cylindracea through the production and dispersal of fragments (Piazzi et al. 2016), sustaining the next generation. On the other hand, high hydrodynamic stress could negatively affect C. cylindracea cover at the peak of the growing season, as this species generally achieves higher abundance at sheltered sites (Ruitton et al. 2005;Vaselli et al. 2008;Iveša et al. 2015). Thus, although seawater temperature and wave action are potential determinants of the peak abundance of this seaweed, we cannot devise a priori a directionality of their effects. Finally, previous studies have shown that enhanced nutrient loading sustains the growth and abundance of C. cylindracea, irrespective of the time of the year at which it happens (Gennaro et al. 2015;Uyà et al. 2017). Thus, we expected the summer peak in abundance of the seaweed to be greater in years characterized by heavy rainfalls, enriching seawater through runoff, either before (fall-spring) or during the active growing phase.

Study sites and data collection
Data on C. cylindracea abundance were collected within the framework of experimental and descriptive studies, carried out between 2005 and 2020 (data sources are found in Appendix 1). The percentage cover of this species was estimated in shallow subtidal sites along an exposed rocky shore between the city of Livorno and Calafuria, about * 5 km south (NW Mediterranean; 43°30 0 N, 10°20 0 E). Those sites are subjected to inputs of sediment, due to run-off from the overtopping sandstone cliffs, and organic and inorganic pollutant discharges from the nearby urban city of Livorno (Tamburello et al. 2012). Macroalgal assemblages at these sites are predominantly colonized by turf-forming algae Bulleri et al. 2010). Data on C. cylindracea cover were collected every 1-2 years between 2005 and 2013 and again in 2020 (8 sampling years in total), during late summer-early fall, corresponding to the period of the year in which C. cylindracea peaks in abundance. The percentage cover of this species was estimated not destructively, at depths varying from 6 to 8 m, by means of a plastic frame subdivided into 25 sub-quadrats. A score from 0 to 4 was given to C. cylindracea in each sub-quadrat and the percentage cover was obtained by summing over the entire set of sub-quadrats. In 2008, benthic assemblages (including C. cylindracea) were photographed underwater, and images were processed in the laboratory to extract from each quadrat the percentage cover of each species, using the same method described above.
Along the coast of Calafuria, C. cylindracea cover was also assessed in rockpools, between 2005 and 2020. We used data collected every 1-2 years between 2005 and 2010 and again in 2016 and 2020 (6 sampling years in total). Rockpool assemblages, characterized by turf-forming algae, barren patches and the presence of canopy-forming algae, such as Cystoseira spp. and Halopteris scoparia, have been invaded by C. cylindracea since the early 2000s (Incera et al. 2010;Maggi et al. 2012). Caulerpa cylindracea cover was visually estimated in the field, during late summer-early fall, by means of a plastic frame subdivided into 25 sub-quadrats, as described above. Details of size and number of quadrats used in each sampling event, for both subtidal and intertidal habitats, are reported in the Table A1 of Appendix 1.

Environmental variables
Data on hourly seawater temperature (°C) and wind speed (m/s) were collected from the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA, www.mareografico.it; station located in Livorno). Onshore winds (between 180 and 360°for the study sites) were used as representative of wave height, which was assumed to be positively correlated to wind speed intensity (Ellis et al. 2019). Data on daily rain (mm) were collected from the Servizio Idrologico Regionale (SIR, www.sir.toscana.it; station located in Livorno). Data were, then, subdivided into four periods throughout the year, each corresponding to a different phase of the life cycle of C. cylindracea (declining phase = October -November; resting phase = December -March; early-growing phase = April -June; fast-growing phase = July -August).
We calculated mean and maximum daily values for seawater temperature and wind speed, while, for rainfall, we calculated maximum monthly values from total daily precipitation, as available from the SIR website. Within each year, an average value of each metric (mean and maximum daily seawater temperature and wind speed and mean daily and maximum monthly precipitations) was calculated for each period of C. cylindracea growth.

Statistical analyses
The full set of environmental variables (mean and maximum seawater temperature, wind speed and rainfall), over the cycle of C. cylindracea growth, was compared in a matrix of pairwise correlations, for both subtidal habitat and rockpools. Since the mean and maximum metrics of each environmental variable, within each period of C. cylindracea growth, were highly correlated ( Fig. S1 and S2 in Supporting information), we performed statistical analyses separately for the mean and maximum metrics, for both the habitats.
For each metric, we then tested the collinearity of environmental variables by means of the variance inflation analysis (VIF) and discarded those variables with a VIF value[5 (Zuur et al. 2009). The final sets of uncorrelated environmental variables used in our analyses are reported in Table 1.
For subtidal habitats, generalized linear mixed models (GLMM), assuming negative binomial distribution, were used to model the relationships between temporal patterns of abundance (estimated as percentage cover) of C. cylindracea and explanatory variables, accounted for the effects of selected environmental variables, separately for the mean and maximum metrics. For both metrics, the environmental variables were included in the fixed part of the models. To deal with temporal autocorrelation (a subset of plots was repeatedly sampled through years), the plot was also included as a random effect in the model (Table S1 in Supporting information).
In the analysis of data from rockpools, zero-inflated GLMMs, assuming a gaussian distribution, were used to model the relationships between temporal patterns of C. cylindracea cover and the selected environmental variables, separately for the mean and the maximum metrics. The zero-inflated model allows to handle the high number of quadrats with 0% cover of C. cylindracea. For both the metrics, the selected environmental variables were included in the fixed part of the models. Since a subset of pools was repeatedly sampled through years, the pool was included in the random part of the model (Table S1 in Supporting information). In order to take into account variations in the size of quadrats used to sample C. cylindracea cover among years (Table A1 in Appendix 1), the size of each quadrat was used as the offset in the analysis. The GLMMs for both habitats were run using the glmmTMB R package (Brooks et al. 2017).
The minimum adequate models (MAM) were generated by means of a step-backward selection procedure (buildglmmTMB function of the buildmer package), using the Akaike information criteria. Residuals were checked with QQ-plots and plots of standardized residuals versus expected values, using the R package DHARMA (Hartig 2020). The same package was used to run a Kolmorov-Smirnoff test to assess heteroscedasticity and goodness-of-fit tests on simulated residuals to check for over-dispersion and zero inflation (Fig. S3 and S4 for subtidal and rockpool analyses, respectively, in Supplementary information). All statistical analyses were performed in R version 3.6.1 (R Development Core Team 2019).

Results
Temporal patterns of C. cylindracea cover in the subtidal habitat Using mean values of environmental variables, the MAM retained two explanatory variables, accounting for the effects of rainfall and wind speed, from the early-growing (April-June) and declining (October-November) phases, respectively (Table 2A). The peak abundance of C. cylindracea was positively related with increasing mean daily wind speed intensity in the period of its declining-growth and with rainfall in the period of its early-growth (Table 2A, Fig. 1a, b).
The MAM fitted with maximum values of environmental variables, retained seawater temperature from the fast-growing phase (July-August), as explanatory variable (Table 2B). There was a significant negative relation between peak C. cylindracea cover and increasing maximum daily seawater temperature in the fast-growth period (Table 2B, Fig. 1c).
Temporal patterns of C. cylindracea cover in rockpools The MAM fitted with mean values of environmental variables retained seawater temperature from the fastgrowing period (July-August), as explanatory variable (Table 3A). For maximum values, seawater temperature from the early-growing period (April-June) was included as explanatory variable (Table 3B).
The peak abundance of C. cylindracea was negatively related with increasing mean and maximum daily seawater temperature in the fast-and earlygrowth periods, respectively (Table 3A,B; Fig. 2a,b).

Discussion
Long-term monitoring data can allow the identification of factors regulating invasion dynamics and help to predict future patterns of spread (Menke et al. 2018;Chefaoui et al. 2019;Brown et al. 2020a;Keller and Shea 2021). In this study, we evaluated whether yearly peaks in the abundance of C. cylindracea were correlated with interannual fluctuations in seawater temperature, wind speed and rainfall. Our results showed that, in both habitats, higher C. cylindracea cover at summer peak was associated with lower seawater temperature during the early-and fastgrowing periods, possibly as a consequence of acclimation to a reduced temperature range of the receiving environments. In addition, in accordance with our predictions, the peak abundance of subtidal C. cylindracea was positively related to increasing fall wind speed intensity and spring daily total precipitation.
Fluctuations of abiotic conditions, such as temperature, waves, currents and food availability, influence the abundance and distribution of marine species, including invasive ones (Iles et al. 2012;Bell et al. 2015;Bertocci et al. 2015;Chefaoui et al. 2019;Burrows et al. 2020;Rilov et al. 2020;Berry et al. 2021). As the oceans are rapidly changing due to anthropogenic activities, monitoring species abundance along with abiotic variables is essential to predict shifts in species distributions and their ecological effects on ecosystem structure and functioning (Bates et al. 2018;Miloslavich et al. 2018). Although the spread of non-native species represents one of the major ecological and socio-economic threats, the lack of time-series can make it difficult to identify the mechanisms underpinning invasion dynamics and to predict how invader populations will respond to future climates.
Temperature is one of the most important factors shaping species distribution and range limits (Parmesan and Yohe 2003; Cleland et al. 2007;Eggert 2012). In marine systems, seawater temperature can affect physiological performance (e.g., timing and length of reproduction, growth rate, body size) and survival limits of a variety of non-native organisms, including seaweeds (Murphy et al. 2017;Forsström et al. 2018;Chefaoui et al. 2019;Brown et al. 2020a). For instance, the abundance and growth of juveniles of invasive crabs have been shown to be positively correlated to warmer seawater in summer during their planktonic larval stage (Forsström et al. 2018). Likewise, patterns of abundance and distribution of the invasive macroalga, Undaria pinnatifida, across its introduced range in Europe, were related to different responses of different life cycle stage (sporophytes and gametophytes) to seasonal temperature variations (Murphy et al. 2017). Thus, the effects of temperature on different life stages of non-native species may regulate their performance once a new region has been colonized and may, subsequently, impact long-term population dynamics and viability (Murphy et al. 2017;Chefaoui et al. 2019;Keller and Shea 2021).
Caulerpa cylindracea presents strong seasonal fluctuations in abundance in the coldest regions of NW Mediterranean, with a fast-growing period in summer and a drastic biomass regression in winter, mostly driven by temperature (Ruitton et al. 2005). A previous study on spreading patterns of C. cylindracea in the northern Adriatic Sea documented an increase of colonized sites in years with warmer winter seawater temperature, suggesting that milder conditions during the winter resting phase could promote the regrowth of C. cylindracea and, thus, its spreading at the peak growing season (Iveša et al. 2015). In contrast, our analysis revealed a negative influence of increasing maximum daily seawater temperature in summer on the peak abundance of subtidal C. cylindracea. Similarly, the abundance of this seaweed in rockpools was negatively related to higher maximum and mean seawater temperature during the periods of emergence from the winter resting stage and active growth. The northern Adriatic Sea represents the northernmost region of the Mediterranean Sea, and it is characterized by lower winter seawater temperature compared to our study sites in the Ligurian Sea. Thus, the contrast between our results and those of Iveša et al. (2015) would suggest an adaptation of introduced populations to seawater temperature at regional scale. A negative influence of warmer seawater on this alga has been reported in previous studies, showing a peak in the occurrence of C. cylindracea at about 20°C and a sharp decrease at higher seawater temperatures in the invaded regions (Verbruggen et al. 2013). This seaweed displays a broad distribution in its native range, from the temperate waters of the western Australian coasts to the tropical waters of northern Australia (Verlaque et al. 2003;Sauvage et al. 2013).
Thus, the negative influence of warmer seawater on peak C. cylindracea cover could indicate a reduced temperature tolerance range for introduced populations (Verbruggen et al. 2013). Ocean warming due to climate change is predicted to affect seaweed abundance and distribution (Molinos et al. 2016;Chefaoui et al. 2019). Increases of species with warm-water affinity and declines of those with cold-water affinity, as well as the poleward expansion of species at cold range boundaries and retreats at warm edges have been documented along both the Pacific coast of America and the Atlantic coast of Europe (Beas-Luna et al. 2020;Burrows et al. 2020). Hence, our results suggest that the predicted seawater warming under future climate scenarios may change the distribution dynamics of C. cylindracea, with a potential retreat from currently occupied areas and a northward shift, as also predicted for other invasive seaweeds in temperate regions (Chefaoui et al. 2019). Experimental studies are, however, warranted to assess the thermal physiological thresholds for this species in the introduced range and to improve predictions of its future spread. Seasonal and interannual fluctuations in resource availability and hydrodynamic conditions can also play a role in driving seaweed abundance and distribution dynamics (Cavanaugh et al. 2011;Bell et al.  2015; Berry et al. 2021). For instance, along exposed coastlines, physical disturbance due to severe winter storms has been shown to be a major determinant of kelp cover. On the other hand, in relatively sheltered environments, interannual variations in nutrient concentration can influence temporal patterns of kelp abundance, with higher macroalgal recruitment and biomass associated with periods of increased nutrient levels (Cavanaugh et al. 2011;Bell et al. 2015). In our study, variation in the peak abundance of subtidal C. cylindracea was explained by cumulative rainfall during spring-growing phase. This could be the result of increased nutrient inputs to the subtidal habitat, due to run-off from overtopping sandstone cliffs during heavy rain events . In accordance with previous manipulative  (Gennaro et al. 2015;Uyà et al. 2017), our study suggests that C. cylindracea can take advantage of inputs of nutrients throughout different stages of its seasonal growth, ultimately fostering its abundance at the peak of the growing season. Temporal patterns of subtidal C. cylindracea were also positively related to an increase in mean wind intensity during the decay season (fall). Water motion has been suggested to have contrasting effects on C. cylindreacea spread (Piazzi et al. 2016). This species has found to be more abundant in sheltered areas (Iveša et al. 2015) and to preferentially colonize complex habitats (e.g., matrix of algal turfs or P. oceanica dead matte), which allow the retention and anchoring of its fragments, suggesting a certain susceptibility of this species to wave stress (Ruitton et al. 2005;Bulleri et al. 2011;Iveša et al. 2015;Piazzi et al. 2016). On the other hand, drifting fragments of C. cylindracea detached by large waves can easily reestablish during summer months, suggesting that intense hydrodynamic forces may enhance the proliferation and spread of this species before the peak cover period (Ceccherelli and Piazzi 2001;Piazzi et al. 2016). Our results suggest that an increase in wave intensity in autumn, when thalli are fully grown, may enhance the generation and dispersion of high-quality propagules (e.g., fragments with stolon, fronds and rhizomes), sustaining the next season generation. Assessing the extent of this mechanism and its potential contribution to the wide spread of C. cylindracea remains critical to advance our understanding of the invasion dynamics of this species.
In rockpools, the peak abundance of C. cylindracea was influenced neither by rainfall nor by wind speed intensity. A recent experiment showed that nutrient enrichment during the resting phase could promote the spring regrowth of C. cylindracea, but there was no carry over effect on the summer peak in abundance (Uyà et al. 2017). Assemblages inside rockpools are unlikely to be nutrient limited, as a consequence of frequent run-off that might increase nutrient concentration in these confined environments (Uyà et al. 2017). In addition, patterns of invasion of rockpools are likely to be much more dependent on subtidal populations of surrounding areas, which may replenish rockpools with algal propagules, rather than variations in seasonal wind speed condition. Rockpools are characterized by large fluctuations in abiotic conditions (e.g., temperature, salinity, pH) over short temporal scales (Incera et al. 2010;Maggi et al. 2012), and patterns of abundance of species, either invasive or native, in these habitats are generally more dependent upon short-term variation in environmental conditions than past climatic events (Incera et al. 2010;Uyà et al. 2017).
Our study reveals the importance of abiotic factors in shaping the invasion dynamics of the well-established C. cylindracea in the NW Mediterranean. In particular, the subtidal population appears to be influenced positively by an increase in the spring cumulative precipitation and autumn wave intensity and negatively by high seawater temperature in the summer period. Similarly, warmer seawater in the spring-summer period negatively affected peak C. cylindracea cover in rockpools. This pattern can be the result of direct effects on the physiology of C. cylindracea and/or indirect effects on native assemblages. Several correlative and experimental studies have shown that indirect effects of abiotic factors (e.g., nutrient enrichment, sedimentation, mechanical disturbance), that promote the spread of opportunistic turf-forming species at the expense of canopy-forming macroalgae, are key determinants of the success and persistence of C. cylindracea (Bulleri et al. 2010(Bulleri et al. , 2011Incera et al. 2010;Piazzi et al. 2016;Uyà et al. 2017). At our study site, turf-forming algae were dominant in both subtidal and intertidal habitats throughout the investigated period and consistently across seasons (Bulleri et al. 2010(Bulleri et al. , 2011Incera et al. 2010;Tamburello et al. 2012;Uyà et al. 2017). This suggests that yearly fluctuations in the abundance of C. cylindracea are unlikely to reflect changes in the structure of native communities. Long-term ecological studies can provide key insights into the responses of marine communities and ecosystems to environmental change (Hughes et al. 2017;Beas-Luna et al. 2020). To date, most long-term observations of biological variables and environmental conditions are restricted to a limited number of key functional groups, such as corals, phytoplankton and kelp forests, and limited to some geographical regions (Miloslavich et al. 2018). Our study highlights the importance to include biological invasions into long-term monitoring programs to assess how invasive species can respond to environmental drivers and to identify areas most vulnerable to invasions, with the ultimate goal of implementing management and adaptive strategies.
Finally, the identification of the mechanisms underpinning invasion dynamics is key to predict future spread of invaders under climate change (Chefaoui et al. 2019;Keller and Shea 2021). For instance, a recent study has modelled the future range of the invasive seaweed, Sargassum muticum, along the coasts of North America and Europe, under climatic scenarios predicted by the end of the century. The authors showed that ocean warming may cause a substantial retreat of this species from the southernmost regions of its current range and a potential northward shift to sub-arctic areas, although constrained by suitable reproductive temperature window (Chefaoui et al. 2019). Likewise, in a terrestrial system, climate warming has been shown to foster the annual population growth of the invasive musk thistle, Carduus nutans, by increasing the size of reproducing individuals and the size and survival of recruits, with potential negative effects on the productivity of pastures and rangelands (Keller and Shea 2021). Our results indicate a possible reduction in the abundance and distribution of C. cylindracea in both subtidal and intertidal habitats under future climatic conditions, in response to seawater warming. Such negative effects might be partially offset by the increasing amount of rainfall and the intensity and frequency of storms which could provide C. cylindracea with greater nutrient availability and facilitate its dispersal.