Cross-basin and cross-taxa patterns of marine community tropicalization and deborealization in warming European seas

Domain: earth_science
Published: 2024-03-08
Editorial outcome: minor_revision

Abstract
Ocean warming and acidification, decreases in dissolved oxygen concentrations, and changes in primary production are causing an unprecedented global redistribution of marine life. The identification of underlying ecological processes underpinning marine species turnover, particularly the prevalence of increases of warm-water species or declines of cold-water species, has been recently debated in the context of ocean warming. Here, we track changes in the mean thermal affinity of marine communities across European seas by calculating the Community Temperature Index for 65 biodiversity time series collected over four decades and containing 1,817 species from different communities (zooplankton, coastal benthos, pelagic and demersal invertebrates and fish). We show that most communities and sites have clearly responded to ongoing ocean warming via abundance increases of warm-water species (tropicalization, 54%) and decreases of cold-water species (deborealization, 18%). Tropicalization dominated Atlantic sites compared to semi-enclosed basins such as the Mediterranean and Baltic Seas, probably due to physical barrier constraints to connectivity and species colonization. Semi-enclosed basins appeared to be particularly vulnerable to ocean warming, experiencing the fastest rates of warming and biodiversity loss through deborealization.

Introduction
Global climate change is fundamentally altering life on Earth1–3. This alteration is primarily due to three universal ecological responses of species to global warming including4 poleward distribution shifts, adjustments in phenology, and reductions in mean body size. Shifts in distribution or seasonal timing help species track their specific thermal niches5–7, whilst body size reduction is a complex process involving physiological and/or eco-evolutionary responses8–10. A poleward distribution shift is the displacement (active or passive) of populations that track species’ thermal preferences11, which is mainly driven by the general trend of decreasing water temperature with latitude12.

Although our knowledge of marine species responses to ocean warming in North Atlantic is substantial13,14, comparative assessments across biotic groups and regions to identify the underlying ecological processes15 associated with marine community expansion, retraction and dispersal constraints are limited to the uneven distribution of monitoring programmes16–20. The role of dispersal limitation and habitat fragmentation in community responsiveness to climate variation has been addressed more thoroughly in terrestrial groups such as plants21,22 compared to marine biodiversity20,23,24. However, organisms’ dispersal ability and ocean connectivity are also important in shaping marine species distribution25,26. For instance, the presence of the European land mass north of the Mediterranean basin poses a considerable limitation to the response of species to sea warming, as it constrains the poleward movement of adults or propagules to more thermally favourable conditions, especially jeopardising endemic species27.

Physiological constraints limit the range of suitable environmental conditions for populations and species28. To maintain viable populations under ocean warming, species might respond (actively or passively) by shifting latitudinally (generally poleward) or by changing their phenology (e.g. earlier onset of spawning, migration or return from dormancy)11. Latitudinal shifts will change the relative abundance of species at specific locations (Fig. 1). Long-term monitoring programmes based on permanent stations have detected such changes in North Atlantic fish communities associated with ocean warming15, based on temporal changes of the Community Temperature Index (CTI). The CTI is a measure of the average thermal affinity of ecological communities weighted by the relative abundance of each species29,30. Thus, temporal change in CTI informs on the turnover of the relative abundance of species (sensu Vellend31) according to their thermal affinity and the response of the community to temperature change.

Changes in the CTI can be decomposed into the four underlying ecological processes15. Positive CTI changes correspond to an increased prevalence in warm-affinity species (i.e., tropicalization), and/or a decrease in cold-affinity species (i.e., deborealization). Tropicalization is hence interlinked to leading (warm range) edges, whilst deborealization to trailing (cold range) edges (Fig. 1). Negative CTI changes correspond to an increase in cold-affinity species (i.e., borealization) and/or a decrease in warm-affinity species (i.e., detropicalization). The reason why tropicalization may prevail (or not) over deborealization has been recently debated32–34. Assuming an ecophysiological equilibrium between habitat suitability and the occurrence of species, we can expect that species track temperature change equally33,35, hence tropicalization would be similar to deborealization. However, in species with slow demography and limited dispersal, lags between climate change and distribution shifts can result in ‘extinction debts’36, where populations temporarily persist under unsuitable conditions, and ‘colonization credits’, where suitable locations are not occupied32. This latter pattern has been observed for instance in some studies on trees32 and demersal fish33, whilst other studies found equally responsive shifts at both range boundaries in marine ectotherms35, and prevalence of the tropicalization in demersal fish15. The prevalence in tropicalization vs deborealization might depend on biological traits of the community (size, dispersal capacity), phenotypic plasticity enabling populations to remain in sub-optimal habitats37, or due to seascape aspects (e.g., marine geomorphological features) limiting species dispersal and ocean connectivity. Thus, it can be expected that seas with limited connectivity to open oceans would constrain species abundance increases (tropicalization+borealization) over species abundance decreases (detropicalization+deborealization) processes.

The large geomorphological and ecological differences lead to large differences in the physical, biogeochemical and ecosystem responses of European Seas to climate change38. All European Seas are warming at least since the 1980’s, but faster in the semi-enclosed basins of the Mediterranean Sea39, Baltic Sea40, and Black Sea (1982–2018)41, compared to the NE Atlantic39,42. Here, we analysed the extent to which long-term trends in marine communities were related to ocean warming across biological groups and European regional seas. To do so, we used the CTI, which tracks the mean thermal affinity of a community29, to quantify temporal rates of turnover of the relative abundance of species. We applied CTI to long-term time series from 65 biodiversity monitoring programmes spanning the last four decades (in the longest case), which accounts for 1,817 species, including zooplankton, benthos, demersal and pelagic assemblages. Across the entire range of organisms and habitats, results show an average rate of increase in CTI of 0.23°C decade−1, meaning a consistent response of the marine communities to ocean warming in European seas. Furthermore, we explored the main underlying ecological processes driving temporal variation in CTI including: tropicalization, detropicalization, borealization, and deborealization15. We hypothesised that an increase in CTI with time linked to ocean warming may vary according to the biological groups, degree of ocean connectivity, and habitat type. Specifically, semi-enclosed sea basins with lower ocean connectivity to warmer waters are expected to experience less tropicalization compared to deborealization than is the case in the well-connected northeast Atlantic region. Our findings indicate that tropicalization dominated in Atlantic sites compared to semi-enclosed basins, supporting the expectations. Cross-region and cross-taxa comparisons may identify climate refugia as important reservoirs of biodiversity43, hotspots of high biodiversity climate velocity44 and sentinel systems45.

Results
We compared the rate of temporal change in CTI (CTIr) and its underlying ecological processes in biodiversity time series collected over the last four decades among six biological groups (hard-bottom and soft-bottom coastal benthic communities, zooplankton, demersal crustaceans, cephalopods, and fish). Biodiversity time series included three European Seas (NE Atlantic, Mediterranean Sea, and Baltic Sea), three habitats (marine benthic or demersal, marine pelagic, and estuarine), and two basin types based on the presence or absence of dispersal barriers to ocean connectivity (non-enclosed sea, i.e., Atlantic Ocean, and semi-enclosed seas).

In almost all biodiversity sampling sites, sea surface temperature increased from 1980 to 2021 (Fig. 2a). On average, the warming rate of sea surface temperature (SST) was 0.32°C decade−1 (t=28.24, n=1730, p<0.0001) (Fig. 2c) and the change of sea temperature from surface down to 100m (ST 100m) was 0.15°C decade−1 (t=9.70, n=1720, p<0.0001). CTIr trends mirrored those in ocean warming with an increase over time at most sites (80.0% of sites were positive, and 47.7% were significantly positive) (Fig. 2b; Supplementary Fig. 1, Supplementary Data 1). On average, the CTI also increased with time (Fig. 2c) at a rate of 0.23°C decade−1 (t=14.69, n=1730, p<0.0001, using a linear mixed model) (Fig. 2c). The relation between CTIr and sea temperature change across sites was positive and significant on average for both SST (t=7.47, n=1730, p<0.0001) and ST 100m (t=4.33, n=1720, p<0.0001), indicating that changes in marine communities can be associated with ocean warming.

We compared CTIr of marine communities across several factors (biological group, habitat, region, and basin type) (Supplementary Fig. 2). By biological group, CTIr trends were significantly positive for coastal soft-bottom benthos, fish, and zooplankton communities, with no significant change for demersal crustaceans and cephalopods (Fig. 3a) (Table 1). By habitat type, CTIr trends were significantly positive in all habitat types (Fig. 3b and Table 1). Estuaries showed the highest CTI value. CTIr was significantly positive in the NE Atlantic and Baltic, but not significant in the Mediterranean Sea (Fig. 3c and Table 1). The mean CTIr was also significant in both categories of basins (connectivity), with slightly higher values in non-enclosed seas compared to semi-enclosed seas (including estuarine types) (Fig. 3d and Table 1).

The selection of factors identified the region and habitat as the best model (AICc=1756.8). Diagnostic plots for the residuals of the selected factors were checked, indicating model reliability (Supplementary Fig. 3).

