New climatic targets against global warming: will the maximum 2 °C temperature rise affect estuarine benthic communities?

1 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 www.nature.com/scientificreports New climatic targets against global warming: will the maximum 2 °C temperature rise affect estuarine benthic communities? Daniel Crespo 1, Tiago Fernandes Grilo2, Joana Baptista1, João Pedro Coelho 3,4, Ana Isabel Lillebø5, Fernanda Cássio6,7, Isabel Fernandes 6,7, Cláudia Pascoal6,7, Miguel Ângelo Pardal1 & Marina Dolbeth3,5 The Paris Agreement signed by 195 countries in 2015 sets out a global action plan to avoid dangerous climate change by limiting global warming to remain below 2 °C. Under that premise, in situ experiments were run to test the effects of 2 °C temperature increase on the benthic communities in a seagrass bed and adjacent bare sediment, from a temperate European estuary. Temperature was artificially increased in situ and diversity and ecosystem functioning components measured after 10 and 30 days. Despite some warmness effects on the analysed components, significant impacts were not verified on macro and microfauna structure, bioturbation or in the fluxes of nutrients. The effect of site/ habitat seemed more important than the effects of the warmness, with the seagrass habitat providing more homogenous results and being less impacted by warmness than the adjacent bare sediment. The results reinforce that most ecological responses to global changes are context dependent and that ecosystem stability depends not only on biological diversity but also on the availability of different habitats and niches, highlighting the role of coastal wetlands. In the context of the Paris Agreement it seems that estuarine benthic ecosystems will be able to cope if global warming remains below 2 °C. Despite specific adaptations to highly dynamic habitats such as estuaries, estuarine organisms can only cope with environmental variation within their specific range of tolerance. As such, species richness is generally low in estuaries, as few taxa evolved towards the required broad tolerance1, 2. Yet, estuaries are highly productive areas3, 4, due to continuous organic matter input from river basins, nutrient dynamics through decomposition and organic matter mineralization and habitat heterogeneity4, 5. It has been accepted that biodiversity with all its components, from the microscopic to the landscape level, is fundamental for the ecosystem functioning6–8. For instance, biological and functional diversity are essential for a complementarity effect and functional redundancy in a system, i.e. different species with similar functions, which usually delivers a better functional performance towards environmental change9–11. Nevertheless, species diversity does not always cope with rapid ecosystem modification under external pressures12, 13, which are becoming more frequent in face of global ecological changes. In the short term, processes and functions can also depend on behavioural shifts14. Therefore, even within rich biological systems, ecosystem functions such as primary and secondary production, decomposition and nutrient cycling can be affected even without drastic changes in the 1Centre for Functional Ecology - CFE, Department of Life Sciences, University of Coimbra, Calçada Martim de Freitas, 3000-456, Coimbra, Portugal. 2MARE - Marine and Environmental Sciences Centre, Laboratório Marítimo da Guia - Faculdade de Ciências da Universidade de Lisboa, Av. Nossa Senhora do Cabo, 939, 2750-374, Cascais, Portugal. 3CIIMAR - Interdisciplinary Centre of Marine and Environmental Research, Novo Edifício do Terminal de Cruzeiros do Porto de Leixões, Avenida General Norton de Matos s/n, 4450-208, Matosinhos, Portugal. 4Department of Chemistry & CESAM, University of Aveiro, Campus Universitário de Santiago, 3810-193, Aveiro, Portugal. 5Department of Biology & CESAM, University of Aveiro, Campus Universitário de Santiago, 3810-193, Aveiro, Portugal. 6Centre of Molecular and Environmental Biology (CBMA), Department of Biology, University of Minho, Campus de Gualtar, 4710-057, Braga, Portugal. 7Institute of Science and Innovation for Bio-sustainability (IB-S), University of Minho, Campus de Gualtar, 4710-057, Braga, Portugal. Correspondence and requests for materials should be addressed to D.C. (email: [email protected]) Received: 14 December 2016 Accepted: 12 May 2017 Published: xx xx xxxx OPEN www.nature.com/scientificreports/ 2 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 community demographic figures. Multilevel faunal interactions (microfauna, meiofauna and macrofauna) represent an additional challenge for tracking flows of energy/biomass and nutrients. Ex situ mesocosm experiments have proved to be efficient in measuring ecosystem functions and assessing the effect of changing diversity in those functions14–16 because control parameters can be easily tuned and models are straightforward to assess, with direct cause-effect responses16–19, although not without limitations. The interactive effects of multiple stressors on ecosystem functioning are difficult to predict, due to synergistic/antagonistic variation of complex natural systems16 or the effects of multifunctionality9. An experimental setup in field conditions is one way to account for more realistic responses and to understand how natural context may control biological responses (e.g. refs 19 and 20), by manipulating control/test variables, but keeping all the remaining natural variability (e.g. ref. 21). Nevertheless, it poses additional challenges in the experimental design and for the stabilization of control parameters17, 22. Global changes induced by anthropogenic impact have increased severely during the 20th century23, 24 and their outcome for biological systems are still unpredictable. However, in spite of the associated uncertainty, there is a general consensus regarding environmental and socio-economic implications derived from such changes25, 26. Several international agreements on greenhouse gases and other climate policies were assumed by major stakeholders in recent decade27, 28. However, an ambitious common goal regarding global warming was only assumed during the recent United Nations Conference on Climate Change, held in Paris, France, in 2015: to keep temperature rise below 2 °C, preferably 1.5 °C, above pre-industrial levels29. Therefore, it is highly relevant to measure responses of different functional groups within this new international framework. In fact, the Paris Conference on Climate Change appeals to the increase of scientific knowledge on the effects of climate, in order to support decision-making (article 7, paragraph 7 ref. 29). As such, the proposed research aims to contribute with science-based knowledge, focusing on the effect of a 2 °C temperature rise on benthic estuarine communities. For this, we performed a 30-day in situ experiment with benthic intertidal communities from a southern European temperate estuary, where a temperature increase was induced artificially. Data on micro and macrofauna diversity and important ecological processes (bioturbation) and functions (nutrient balance) were measured in order to contribute for a better knowledge about the effects of a mild temperature rise in the functioning of estuarine ecosystems. We expect changes in the benthic communities, namely decreased levels of diversity for both macrofauna and microfauna and changes in the dominant species under the effect of warmness, taking into account results from previous records on the effects of temperature increase in macrobenthic communities (e.g. refs 3 and 13). As for the ecosystem processes and functions, we also assume changes owing the differential communities expected under the warmness scenario. However, we cannot anticipate the change trajectory. Results Efficiency of the temperature increase. The warmness treatments were able to produce a differential temperature between the open boxes – “no warmness effect”, and the closed boxes – “warmness effect” (Supplementary Fig.S1). The temperature difference measured in the sediment at the daily peak was 2.35 °C for the sandflat and 1.52 °C for the Zostera bed (Supplementary Fig.S1), whose values were within the ones considered acceptable under the framework of the Paris United Nations Conference on Climate Change. Also, values for the continuous record (which includes high tides and night periods) denoted an average temperature rise on both sites (0.36 °C in the sandflat and 0.21 °C in the Zostera bed, Supplementary Fig.S1). Effect of the warmness treatment, habitat and time. Macrobenthic communities. For the macrobenthic communities’ species richness, we found two significant interactions with site, in combination with time and in combination with the temperature treatment effect (Table1). At T10 species richness was higher in the sandflat than in the Zostera, particularly for the “warmness” treatment (Fig.1a, Table1). However, at T30 the values were very similar between sites and treatments (Fig.1a). The Shannon-Wiener index (H’) for macrobenthic density was higher in the sandflat for both times (Fig.1b), with a significant interaction detected between site and time (Table1). For the H’ estimated with biomass, only site became significant (Table1), with the Zostera showing higher H’ values, except in the “no warmness” treatment, in T10 (Fig.1c). Regarding the macrobenthic community data, significant interactions between site and treatment and site and time were found (Table1). However, differences were clearer per site, with Peringia ulvae and Hediste diversicolor associated with the seagrass, while small opportunist polychaetes and Cerastoderma edule were associated with the sandflat (Fig.2a). Despite the significant interaction found between site and time, differences were significant for both time periods and sites, while for the treatment effect statistical differences were only found between control and the warmness treatment for the sandflat alone (Table1). For biomass data, we had similar results as the ones for density (two significant interactions, Table1). Again, differences were very clear per site (Fig.2b) and there were significant differences between T10 and T30 for the sandflat samples and when comparing sites within each time period (Table1). Regarding the effect of temperature rise on macrobenthic fauna, results were clearer than with the density (Fig.2b). However, PERMANOVA provided the same results: no statistical differences were found between control and the “no warmness effect” treatment, used here as a control for the box device, but there were significant differences between control vs “warmness effect” (Table1, Fig.2b). These differences seemed more evident for the sandflat (Fig.2b), despite being significant for both sites (Table1). Microbial communities. Microfauna diversity and composition was only assessed in T10. A similar trend variation was observed for the fungal richness and Shannon-Wiener index, but with clearer differences among areas and treatments for the richness (Fig.3a,b). Both indices were significantly higher in the sandflat area than in Zostera, with the highest values found on the control treatment associated with sandflat (Table2, Fig.3). Differences were also evident with regard to the warmness treatment effect (Table2), particularly for the richness www.nature.com/scientificreports/ 3 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 (Fig.3a,b). Richness values were significantly lower in the “warmness” compared to the “no warmness” and control treatments (Tukey’s, p < 0.05), while Shannon-Wiener values were significantly lower in the “warmness” compared to the control (Tukey’s, p = 0.0182). The trend variation in the diversity indices of the ciliate communities was also similar, but again, differences were clearer with richness (Fig.3c,d). For both indices, the interaction between site and temperature was significant (Table2). For the sandflat, both indices were highest in the control treatment and lowest in the “warmness” treatment (Fig.3c,d). However, this trend was the opposite for the Zostera, which had the lowest values in the control treatment compared to the others (Fig.3e,f). The same pattern was observed for the Shannon index (2-way ANOVA, p = 0.003; Table2; Fig.3d), with temperature treatment decreasing ciliates communities’ evenness in sandflat, but increasing ciliates communities’ evenness in Zostera. Bacterial richness and Shannon-Wiener index varied within similar values for both sites and treatments (Fig.3e,f), with no significant differences detected among factors (Table2). Bacterial communities contained on average 30 species (Fig.3e). Regarding the composition and structure of the fungal, ciliate and bacterial communities, significant interactions were found between site and temperature for all communities (Table3). However, the pairwise comparisons were not significant for any of the terms/factors pairs (Table3). As such, we have considered the significance among main terms, which were always significant (Table3). Despite this result, the variability explained in the PCO plots for the fungal and ciliate communities was relatively low (<25.1%) and the discrimination among factors was not that clear (Fig.4a,b). In both these communities, significant differences were found per site and among the control with the other treatments (Table3). For the bacterial communities, statistical differences were found per site and among all temperature treatments (Table3). The variability explained in the PCO plot was higher than for the other community (about 60%), with the control samples clearly separated from the other treatments (Fig.4c). For the Zostera, differences among treatments were also clearer (Fig.4c). Biological process – particle reworking. Evidence that bioturbation was affected by the different levels of the tested factors was found. The minimal adequate models were GLS regressions (Table4), with a variance structure of different spreads for each independent fixed factor (full structure of the models, correction of the random part, pairwise comparison matrix for significant differences among levels and graphic representation of model predictions in Supplementary Material - Model 1 to Model 4; mean values and standard errors for each measurement in Supplementary TableS1). Surface Boundary Roughness (SBR), f-SPILmean and f-SPILmax were significantly affected by the interaction of all factors (site × temperature × time), while f-SPILmedian only varied with site (Table4, models structure for the four bioturbation components described in Supplementary Material). SBR values ranged from 0.358 to 2.470 cm. There was always a decrease in SBR when comparing the effect of the warmness treatment, a result that was consistent for both sites and T10 (Fig.5a, Supplementary TableS1). For T30, SBR increased with the “warmness effect” (Fig.5a, Supplementary TableS1). SBR increased with time for both temperature treatments in the Zostera site. However, in the sandflat SBR decreased with time for the “no warmness effect” treatment and increased with time for the “warmness effect” treatment, in agreement with the significant interaction detected. Despite the observed tendencies, when comparing relevant significant differences (i.e. comparisons of treatments that share the same levels of two different factors), those were only found between the “no warmness” vs “warmness” effects in the sandflat at T10 (p = 0.0072), between T10 and T30 in the sandflat under the “no warmness effect” (p = 0.0009) and between sites at T30 under the “no warmness effect” (p = 0.039) (Supplementary TableS2). f-SPILmean values ranged from 0.689 to 6.476 cm and reflected the complex 3-way interaction. At T10, f-SPILmean values were lower with “no warmness effect” treatment and increased with the “warmness effect” for the Dependent variable Significant terms d.f. Pseudo-F p-perm Terms/levels of factor p-perm Species Richness site × temperature 2 4.2065 0.021 warmness [sandflat vs Zostera] 0.021 site × time 1 22.011 0.002 T10 [sandflat vs Zostera] 0.001 Sandflat, Zostera: [T10 vs T30] <0.015 Shannon (density) site × time 1 12.761 0.002 T10, T30: [sandflat vs Zostera]<0.004 Sandflat, Zostera: [T10 vs T30] <0.05 Shannon (biomass) site 1 4.12 0.05 sandflat vs Zostera 0.05 Benthic community density site × temperature 2 1.8469 0.037 Sandflat: [Control vs warmness] 0.056 All treatments: [sandflat vs Zostera] <0.05 site × time 1 2.8878 0.01 For T10, T30: Zostera vs sandflat 0.001 Sandflat, Zostera: T10 vs T30 <0.005 Benthic community biomass site × temperature 2 2.171 0.032 Sandflat: [Control vs warmness] 0.017 All treatments: [sandflat vs Zostera] <0.05 site × time 1 3.2585 0.019 For T10, T30: Zostera vs sandflat <0.002 Sandflat, Zostera: [T10 vs T30] <0.04 Table 1. Summary of significant terms from the 3-way PERMANOVA analyses, with macrobenthic community and respective diversity indices as dependent variables, and site, treatment and time as explanatory variables, with indication of the significant pairwise comparisons. www.nature.com/scientificreports/ 4 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 sandflat (p = 0.0339, Supplementary TableS3), while the inverse tendency was found for the Zostera (Fig.5b, Supplementary TableS1). At T30, in the Zostera bed, f-SPILmean increased for the “warmness effect” compared to the “no warmness” one (Fig.5b, Supplementary TableS1). In the sandflat, f-SPILmean values decreased with time for “warmness effect” (p = 0.0199, Supplementary TableS3), while for the Zostera bed f-SPILmean decreased with time for the “no warmness effect” treatment and increased with the “warmness effect”. Generally, Zostera bed had higher values of f-SPILmean than the sandflat at T30, for both “no warmness effect” (p = 0.0175) and “warmness effect” (p = 0.0072, Supplementary TableS3). f-SPILmedian only varied with site and values ranged from 0 to 6.949 cm (Fig.5c). These were considerably higher for the Zostera bed than for the sandflat (p = 0.004, Supplementary TableS1). Values of f-SPILmax ranged from 8.561 to 12.490 cm and again reflected the complex 3-way interaction. Overall, there was a tendency for f-SPILmax to increase with time, both in the sandflat and in the Zostera bed (Fig.5d). In the sandflat the increase was higher in the “no warmness effect” treatment than in the “warmness effect” treatment that showed a circumstantial decrease (Fig.5d, Supplementary TableS1). In the Zostera bed the larger increase over time was found in the “warmness effect” treatment (relevant significant differences only detected for this scenario, p = 0.0027, Supplementary TableS4), with a small decrease in the “no warmness effect” treatment (Fig.5d, Supplementary TableS1). When comparing temperature treatments, the “warmness effect” showed higher or similar values to the “no warmness effect”, except at T10, for the Zostera bed (Fig.5d). Nevertheless, relevant significant differences were only found between the sandflat and Zostera bed, at T10, under the “warmness effect” (p = 0.0001, Supplementary TableS4). Sediment OM, C and N pools and intertidal water nutrient concentrations. The average organic matter content in the sediments was 1.96 ± 0.02% (mean ± SD) in sandflat and 2.58 ± 0.02% in Zostera, respectively, with significant differences between sites (matched-pairs t-test = −206.621, df = 1, p-value = 0.003). Regarding C and N pools in the sediment at the beginning of the experiment, differences were observed between sites, with higher carbon and nitrogen concentrations observed in the Zostera meadows, as a result of the higher organic matter content (C [mean ± s.d]: 1.138 ± 0.25 mg L−1 in sandflat and 1.74 ± 0.36 mg L−1 in Zostera; N: 0.086 ± 0.027 mg L−1 in sandflat and 0.135 ± 0.027 mg L−1 in Zostera). However, these distinct sediment characteristics and distinct plant coverage were not reflected in the dissolved nutrient concentrations in the shallow low water pools. Overall, the Figure 1. Mean + s.d. (n = 3) species richness (a) and Shannon-Wiener index of the macrobenthic communities measured with density (b) and biomass (c) in the sandflat and Zostera areas, under control, no warmness and warmness treatment, for the two time periods. www.nature.com/scientificreports/ 5 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 range of nutrient concentrations at the sampling sites and time were within the described for the system30, taking into account tidal and circadian cycles. However, no distinct pattern was observed in the nitrogen and phosphorus concentrations of the shallow intertidal pools for all treatments, sites and times. Concentrations of dissolved inorganic nitrogen ranged from 0.005 to 0.58 mg L−1, with most values <0.25 mg L−1, while PO4-P concentrations ranged between 0.02 and 0.15 mg L−1, but most values were <0.05 mg L−1. Discussion The importance of in situ experiments in the context of global warming studies. The use of field experiments to assess the effects of biodiversity on ecosystem functioning could be complemented with laboratorial work31, 32 in order to fine tune the potential causal-effect relationship. Aquatic ecosystems are highly complex systems where organisms are challenged by other biological components, as well as by abiotic influences, with several levels of interactions. This is particularly true for estuaries, where environmental variables shift over several cycles33, 34. Therefore, the study of transitional systems in laboratory has challenges that are inherent to that complexity. Despite the high level of accuracy allowed by laboratorial work, in situ experiments provide a more realistic interpretation of natural processes that actually occur in the system. Nevertheless, in situ experiments face limitations concerning universality (e.g. refs 21 and 35): the outcomes can be site specific as context modulates biotic responses, either due to physiological events or behavioural alteration. Simultaneously, climate changes, including global warming, may have very specific local responses. Therefore, it is of great interest to compile data from a wide geographic range in order to better understand climate change consequences. The present work brings some insights on these possible consequences on an estuarine system, in particular regarding a temperature increase, that is within the range recently preconized by the Paris United Nations Conference on Climate Change −1.5 to 2 °C. Nevertheless, this temperature increase should be regarded as an ambitious goal, fairly more optimistic than most of global climatic models predictions26. The temperature increase effect on different habitats and in benthic communities. A more complete understanding of the consequences of warmness in benthic estuarine communities could be provided by looking upon the meiofauna, as they play a key role linking microand macrofauna. Nevertheless, in such a complex and detailed in situ approach, already involving a large multidisciplinary team and limited funding, we realized that all fauna communities could not be examined. Moreover, a previous study in the same location (the Mondego estuary)36, which followed the organic matter mineralization and nutrient dynamics (including denitrification) mediated by the three faunal groups found that microfauna was responsible for the majority of the ecosystem processes (67%), followed by macrofauna (18%) and to a lesser extent by meiofauna (7%). Consequentially, we addressed the microfauna, as it is known to be the most reactive faunal community to environmental change, and the macrofauna, which shows a slower response to external pressures. Figure 2. PCO ordination for macrobenthic community density and biomass for each area and treatment. The length and direction of each vector indicates the strength and sign of the relationship between species abundance or biomass and PCO axes, based on a Spearman correlation (only vector with length >0.5 were represented). www.nature.com/scientificreports/ 6 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 In our study, temperature was measured on the top layers of the sediment, suggesting that the air inside the closed boxes (warmness effect) could be warmer than in the open boxes. This temperature differential could emulate plausible values found during heat wave events, despite the fact that the experimental procedures were conducted during fairly routine summer period. The temperature increase in the sediment was within the range of the acceptable temperature increase agreed at the Paris Conference29. Although the temperature difference between the treatments was not constant, due to night temperature and high tides periods, we found a general 1.5 to 2 °C increase for the warmness effect treatment during the daily peak measurements. Even so, the Zostera bed seems to be acting as a temperature buffer, since temperature range was generally lower than in the bare sandflat. Three different mechanisms could explain this result: 1) a shading effect that the Zostera’ leaves produce on the underlying sediment, which could reduce the amount of sunlight that reaches the sediment and thus warming; 2) a different colour between the Zostera bed and the sandflat areas, since the sandflat has darker shade and, therefore, better efficiency in light absorption, becoming warmer than the Zostera bed; 3) a differential thermal buffer due to distinct water retention capacities provided by the seagrass leaves and rhizomes. Regarding the benthic communities, a differential temperature regime could limit or stimulate some natural processes and cause shifts on the community structure and functioning associated with those communities3, 37, 38. However, our 1.5 to 2 °C temperature increase seemed less influential than initially expected. In fact, the differences between sites (Zostera bed and sandflat) were more evident than the differences between temperature treatments, for both the macrofauna and microbiological communities. Against expectations, our macrofauna diversity results were inconsistent with previous studies dealing with the effects of the temperature rise, which have reported increased mortalities37, 38 and decreased diversity13 under such scenario. Our temperature increase Figure 3. Mean + s.d. (n = 3) richness of fungal (a), ciliates (c) and bacterial (e) communities and ShannonWiener index of fungal (b), ciliates (d) and bacterial (f) communities on sediment samples from sandflat or Zostera under control, no warmness and warmness treatment. www.nature.com/scientificreports/ 7 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 was lower than the considered in those studies, and diversity was sometimes higher in the warmness effect treatment than in the control and non-warmness ones. Similarly, the period during which the communities were exposed to the heating effect (10 days or 30 days) did not have a consistent/clear effect, except for a higher variability among treatments in the sandflat, particularly at T10. In other words, there was less variation in diversity indices among treatments for the Zostera bed. Still, when examining the composition and structure of the macrofauna communities, some differences emerged due to the temperature increase (warmness effect), particularly for the biomass levels in the sandflat, the area where the temperature increase was higher. Again, these results suggest that the Zostera bed may mitigate the temperature increase effect on the macrobenthic community, and that time may have allowed communities to adjust to the disturbance (e.g. under an acclimation process39). Dependent variable Effect S.S. d.f. M.S. F-value p-value Fungal richness temperature 94.33 2 47.17 7.13 0.009087 site 32.00 1 32.00 4.84 0.048130 temperature × site 14.33 2 7.17 1.08 0.369168 error 79.33 12 6.61 Fungi ShannonWiener index temperature 0.32 2 0.16 5.63 0.018869 site 0.153 1 0.15 5.30 0.040004 temperature × site 0.06 2 0.03 1.11 0.360000 error 0.35 12 0.03 Ciliate richness temperature 12.33 2 6.17 1.82 0.204069 site 6.72 1 6.72 1.98 0.184389 temperature × site 106.78 2 53.39 15.75 0.000440 error 40.67 12 3.39 Ciliate ShannonWiener index temperature 0.14 2 0.07 1.70 0.225634 site 0.01 1 0.01 0.31 0.588654 temperature × site 0.80 2 0.40 9.50 0.003361 error 0.50 12 0.04 Bacterial richness temperature 4.25 × 1012 2 2.13 × 1012 1.01 0.393082 site 8.71 × 1011 1 8.71 × 1011 0.41 0.532230 temperature × site 1.02 × 1012 2 5.12 × 1011 0.24 0.788013 error 2.53 × 1013 12 2.11 × 1012 Bacteria ShannonWiener index temperature 2.11E-03 1 2.11E-03 0.46085a0.558b site 6.69E-03 2 3.34E-03 0.73183a0.524b temperature × site 1.04E-02 2 5.19E-03 1.1361a0.367b Residual 5.48E-02 12 4.57E-03 Table 2. Results from the 2-way ANOVA for the effect of site and temperature on fungal, ciliate and bacterial communities’ richness and Shannon-Wiener indices. For the bacterial Shannon-index a 2-way PERMANOVA was applied, as data did not pass the assumptions of parametric tests; apseudo-F, bP(perm). Dependent variable Significant terms d.f. Pseudo-F p-perm Terms/levels of factor p-perm Fungi site × temperature 2 1.9417 0.02 no significant terms ns site 1 2.1485 0.02 sandflat vs Zostera 0.02 temperature 2 1.7761 0.01 control vs no warmness 0.04 control vs warmness 0.03 Ciliate site × temperature 2 1.7003 0.03 no significant terms ns site 1 2.956 0.002 sandflat vs Zostera 0.004 temperature 2 1.6075 0.03 control vs no warmness 0.03 control vs warmness 0.03 Bacteria site × temperature 2 5.4363 0.001 no significant terms ns site 1 2.8995 0.003 sandflat vs Zostera 0.004 temperature 2 9.0145 0.001 control vs no warmness 0.003 control vs warmness 0.002 no warmness vs warmness 0.002 Table 3. Summary of significant terms from the 2-way PERMANOVA analyses with fungal, ciliate or bacterial communities as dependent variables, and site and temperature as explanatory variables, with indication of the significant pairwise comparisons. www.nature.com/scientificreports/ 8 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 The microfauna communities, analyzed only at T10, were slightly affected by the warmness treatments, similarly to the results observed for the macrofauna. In general, fungal and ciliate richness was higher in sandflat than in the Zostera. Under the temperature increase scenario, fungal communities declined in both areas while the ciliate ones decreased only in the sandflat. Significant differences were also detected for the structure of those communities regarding site and temperature. However, these differences were not clear in PCO plots, and whose variability explained was considerably low (<25.1%). On the other hand, bacterial diversity did not respond clearly to the treatments or site. However, the differences on the bacterial community structure were more expressive than fungal and ciliate communities: there were also significant differences with site and temperature and samples clustered taking into account these differences, particularly regarding the control and the other treatments. Overall, it seems that the warmness treatment had a larger impact in the microbial diversity than in Figure 4. PCO ordination for fungal (a), ciliates (b) and bacterial (c) communities for each site and temperature treatment. Dependent variable Significant terms d.f. L-ratio p-value SBR site × temperature × time 7 16.813 0.0186 f-SPILmean site × temperature × time 7 20.003 0.0056 f-SPILmedian site 1 8.153 0.0043 f-SPILmax site × temperature × time 7 20.316 0.0049 Table 4. Summary of significant terms from the GLS for particle reworking measurements, with bioturbation components as dependent variables and site, treatment and time as explanatory variables. www.nature.com/scientificreports/ 9 Scientific RepoRts | 7: 3918 | DOI:10.1038/s41598-017-04309-0 the macrobenthic diversity. This was expected, since microfauna communities respond much faster to external factors than the macrofauna ones40, 41, especially due to faster generation times42. There is also the possibility that the macrofauna under the warmness stimulus could interfere with microbial communities by competitive interaction43. Despite the observed differences, in light of our initial question, how a 1.5 to 2 °C temperature increase would affect benthic communities, the variation trends obtained were not very expressive. Also, daily environmental fluctuations could conceal the influence of the warmness treatment on microbial communities, by diluting the differential temperature between treatments. Consequences on ecosystem processes and functions. Our next question was whether the changes in the benthic communities due to the warming would reflect in changes on the processes and functions that are usually associated with those communities. However, since differences in the benthic communities with regard to the temperature rise effect were not particularly consistent, we expect the same pattern to occur with these processes and functions. For almost all components of bioturbation, there were significant interactions between all factors, suggesting that temperature had an effect on bioturbation that depended on the site and sampling time. We used smaller sized corers when compared with other studies (usually width >10 cm: e.g. refs 19, 22 and 44). However, we were able to measure the general tendency in fauna movements because we evaluated the luminophores’ statistical distribution along an unidimensional direction and not taking into account particular features of biogenic structures19, 45, 46. Also, bioturbation is mediated by organisms within a large range of sizes47, 48, which implies that even if larger animals are kept outside the corers, particle movement was still detected. Surface boundary roughness (SBR) increased with temperature only for the sandflat and for the longer time (T30), which could be related with a temperature gradient decreasing from the sediment top downwards that stimulates some activity in the top layer of the sediment. The temperature increase on the top sediment layers was expected to trigger the downwards movement of animals. In a scenario of warmness, the same would happen, even if deeper sediment layers were able to buffer the overlaying temperature increase. The warmness treatment had no significant differences in f-SPILmed (Median Luminophore depth), which reflects the short-term depth of mixing20 and could be a good proxy for the avoidance/downwards movement behaviour. Differential values of f-SPILmed between treatments may be a consequence of differential heat avoidance behaviour by the infaunal organisms, although this was not the case in our experiment. Nevertheless, for the other conditions, the variation was the inverse and not conclusive with regard to the temperature rise effects. For the other components of bioturbation, the variation tendency was also not clear or conclusive with regard to the temperature increase. We expected an overall bioturbation modification with temperature, as species respond differently to thermal pressures14, 44. However, the interaction of several factors (including those not assessed in the experiment, e.g. salinity, water flow, turbidity) resulted in indistinct tendencies, with large variation in the sediment mixing measurements within each Figure 5. The significant effects of warmness treatment, site and time on surface boundary roughness (SBR) (a), f-SPILmean (b), f-SPILmedian (c) and f-SPILmean (d) (cm, mean ± s.e.). For clarity, jitter has been applied to the x = argument of the plot function to avoid over-plotting.