The underlying ecological processes explaining the increase in CTI can generally be attributed to the prevalence of tropicalization or deborealization in most of sites (76.9% of sites), whilst detropicalization and borealization dominated at fewer sites (23.1%) (Fig. 2d). On a per-species basis (encompassing 5324 cases of 1817 species), the tropicalization and deborealization account for 59.3% of the overall process intensity (Fig. 4). Among the four processes considered, tropicalization was the most frequent (53.8%) on a per-site basis, followed by deborealization (18.5%) (Fig. 2d). In those biodiversity time series where the CTI increased (i.e., 52 sites, 80% of sites), we analysed tropicalization relative to deborealization across biological groups, habitat, regions, and basin types (Table 2). By biological group, tropicalization relative to deborealization was most frequent for coastal benthos (soft- and hard-bottom), whilst the other groups showed more equal process intensities (Table 2). By habitat type, tropicalization only dominated in benthic/demersal habitats (Table 2).

By region, tropicalization dominated in NE Atlantic, whilst Baltic Sea and specially the Mediterranean presented more equal intensity (Table 2).

By basin type (i.e., ocean connectivity), tropicalization relative to deborealization strongly dominated the community responses in the non-closed seas (i.e., Atlantic Ocean), whilst semi-enclosed seas presented more equal process intensities (Table 2 and Fig. 5a). Comparing factors individually, basin type explained the greatest amount of variation of the process intensity according to the lowest AICc. Two models were selected as the best ones with the same AICc (i.e., 33.9): (1) habitat and region, and (2) habitat, region, and basin type. The results of the analysis of process intensity at per-species basis were qualitatively similar to those at per-site basis (Supplementary Table 1).

Species’ abundance increases (i.e., tropicalization+borealization) relative to species’ abundance decreases (i.e., detropicalization+deborealization) dominated in non-closed compared to semi-enclosed seas (p=0.00018; Fig. 5b).

Methods
Temperature data
To provide a common long-term dataset of sea temperatures through water column (surface and down to 100m), data from the NCEP Global Ocean Data Assimilation System (GODAS) (www.cpc.ncep.noaa.gov/products/GODAS/) were extracted for each site for the 1980−2020 period. GODAS provides monthly water temperature data to a 0.333° × 1° latitude-longitude grid. Annual data from the grid centred to each biodiversity sampling site was extracted, averaging a spatial window according to the area of the site. We analysed the slope of the temperature trends and its significance using a linear mixed model80 of sea temperature with year as fixed effect, the sampling sites as random effect, and adding a temporal autoregressive function. t- and p-values of the two-sided Wald test for the estimated slope coefficients are provided.

Biodiversity time series
We used 65 long-term time series corresponding to 1817 species taxonomically identified at species level (71 zooplankton species79,81, 238 coastal hard-bottom benthic species18,77,82,83, 923 coastal soft-bottom benthic species83–85, 63 cephalopod species86,87, 104 demersal crustacean species86,87, 418 fish species86–92). Fish includes pelagic, demersal, and estuarine fish species. Demersal crustaceans were collected offshore from bottom trawl fishing and surveys. Cephalopods includes demersal and pelagic species collected offshore from bottom trawl and purse-seine fishing and surveys. Coastal soft-bottom benthic macroinvertebrates includes species sampled in coastal intertidal or subtidal soft-bottom sites, using grab or box corer samples. Coastal hard-bottom benthos includes macroinvertebrates, macroalgae, lichens, and coralligenous assemblages in intertidal or subtidal, coastal, hard-bottom sites. Zooplankton includes mainly copepods and other species collected with zooplankton net hauls. Time series ranges over varying durations during the last four decades (from 1980 to 2022 in the longest case). Most of time series span more than 15 years (92%), with median of 25 years, minimum of 9 years, and 20 and 39 years the 25% and 75% quartiles, respectively. Detailed information on each time series such as data sources, sampling procedures, time series length, biological community, and location is provided in Supplementary Information 1.

Species relative abundance turnover: Community Temperature Index through time
To test whether temporal changes in the composition of species in a community respond to warming according to their thermal preferences, we used the CTI, which is a measure of the average thermal affinity of ecological communities, weighted by the relative abundance29,30:

CTI = \sum_{s=1}^{n} T_s \times \log(A_s + 1) (1)

where n is the number of species in the community, Ts is the temperature preference of each species (s) and As is the relative abundance of species s (i.e., the abundance of species s divided by the total number of individuals in the community at a site). The thermal preferences were determined for each species by matching occurrence records collected from OBIS (Ocean Biogeographic Information System; www.iobis.org) with annual means of GODAS local SST during the 1980−2021 period67. OBIS records were quality-checked removing duplicates. To characterise the thermal preferences of the species local temperatures derived from GODAS were used, which was available from OBIS for each observation. The midpoint between the 5th and 95th percentile of the temperature distribution occupied by each species was calculated as a measure of central tendency of their realised thermal distribution30. Subsequently, CTI was calculated at each station using the thermal midpoint values for each species recorded weighted by their log(abundance+1)30. CTI computation have been coded in R language93 and available for one of the time series in a public repository (https://doi.org/10.5281/zenodo.10708267, https://doi.org/10.5281/zenodo.10708267).

Inter-annual temporal change rates in CTI (CTIr) were estimated using two approaches: (1) for each site independently, and (2) jointly for all sites to provide an overall estimate of CTI change. For the first approach, we estimated the CTIr for each site. To avoid potential biases in the estimation of CTI due to temporal autocorrelation94, partial autocorrelation has been checked for each yearly data time series. Subsequently, CTIr was estimated fitting linear models using generalised least squares, and adding an autoregressive function when autocorrelation was detected. For the second approach, we tested whether the CTI change through time for all sites on average is different from zero using a linear mixed model of CTI with year as the fixed effect, the sampling sites as the random effect, and adding a temporal autoregressive function. The relation between CTI and sea temperature (surface and 100m integrated column) through time and across sites were tested using linear mixed models with site as random effect, and adding a temporal autoregressive function. t- and p-values of the two-sided Wald test for the estimated slope coefficients are provided.

Tropicalization and borealization of communities
CTI changes were decomposed into four underlying process following McLean et al. 15: tropicalization (increasing warm-affinity species), deborealization (decreasing cold-affinity species), borealization (increasing cold-affinity species), and detropicalization (decreasing warm-affinity species). This categorisation is computed at per-species basis by (i) calculating species’ thermal bias (thermal preference – CTIstation) and species’ abundance change, and (ii) assign the corresponding process to each species, e.g., tropicalization if both species’ thermal bias and abundance change are positive, deborealization if both species’ thermal bias and abundance change are negative. The intensity of the ecological processes underlying the temporal change in CTI at each station and over all datasets were also examined by (i) calculating the difference between each species’ thermal preference and the mean of the community, (ii) multiplying this value by each species’ change in abundance, and (iii) taking the sum of the resulting values for all species within each process15. Thus, underlying process at per-site basis refer to the community, e.g. tropicalization intensity of the community corresponds to an increase in the number of warm-affinity species, and deborealization to a decrease in cold-affinity species. According to the latitudinal position of the sampling, a given species can be a cold-affinity species at a certain latitude and a warm-affinity species at another latitude. Increases in CTI are expected to occur when the combination of tropicalization and deborealization is stronger than the combination of borealization and de-tropicalization.

CTI analysis and underlying process (i.e., tropicalization and borealization) have been coded in R language93 and available in a public repository (https://doi.org/10.5281/zenodo.10708267, https://doi.org/10.5281/zenodo.10708267). As an overall statistic, the percentages of the prevailing underlying process ((de)tropicalization, and (de)borealization) over all biodiversity time series were also calculated.

Analysis across regions and biotic groups
We compared the CTIr of marine communities across different factors (biological groups, habitat, sea, and basin type) and their levels (Biological groups: coastal hard-bottom benthos, coastal soft-bottom benthos, zooplankton, demersal crustaceans, cephalopods, fish; Habitat: benthic / demersal, estuarine, and pelagic; Sea: Baltic (including Kattegat), Mediterranean, Northeast Atlantic, and Basin type: Non-enclosed (Atlantic Ocean), Semi-enclosed (Mediterranean, Baltic Sea, and estuarine types). The comparison was based on linear mixed models of CTI with the interaction of year and factor as the fixed effect, the sampling sites as the random effect, and adding a temporal autoregressive function to test if the CTIr mean of a given level differs from zero. t- and p-values of the two-sided Wald test for the estimated coefficients are provided. To identify the most representative factors, we selected the best model using the Akaike’s Information Criterion corrected (AICc) by comparing all combinations95. Diagnostic plots for the residuals of the selected factors were checked to ensure model reliability.

As positive CTI changes correspond to an increased prevalence in warm-affinity species (i.e., tropicalization), and/or a decrease in cold-affinity species (i.e., deborealization), the aim here is to test if the intensity of tropicalization minus deborealization at per-site basis varies across levels of the different factors. The reason why tropicalization may prevail (or not) over deborealization might depend on biological traits (size, dispersal capacity), and seascape aspects (e.g., marine geomorphological features) limiting species dispersal and ocean connectivity, which has been recently debated. For positive CTIr, we tested different factors. Therefore, we compared the means of intensity of underlying processes of increased CTI (in particular, tropicalization minus borealization) of marine communities across different factors and their levels using linear models96 and selected the best model using the AICc by comparing all combinations, or with forward stepwise selection in case of singularities caused from interaction among factors95). For comparison purposes, we also modelled the intensity of underlying processes at per-species basis with a linear mixed model for CTIr > 0 as a function of factors, with sites as random effect.

It can be also expected that seas with limited connectivity to open oceans would constrain species abundance increases (i.e., tropicalization and borealization) over abundance decreases (i.e., detropicalization and deborealization) processes in both CTI increase or decrease. Therefore, for basin type, we also tested whether species abundance increase (tropicalization+borealization) minus abundance decrease (detropicalization+deborealization) processes varied across levels at per-site basis.

All these analysis have been coded in R language93 and available in a public repository (https://doi.org/10.5281/zenodo.10708267, https://doi.org/10.5281/zenodo.10708267).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a neat paper that attempts to address an important ecological phenomenon. I\nparticularly like the authors' intent to explore CTI change in the light its constituent\ncomponents (topicalisation, demoralisation, etc.). As such, I would really like to see this\nwork published. However, to properly assess the robustness of its findings, more detail\nabout the analytical approaches is required—and depending on the outcome of that\nprocess, potentially also reanalysis of the data.\nI have several major comments, as well as minor editorial suggestions.\nMajor comments:\n1 – There were serval unanswered questions about how the various groupings were\nconstituted and what the group names really meant. Since this is a central focus of the\nmanuscript, more detail/rationale is required here. Line-by-line examples follow.\n*Line 183: Fig. 1b. Were crustaceans and molluscs included among the benthos, if not, why\nnot (I don’t see an answer in the Methods)?\n*Lines 199–202: These results raise a further issue with separating molluscs and crustaceans\nfrom other benthos, unless the other benthos was sessile/burrowing, in which case, better\nnames might be useful.\n2 – There were several issues with the analyses. Although the authors demonstrate\nknowledge of most of the issues (e.g., temporal [and to some degree, spatial]\nautocorrelation, the need for mixed models, etc.), the application of techniques was\ninconsistent, and in some cases, (given the information provided in the manuscript)\nquestionable. The most important of these is the apparent abandonment of the results from\nthe mixed-effects model to produce test statistics in the main results table (Table 1),\nrendering their interpretation open to questions about the actual alpha-level (probability of\ncommitting a type I error) used in the analyses. There are several other issues listed line-by-\nline, below. Some of these issues might be resolved by more careful explanation of the\nrationale and philosophy of the individual analyses (models), but in other cases,\nreconsidering the approaches to analyses might be warranted. Irrespective, it would be very\nuseful to include diagnostic plots of the final fitted models in the Supplement, as well as a\nshort section on caveats in the Discussion.\n*Lines 203–205: Given that there is only one observation from estuaries, would it not be\nmore rigorous to simply omit that datum than to point out that a single point had the\nstrongest trend, albeit non-significant? I am surprised that the algorithm managed to\nestimate a p-value from a single point…although this would have depended on the structure\nof the model, I guess. Nevertheless, the inclusion of this point is somewhat problematic.\n*Line 216–217: Fig. 3 – How were pies computed? On a per-species basis? How were these\nweighted when there were multiple locations/sites? This question further challenges the\nwisdom of including “estuary” as a habitat with only a single observation.\n*Line 244–247: Table 2 and associated results – I’m not convinced that an ANOVA is the\nright way to go here. These are counts, so a log-linear analysis (saturated Poisson glm that\nincludes interactions between grouping factors and process to test for independence of\nfrequencies) seems more defensible…or perhaps a Poisson generalised mixed model that\nalso includes the random effects?\n*Lines 428–429: No consideration of temporal autocorrelation here? This needs to be\naddressed or at least listed as a caveat, with a brief mention of the likely consequences. You\ndo this in Lines 465–467 for CTIr, so why not temperature?\n*Lines 433–439: A more detailed explanation for the soft-bottom/hard-bottom benthos vs\nmollusca/crustacea classifications are needed here, as is more detail on the durations of the\ntime series (at least median, quartiles and range). Data are critical in this analysis, so such\nbrief coverage is problematic.\n*Lines 443–461: Move all of the technical information into a single paragraph before\nexplaining how the metric is used. Otherwise, the reader is left wondering how Ts and A are\ndetermined.\n*Lines 468–471: First, temporal autocorrelation seems to be missing from this analysis.\nSecond, the description here clashes with information in the caption for Figure 2: if\nsite/taxonomic group is the random factor, how can Fig. 2c be showing random effects of\nthe mixed model? For models of biological group, surely the random effect should be site\nnested within group (i.e., the effect of site doesn’t disappear because you are more\ninterested in taxa)? Finally, what is the effect of multiple tests, here? It is true that\ncorrections for multiple tests are usually deployed for t-tests and the like, but the principle\nholds for any sequence of hypothesis tests. Here, you are refitting different models to the\nsame response data. This increases the risk of a Type I error. Given the marginal significance\nof several of the results in Table 1, this becomes more problematic. One solution might be\nto fit a model with all predictors included, then simplify on the basis of AIC or some other\ninformation-theoretic criterion, and interpret only the predictors that remain? Irrespective,\nthis is a serious issue here. Proper diagnostic plots of the final fitted models would also be\nuseful. Finally, since the results depend on these analyses, a clearer and more detailed\nexplanation of the underlying models (including their underlying rationale/philosophy)\nwould really help to clarify what is being done.\n*Lines 475–484: This does not provide any real insight into how the frequency(?) pies were\ncomputed for presentation in Figures 2 & 3. Again, without understanding what these plots\nactually represent, it is impossible to properly assess their meaning.\n*Lines 493–496: Here, it seems as if the mixed-effects models are simply ignored in favour\nof general linear models. Why? What happened to the random effects? Surely it would be\nbetter to interpret the coefficients from the mixed-effects models, which are likely more\nrobust?\n*Lines 498–500: What is the hypothesis here? That CTIr is the same across all levels of the\npredictor? If so, what levels vary from what other levels, and what is the importance of this\ninsight? Also, given the massive differences in sample sizes among levels for most groups, it\nseems unlikely that assumptions of homoscedasticity of variances is met. As mentioned\npreviously, if these numbers are based on counts, it might be better to explore some form\nof Poisson GLM (or better still, GLMM—so that the random effects can be included).\nMinor editorial comments:\nOverall, the manuscript was well constructed, but it would benefit from a thorough copy\nedit. Some of the suggestions below are a matter of taste, but others reflect a lack of\nconsistency in writing.\nLine 1 (and elsewhere): “cross-taxon”, not “cross-taxa”; adjectives like this are singular, not\nplural. A good example is available at the front of your title.\nLine 65 (and elsewhere): add “concentrations” after “oxygen”.\nLines 75–76: Insert commas: “…semi-enclosed basins, such as the Mediterranean and Baltic\nSeas, where physical…”\nLines 106–107: I suggest a slight revision to the sentence “All European Seas HAVE BEEN\nwarming since AT LEAST the 1980’s, but WITH faster RATES in the semi-enclosed basins…”\nLines 112–114: This sentence is tricky. Time series cannot “correspond” to species. And it is\nunclear whether all time series cover the last 4 decades or whether they are simply located\nwithin this timeframe (the same issue exists in the Abstract). Consider revising for clarity.\nLine 119: The phrase “European Seas basins“ reads awkwardly because of the pluralisation\nof “Seas”. Perhaps “…basins within semi-enclosed European Seas…”?\nLine 121: “…than IS THE CASE in the well-connected…”\nLine 123: The concept of a “hotspot of high biodiversity climate velocity” eludes me.\nWording requires some thought here.\nLine 130 and elsewhere: “i.e.” and “e.g.” should preferably be followed by a comma, but as\nit stands, this notion is inconsistently applied, throughout.\nLines 132–134: Are the signals or the communities associated with warming? As written, it\nseems to be the latter?\nLine 150: The use of “their” seems strange, here.\nLines 171–172: Consider “…increased from 1980–2021, both AT the surface…and IN the\nwater column…”\nLines 173: Consider “rate of ocean warming at the surface was…”. Also, note that SST and\nST100 are undefined abbreviations at this point in the MS. I suspect that you have defined\nthese in the Methods, but the arrangement of this MS means that the reader might not\nhave got there, yet.\nLine 186: GODAS is undefined at this point.\nLine 188: Do you mean “…rate of change in CTI…?”\nLine 189: Consider replacing “their” with “corresponding”. Also, “site labelling”, not “sites\nlabelling”\nLines 191–192: This wording seems to suggest that year was included as both a random and\na fixed effect? This needs clarification.\nLines 198–199: Why a mix in singular and plural descriptors in parentheses?\nLines 202–203: “…weaker trends…”? I might be better to frame these as “…trends\nstatistically indistinguishable from zero…”?\nLine 214–215: Table 1 – what are the asterisks? It would also be good to explain the\nbolding.\nLine 226: Most frequent on a per-species basis? Or per location? Or some weighting?\nLine 245: “whether”, not “if”\nLines 262–271: This material might better suit a summary table?\nLines 255–416: I read these, but cannot assess the robustness of interpretation until\nquestions around the analyses are resolved.\nLine 425: “…FOR the PERIOD 1980−2020…”\nLine 426: “on” not “to”\nLines 449–450: What is meant by “A, the species’ abundance divided by the abundances of\nall species”? Do you mean “A, the species’ abundance divided by the total number of\nindividuals at a site”? Would it not be simpler to refer to this as “relative abundance” or\nsimilar?\nLines 460–461: Delete “each”.\nLine 463: A “temporal trend” (= change in a variable per unit time) is by definition a “rate”.\nLines 463–471: After having read the paragraph remain unclear about what is meant by\n“independently” and “jointly”.\nLines 467–468: The meaning of “CTIr mean change was tested using: 1) a t-test weighted by\nthe number of species present in each site…” is unclear to me. Statistical tests require an a\npriori hypothesis. What is the hypothesis about the trend being tested? Is this a t-test of the\nhypothesis that the slope of the relationship is zero?\nLine 486: Consider “Analysis across regions and biotic groups”\nLine 494: “whether” not “if”\nLines 500–502: I might have missed it, but I didn’t see this in the Results?", "reviewer_2": "This manuscript calculates CTI change in European seas and decomposes CTI to detect the\nrelative influence of deborealization vs tropicalization in long-term monitoring datasets. This\nlatter element is the most exciting contribution, given the extensive literature that already\nexists on changing CTI (Burrows et al. 2019), species richness (Blowes et al. 2019),\nabundance change relative to range position (Hastings et al. 2020), and community\ncomposition (Rutterford et al. 2023)—to name a few—in this well-studied region. (Given\nthat all these papers, much like the present manuscript, aim to indirectly gain inference into\necological processes by analyzing long-term surveys, I was surprised to read L92-94 in this\nmanuscript which states that such studies are “scarce”.)\nThe manuscript is well-written and I agree that its goals—trying to understand the degree to\nwhich tropicalization vs. deborealization are driving community change in warming seas—\nare important. However, I have several issues that make it hard to assess the results of the\npaper as is. The first is the conflation of tropicalization with range expansion and\ndeborealization with range contraction, which is captured in Fig. 1 and motivates much of\nthe analysis. This assumes that the abundant-center hypothesis is true, which has very\nlimited empirical support (Dallas et al. 2020) and the authors do not test for the species\nstudied. For this reason, I was very confused by the attempt to back out retraction and\nexpansion processes from CTI change (L498-502). The methods describing this are brief and\nnot referenced (and it’s not clear if this is expansion/retraction of populations, or the whole\nrange, in which case the authors should ensure that they are actually capturing range edges\nin the surveys), so maybe there is a very clear rationale that I just don’t know. Without that\nrationale, my understanding is: when we see tropicalization and deborealization, it is just\nthat—abundance change of species with relatively higher or lower thermal affinities—and\nnot necessarily range change unless we actually go measure range shifts.\nRelatedly, while I understand the authors may be up against a limit to references, I\nencourage you to unpack and support the hypotheses a bit more (e.g., L143-148). I disagree\nthat tropicalization and borealization processes are “poorly understood”; there are lots of\nlarge-scale studies and predictions to draw from here. For example, it’s possible that cold\nrange edges of marine fishes are tracking temperature better than warm edges (a long-\nstanding prediction of terrestrial biogeography theory), which would predominantly lead to\ntropicalization (Fredston‐Hermann et al. 2020). On the other hand, some fish habitat models\nproject faster climate velocities at trailing edges than at leading edges (Robinson et al.\n2015), and some terrestrial species have even greater “colonization credit” than they have\n“extinction debt” (Talluto et al. 2017), leading to range contraction and therefore suggesting\npredominant deborealization.\nMy second concern is that there are extensive deficiencies in this manuscript’s reporting\nsummary and data/code availability relative to the Nature Portfolio standards and our field’s\nbest practices for open and transparent science which make it challenging to assess the\nmethods and impossible to reproduce them. Unless I missed it, sample sizes are not\nreported in the main text. The reporting summary points to Supp. Tab. 1 but I can’t find\nsample sizes there either. They should be reported as the number of records actually\nanalyzed, not just number of species per dataset, and broken down by system and clade so\nreaders can understand which groups and surveys drove the results. The code and raw data\nare also missing. I understand that not all the data are public but at a minimum the authors\nshould share all of the code as well as the subset of the raw and processed datasets that\nthey are able to publish. A data availability statement is also missing from the manuscript. In\nthe reporting summary section on code availability, the authors pasted a citation for R\nrather than access information for their code. In the reporting summary section on data\navailability, the authors stated that raw data for some surveys are unavailable and that the\nprocessed data is described in a table, rather than providing access information for the raw\ndata that are available.\nI have a number of questions about the data analysis and statistics—which are really the\nheart of this paper—which are not answered in the brief section summarizing them (L452-\n471). Were the OBIS records filtered or quality-checked at all? How exactly were they\nmerged with GODAS (the text says “annual means”—is that the year of the OBIS record?)\nWere models fitted to individual surveys or were surveys pooled? If the former, how did the\nauthors correct for multiple testing and confirm that each individual time-series was not\nunderpowered? How precisely was partial autocorrelation “checked for each yearly data\ntime series”? Where are the model results reported and how many models were there?\nWhat were the two groups being compared in the t-test? Did the authors analyze the effect\nof sampling intensity / time-series length / sample size on the significance of each individual\nCTI trend?\nFinally, I’m not convinced that this post-hoc ANOVA approach and the pie and violin charts\nin Fig. 3-4 are the best way to explore which species drove CTI change. More plots of the\nraw data, like a main text figure showing species’ thermal bias (thermal affinity – CTI) vs.\nabundance change, would really help to show this visually. A single species’ influence on CTI\nis a combination of its abundance (since CTI in this analysis is abundance-weighted, I think)\nand the magnitude of its thermal bias. Ranking species by their influence on the CTI score\nand then reporting whether the changes in the most influential species are consistent with\ndeborealization / tropicalization / etc and whether abundance change or differential\nthermal bias seem to be most important would keep the statistics “closer to the data” than\ndoing an ANOVA on coefficients from a regression, and would enable the authors to unpack\nthe underlying processes they focus on more than the current methods allow. I’m sure there\nare lots of other ways to think about this but I definitely encourage the authors to take a\nstep back, connect the statistical approach with the processes they want to measure, focus\nmore on effect sizes and which species are driving results, and less on p-values.\nBecause of these fundamental questions, I haven’t really reviewed the results and\ndiscussion here. As I said earlier, I think this is an important research direction and an\nexciting dataset, and I would be happy to review a revised manuscript that resolves these\nmethodological issues. I don’t think these issues are insurmountable: expanding on the\nmethods, sharing the data and code, and refocusing the text and results on how shifts in\nabundance of different species are driving different patterns of community turnover in\nthese regions will turn this into a great paper."}, "author_response": "AR: The terminology employed was not clear as the reviewer pointed out. In benthos, there\nare crustaceans and molluscs among other macroinvertebrates, but these organisms were\nmainly sessile (in hard-bottom) or inhabiting the sediment and located in coastal areas (most\nof them at intertidal areas), while the groups of crustaceans and molluscs (now named\ncephalopods more precisely) were collected offshore using fishing or bottom trawl surveys,\nhence representing different benthic communities. We therefore renamed the biological\ngroups and defined them in Methods section (lines 466-479):\n1. Fish: pelagic, demersal and estuarine fish species.\n2. Crustacea was renamed to demersal crustaceans: Demersal crustaceans collected\noffshore using bottom trawl fishing and surveys.\n3. Mollusca was renamed to Cephalopods: Demersal and pelagic cephalopods collected\noffshore using bottom trawl and purse-seine fishing and surveys.\n4. SB benthos was renamed to Coastal soft-bottom benthic macroinvertebrates: benthic\nmacroinvertebrates species in coastal intertidal or subtidal soft-bottom sites, using\ngrab or box corer samples.\n5. HB benthos was renamed to Coastal hard-bottom benthos: benthic\nmacroinvertebrates, macroalgae, lichens, and coralligenous assemblages in intertidal\nor subtidal, coastal, hard-bottom sites\n6. Zooplankton: includes copepods and other zooplankton species collected by vertical\nnet hauls\n2 – There were several issues with the analyses. Although the authors demonstrate knowledge\nof most of the issues (e.g., temporal [and to some degree, spatial] autocorrelation, the need\nfor mixed models, etc.), the application of techniques was inconsistent, and in some cases,\n(given the information provided in the manuscript) questionable. The most important of these\nis the apparent abandonment of the results from the mixed-effects model to produce test\nstatistics in the main results table (Table 1), rendering their interpretation open to questions\nabout the actual alpha-level (probability of committing a type I error) used in the analyses.\nThere are several other issues listed line-by-line, below. Some of these issues might be\nresolved by more careful explanation of the rationale and philosophy of the individual analyses\n(models), but in other cases, reconsidering the approaches to analyses might be warranted.\nIrrespective, it would be very useful to include diagnostic plots of the final fitted models in the\nSupplement, as well as a short section on caveats in the Discussion.\nAR: In both temperature and CTI change over time analysis we now used linear mixed models\n(Zuur et al. 2009) consistently with sampling site as the random effect and year as the fixed\neffect, and adding a temporal autoregressive function (Mudelsee, 2019). The analysis of CTIr\nby factors (biological group, basin type, habitat, region) (Table 1) was now undertaken using\nalso consistent linear mixed models with temporal autoregressive function. Further we added\na model selection analysis to identify the most representative factors underpinning CTI\nchanges over time and addressing type I error, in particular, using the Akaike’s Information\nCriterion corrected (AICc) by comparing all model combinations (Burnham & Anderson, 2002)\n(see Methods section, lines 504-516 and 547-558).\nDiagnostic plot: As suggested, we performed diagnostic plots of residuals for the two selected\nfactors (See new Supplementary Figure 3). Diagnostic plots did not show any major pattern,\nindicating that the model is reliable and adequate for this analysis.\n*Lines 203–205: Given that there is only one observation from estuaries, would it not be more\nrigorous to simply omit that datum than to point out that a single point had the strongest\ntrend, albeit non-significant? I am surprised that the algorithm managed to estimate a p-value\nfrom a single point…although this would have depended on the structure of the model, I\nguess. Nevertheless, the inclusion of this point is somewhat problematic.\nAR: We included two new estuarine time series (“INRAE fish Pertuis” and “INRAE fish Basque”,\nsee Supplementary Table 1) thanks to efforts in compiling new data and research\ncollaboration. Now, three estuarine time series have been analysed with the improved and\nmore consistent time series analysis. The two included estuarine time series showed also\npositive increase in CTI over time, as the first original one.\n*Line 216–217: Fig. 3 – How were pies computed? On a per-species basis? How were these\nweighted when there were multiple locations/sites? This question further challenges the\nwisdom of including “estuary” as a habitat with only a single observation.\nAR: Pie charts were computed estimating the percentage of the prevailing underlying\necological process ((de)tropicalization, and (de)borealization) for each case study and factor\n(Biological group, habitat, region, basin type). (see caption of Figure 3). As mentioned before,\nnow there are three estuarine time series.\n*Line 244–247: Table 2 and associated results – I’m not convinced that an ANOVA is the right\nway to go here. These are counts, so a log-linear analysis (saturated Poisson glm that includes\ninteractions between grouping factors and process to test for independence of frequencies)\nseems more defensible…or perhaps a Poisson generalised mixed model that also includes the\nrandom effects?\nAR: In Table 2, the variable “Tropicalization minus Deborealization” are not counts, it is a\ndifference between percentages. Also, the underlying ecological processes (tropicalization,\ndeborealization, …) are computed for each site and over the entire time period, hence, each\ntime series yield a single value of each underlying process. For these reasons, we compared\nthe means of theunderlying ecological processes intensities when CTIr increased (in particular,\ntropicalization minus borealization) across different factors and their levels. To do so, we used\nlinear models (Chambers, 1992) and added a selection of the best model using the AICc by\ncomparing all factor combinations, or with forward stepwise selection in case of singularities\namong factors (Burnham & Anderson, 2002) (see section Methods lines 560-574). For\ncomparison purposes, we also modelled the intensity of underlying processes at per-species\nbasis with a linear mixed model for CTIr>0 as a function of factors, with sites as random effect\n(lines 572-574); both analyses provide similar results (see new Supplementary Table 2, shown\nbelow).\n*Lines 428–429: No consideration of temporal autocorrelation here? This needs to be\naddressed or at least listed as a caveat, with a brief mention of the likely consequences. You do\nthis in Lines 465–467 for CTIr, so why not temperature?\nAR: Accordingly with this comment and previous ones in relation to the consistency of analysis,\nwe analysed the slope of the temperature trends and its significance using a linear mixed\nmodel (Zuur et al. 2009) of sea temperature with year as fixed effect, the sampling sites as\nrandom effect, and adding an autoregressive function to take into consideration the temporal\nautocorrelation (Mudelsee, 2019) (See Method section lines 551-559).\n*Lines 433–439: A more detailed explanation for the soft-bottom/hard-bottom benthos vs\nmollusca/crustacea classifications are needed here, as is more detail on the durations of the\ntime series (at least median, quartiles and range). Data are critical in this analysis, so such brief\ncoverage is problematic.\nAR: as mentioned above, we renamed the biological groups and defined in more detail in\nMethods section (lines 466-479):\n1. Fish: pelagic, demersal and estuarine fish species.\n2. Crustacea was named to demersal crustaceans: Demersal crustaceans collected\noffshore using bottom trawl fishing and surveys.\n3. Mollusca was named Cephalopods: Demersal and pelagic cephalopods collected\noffshore using bottom trawl and purse-seine fishing and surveys.\n4. SB benthos was named Coastal soft-bottom benthic macroinvertebrates: benthic\nmacroinvertebrates species in coastal intertidal or subtidal soft-bottom sites, using\ngrab or box corer samples.\n5. HB benthos was named Coastal hard-bottom benthos: benthic\nmacroinvertebrates, macroalgae, lichens, and coralligenous assemblages in intertidal or subtidal,\ncoastal, hard-bottom sites.\n6. Zooplankton: includes copepods and other zooplankton species collected by vertical\nnet hauls.\nThe time series analysed here have varying durations over the last four decades (from 1980 to\n2022 in the longest case). Most of time series span more than 15 years (92%), with median of\n25 years, minimum of 9 years, and 20 and 39 years the 25% and 75% quartiles, respectively\n(see Methods section (lines 474-479).\n*Lines 443–461: Move all of the technical information into a single paragraph before\nexplaining how the metric is used. Otherwise, the reader is left wondering how Ts and A are\ndetermined.\nAR: This paragraph was rewritten accordingly (see lines 483-502).\n*Lines 468–471: First, temporal autocorrelation seems to be missing from this analysis.\nAR: Temporal autocorrelation has been now taken into consideration in the new linear mixed\nmodel analysis (see above and section Methods).\nSecond, the description here clashes with information in the caption for Figure 2: if\nsite/taxonomic group is the random factor, how can Fig. 2c be showing random effects of the\nmixed model? For models of biological group, surely the random effect should be site nested\nwithin group (i.e., the effect of site doesn’t disappear because you are more interested in\ntaxa)?\nAR: CTI change over time analysis was now based on linear mixed models (Zuur et al. 2009)\nconsistently with sampling site as the random effect and year as the fixed effect, and adding a\ntemporal autoregressive function (Mudelsee, 2019) (see Methods section, lines 547-558). In\nFig. 2c, partial residuals of CTI overtime based on the mentioned model are showed (axis title\nclarified) and legend explanation was rewritten to: “partial residuals of CTI were calculated as\nCTI minus the random effect of sampling site of the linear mixed model with year as fixed\neffect and site as random effect”, see lines 207 (Fig 2 caption). Hopefully, the figure is clearer\nnow.\nFinally, what is the effect of multiple tests, here? It is true that corrections for multiple tests\nare usually deployed for t-tests and the like, but the principle holds for any sequence of\nhypothesis tests. Here, you are refitting different models to the same response data. This\nincreases the risk of a Type I error. Given the marginal significance of several of the results in\nTable 1, this becomes more problematic. One solution might be to fit a model with all\npredictors included, then simplify on the basis of AIC or some other information-theoretic\ncriterion, and interpret only the predictors that remain? Irrespective, this is a serious issue\nhere. Proper diagnostic plots of the final fitted models would also be useful. Finally, since the\nresults depend on these analyses, a clearer and more detailed explanation of the underlying\nmodels (including their underlying rationale/philosophy) would really help to clarify what is\nbeing done.\nAR: As mentioned above, we now used linear mixed models (Zuur et al. 2009) to analyse CTI\nchange consistently with sampling site as the random effect and year as the fixed effect, and\nadding a temporal autoregressive function (Mudelsee, 2019). The analysis of CTIr by groups\n(Table 1) was now undertaken using also consistent linear mixed models with temporal\nautoregressive function, and we added model selection to identify the most representative\nfactors and addressing type I error, in particular, using the Akaike’s Information Criterion\ncorrected (AICc) by comparing all model combinations (Burnham & Anderson, 2002) (see\nMethods section with more detailed explanations, lines 547-558).\nAs suggested, we performed diagnostic plots of residuals for the two selected factors of the\nlinear mixed model (See new Supplementary Figure 3, showed below). Diagnostic plots did not\nshow any major pattern, indicating that the model is reliable and adequate.\n*Lines 475–484: This does not provide any real insight into how the frequency(?) pies were\ncomputed for presentation in Figures 2 & 3. Again, without understanding what these plots\nactually represent, it is impossible to properly assess their meaning.\nAR: For Fig 2d pie chart, we computed the percentages of the prevailing underlying process\n((de)tropicalization, and (de)borealization) over all biodiversity time series (text added in lines\n217-219, Fig 2 caption). Pie charts of Figure 3 were computed estimating the percentage of the\nprevailing underlying process ((de)tropicalization, and (de)borealization) in each case study\nand by factor (Biological group, habitat, region, basin type) (see legend of Figure 3).\n*Lines 493–496: Here, it seems as if the mixed-effects models are simply ignored in favour of\ngeneral linear models. Why? What happened to the random effects? Surely it would be better\nto interpret the coefficients from the mixed-effects models, which are likely more robust?\nAR: As mentioned above, we now used linear mixed models (Zuur et al. 2009) to analyse CTI\nchange consistently with sampling site as the random effect and year as the fixed effect, and\nadding a temporal autoregressive function (Mudelsee, 2019). The analysis of CTIr by groups\n(Table 1) was now undertaken using also consistent linear mixed models with temporal\nautoregressive function, and we added model selection to identify the most representative\nfactors and addressing type I error, in particular, using the Akaike’s Information Criterion\ncorrected (AICc) by comparing all model combinations (Burnham & Anderson, 2002) (see\nMethods section, lines 547-558).\n*Lines 498–500: What is the hypothesis here? That CTIr is the same across all levels of the\npredictor? If so, what levels vary from what other levels, and what is the importance of this\ninsight? Also, given the massive differences in sample sizes among levels for most groups, it\nseems unlikely that assumptions of homoscedasticity of variances is met. As mentioned\npreviously, if these numbers are based on counts, it might be better to explore some form of\nPoisson GLM (or better still, GLMM—so that the random effects can be included).\nAR: As positive CTI changes correspond to an increased prevalence in warm-affinity species\n(i.e., tropicalization), and/or a decrease in cold-affinity species (i.e., deborealization), here, the\naim is to test if tropicalization minus deborealization varies across levels of the different\nfactors (text now added in Methods lines 560-565). As mentioned above, the variable\n“Tropicalization minus Deborealization” are not counts, it is a difference between percentages.\nAlso, the underlying processes (tropicalization, deborealization, …) are computed for each site\nand over the time period, hence, each time series yield a single value for each underlying\nprocess. For these reasons, the analysis was based on the comparison of the means of\nintensity of underlying processes (tropicalization minus borealization), for the set of time\nseries where CTIr increased, across different factors and their levels using linear models\n(Chambers, 1992), and we selected the best model using the AICc by comparing all factor\ncombinations, or with forward stepwise selection in case of singularities among factors\n(Burnham & Anderson, 2002) (see section Methods lines 560-574).\nConcerning homoscedasticity of variances, we undertook a Levene test, all groups were not\nsignificant (p>0.05) indicating equal variances, except Basin group with p=0.0217. However, we\nperformed a non-parametric test (Kruskal Wallis test) that indicates, in agreement with the\nlinear model and AICc, that Bain type is the most significant factor (p<0.0001) and that\ntropicalization prevailed deborealization in non-enclosed seas.\nMinor editorial comments:\nOverall, the manuscript was well constructed, but it would benefit from a thorough copy edit.\nSome of the suggestions below are a matter of taste, but others reflect a lack of consistency in\nwriting.\nLine 1 (and elsewhere): “cross-taxon”, not “cross-taxa”; adjectives like this are singular, not\nplural. A good example is available at the front of your title.\nAR: corrected, “Cross-basins and cross-taxon patterns in biodiversity turnover in warming\nseas”\nLine 65 (and elsewhere): add “concentrations” after “oxygen”.\nAR: added\nLines 75–76: Insert commas: “…semi-enclosed basins, such as the Mediterranean and Baltic\nSeas, where physical…”\nAR: added\nLines 106–107: I suggest a slight revision to the sentence “All European Seas HAVE BEEN\nwarming since AT LEAST the 1980’s, but WITH faster RATES in the semi-enclosed basins…”\nAR: corrected\nLines 112–114: This sentence is tricky. Time series cannot “correspond” to species. And it is\nunclear whether all time series cover the last 4 decades or whether they are simply located\nwithin this timeframe (the same issue exists in the Abstract). Consider revising for clarity.\nAR: CTI was applied to long-term time series from 63 biodiversity programs spanning the last\nfour decades (in the longest case), which accounts for 1806 species, including zooplankton,\nbenthos, demersal and pelagic assemblages\nLine 119: The phrase “European Seas basins“ reads awkwardly because of the pluralisation of\n“Seas”. Perhaps “…basins within semi-enclosed European Seas…”?\nAR: rewritten: “semi-enclosed sea basins with lower ocean connectivity to warmer waters are\nhypothesized to experience less tropicalization compared to deborealization than is the case in\nthe well-connected northeast Atlantic region”\nLine 121: “…than IS THE CASE in the well-connected…”\nAR: corrected, see previous comment\nLine 123: The concept of a “hotspot of high biodiversity climate velocity” eludes me. Wording\nrequires some thought here.\nAR: we rewrote the sentence to clarify “biodiversity hotspots associated to climate velocity”\nLine 130 and elsewhere: “i.e.” and “e.g.” should preferably be followed by a comma, but as it\nstands, this notion is inconsistently applied, throughout.\nAR: revised commas and use of terms; e.g. as “for example” and i.e. as “that is”\nLines 132–134: Are the signals or the communities associated with warming? As written, it\nseems to be the latter?\nAR: revised: “Long-term monitoring programs based on permanent stations have detected\nsuch changes in North Atlantic fish communities, based on temporal changes of the CTI”\nLine 150: The use of “their” seems strange, here.\nAR: changed to “its”\nLines 171–172: Consider “…increased from 1980–2021, both AT the surface…and IN the water\ncolumn…”\nAR: corrected\nLines 173: Consider “rate of ocean warming at the surface was…”. Also, note that SST and\nST100 are undefined abbreviations at this point in the MS. I suspect that you have defined\nthese in the Methods, but the arrangement of this MS means that the reader might not have\ngot there, yet.\nAR: changed and defined here\nLine 186: GODAS is undefined at this point.\nAR: defined\nLine 188: Do you mean “…rate of change in CTI…?”\nAR: corrected\nLine 189: Consider replacing “their” with “corresponding”. Also, “site labelling”, not “sites\nlabelling”\nAR: corrected\nLines 191–192: This wording seems to suggest that year was included as both a random and a\nfixed effect? This needs clarification.\nAR: The linear mixed model includes year as fixed effect and site as random effect. We rewrote\nthe sentence to clarify this as follows: “Partial residuals of CTI across time, calculated as CTI\nminus the random effect of sampling site of the mixed model, see Methods and Fig 2 caption”\n(lines 207-2018)\nLines 198–199: Why a mix in singular and plural descriptors in parentheses?\nAR: corrected, all in singular\nLines 202–203: “…weaker trends…”? I might be better to frame these as “…trends\nstatistically indistinguishable from zero…”?\nAR: changed\nLine 214–215: Table 1 – what are the asterisks? It would also be good to explain the bolding.\nAR: asterisks are removed, bolding (significant values, p<0.05) is indicated now.\nLine 226: Most frequent on a per-species basis? Or per location? Or some weighting?\nAR: on a per-site basis, now indicated\nLine 245: “whether”, not “if”\nAR: changed\nLines 262–271: This material might better suit a summary table?\nAR: This information covers some of our results and it is compared with specific studies of sea\nwarming\nLines 255–416: I read these, but cannot assess the robustness of interpretation until\nquestions around the analyses are resolved.\nAR: Discussion has now been slightly changed according to the new results\nLine 425: “…FOR the PERIOD 1980−2020…”\nAR: corrected\nLine 426: “on” not “to”\nAR: corrected\nLines 449–450: What is meant by “A, the species’ abundance divided by the abundances of all\nspecies”? Do you mean “A, the species’ abundance divided by the total number of individuals\nat a site”? Would it not be simpler to refer to this as “relative abundance” or similar?\nAR: changed: “where n is the number of species in the community, T is the temperature\ns\npreference of each species (s) and A is the relative abundance of species s (i.e., the abundance\ns\nof species s divided by the total number of individuals in the community at a site)”\nLines 460–461: Delete “each”.\nAR: deleted\nLine 463: A “temporal trend” (= change in a variable per unit time) is by definition a “rate”.\nAR: “rate” deleted\nLines 463–471: After having read the paragraph remain unclear about what is meant by\n“independently” and “jointly”.\nAR: hopefully now is clarified: “Temporal trend of CTI were estimated: 1) for each site\nindependently, and 2) jointly for all sites to provide an overall estimate of CTI change.”\nLines 467–468: The meaning of “CTIr mean change was tested using: 1) a t-test weighted by\nthe number of species present in each site…” is unclear to me. Statistical tests require an a\npriori hypothesis. What is the hypothesis about the trend being tested? Is this a t-test of the\nhypothesis that the slope of the relationship is zero?\nAR: We have rewritten the whole paragraph to clarify the methods:\n“Temporal trends of CTI were estimated using two approaches: 1) for each site independently,\nand 2) jointly for all sites to provide an overall estimate of CTI change. 1) We estimated the CTI\nr\nfor each site. To avoid potential biases in the estimation of CTI due to temporal\nr\nautocorrelation, partial autocorrelation has been checked for each yearly data time series.\nSubsequently, CTI was estimated fitting linear models using generalized least squares, and\nr\nadding an autoregressive function when autocorrelation was detected. 2) We tested whether\nthe CTI change for all sites on average is different from zero using two statistical methods: a) a\nt-test weighted by the number of species present in each site, and b) a fit of a linear mixed\nmodel of CTI with year as fixed effect and sampling sites as random effect. Correlation\nbetween CTI and sea temperature (surface and 100 m integrated column) at each site were\nr\ntested using linear mixed model with biological group as random effect. ”\nLine 486: Consider “Analysis across regions and biotic groups”\nAR: changed\nLine 494: “whether” not “if”\nAR: changed\nLines 500–502: I might have missed it, but I didn’t see this in the Results?\nAR: shown in new Figure 5\nReviewer #2 (Remarks to the Author):\nThis manuscript calculates CTI change in European seas and decomposes CTI to detect the\nrelative influence of deborealization vs tropicalization in long-term monitoring datasets. This\nlatter element is the most exciting contribution, given the extensive literature that already\nexists on changing CTI (Burrows et al. 2019), species richness (Blowes et al. 2019), abundance\nchange relative to range position (Hastings et al. 2020), and community composition\n(Rutterford et al. 2023)—to name a few—in this well-studied region. (Given that all these\npapers, much like the present manuscript, aim to indirectly gain inference into ecological\nprocesses by analyzing long-term surveys, I was surprised to read L92-94 in this manuscript\nwhich states that such studies are “scarce”.)\nAR: Thanks for the references. Those references not cited in the original version (Blowes et al.\n2019, Hastings et al. 2020) were now added (together with others (Pinsky et al. 2019, Smale et\nal. 2019) in the new version to improve the context and discussion. We agree that substantial\nnumber of studies exist exploring the ecological processes explaining community response to\nclimate change. We are specifically interested in explaining deborealization and tropicalization,\nwhich is not studied so extensively and comparatively. We changed the sentence: “Although\nour knowledge of marine species responses to climate change is substantial, at least in the\nNorth Atlantic Ocean (Philippart et al. 2011, Peck and Pinnegar 2018), comparative\nassessments to identify the underlying ecological processes (McLean et al. 2021) associated\nwith marine community expansion, retraction and dispersal constraints are limited to unevenly\ndistribution of monitoring programs (Poloczanska et al. 2016, Bowler et al. 2017, Blowes et al.\n2019, Burrows et al. 2019, Burrows et al. 2020)”.\nThe manuscript is well-written and I agree that its goals—trying to understand the degree to\nwhich tropicalization vs. deborealization are driving community change in warming seas—\nare important. However, I have several issues that make it hard to assess the results of the\npaper as is. The first is the conflation of tropicalization with range expansion and\ndeborealization with range contraction, which is captured in Fig. 1 and motivates much of\nthe analysis. This assumes that the abundant-center hypothesis is true, which has very\nlimited empirical support (Dallas et al. 2020) and the authors do not test for the species\nstudied. For this reason, I was very confused by the attempt to back out retraction and\nexpansion processes from CTI change (L498-502). The methods describing this are brief and\nnot referenced (and it’s not clear if this is expansion/retraction of populations, or the whole\nrange, in which case the authors should ensure that they are actually capturing range edges\nin the surveys), so maybe there is a very clear rationale that I just don’t know. Without that\nrationale, my understanding is: when we see tropicalization and deborealization, it is just\nthat—abundance change of species with relatively higher or lower thermal affinities—and\nnot necessarily range change unless we actually go measure range shifts.\nAR: We agree with the reviewer on that we are not actually testing geographical range change,\nand focusing our analysis in abundance change of species relative to their thermal affinities.\nNow we changed the use of “expansion” term to simply “tropicalization + borealization” and\n“retraction” to “detropicalization + deborealization”. We also expanded the hypothesis in\nMethods (lines 560-574): “As positive CTI changes correspond to an increased prevalence in\nwarm-affinity species (i.e., tropicalization), and/or a decrease in cold-affinity species (i.e.,\ndeborealization), the aim here is to test if tropicalization minus deborealization varies across\nlevels of the different factors. The reason why tropicalization may prevail (or not) over\ndeborealization might depend on biological traits (size, dispersal capacity), and seascape\naspects (e.g., marine geomorphological features) limiting species dispersal and ocean\nconnectivity, which has been scarcely tested with observations. For positive CTIr, we hence\ntested different factors. […] It can be also expected that seas with limited connectivity to open\noceans would constrain species abundance increases (i.e., tropicalization and borealization)\nover abundance decreases (i.e., detropicalization and deborealization) processes in both CTI\nincrease or decrease. Hence, for basin type, we tested whether (tropicalization + borealization)\nminus (detropicalization + deborealization) processes varied across levels for all dataset (i.e.\nCTIr positive and negative)”.\nRelatedly, while I understand the authors may be up against a limit to references, I\nencourage you to unpack and support the hypotheses a bit more (e.g., L143-148). I disagree\nthat tropicalization and borealization processes are “poorly understood”; there are lots of\nlarge-scale studies and predictions to draw from here. For example, it’s possible that cold\nrange edges of marine fishes are tracking temperature better than warm edges (a long-standing\nprediction of terrestrial biogeography theory), which would predominantly lead to\ntropicalization (Fredston‐Hermann et al. 2020). On the other hand, some fish habitat models\nproject faster climate velocities at trailing edges than at leading edges (Robinson et al. 2015),\nand some terrestrial species have even greater “colonization credit” than they have\n“extinction debt” (Talluto et al. 2017), leading to range contraction and therefore suggesting\npredominant deborealization.\nAR: We agree on the comments and accordingly we expanded and rephrased this paragraph:\n“The reason why tropicalization may prevail (or not) over deborealization has been recently\ndebated (Robinson et al. 2015, Talluto et al. 2017, Fredston-Hermann et al. 2020). Assuming an\necophysiological equilibrium between habitat suitability and the occurrence of species, we can\nexpect that marine species track temperature change equally (Fredston-Hermann et al. 2020;\nSunday et al. 2012), hence tropicalization would be similar to deborealization. However, in\nspecies with slow demography and limited dispersal, lags between climate change and\ndistribution shifts can result in ‘extinction debts’ (Jablonski 2001), where populations\ntemporarily persist under unsuitable conditions, and ‘colonization credits’, where suitable\nlocations are not occupied (Talluto et al. 2017). This latter pattern has been observed for\ninstance in trees (Talluto et al. 2017), whilst other studies found equally responsive shifts at\nboth range boundaries in marine ectotherms (Sunday et al. 2012), and prevalence in\ntropicalization in demersal fish (McLean et al. 2021). The prevalence in tropicalization vs\ndeborealization might depend on biological traits of the community (size, dispersal capacity),\nand also seascape aspects (e.g., marine geomorphological features) limiting species dispersal\nand ocean connectivity. Thus, it can be expected that seas with limited connectivity to open\noceans would constrain species abundance increases (tropicalization and borealization) over\nspecies abundance decreases (detropicalization and deborealization).”\nMy second concern is that there are extensive deficiencies in this manuscript’s reporting\nsummary and data/code availability relative to the Nature Portfolio standards and our field’s\nbest practices for open and transparent science which make it challenging to assess the\nmethods and impossible to reproduce them. Unless I missed it, sample sizes are not reported\nin the main text. The reporting summary points to Supp. Tab. 1 but I can’t find sample sizes\nthere either. They should be reported as the number of records actually analyzed, not just\nnumber of species per dataset, and broken down by system and clade so readers can\nunderstand which groups and surveys drove the results. The code and raw data are also\nmissing. I understand that not all the data are public but at a minimum the authors should\nshare all of the code as well as the subset of the raw and processed datasets that they are able\nto publish. A data availability statement is also missing from the manuscript. In the\nreporting summary section on code availability, the authors pasted a citation for R\nrather than access information for their code. In the reporting summary section on data\navailability, the authors stated that raw data for some surveys are unavailable and that the\nprocessed data is described in a table, rather than providing access information for the raw\ndata that are available.\nAR: All codes are available in a public repository (https://zenodo.org/records/10149019, DOI:\n10.5281/zenodo.10149018). Concerning data availability, Data generated during the study and\nthat support their findings (all CTI time series, i.e. CTI per year for each site, and all underlying\nprocess scores ((de)tropicalization, (de)borealization) at per-species and per-site basis) are\navailable in a public repository (https://zenodo.org/records/10149019, DOI:\n10.5281/zenodo.10149018), and in Supplementary Table 1. Biodiversity original data (i.e.,\nspecies abundance at each year for each site survey) is subject to restrictions as it pertains to\nthe corresponding institution. Certain data is, however, available from the corresponding\nauthor upon request.\nSample size encompasses the number of species to calculate CTI, and time period span to\nestimate CTI change. With the new analysis, this means that 65 sites account for 1730\n“samples”. In most of datasets, the abundance is estimated with number of individuals, but in\nsome datasets the abundance is determined by semi-quantitative sampling of surface\ncoverage similar to Braun-Blanquet. However, the reader can now check species list in\navailable dataset.\nI have a number of questions about the data analysis and statistics—which are really the heart\nof this paper—which are not answered in the brief section summarizing them (L452-471).\nWere the OBIS records filtered or quality-checked at all? How exactly were they merged with\nGODAS (the text says “annual means”—is that the year of the OBIS record?) Were models\nfitted to individual surveys or were surveys pooled? If the former, how did the authors correct\nfor multiple testing and confirm that each individual time-series was not underpowered? How\nprecisely was partial autocorrelation “checked for each yearly data time series”? Where are\nthe model results reported and how many models were there? What were the two groups\nbeing compared in the t-test? Did the authors analyze the effect of sampling intensity / time-\nseries length / sample size on the significance of each individual CTI trend?\nAR: OBIS records were quality-checked removing duplicates. To characterize the thermal\npreferences of the species local temperatures derived from GODAS were used, which was\navailable from OBIS for each observation. This information has been explicitly included (line\n490-502).\nAR: As explained in previous comments, we now used linear mixed models (Zuur et al. 2009) to\nanalyse CTI change consistently with sampling site as the random effect (i.e. all surveys\nanalysed together) and year as the fixed effect, and including a temporal autoregressive\nfunction (Mudelsee, 2019), and adding model selection to identify the most representative\nfactors and addressing type I error using the Akaike’s Information Criterion corrected (AICc) by\ncomparing all model combinations (Burnham & Anderson, 2002) (see Methods section, lines\n504-516 and 547-558). In this way, we avoid multiple testing now.\nFinally, I’m not convinced that this post-hoc ANOVA approach and the pie and violin charts\nin Fig. 3-4 are the best way to explore which species drove CTI change. More plots of the\nraw data, like a main text figure showing species’ thermal bias (thermal affinity – CTI) vs.\nabundance change, would really help to show this visually. A single species’ influence on CTI is\na combination of its abundance (since CTI in this analysis is abundance-weighted, I think) and\nthe magnitude of its thermal bias. Ranking species by their influence on the CTI score and then\nreporting whether the changes in the most influential species are consistent with\ndeborealization / tropicalization / etc and whether abundance change or differential thermal\nbias seem to be most important would keep the statistics “closer to the data” than doing an\nANOVA on coefficients from a regression, and would enable the authors to unpack the\nunderlying processes they focus on more than the current methods allow. I’m sure there are\nlots of other ways to think about this but I definitely encourage the authors to take a step\nback, connect the statistical approach with the processes they want to measure, focus more\non effect sizes and which species are driving results, and less on p-values.\nAR: The analysis at per-species suggested by the reviewer is indeed interesting and confirmed\nthe original analysis at per-site basis. We addressed this idea and plot species’ thermal bias\n(thermal affinity – CTI) vs species’ abundance change (see new Figure 4, showed below). The\nextreme tropicalized and deborealized species have been indicated in Figure 4 and discussed in\nDiscussion section (lines 383-401). From underlying processes calculated at species level, we\ncan now calculate the intensity of process at overall and per-site basis in order to provide\nestimates of the process at community level and analysis the potential effect of factors. To do\nthat, we compared the means of intensity of underlying processes of increased CTIr (in\nparticular, tropicalization minus borealization) of marine communities across different factors\nand their levels using linear models and selected the best model using the AICc by comparing\nall combinations, or with forward stepwise selection in case of singularities caused from\ninteraction among factors. For comparison purposes to check if packing or unpacking data (as\nsuggested by the reviewer) yield similar or different results, we also modelled the intensity of\nunderlying processes at per-species basis with a linear mixed model for CTIr>0 as a function of\nfactors, with sites as random effect (lines 276, 572). As mentioned, both analyses provide\nsimilar results (see new Supplementary Table 2)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done a good job in addressing my initial queries. I have only a few minor\ncomments, listed below.\nTitle: “Cross-basin” singular, like, “cross-taxon”.\nLine 79: “…in the CONTEXT OF ocean warming…”. Throughout the MS the tendency to\ncompound nouns should be avoided; alternatively, please hyphenate correctly to avoid\nambiguity.\nRemainder of the MS: I will refrain from further language editing and leave this to the\nEditorial staff to deal with. In my opinion, many grammatical errors and complications\nremain or have been introduced in revision.\nLines 226–227: Contrary to the wording, the confidence intervals seem to be appropriately\ndisplayed?\nLines 279–299: These descriptions to me don’t really match what the pie charts are\nshowing. I guess it depends on what is meant by “equal” and “dominate”. For example, in\npanel a, the only group with what looks like “equal(ish) intensities” to me is crustacea.\nFigure 4: It’s not clear to me what the crossed arrows in the middle of the plot are meant to\nillustrate. Please either explain in the caption or remove.\nLine 371–373: Given information in Figure 4, dominance is ~60:40. Is that really “highly\ndominant”? This ratio is a lot lower than I would have expected.\nLine 403: Figure 5b doesn’t seem to show changes in abundance?\nOverall, I believe that the work in this manuscript is robust, but the writing remains an issue,\nmaking it difficult to read/understand in places.", "reviewer_2": "The manuscript is vastly improved and I think the new analysis is much more statistically\nrobust. Thank you for including Figure 4, which addresses many of my questions (although\nnote that the figure legend and figure don't seem to match--the figure has arrows that\naren't explained in the legend, and the legend mentions extreme values that aren't labeled\non the figure.) I have no additional comments."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThe authors have done a good job in addressing my initial queries. I have only a few\nminor comments, listed below.\nTitle: “Cross-basin” singular, like, “cross-taxon”.\nAR: done\nLine 79: “…in the CONTEXT OF ocean warming…”. Throughout the MS the\ntendency to compound nouns should be avoided; alternatively, please hyphenate\ncorrectly to avoid ambiguity.\nAR: done\nRemainder of the MS: I will refrain from further language editing and leave this to\nthe Editorial staff to deal with. In my opinion, many grammatical errors and\ncomplications remain or have been introduced in revision.\nAR: I leave Editorial editing to deal with.\nLines 226–227: Contrary to the wording, the confidence intervals seem to be\nappropriately displayed?\nAR: Corrected. Sentence removed.\nLines 279–299: These descriptions to me don’t really match what the pie charts\nare showing. I guess it depends on what is meant by “equal” and “dominate”. For\nexample, in panel a, the only group with what looks like “equal(ish) intensities” to\nme is crustacea.\nAR: Reviewer is correct. Descriptions refer to the model analysis on the relative\nimportance of tropicalization with respect to deborealization (shown in Table 2)\nwhich is now indicated instead of Figure 3 pie charts.\nFigure 4: It’s not clear to me what the crossed arrows in the middle of the plot are\nmeant to illustrate. Please either explain in the caption or remove.\nAR: For all species in all case studies, arrows represent the mean value of Species’\nthermal bias with respect to Species’ abundance change for each underlying\nprocess (i.e. tropicalization, borealization, detropicalization, deborealization). Now\nincluded in Figure 4 caption.\nLine 371–373: Given information in Figure 4, dominance is ~60:40. Is that really\n“highly dominant”? This ratio is a lot lower than I would have expected.\nAR: The analysis at per-site basis (Fig 2d) and at per-species basis (Figure 4) agree\nin that tropicalization and deborealization dominates detropicalization and\nborealization, but proportions are not the same. The underlying ecological\nprocesses explaining the increase in CTI can generally be attributed to the\nprevalence of tropicalization or deborealization in most of sites (76.9% of sites),\nwhilst detropicalization and borealization dominated at fewer sites (23.1%) (Figure\n2d). We removed “highly”.\nLine 403: Figure 5b doesn’t seem to show changes in abundance?\nAR: Figure 5b show Species’ abundance increases (i.e.,\ntropicalization+borealization) relative to species’ abundance decreases (i.e.,\ndetropicalization+deborealization). This is now included in the text.\nOverall, I believe that the work in this manuscript is robust, but the writing remains\nan issue, making it difficult to read/understand in places.\nReviewer #2 (Remarks to the Author):\nThe manuscript is vastly improved and I think the new analysis is much more\nstatistically robust. Thank you for including Figure 4, which addresses many of my\nquestions (although note that the figure legend and figure don't seem to match--\nthe figure has arrows that aren't explained in the legend, and the legend mentions\nextreme values that aren't labeled on the figure.) I have no additional comments.\nAR: For all species in all case studies, arrows represent the mean value of Species’\nthermal bias with respect to Species’ abundance change for each underlying\nprocess (i.e. tropicalization, borealization, detropicalization, deborealization). Now\nincluded in Figure 4 caption.\nExtreme values are now labelled with species names in new Figure 4."}]

Decision letter