Philippa Kaur

Funicular bodies: structure, function, and developmentThe ultrastructural and functional complexity of funicular bodies in bryozoans25,46,47 (our data), as well as the reduction of the genome in their symbiotic bacteria (e.g., Bugula neritina48), point at a long-term co-evolution between these organisms, also suggesting that infected bryozoan colonies spend a significant part of their energy budget supporting numerous bacteria inside FBs. All known prokaryote symbionts are apparently non-pathogenic for bryozoans49,62. Instead, the overall evidence indicates a mutualistic relationship. Symbiotic bacteria are known to produce toxic substances (bryostatins in Bugula63; bryoanthrathiophene in Watersipora64) that protect bryozoan larvae from predators42,44,45. Another potential role of bacterial secondary metabolites is the chemical defense of early developmental stages during larval settlement and metamorphosis40,65. Such chemical defence may also prevent epibiotic overgrowth of the bryozoan colonies by other bacteria and algae (reviewed in66,67). Similar functions can be assumed for symbiotic prokaryotes of Dendrobeania fruticosa, although experimental evidence is still required.FBs with bacterial symbionts in colonies of D. fruticosa show signs of high specialization. The walls of the funicular body completely isolate its internal cavity from the surrounding zooidal cavity: the cells of the outer layer overlap each other, whereas cells of the inner layer have tight Z-shaped contacts. Such isolation probably creates a specific environment inside the symbiont-containing space that helps maintain a growing population of bacterial symbionts. A massive protein-synthesizing apparatus was observed in the cells of the inner layer. In addition, the ‘pocket’-like structure of the internal cavity of FB and the abundance of cytoplasmic processes, some of which protruded deep into this cavity, contribute to increasing both the general inner surface area of FB and the contact area between its inner cells and bacteria. From their side, bacteria provide this contact by numerous pili. Numerous pits and microvesicles associated with the apical membrane of the inner cells imply an active exchange of substances between the host and its symbionts (e.g., nutrition provided to growing and multiplying bacteria and absorption of the potential wastes they produce).The ultrastructure and functional morphology of FBs in Aquiloniella scabra, the only species in which the ultrastructure of these organs was described, are very similar to those in D. fruticosa. Its FBs also consist of two cell types, but the inner cell layer is formed by one or a few cells, whereas external cells form a multilayered envelope25. In both species, the main function of FBs was considered to be the chambers/organs for symbiont incubation and nourishment.FBs with symbiotic bacteria have been found in several bryozoan species from four different families20,25,47,49,51 (our data). Remarkably, all these families (Bugulidae, Beaniidae, Candidae, Epistomiidae) belong to the clade Buguloidea. Encapsulated aggregations of bacteria within the zooidal cavity were also described in two species of Watersipora from the phylogenetically distant family Watersiporidae (Smittinoidea)62,68. In that case, prokaryotes (“mollicutes”) were also enveloped by flattened bryozoan cells, although the current scarcity of data makes it difficult to compare these cell aggregations with funicular bodies.Among Buguloidea, FBs show the same basic structure, although species described by Lutaud49, Dyrynda and King51, and Mathew with co-authors46 based on light microscopy require re-examination by TEM. Moreover, FBs are in all instances associated with the funicular system of the colony. This similarity could indicate the single origin of the bacterial symbiosis and FBs within Buguloidea. Still, it remains unknown (although rather probable) whether that system transports nutrients from feeding polypides to FBs because no communication between the lacunae of the funicular cords and the FB internal cavity was detected. Another, more likely option, however, is the independent acquisition of bacteria (see also below) and a similar ‘reaction’ of host tissues to invaders, resulting in the formation of bacterial organs (FB) with a similar bauplan. The presence of incapsulated bacteria in non-related Watersipora supports the second interpretation. Answering this question will require both ultrastructural studies and molecular identification of bacteria.As for the initial source of somatic cells and development of FBs, two variants have been suggested. Describing the development of FBs in Bugulina turbinata, Lutaud49 stated that the epithelial cells of the cystid wall were transformed into the inner cell layer of the FB, whereas peritoneal cells formed its external lining. By contrast, while studying Aquiloniella scabra, Karagodina with co-authors25 suggested that bacteria are engulfed by one of the funicular cells, which becomes a ‘bacteriocyte’ that is later enveloped by neighboring funicular cells. In the latter case, FBs are considered to be modified expanded parts of the funicular system. This is consistent with experiments by Sharp and co-authors40, who detected groups of labeled bacteria in the funicular cords of Bugula neritina, and also multiple bacteria developing inside enlarged funicular cords (in fact, very large FBs) in the related Paralicornia sinuosa20.Our TEM study of FBs in D. fruticosa showed that they are not swollen parts of the funicular cords, as was stated by Vishnyakov with co-authors for B. neritina20. It is more likely that the funicular cords and processes of their cells contact the external cell layer of FBs. According to the third scenario that we present here, the inner cell layer of FBs in D. fruticosa, as well as in other studied bugulids, most likely originates from the coelomocytes which accumulate bacteria via phagocytosis. Such solitary cells were described inside the zooidal cavity in B. neritina20. Instead of being digested, engulfed bacteria could trigger the coelomocyte divisions resulting in the formation of the inner cell layer. In contrast, in A. scabra, the coelomocyte can remain single or undergo only a few divisions. We propose that the external cell layer of FB in D. fruticosa originates from the funicular cells because these cell types are ultrastructurally similar. Finally, it is also possible that the exact process and sources of FB formation differ in different species.Multiple and lobed FBs found in two zooids of D. fruticosa could indicate a potential mode of their multiplication. The case of P. sinuosa requires additional study, but currently we believe that its bacteria-bearing ‘funicular cords’ are very large, elongated FBs, as well (see20).Symbiont circulation in the bryozoan life cycleThe taxonomic diversity of bryozoan hosts and their symbiotic bacteria—supported by a variety of sites in the bryozoan zooids and larvae where symbionts have been found—unambiguously point to multiple independent origins of symbiotic associations between bacteria and cheilostome Bryozoa20,25.Rod-shaped bacteria are the most common symbionts in the superfamily Buguloidea. Although superficially similar, these bacteria strongly differ in their maximum size, suggesting the presence of different procaryote species. Thus, bacteria detected in the larvae of Bugulina simplex and in FBs of Aquiloniella scabra can reach 10 μm in length25,69, while healthy symbionts (see below) inside coelomocytes and presumably peritoneal cells of Bugula neritina were only 2.5 μm long20. The maximum length of bacteria in FBs of D. fruticosa never exceeded 5 μm. Coccoid bacteria in the tentacles of B. neritina were 0.5–0.7 μm in diameter20. Else, oval or irregularly-shaped mycoplasma-like α-Proteobacteria were detected in the genus Watersipora62,68,70. By contrast, the symbionts identified in B. neritina and B. simplex belong to γ-Proteobacteria69,71,72.Apart from FBs, prokaryote symbionts were described extracellularly in colonies of different cheilostome species: in vestibular glands of autozooids, inside polymorphic zooids (avicularia)38,39,50, in tentacles and funicular cords20. They have also been found intracellularly: inside coelomocytes, epithelial and peritoneal cells of the body wall, and pharyngeal cells20,50. In addition, Woollacott and Zimmer73 described bacteria in the ‘channels’ of the funicular cords associated with brood chambers. However, the TEM image they published shows bacteria inside large vacuoles of the funicular cells—seemingly not in the lacunae between these cells, recalling the aforementioned idea of a ‘bacteriocyte’. Finally, bacteria have also been found in the pallial sinus of bryozoan larvae (62,68,74; reviewed in72). All these diverse data have led to two opposite views on the acquisition and circulation of symbionts in the bryozoan life cycle.Discovery of bacteria in both colonies and larvae of B. neritina was regarded as possible evidence of their transmission from larvae to adults74. This assumption was experimentally proven using both labeled bacteria and their metabolites (bryostatins) in the larvae and preancestrulae developing from them40. Moreover, the presence of bacteria within brood chambers (ovicells) in this species was considered as proof for the next step—the transition of symbionts from the colony to the incubated larvae, i.e., the vertical transfer of symbionts. Symbiotic bacteria populating larvae (and making them unpalatable for predators) are incorporated into the preancestrular tissues during larval metamorphosis, and then found inside zooidal buds in early colonies (a symbiont association with the host cells was not specified) and funicular cords of rhizoids in adult colonies40.The next step of the bacterial development could be the formation of FBs as a locus of symbiont reproduction. Mature FBs, full of bacteria, were considered to be the starting point for the transfer of prokaryotes (by an unknown mechanism) from the zooid to the brood cavity (via funicular cords associated with both FB and ovicells), and then to the incubated larvae47. Light microscopic data demonstrated: (1) the association of FBs with tube-like funicular cords47, and (2) the presence of groups of ‘bacterial bodies’ (small aggregations of bacteria) inside the ooecial vesicle (membranous-epithelial ‘plug’ that closes the entrance to the ovicell; mentioned in B. neritina47, and shown in images of the related Bugulina flabellata75,76,77). In addition, ultrastructural data proved the presence of bacteria inside funicular cords, more precisely—inside their funicular cells (see above), “extending to the ooecial vesicle” (73, p. 362). Elsewhere, Sharp and co-authors (40, p. 697) used fluorescence microscopy to demonstrate the presence of bacteria “within the ovicells”, and suggested that they are transported there across the colony via funicular cords that also house bacteria. Combined, all these data imply that bacteria move from FBs to the ooecial vesicle, accumulate there, and then somehow enter the brood cavity, which contains a larva, either through or in-between the epithelial cells and the cuticle of the ooecial vesicle. Findings of bacteria inside larvae and adult colonies of two species from the non-related genus Watersipora62,68 further strengthened the hypothesis of vertical transfer, which has subsequently been widely accepted by many authors20,30,47,78. Despite extensive TEM studies, no bacteria have been found inside the funicular cords in B. neritina (Vishnyakov & Ostrovsky, unpublished data), which contradicts the data of Sharp and coauthors40 obtained by fluorescent microscopy. Accordingly, it was suggested that coelomocytes carry symbionts to the ooecial vesicle instead20.Nonetheless, the hypothesis of vertical transfer faces serious objections based on life history, molecular and morphological data. In Dendrobeania fruticosa in the White Sea, for example, larval production occurs predominantly in autumn (mainly in the distal parts of branches, Fig. 1A), and no bacteria are present in colonies during this period. In addition, molecular population studies revealed that B. neritina is a complex of sibling species, both symbiotic and aposymbiotic, some of which live in sympatry, with the horizontal transfer between colonies being the most parsimonious explanation for the distribution of bacteria between siblings79. A study of the genome of symbiotic bacteria showed that they may be able to live outside the host48, which is consistent with the hypothesis of horizontal transfer (which is not the same as environmental transmission, see below).TEM data showed no communication between the FB cavity and lacunae of the funicular cords in the studied species, in particular in B. neritina (Vishnyakov & Ostrovsky, unpublished data), Aquiloniella scabra25 and D. fruticosa (this study). Coelomocytes (and presumed peritoneal cells) with bacteria embedded in their cytoplasm were indeed recorded inside the ooecial vesicle in B. neritina20. Nevertheless, extensive TEM studies of B. neritina ovicells at various stages of placental development (Vishnyakov and Ostrovsky, unpublished data) have not revealed bacteria between placental cells adjoining a developing embryo. The fact that these cells are provided with both tight and adherens junctions (e.g., B. neritina and Bicellariella ciliata73,80), and additionally are covered by a cuticle (albeit thin) raises the question of whether coelomocytes with bacteria and/or bacteria alone can move through the very thick hypertrophied placental epithelium. Interestingly, Miller with co-authors48 detected a gene encoding chitinase in the genome of the symbiont of B. neritina that could potentially be used for cuticle piercing.Another opportunity for the vertical transfer of symbionts is their transport via the supraneural coelomopore during oviposition (see20,81). In this case, free bacteria in the cavity of the maternal zooid could stick to the ovulated oocyte before its transfer into the brood cavity via the coelomopore. However, free bacteria were never recorded in the zooidal cavity. For the Watersipora species, an assumed variant of symbiont transmission is through a strand of mucus extending from the maternal zooid to the released larva and tethering it for a few minutes68.Environmental transmission, when bacteria are acquired from the surrounding seawater, is an alternative option for symbiont acquisition. It may potentially occur either via infection of brooded larvae inside the ovicell by bacteria entering the brood cavity from the external environment or via infection of larvae during the free-swimming period by bacteria from the water column. Published images by Sharp and co-authors40 showed the presence of both symbiotic bacteria and bryostatins both inside the ooecial vesicle and in the peripheral part of the brood cavity, close to the entrance of the ovicell. Although these authors stated that such close “locations of both the bacterial symbionts and the bryostatins demonstrate that the B. neritina–‘E. sertula’ association has a delivery system for both the symbionts and the bryostatins to embryos within the ovicell” (40, p. 699), we argue, based on the above-mentioned data, that this statement remains a probable yet unproven speculation. Our numerous unpublished TEM images indicate the presence of a large number of bacteria filling the brood space between the embryos and the ovicell wall in the brood chambers with and without developing larvae. These bacteria, attracted by some chemical signal(s), can enter the ovicells from outside and infect larvae. Until the transfer of bacteria through the wall of the ooecial vesicle or during oviposition is documented, environmental transmission remains the more probable method. Notably, recent studies on sponge microbiotas showed that the environmental transmission is widespread in this group of suspension feeders82.The hypothesis of external acquisition of bacteria by the bryozoan hosts leaves different infection pathways open. Some of these could potentially develop into vertical transfer. Beyond the infection of larvae, prokaryotes could enter feeding autozooids via the mouth (and further through the intestinal epithelium into the zooidal cavity), through the coelomopore—a presumed entrance for alien sperm76,81,83, or by direct infection of the tentacles (probably by penetration through the outer epithelium).We should stress that FBs were absent in zooidal buds and the youngest zooids with functional polypides in the growing branch tip of one colony of D. fruticosa collected in June. This suggests that bacteria are not transmitted from the older colony parts (and, thus, are not inherited from the founding larvae), but obtained from the external medium since older (and more proximal) zooids had FBs. This idea is supported by the lack of signs of transfer of bacteria between zooids via communication pores and their pore-cell complexes in Dendrobeania fruticosa and Aquiloniella scabra in our TEMs. In contrast, fluorescence microscopy showed symbionts in the non-feeding zooidal bud of the newly-formed small colony of B. neritina40. Bacteria were also present in the preancestrula formed during larval metamorphosis. It remains unknown whether they can move from the preancestrula to the bud along with coelomocytes, with growing funicular cords or both before the formation of transverse walls that isolate newly budded zooids. So, interzooidal transfer to budding sites is possible, and the youngest zooids with functional polypides in the growing tips of D. fruticosa could already receive bacteria too, but FBs were not yet developed. Thus the question of interzooidal/intracolonial transport of bacteria remains open.Two infection pathways—via larvae and by direct penetration through the tissues of the functional polypide, potentially exist in the same species. This is the case in B. neritina, which has morphologically different symbionts in FBs and in the tentacles20. Bacteria in the epithelial wall of the tentacle sheath and ooecial vesicle (see above) could potentially get there via both pathways, or enter the zooid via the coelomopore (third way), subsequently becoming entrapped by coelomocytes or cells of the cystid wall.Finally, the presence of bryozoan sibling species, some of which have symbionts while others do not, and the presence of symbiotic and aposymbiotic colonies within the same species79, suggests that bacteria can be lost and acquired anew at both short- and long-term time scales, as occurs in hermatypic corals and their symbiotic zooxanthellae84,85,86. In this light, it would be important to know whether FBs can develop anew in D. fruticosa after overwintering in the same zooids, or whether they appear only in newly budding zooids. What is the source of bacteria in overwintered colonies? Is it an external infection, or some ‘survivor’-cells (descendants of the bacterial pool from the larva that overwintered inside epithelial cells and/or coelomocytes), or both? This will require further study.Symbiont population dynamics in Dendrobeania fruticosa and its potential driversWhatever the route used by bacteria to enter zooids, they are apparently immediately ‘trapped’ by somatic cells. Free bacteria have never been observed inside the zooidal cavity, another argument against their passage through the wall of the ooecial vesicle to the ovicell.We have shown that regardless of the as-yet-unknown mode of FB development, these temporary organs and the symbionts inside them undergo seasonal changes. Early and mature FBs with non-modified morphology and ‘healthy’ bacteria were found in young zooids only in the colonies collected in June. In one of these colonies, older zooids contained FBs at the initial stage of degradation. In the same month, one colony possessed FBs either at the early-advanced degradation stage (in young zooids) still containing numerous bacteria or at the late-advanced and even terminal stages in old, presumably overwintered zooids. At the initial stage of modification, the slightly developed ‘interlayer’ space between the inner and outer FB cell layers contained fibrils (presumed virus-like particles). Interestingly, the cells of the inner lining apparently engulfed some of the bacterial cells by phagocytosis, supporting our interpretation of the origin of these cells from coelomocytes in this species.One to two months later (August–September), all examined FBs were either at the middle-, or the late-advanced stages of degradation. The number of bacteria in the FB internal cavity distinctly decreased, the inner layer of cells became thinner, and in some regions remained only as a double membrane. The protein-synthesis apparatus was seen only occasionally, and engulfed bacteria were no longer visible inside the inner cells. A wide ILS, formed between the cells of the outer and inner layers, contained abundant putative virus-like particles (Figs. 2, 10). All these FBs were recorded in feeding and non-feeding zooids in the non-growing distal parts of colony branches.Polypide recycling and a seasonal drop of planktonic food alone cannot explain these changes. Firstly, modified FBs were found in zooids with both degenerated and functioning polypides. Secondly, the initial stages of FBs degradation were detected in June when phytoplankton is abundant in the White Sea (e.g.87). We therefore propose the following scenario for the sequence of changes in FBs and their possible causes. In June, newly-formed zooids build funicular bodies containing bacteria that were acquired either from outside or via internal transfer from older colony parts. During that month, FBs begin to degrade. This process continues throughout the rest of the summer. The mid-advanced stage of FB degradation, with few modified bacteria surviving and distributed on the periphery of the FB cavity, was recorded in August. This stage is reminiscent of the final stage in Lutaud’s49 descriptions of the gradual destruction of FBs accompanied by the disappearance of bacterial symbionts in Bugulina turbinata. In late September, FBs change and bacteria disappear, probably through viral lysis (the induction of prophages) (Fig. 2). Young and non-modified FBs were never encountered in August and September, indicating that development of FBs occurs only in young zooids at the growing tips of colony branches in June.Vishnyakov et al.20 recently described the degradation of symbiotic prokaryotes in Bugula neritina and Paralicornia sinuosa accompanied by a change in bacterial morphology similar to bacteriophage-mediated lysis. Degradation process was accompanied by the appearance of polyhedral VLPs in B. neritina and by the formation of structures similar to the so-called metamorphosis-associated contractile complexes (MACs) in P. sinuosa. These complexes are phage-related structures whose activity eventually results in the cell lysis, see88. Although the fate of FBs was not analyzed in their paper, these two VLP variants were observed both inside the bacteria and in the FB internal space by Vishnyakov and co-authors20.In D. fruticosa, presumed VLPs in the form of spherical complexes (as clusters of straight filaments) are present inside bacteria and, together with their fragments, inside the ILS. The filaments were also observed in the free state, frequently curved (apparently flexible) in the internal cavity of FBs. ILS was mostly filled with ‘fibrils’ (potentially representing modified/corrupted filaments) and ‘globules’, although filaments were incidentally recorded inside ILS too. We suggest that, in D. fruticosa, filaments may represent bacteriophage virions. This interpretation is supported by their appearance being associated with the degradation of bacterial cells, as in the case of VLP in B. neritina and P. sinuosa. Nonetheless, the morphology of the spherical complexes built from the filaments in D. fruticosa is unique: they do not resemble any known group of bacterial viruses. We found these putative VLPs inside ILS between the outer and inner FB cell layers. It remains unclear whether they travel there from the FB internal cavity or self-assemble inside ILS from individual filaments that were incidentally met there too. We add that spherical complexes, complete or partial, were recorded inside ILS of the funicular bodies in non-overwintered (collected on 31 September 2019) and presumably overwintered (14 June 2021) zooids.A filamentous morphology is known from only one bacteriophage from the order Tubulavirales, which includes two families Inoviridae and Plectoviridae89. Although inoviruses or plectoviruses have never been reported to assemble in any regular macrocomplexes, the ability of the filamentous Pf phages (inoviruses) of Pseudomonas aeruginosa to form nearly regular liquid crystalline assemblages was recently demonstrated90 (see also review91). Interestingly the formation of such crystals required interactions with bacterial or eukaryotic polymeric molecules such as polysaccharides, DNA and probably mucin90,91.The development of filamentous phages in bacterial cells usually does not kill the cells because these viruses assemble, along with extrusion from the infected bacterium, without disrupting its cell wall (reviewed in92). However, cell death mediated by filamentous prophage induction has been reported in P. aeruginosa due to the emergence of so-called superinfective phage variants93,94,95. Accordingly, the degradation of bacteria in D. fruticosa associated with bacterial viruses is possible.Since the assembly of known filamentous phages is associated with their extrusion from the cell, virions should not accumulate inside bacteria. Although filamentous assembly intermediates may be present (see92), they are not expected to accumulate in such large quantities and/or form superstructures in the bacterial cytoplasm like the ones we found in D. fruticosa (Fig. 9). Our observations revealed nothing resembling the extrusion of a filamentous phage from the surface of a bacterial cell. Therefore, if the described filaments are indeed VLPs, they may represent a new type of bacterial virus.Even though many details remain unknown, we assume that the filaments, ‘globules’, ’fibrils’, and spherical structures in D. fruticosa are of viral origin. Their development following the total disappearance of bacteria in FBs indicates their bacteriophage nature. If so, our observations support the idea that viruses control the number of symbionts in their bryozoan host20, as has been reported in some insects85,86,96. More

Identification of windthrow eventsLandsat images from January 1st 2018 to December 31st 2019 were filtered on 20% or less of cloud coverage, and only the least cloudy image at each location was selected to make an image composite covering the entire Amazon region. In total, 395 least cloudy Landsat 8 images within the Amazon boundary during 2018–2019 were selected and displayed in false color (red: shortwave infrared band, green: near-infrared band, blue: red band) on Google Earth Engine for windthrow events identification (Supplementary Fig. 6). Hollow regions on Supplementary Fig. 6 (2.8% of the total area of the Amazon region) indicated that no clear images with 1 year before the identification were displayed in bright green colors (due to reflectance in near-infrared band from the pioneer species). “Old” windthrows account for ~80% of total identified windthrows, and they were verified using historical Landsat images that can go as far as 1984 (when Landsat 5 was launched). “Old” windthrows were validated once they were found with clear shape and more distinguish color on the historical Landsat images (Supplementary Fig. 7c). 10–15% of “old” windthrows without fan-shape were eliminated from this study because it was hard to identify if they were windthrows or other types of forest disturbance. The minimum size of windthrows identified in this study was 25,000 m2. This process generated the location and rough size of 1012 visible (both “old” and “new”) windthrow scars with fan-shaped patch, scattered small disturbance pixels tails, and an area of over 25,000 m2 (Supplementary Fig. 8). Based on a gap-size probability distribution function that simulates the entire disturbance gradient from all sizes of windthrows19, the proportion of total tree mortality represented by large windthrows ( >25,000 m2) identified in this study is 0.5–1.1%.Among 1012 visual identified windthrows, the occurrence year of 125 windthrows were identified using Landsat 5,7,8, MODIS, and TRMM dataset (Supplementary Table 2), and 38 windthrows from these 125 windthrows had clear remote sensing evidence to validate their occurring date (Supplementary Table 3). It is difficult to get the accurate year and date of occurrence of all identified windthrows. Previous studies showed that windthrows in the northwestern Amazon took ~20 years to recover to 90% of “pre-disturbance” biomass from all damage classes while forests in the central Amazon took ~40 years to recover40. The biomass recovery depends on the windthrow severity and time since disturbance33. Based on the recovery time (20–40 years) and the time of windthrow identification (2018–2019), we estimated that these 1012 windthrows most likely occurred within 30 years (between 1990 and 2019), and the estimated occurrence period was validated using the range of the occurrence year (1986–2019) of 125 windthrow cases.Windthrow density dataThe windthrow density shown in Fig. 1b was generated using 1012 windthrow points in QGIS45. We created a 2.5° by 2.5° grid map, and the windthrow density was calculated by counting the number of windthrows in each grid. These values were then converted to a density with units of counts of windthrows per 10,000 km2. We chose 2.5 degrees to aggregate the data to make sure that over 50% of grids have at least 1 windthrow event while still preserving the spatial distribution of mean afternoon CAPE over the Amazon. The contour lines displayed in Fig. 1c were generated using the “Contour” function on the windthrow density map in QGIS.Meteorological dataTo derive the correlation between windthrow density and meteorological variables, we used ERA 5 global reanalysis hourly CAPE on single levels from 1979 to present at 0.25° × 0.25° resolution provided by the European Center for Medium-Range Weather Forecasts. ERA 5 CAPE was computed by considering parcels of air departing at different pressure levels below the 35 kPa level, with maximum–unstable algorithm under a pseudo-adiabatic assumption46. Afternoon mean CAPE map was calculated as the average of hourly CAPE data from 17:00–23:00 UTC (13:00–19:00 local time in Amazon) over all the months between 1990 and 2019. We chose to average CAPE over 30 years because these windthrow events occurred in these 30 years and calculating the average can help capture the overall spatial pattern of CAPE and minimize the influence of interannual climate variability on windthrow events.To project future windthrow density in the Amazon for the end of the 21st century, we analyzed meteorological output from 10 ESMs that participated in CMIP6 (https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6). The models used in this research were listed in Supplementary Table 1. We extracted daily surface temperature (tas), specific humidity (huss), surface pressure (ps), temperature (ta) from these models to calculate daily nondilute, near-surface-based, adiabatic CAPE. CMIP6 CAPE was calculated by considering the buoyancy of a near-surface parcel lifted adiabatically to a series of discrete pressure levels (100 kPa to 10 kPa in increments of 10 kPa). CMIP6 CAPE is calculated as follows:$${CAPE}=mathop{sum }limits_{i=1}^{10}{{{{{rm{d}}}}}}p{{{{{rm{H}}}}}}({b}_{i}){b}_{i}$$
(1)
Where ({{{{{rm{d}}}}}}p) = 10 kPa, H is the Heaviside unit step function, and ({b}_{i}=frac{1}{{rho }_{i}}-frac{1}{{rho }_{e,i}}), with ({rho }_{i}) being the parcel density at pressure level i and ({rho }_{e,i}) being the environmental density at pressure level i.The future projections in our analysis were based on SSP585, a high-emission scenario with high radiative forcing by the end of the century. We calculated mean daily CAPE over 1990–2015 as current CMIP6 CAPE and mean daily CAPE over 2070–2099 as future CMIP6 CAPE. Since different approaches were used to calculate ERA 5 CAPE and CMIP6 CAPE47, the absolute CAPE values of the two datasets are not comparable. Therefore, for each ESM model, we scaled future CMIP6 CAPE by multiplying, grid-wise, the delta CAPE generated from an individual model in CMIP6 with the ERA 5 current mean afternoon CAPE (Fig. 1c) as follows:$${delta},{CAPE}=(CAPE_{CMIP6_{,future}},-CAPE_{CMIP6_current})/CAPE_{CMIP6_current}$$
(2)
$$CAP{E}_{scaled_CMIP6_,future}=(1+delta,CAPE)times CAP{E}_{ERA5}{_}_{current}$$
(3)
The delta CAPE indicated the projected increase in CAPE from 1990–2015 to the end of the 21st century. In this way, a scaled CMIP6 future CAPE map was generated for each model, and an ensemble-mean scaled CMIP6 CAPE map over 10 ESM models can be found in Supplementary Fig. 5b. The scaled CMIP6 future CAPE values were within plausible range compared to the ERA 5 current mean afternoon CAPE values, and both current and future CAPE maps were used to produce the increase in area with high CAPE values ( >1023 J kg−1) in Table 1. However, it is worth noting that the scaling with relative changes in delta CAPE (%) is more sensitive to CMIP historical baseline conditions than absolute changes of CAPE (J kg−1), which will likely introduce a larger scaled spread (min/max CAPE changes).The increase in area with storm-favorable environments was calculated as follows:$$Increase=(are{a}_{future}-area_{current})/are{a}_{current}$$
(4)
Where areacurrent is the area of CAPE  > 1023 J kg−1 for current ERA 5 CAPE, and areafuture is the area of CAPE  > 1023 J kg−1 for the scaled CMIP6 future CAPE.A model of windthrow densityWe developed a model based on the relationship between satellite-derived windthrow density and mean afternoon CAPE from the ERA 5 reanalysis over 1990–2019. The non-parametric model provides a look-up table of windthrow density as a function of CAPE within the range of observations. Counts of observed windthrow events and Amazon’s area were separately binned by CAPE using the same bins, producing two histograms of CAPE. The ratio of the former to the latter gives the density of windthrow events (windthrow events per area) as a function of CAPE. To avoid noise at the tails of the histograms, the six CAPE bins were chosen such that each bin would have about the same number of windthrow events (either 168 or 169). The total number of windthrow events is given by the sum over bins of the product of windthrow density and area. The minimum and maximum of current ERA 5 mean afternoon CAPE was 42 and 1549. The minimum CAPE value of the first bin was extended to 0 and the maximum CAPE value of the last bin was extended to infinity under the assumption that the windthrow density is similar for neighboring values. Based on the windthrow density and CAPE relationship used in the model, it is the increase in the area with high CAPE that then leads to an increase in the number of windthrow events.It is worth noting that the future windthrow density produced by models may be underestimated because the windthrow observations within regions with high CAPE were incomplete due to high cloud coverage. Moreover, the non-parametric model makes the conservative assumption that the windthrow density does not increase at higher, as-yet unobserved values of mean afternoon CAPE.Future projections of windthrow densityWe combined the non-parametric relationship (Fig. 2a) with the future CAPE map generated from the ten CMIP6 ESMs (adjusted by ERA 5 mean CAPE values) to estimate the changes in windthrow density at the end of the 21st century. We estimated uncertainties for windthrow density projections by combining information about model-to-model differences. The analysis yielded a set of 10 estimates. The overall windthrow density increase and uncertainty were estimated using the mean increase and one standard deviation from the ensemble of the 10 models. More

MAWF. National Rangeland Management Policy & Strategy. Restoring Namibia’s Rangelands (2012).SAIEA. Strategic environmental assessment of large-scale bush thinning and value-addition activities in Namibia. (2016).de Klerk, J. N. Bush Encroachment in Namibia. Report on Phase 1 of the Bush Encroachment Research, Monitoring and Management Project. (2004).Muntifering, J. R. et al. Managing the matrix for large carnivores: A novel approach and perspective from cheetah (Acinonyx jubatus) habitat suitability modelling. Anim. Conserv. 9, 103–112 (2006).Article 

Google Scholar 
Marker, L. L., Dickman, A. J., Mills, M. G. L., Jeo, R. M. & Macdonald, D. W. Spatial ecology of cheetahs on north-central Namibian farmlands. J. Zool. 274, 226–238 (2008).Article 

Google Scholar 
Wang, J. et al. Impacts of juniper woody plant encroachment into grasslands on local climate. Agric. For. Meteorol. 307, 108508 (2021).Article 
ADS 

Google Scholar 
Shen, X. et al. Effect of shrub encroachment on land surface temperature in semi-arid areas of temperate regions of the Northern Hemisphere. Agric. For. Meteorol. 320, 108943 (2022).Article 
ADS 

Google Scholar 
Shen, X. et al. Vegetation greening, extended growing seasons, and temperature feedbacks in warming temperate grasslands of China. J. Clim. 35, 5103–5117 (2022).Article 
ADS 

Google Scholar 
Martins, A. R. O. & Shackleton, C. M. Population structure and harvesting selection of two palm species ( Hyphaene coriacea and Phoenix reclinata ) in Zitundo area, southern Mozambique. For. Ecol. Manage. 398, 64–74 (2017).Article 

Google Scholar 
Brown, G. W., Murphy, A., Fanson, B. & Tolsma, A. The influence of different restoration thinning treatments on tree growth in a depleted forest system. For. Ecol. Manage. 437, 10–16 (2019).Article 

Google Scholar 
Belsky, A. Influences of trees on savanna productivity: Tests of shade, nutrients and grass–tree competition. Ecology 75, 922–932 (1994).Article 

Google Scholar 
Hagos, M. G. & Smit, G. N. Soil enrichment by Acacia mellifera subsp. detinens on nutrient poor sandy soil in a semi-arid southern African savanna. J. Arid Environ. 61, 47–59 (2005).Ludwig, F., Kroon, H. D., Berendse, F. & Prins, H. H. T. The influence of savanna trees on nutrient, water and light availability and the understorey vegetation. Plant Ecol. 170, 93–105 (2004).Article 

Google Scholar 
Ridolfi, L., Laio, F., D’Odorico, P. & D’Odorico, P. Fertility island formation and evolution in dryland ecosystems. Ecol. Soc. 13, 13 (2008).Article 
MATH 

Google Scholar 
Wiegand, K., Ward, D. & Saltz, D. Multi-scale patterns and bush encroachment in an arid savanna with a shallow soil layer. J. Veg. Sci. 16, 311–320 (2005).Article 

Google Scholar 
Burke, A. Savanna trees in Namibia – Factors controlling their distribution at the arid end of the spectrum. Flora Morphol. Distrib. Funct. Ecol. Plants 201, 181–201 (2006).
Google Scholar 
Buyer, J. S., Schmidt-Küntzel, A., Nghikembua, M., Maul, J. E. & Marker, L. Soil microbial communities following bush removal in a Namibian savanna. Soil 2, 101–110 (2016).Article 
ADS 
CAS 

Google Scholar 
Dwivedi, V. & Soni, P. A review on the role of soil microbial biomass in eco-restoration of degraded ecosystem with special reference to mining areas. J. Appl. Nat. Sci. 151–158 (2011). https://doi.org/10.31018/jans.v3i1.173.Smit, G. N. An approach to tree thinning to structure southern African savannas’ for long-term restoration from bush encroachment. J. Environ. Manage. 71, 179–191 (2004).Article 
CAS 

Google Scholar 
MAWF. Forestry and environmental authorisation process for bush harvesting projects. 34 (2017).Smit, G. N., de Klerk, J. N., Schneider, M. B. & van Eck, J. Detailed assessment of the biomass resource and potential yield in a selected bush encroached area of Namibia. 141 (2015).Marker, L. et al. and Distribution. in Cheetahs: Biology and Conservation: Biodiversity of the World: Conservation from Genes to Landscapes (eds. Marker, L., Boast, L. k. & Schmidt-Küntzel, A.) 9–20 (John Fedor, 2018). https://doi.org/10.1016/B978-0-12-804088-1.00004-6.NAPHA. Namibia Professional Hunting Association. http://www.napha-namibia.com/conservation/huntable-species/carnivora/ (2015).SWA. Nature Conservation Ordinance, 1975 (No. 4 of 1975). vol. 1975 (1975).Marker, L. et al. The status of Key pre species and the Consequences of Prey Loss for Cheetah Conservation in North and West Africa. in Cheetahs: Biology and Conservation: Biodiversity of the World: Conservation from Genes to Landscapes. (eds. Marker, L., Boast, L. k. & Schmidt-Küntzel, A.) 151–161 (John Fedor, 2018). https://doi.org/10.1016/B978-0-12-804088-1.00004-6.Kiruki, H. M., van der Zanden, E. H., Gikuma-Njuru, P. & Verburg, P. H. The effect of charcoal production and other land uses on diversity, structure and regeneration of woodlands in a semi-arid area in Kenya. For. Ecol. Manage. 391, 282–295 (2017).Article 

Google Scholar 
Harmse, C. J., Kellner, K. & Dreber, N. Restoring productive rangelands: A comparative assessment of selective and non-selective chemical bush control in a semi-arid Kalahari savanna. J. Arid Environ. 135, 39–49 (2016).Article 
ADS 

Google Scholar 
Nghikembua, M. T. et al. Response of wildlife to bush thinning on the north central freehold farmlands of Namibia. For. Ecol. Manage. 473, 118330 (2020).Article 

Google Scholar 
Soto-Shoender, J. R., McCleery, R. A., Monadjem, A. & Gwinn, D. C. The importance of grass cover for mammalian diversity and habitat associations in a bush encroached savanna. Biol. Conserv. 221, 127–136 (2018).Article 

Google Scholar 
Strohbach, B. J. Environmental information service, Namibia for the Ministry of Environment and Tourism, the Namibian Chamber of Environment and the Namibia University of Contribution to the knowledge of southern African Lepismatidae. Namibian J. Environ. 8327, 14–33 (2017).
Google Scholar 
Smit, N. BECVOL 3: An expansion of the aboveground biomass quantification model for trees and shrubs to include the wood component. Afr. J. Range Forage Sci. 31, 179–186 (2014).Article 

Google Scholar 
Zimmerman, I. Causes and Consequences of Fenceline Contrasts in Namibia. (Free State, Bloemfontein, South Africa, 2009).Dwyer, J. M. & Mason, R. Plant community responses to thinning in densely regenerating Acacia harpophylla forest. Restor. Ecol. 26, 97–105 (2018).Article 

Google Scholar 
Thomas, S. C. & Martin, A. R. Carbon content of tree tissues: A synthesis. Forests 3, 332–352 (2012).Article 

Google Scholar 
Chave, J. et al. Improved allometric models to estimate the aboveground biomass of tropical trees. Glob. Chang. Biol. 20, 3177–3190 (2014).Article 
ADS 

Google Scholar 
Djomo, A. N. & Chimi, C. D. Tree allometric equations for estimation of above, below and total biomass in a tropical moist forest: Case study with application to remote sensing. For. Ecol. Manage. 391, 184–193 (2017).Article 

Google Scholar 
Ganamé, M., Bayen, P., Dimobe, K., Ouédraogo, I. & Thiombiano, A. Aboveground biomass allocation, additive biomass and carbon sequestration models for Pterocarpus erinaceus Poir. in Burkina Faso. Heliyon 6, (2020).Picard, N., Saint-André, L. & Henry, M. Manual for building tree volume and biomass allometric equations: From field measurement to prediction. Food and Agricultural Organization of the United Nations, Rome, and Centre de Coopération Internationale en Recherche Agronomique pour le Développement. (2012).Feyisa, K. et al. Allometric equations for predicting above-ground biomass of selected woody species to estimate carbon in East African rangelands. Agrofor. Syst. 92, 599–621 (2018).Article 

Google Scholar 
Boys, J. M. & Smit, N. G. Development of an Excel Based Bush Biomass Quantification Tool. (2020).Ngomanda, A. et al. Site-specific versus pantropical allometric equations: Which option to estimate the biomass of a moist central African forest?. For. Ecol. Manage. 312, 1–9 (2014).Article 

Google Scholar 
CCF. Cheetah Conservation Fund Bush PTY (Ltd). https://bushblok.com/management-plan (2019).Wykstra, M. et al. Improved and Alternative Livelihoods: Links Between Poverty Alleviation, Biodiversity, and Cheetah Conservation. in Cheetahs: Biology and Conservation: Biodiversity of the World: Conservation from Genes to Landscapes (eds. Marker, L., Boast, L. k. & Schmidt-Küntzel, A.) 223–237 (John Fedor, 2018).Zimmermann, I. et al. The influence of two levels of debushing in Namibia’s Thornbush Savanna on overall soil fertility, measured through bioassays. Namibian J. Environ. 1, 52–59 (2017).
Google Scholar 
Nghikembua, M. T. et al. Restoration thinning reduces bush encroachment on freehold farmlands in north-central Namibia. For. An Int. J. For. Res. 1–14 (2021). https://doi.org/10.1093/forestry/cpab009.Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Mendelsohn, J., Jarvis, A., Roberts, C. & Robertson, T. Atlas of Namibia: A portrait of the land and its people. (2003).Curtis, B. A. & Mannheimer, C. A. Tree Atlas of Namibia; National Botanical Research Institute (NBRI). (National Botanical Research Institute, 2005).Mannheimer, C. & Curtis, B. Le Roux and Muller’s Field Guide to the Trees & Shrubs of Namibia. (Macmillan Education Namibia (PTY) LTD, 2009).Honsbein, D., Shiningavamwe, K., Iikela, J. & de la Puerta Fernandez, Maria, L. Animal Feed from Namibian Encroacher Bush. (2017).Coates Palgrave, K. Trees of Southern Africa. (Struik publishers, 1993).Nghikembua, M., Harris, J., Tregenza, T. & Marker, L. Spatial and temporal habitat use by GPS collared male cheetahs in modified bushland habitat. Open J. For. 06, 269–280 (2016).
Google Scholar 
Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S. (Springer, 2002).Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
Mehtatalo, L. & Lappi, J. Forest biometrics with examples in R. (Taylor & Francis Inc, 2020).R Core Team. R: A language and environment for statistical computing (2021).Birch, C. & Middleton, A. Economics of Land Degradation Related To Bush Encroachment in Namibia. (2017).Neke, K. S., Owen-Smith, N. & Witkowski, E. T. F. Comparative resprouting response of Savanna woody plant species following harvesting: the value of persistence. For. Ecol. Manage. 232, 114–123 (2006).Article 

Google Scholar 
Smit, N. Response of Colophospermum mopane to different intensities of tree thinning in the Mopane Bushveld of southern Africa. Afr. J. Range Forage Sci. 31, 173–177 (2014).Article 

Google Scholar 
DAS. Bush Control Manual. (John Meinert Printing, 2017).Leinonen, A. Wood chip production technology and costs for fuel in Namibia. VTT Tiedotteita – Valtion Teknillinen Tutkimuskeskus (2007).Chakanga, M. A preliminary analysis of the Economic Plots Research Data of the CCF Bush Project. (2003).West, P. W. Thinning. in Growing Plantation Forests vol. 9783319018 115–129 (Springer International Publishing, 2014).Dwyer, J. M., Fensham, R. & Buckley, Y. M. Restoration thinning accelerates structural development and carbon sequestration in an endangered Australian ecosystem. J. Appl. Ecol. 47, 681–691 (2010).Article 

Google Scholar 
Groengroeft, A., de Blécourt, M., Classen, N., Landschreiber, L. & Eschenbach, A. Acacia trees modify soil water dynamics and the potential groundwater recharge in savanna ecosystems. in Climate change and adaptive land management in southern Africa – assessments, changes, challenges, and solutions (ed. (eds. Revermann, R. et al.) 177–186 (Klaus Hess Publishers, 2018).Richter, C. G. F., Snyman, H. A. & Smit, G. N. The influence of tree density on the grass layer of three semi-arid savanna types of southern africa. Afr. J. Range Forage Sci. 18, 103–109 (2001).Article 

Google Scholar 
Stafford, W. et al. The economics of landscape restoration: Benefits of controlling bush encroachment and invasive plant species in South Africa and Namibia. Ecosyst. Serv. 27, 193–202 (2017).Article 

Google Scholar 
Joubert, D. F., Smit, G. N. & Hoffman, M. T. The influence of rainfall, competition and predation on seed production, germination and establishment of an encroaching Acacia in an arid Namibian savanna. J. Arid Environ. 91, 7–13 (2013).Article 
ADS 

Google Scholar  More

Quantifying heating of coastal ecosystemsA variety of metrics have been developed to quantify MHWs, but to date, they have mainly focused on surface heating evident from SST. Surface MHWs have been defined based on exceedance of 90% confidence intervals calculated seasonally from historical SST10,11. Defined in this way, the threshold temperature for quantifying a MHW is likely to vary seasonally from, for example, winter to summer. However, the physiological and ecological relevance of seasonally changing thresholds remains unclear, especially for tropical biota such as corals that typically inhabit a relatively narrow range of temperatures close to their physiological thermal limits12,13,14. Assessing heating in the specific context of temperature-induced coral bleaching has instead focused on calculating cumulative degree heating, sensitive to both the magnitude and duration of heating, above a fixed, putative ‘bleaching threshold’ defined by the local Maximum Monthly Mean SST (MMM); i.e., the mean summer-time peak temperature predicted to initiate coral stress and bleaching. Such heat accumulation has most commonly been expressed as Degree Heating Weeks (DHW in °C-weeks) accumulated over a 12-week period15,16, or sometimes using monthly SST data to compute Degree Heating Months (DHM in °C-months) over 3 months17,18. Numerous studies have documented coral bleaching in shallow water linked to periods of anomalously high SST and accumulated DHW, especially during El Niño events1,2,4,19,20,21. However, DHW generally only explain a limited proportion of the observed variation in bleaching, even for communities in very shallow water (e.g., 50% of bleaching variation at 2 m depth7). A number of studies have documented bleaching that was less than that predicted based on contemporaneous SST and DHW. This discrepancy has sometimes been hypothesised to reflect ongoing coral acclimatisation to increasing temperatures and/or shifts in community composition towards more heat-tolerant genotypes and species22,23,24,25,26. Bleaching rates higher than predicted by SST, under limited DHW, have also been documented, and can be species-specific and more pronounced in heat-intolerant cryptic species27. It is unclear to what extent the relatively limited power of SST metrics such as DHW to predict coral bleaching results from an incomplete description of environmental conditions, a lack of nuance in describing biological thresholds and organisms’ reactions to elevated temperatures, genetic variation and cryptic species27,28, or, most likely, a combination of factors.Estimates of surface MHW severity are also highly sensitive to both the spatial and temporal scales of SST data considered. Present-day satellite products allow degree heating to be calculated at relatively fine spatial and temporal resolution using SST data measured daily over pixels of ~25 km2 size (e.g., NOAA Coral Reef Watch29). While many heating assessments in the context of coral bleaching continue to focus on long-term heating based on temperature anomalies accumulated across months (12 weeks or 3 months for DHW or DHM, respectively)18,20,30, here we focus on higher-resolution heating calculated as Degree Heating Days (DHD in °C-days) using daily SST data. Calculation of DHD over 12-day windows is analogous to the more coarse resolution DHW (i.e., weekly data over 12 weeks)31 but with different units, finer temporal resolution, and shorter time lags between the actual elevation of environmental temperatures and the resulting accumulated heating metric (see the DHD and DHW comparison in Fig. S5).To investigate the role of spatial scale in characterising MHWs, we analysed SST at a range of scales around Moorea between 1985 and 2019: 2° × 2° (~50,000 km2), 1° × 1° (~12,300 km2) and 0.1° × 0.1° (~100 km2) (see Fig. S1). The importance of spatial scales in heating severity apparent at the surface is reinforced by the regional heterogeneity in SST (see Fig. 2 and data animations in the Supplementary Information). Considering heating at local scales and finer temporal resolutions maximises the potential to detect, characterise and compare heating events at the surface in a way that is relevant to in situ MHW conditions (see also Guo et al.32. for the importance of scales in MHW assessments). For instance, when assessed using NOAA’s ‘Regional Bleaching Heat Stress Gauges’ for the Society Islands, DHW reached 4.54 °C-weeks in April 2016 and 5.35 °C-weeks in May 201933, indicating moderate likelihood of bleaching during both events, although inherent to calculations of accumulated DHW, the maximum heating was centred multiple weeks after the in situ heating events actually occurred (e.g., see Fig. S5e, f). Heat accumulation using higher resolution DHD across 12-day windows around Moorea itself (2° × 2°) was more closely aligned temporally to the heating events, suggesting a much higher bleaching risk in 2016 (6.83 °C-days) than 2019 (2.10 °C-days; Table S3). This likely reflects the localised heterogeneity in SST during the 2019 MHW, when hotter surface conditions prevailed north of Moorea (Fig. 1b). At a local scale, assessed using SST within a ~10 × 10 km area north of Moorea (see Fig. S1), heating was similar to regional estimates in 2016 (5.73 °C-days), but in 2019 revealed an intense, localised heatwave over Moorea’s north shore (15.4 °C-days; see Fig. 1e, f; Table S1). Because of these demonstrated advantages of higher temporal-resolution analysis, we rely on DHD rather than DHW for our analysis of MHW patterns among years at Moorea.Fig. 2: Regional sea-surface temperature (SST) variability during the peak of the 2016 and 2019 surface marine heatwaves around Moorea suggested hotter conditions during 2016.Panels show SST over a, b 20° × 20° and c, d 10° × 10° during 2016 and 2019, respectively, focusing on the date of peak SST observed over Moorea’s north shore (8 April 2016 and 4 April 2019). Dashed squares in (a, b) show extent of (c, d) and those in (c, d) the extent of data shown in Fig. 1a, b (2° × 2°). Coastlines based on Wessel and Smith77.Full size imageContrasting surface and subsurface heatingComparing surface and subsurface MHWs is challenging for many coastal ecosystems, even once SST data of appropriately fine spatial and temporal scale are obtained, due to a lack of long-term, in-situ temperature data through which to assess mean climatological patterns below the sea surface. Moorea represents one of the few coral reef systems with consistent in-situ observations over timescales (decades) and depths (sea surface to 40 m) relevant to understanding the oceanographic drivers of subsurface heating and their impacts on coral bleaching. Our analysis of long-term SST records indicates there have been 16 local-scale (0.1° × 0.1°) surface MHWs over the north shore of Moorea relative to a MMM-based bleaching threshold of 29.8 °C (Table S4), compared to 14 regional-scale (2° × 2°) events (Table S3). Localised heating over the north shore was often greater than regional SST would suggest, with the hottest event recorded in 2003 reaching 17.7 °C-days locally (Table S4), compared to 15.5 °C-days regionally (Table S3). Further details on historical events can be found in ‘Surface MHW history around Moorea’ in the Supplementary Information, which provides context for the six more recent surface MHWs that have occurred over the north shore since continuous in-situ, reef-level observations began at Moorea in 2005 (MHWs in 2007, 2012, 2015, 2016, 2017, and 2019; Table S4). Two recent, contrasting events in 2016 and 2019 demonstrate the extent to which thermal environments at depth, and the associated severity of coral bleaching, can vary substantially from predictions based on sea-surface conditions.MHW severity based only on SST can miss important information on the conditions experienced by organisms at depths greater than the surface skin layer quantified through remote sensing, which may only be few millimetres thick34. For example, although the localised peak in sea-surface temperatures were essentially identical between 2016 (30.1 °C; Fig. 1a) and 2019 (30.2 °C; Fig. 1b), and regional surface heating metrics and warnings were similar33, markedly different heat accumulation occurred due to the different duration that temperatures remained above the putative coral bleaching threshold (MMM + 1 = 29.8 °C; up to 2 days in April 2016 compared to up to 11 days in April 2019; Table S1). Yet these results from local SST—of similar SST maximums in 2016 and 2019, but longer durations above the threshold in 2019—only capture some of the significant differences that led to constating MHW severity and ecological outcomes between years and across depths.Daily average temperatures measured in situ at reef level in water depths of 10–40 m over ~15 years are well correlated with daily SST (r2 = 0.94–0.78 at 10–40 m). However, the strength of the relationship between SSTs and in situ temperatures declines with increasing water depth, even when in situ temperatures are averaged to a daily resolution31,35. The potential for subsurface attenuation of heating over coral reefs has previously been demonstrated using high-resolution in-situ water temperature data in the context of both regional upwelling and local internal-wave climates31,36,37, with observations of periodic transport of deeper, cooler water onto reef habitats at a large number of reefs globally31,35,38. The propagation of internal-wave energy is associated with significant vertical displacements of density isopycnals and isotherms; e.g., ~60 m displacements along the Hawaiian Ridge39. Upon encountering a sloping bottom, internal-wave dynamics become complex, and, for habitats on the fore reef slope, typically result in rapid, periodic cooling (rather than oscillations around a mean temperature) as water masses associated with e.g., 24–27 °C isotherms are vertically advected onto the reef and recede again31,40,41. In deeper reef habitats there may also be periodic heating associated with exposure to warmer surface water masses when internal waves lead to downward displacement of isotherms31, but the overall magnitude of any resulting net heating is small across the depth range considered here (i.e., no average heating at depths of 40 m and less; Fig. S4).To separate the effects of low- and high-frequency processes driving heating across the reef slope, we used a filtering approach specifically designed and validated by Wyatt et al.31. to estimate coral reef thermal regimes without internal waves. In situ temperatures were filtered to isolate variability at frequencies higher and lower than the local inertial period (~40.0 h), effectively removing the effects of internal waves from lower frequency processes (i.e., multi-day weather patterns and seasonal effects; see Fig. S3). Contrasting the observed and filtered in-situ temperature variations (black and white lines, respectively, in Fig. 3) highlights differences in the processes driving the 2016 and 2019 subsurface MHWs. During the 2019 MHW around Moorea, the filtered, or ‘non-internal wave’ (NIW), temperatures closely resembled the observed temperatures (Fig. 3e–h, m–p) implying limited internal-wave cooling (IWC). Consistent warming across the water column was evident in 2019 and temperatures remained above the coral bleaching threshold for multiple days during early to late April (Table S1). By contrast, the 2016 MHW was characterized by temperatures remaining generally below the bleaching threshold and significant high-frequency variability indicative of IWC across depths, such that temperatures only exceed the predicted bleaching threshold for hours or less at a time (Fig. 3a–d, i–l; Table S1). The high-resolution temperature observations show that IWC was greatly reduced during 2019 (Fig. 3e–h). The power spectral density of observed temperatures, concentrated at semi-diurnal frequencies and consistent with internal-wave forcing at this location41, was significantly higher in 2016 and lower in 2019 than the average across years at 10 m depth (Fig. 4a; see inset). In deeper water at 20–40 m depths, temperature variance within the semi-diurnal frequency band increased relative to shallow depths and became more similar between the two events, such that at 40 m the semi-diurnal variability was equivalent in 2016 and 2019 (Fig. 4b–d; see insets). However, this similarity in temperature variance does not indicate an equivalent magnitude of IWC, since variability during 2019 (Fig. 3p) was around a warmer background temperature closer to the coral bleaching threshold. Extending the comparison of 2016 and 2019 to other recent local MHWs demonstrates two distinct types of events: greater IWC across reef depths during 2012, 2015, and 2016 MHWs, versus reduced IWC during the 2007, 2017 and 2019 MHWs (Fig. 5).Fig. 3: Contrasting reef-level temperature variations across depths on Moorea’s north shore during the 2016 and 2019 marine heatwaves.Panels on the left show the observed high-frequency water temperature variations (black lines, measured at 2-min intervals) during the hottest months (Apr–May) in a–d 2016 and e–h 2019 at a, e 10, b, f 20, c, g 30 and d, h 40 m depths. Right panels focus on relative variation during the heatwave peaks across the same depths: i–l* 06–12 Apr 2016 and m–p 11–17 Apr 2019. Non-internal-wave temperature variations are shown based on observed temperatures filtered to remove the high-frequency influence of internal waves (white lines). The satellite-derived sea-surface temperatures (SST; grey line) are shown for comparison to in situ temperatures. The horizontal dashed line shows the ‘bleaching threshold’ (maximum monthly mean + 1 °C) and the background shading provides a reference relative to temperatures above (red), equal (yellow) and below (blue) this threshold. The red dashed squares denote the axis limits in the right panels. *Note: 40 m logger during 2016 incorrectly recorded at 2-h interval.Full size imageFig. 4: Reduced semi-diurnal temperature variability during the summer of 2019 in shallower water on the north shore of Moorea.Power spectral density (PSD) plots (logarithmic scale) were computed within a 12-day window at a 10 m, b 20 m, c 30 m, and d 40 m water depths during the summer months (Dec–May) in 2016 (blue), 2019 (red), and 2004–2018 (black; excluding 2016 and 2019). Shading shows the 95% confidence intervals for each PSD. The tidal constituents (dotted lines) show variance consistent with semi-diurnal (M2) forcing across depths, with diurnal (K1) forcing in 10 m of water along with some variability consistent with the shallow water lunar overtide (M4). Insets show details of semi-diurnal differences at each depth.Full size imageFig. 5: Comparison of internal-wave cooling (IWC) across depths during recent surface marine heatwaves (MHWs) around Moorea.Based on average a daily temperature variance (in °C2) and b IWC ((overline{{{{{{rm{IWC}}}}}}}) in °C) during the six recent local MHWs (see Table S4 for dates) that can be grouped into: (1) high IWC events during 2012 (blue), 2015 (green) and 2016 (purple); and, (2) low IWC events that coincided with bleaching events in 2007 (red) and 2019 (orange), along with early 2017 (yellow). The daily variance and (overline{{{{{{rm{IWC}}}}}}}) during Dec–Apr across all years (2005–2019) is shown for reference (black dashed lines). Contours of c–f heat accumulation as degree heating days (DHD in °C days) and g–j degree cooling days due to internal waves (DCDIW in °C days) across depths are shown for 2007, 2012, 2016 and 2019. Due to data gaps in the in situ records, contours are not shown for the 2015 or 2017 events.Full size imageThe magnitude of temperature fluctuations produced by IWC, i.e., occurring at the semi-diurnal frequency, are of similar magnitude to long-term ocean warming and climate change threatening coral reefs globally. The average IWC ((overline{{{{{{rm{IWC}}}}}}})) during the high-IWC MHWs (2012, 2015 and 2016) was between 0.14 and 0.60 °C (Fig. 3a; Table S4) and comparable to the overall SST increase measured over tropical coral reefs during the last four decades (~0.65 °C18). As a result of the subsurface cooling caused by internal waves, the 2016 MHW, which was moderate at the surface (5.7 °C-days), was mild at 10 m ( More

Feit, B. et al. Landscape complexity promotes resilience of biological pest control to climate change. Proc. R. Soc. B. 288, 20210547 (2021).Article 

Google Scholar 
Hall-Spencer, J. M. & Harvey, B. P. Ocean acidification impacts on coastal ecosystem services due to habitat degradation. Emerg. Top. Life Sci. 3, 197–206 (2019).Article 
CAS 

Google Scholar 
Loke, L. H. L. & Todd, P. A. Structural complexity and component type increase intertidal biodiversity independently of area. Ecology 97, 383–393 (2016).Article 

Google Scholar 
Oliver, T. H. et al. Biodiversity and resilience of ecosystem functions. Trends in Ecol. Evol. 30, 673–684 (2015).Article 

Google Scholar 
Bullock, J. M. et al. Future restoration should enhance ecological complexity and emergent properties at multiple scales. Ecography ecog. 4, 05780 (2022).Ortega, J. C. G., Thomaz, S. M. & Bini, L. M. Experiments reveal that environmental heterogeneity increases species richness, but they are rarely designed to detect the underlying mechanisms. Oecologia 188, 11–22 (2018).Article 

Google Scholar 
Griffin, J. N., Byrnes, J. E. K. & Cardinale, B. J. Effects of predator richness on prey suppression: a meta-analysis. Ecology 94, 2180–2187 (2013).Article 

Google Scholar 
Katano, I., Doi, H., Eriksson, B. K. & Hillebrand, H. A cross-system meta-analysis reveals coupled predation effects on prey biomass and diversity. Oikos 124, 1427–1435 (2015).Article 

Google Scholar 
Loke, L. H. L., Ladle, R. J., Bouma, T. J. & Todd, P. A. Creating complex habitats for restoration and reconciliation. Ecol. Eng. 77, 307–313 (2015).Article 

Google Scholar 
Torres-Pulliza, D. et al. A geometric basis for surface habitat complexity and biodiversity. Nat. Ecol. Evol. 4, 1495–1501 (2020).Article 

Google Scholar 
Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366 (2000).Article 

Google Scholar 
Chesson, P. & Kuang, J. J. The interaction between predation and competition. Nature 456, 235–238 (2008).Article 
CAS 

Google Scholar 
Terborgh, J. W. Toward a trophic theory of species diversity. Proc. Natl. Acad. Sci. USA 112, 11415–11422 (2015).Article 
CAS 

Google Scholar 
Pringle, R. M. et al. Predator-induced collapse of niche structure and species coexistence. Nature 570, 58–64 (2019).Article 
CAS 

Google Scholar 
Sandom, C. et al. Mammal predator and prey species richness are strongly linked at macroscales. Ecology 94, 1112–1122 (2013).Article 

Google Scholar 
Grabowski, J. H. Habitat complexity disrupts predator-prey interactions but not the trophic cascade on oyster reefs. Ecology 85, 995–1004 (2004).Article 

Google Scholar 
Crowder, L. B. & Cooper, W. E. Habitat structural complexity and the interaction between bluegills and their prey. Ecology 63, 1802 (1982).Article 

Google Scholar 
Almany, G. R. Does increased habitat complexity reduce predation and competition in coral reef fish assemblages? Oikos 106, 275–284 (2004).Article 

Google Scholar 
Anderson, T. L. & Semlitsch, R. D. Top predators and habitat complexity alter an intraguild predation module in pond communities. J. Anim. Ecol. 85, 548–558 (2016).Article 

Google Scholar 
Brothers, C. A. & Blakeslee, A. M. H. Alien vs predator play hide and seek: How habitat complexity alters parasite mediated host survival. J. Exp. Mar. Biol. Ecol. 535, 151488 (2021).Article 

Google Scholar 
Horinouchi, M. et al. Seagrass habitat complexity does not always decrease foraging efficiencies of piscivorous fishes. Mar. Ecol. Prog. Ser. 377, 43–49 (2009).Article 

Google Scholar 
Ryer, C., Stoner, A. & Titgen, R. Behavioral mechanisms underlying the refuge value of benthic habitat structure for two flatfishes with differing anti-predator strategies. Mar. Ecol. Prog. Ser. 268, 231–243 (2004).Article 

Google Scholar 
Flynn, A. J. & Ritz, D. A. Effect of habitat complexity and predatory style on the capture success of fish feeding on aggregated prey. J. Mar. Biol. Ass. 79, 487–494 (1999).Article 

Google Scholar 
Klecka, J. & Boukal, D. S. The effect of habitat structure on prey mortality depends on predator and prey microhabitat use. Oecologia 176, 183–191 (2014).Article 

Google Scholar 
James, P. L. & Heck, K. L. The effects of habitat complexity and light intensity on ambush predation within a simulated seagrass habitat. J. Exp. Mar. Biol. Ecol. 176, 187–200 (1994).Article 

Google Scholar 
Michel, M. J. & Adams, M. M. Differential effects of structural complexity on predator foraging behavior. Behav. Ecol. 20, 313–317 (2009).Article 

Google Scholar 
Preisser, E. L., Bolnick, D. I. & Benard, M. F. Scared to death? The effects of intimidation and consumption in predator-prey interactions. Ecology 86, 501–509 (2005).Article 

Google Scholar 
Preisser, E. L., Orrock, J. L. & Schmitz, O. J. Predator hunting mode and habitat domain alter nonconsumptive effects in predator-prey interactions. Ecology 88, 2744–2751 (2007).Article 

Google Scholar 
Rypstra, A. L., Schmidt, J. M., Reif, B. D., DeVito, J. & Persons, M. H. Tradeoffs involved in site selection and foraging in a wolf spider: effects of substrate structure and predation risk. Oikos 116, 853–863 (2007).Article 

Google Scholar 
Janssen, A., Sabelis, M. W., Magalhães, S., Montserrat, M. & van der Hammen, T. Habitat structure affects intraguild predation. Ecology 88, 2713–2719 (2007).Article 

Google Scholar 
Grabowski, J. H., Hughes, A. R. & Kimbro, D. L. Habitat complexity influences cascading effects of multiple predators. Ecology 89, 3413–3422 (2008).Article 

Google Scholar 
Hughes, A. R. & Grabowski, J. H. Habitat context influences predator interference interactions and the strength of resource partitioning. Oecologia 149, 256–264 (2006).Article 

Google Scholar 
Bonett, D. G. Meta-analytic interval estimation for standardized and unstandardized mean differences. Psychol. Methods 14, 225–238 (2009).Article 

Google Scholar 
Huey, R. B. & Pianka, E. R. Ecological consequences of foraging mode. Ecology 62, 991–999 (1981).Article 

Google Scholar 
Egger, M., Smith, G. D., Schneider, M. & Minder, C. Bias in meta-analysis detected by a simple, graphical test. BMJ 315, 629–634 (1997).Article 
CAS 

Google Scholar 
Ritchie, E. G. & Johnson, C. N. Predator interactions, mesopredator release and biodiversity conservation. Ecol. Lett. 12, 982–998 (2009).Article 

Google Scholar 
Chaplin-Kramer, R., O’Rourke, M. E., Blitzer, E. J. & Kremen, C. A meta-analysis of crop pest and natural enemy response to landscape complexity: pest and natural enemy response to landscape complexity. Ecol. Lett. 14, 922–932 (2011).Article 

Google Scholar 
Paxton, A. B. et al. Meta-analysis reveals artificial reefs can be effective tools for fish community enhancement but are not one-size-fits-all. Front. Mar. Sci. 7, 282 (2020).Article 

Google Scholar 
Eggleston, D. B., Lipcius, R. N., Miller, D. L. & Coba-Cetina, L. Shelter scaling regulates survival of juvenile Caribbean spiny lobster Panulirus argus. Mar. Ecol. Prog. Ser. 62, 79–88 (1990).Rogers, A., Blanchard, J. L. & Mumby, P. J. Fisheries productivity under progressive coral reef degradation. J. Appl. Ecol. 55, 1041–1049 (2018).Article 

Google Scholar 
Gontijo, L. M. Engineering natural enemy shelters to enhance conservation biological control in field crops. Biol. Control 130, 155–163 (2019).Article 

Google Scholar  More

Pausas, J. G. & Bond, W. J. On the three major recycling pathways in terrestrial ecosystems. Trends Ecol. Evol. 35, 767–775 (2020).
Google Scholar 
Manzano, P. & White, S. R. Intensifying pastoralism may not reduce greenhouse gas emissions: wildlife-dominated landscape scenarios as a baseline in life cycle analysis. Clim. Res. 77, 91–97 (2019).
Google Scholar 
Röös, E. et al. Greedy or needy? Land use and climate impacts of food in 2050 under different livestock futures. Glob. Environ. Change 47, 1–12 (2017).
Google Scholar 
Harwatt, H., Ripple, W. J., Chaudhary, A., Betts, M. G. & Hayek, M. N. Scientists call for renewed Paris pledges to transform agriculture. Lancet Planet. Health 4, E9–E10 (2020).
Google Scholar 
Bar-On, Y. M., Phillips, R. & Milo, R. The biomass distribution on Earth. Proc. Natl Acad. Sci. USA 115, 6506–6511 (2018).CAS 

Google Scholar 
Barnosky, A. D. Megafauna biomass trade-off as a driver of Quaternary and future extinctions. Proc. Natl Acad. Sci. USA 105, 11543–11548 (2008).CAS 

Google Scholar 
Smith, F. A. et al. Exploring the influence of ancient and historic megaherbivore extirpations on the global methane budget. Proc. Natl Acad. Sci. USA 113, 874–879 (2016).CAS 

Google Scholar 
Zimov, S. A., Zimov, N. S., Tikhonov, A. N. & Chapin, F. S. III Mammoth steppe: a high-productivity phenomenon. Quat. Sci. Rev. 57, 26–45 (2012).
Google Scholar 
Adams, J. M., Faure, H., Faure-Denard, L., McGlade, J. M. & Woodward, F. I. Increases in terrestrial carbon storage from the Last Glacial Maximum to the present. Nature 348, 711–714 (1990).CAS 

Google Scholar 
Nogués-Bravo, D., Rodríguez, J., Hortal, J., Batra, P. & Araújo, M. B. Climate change, humans, and the extinction of the woolly mammoth. PLoS Biol. 6, e79 (2008).
Google Scholar 
Bond, W. J. Open Ecosystems: Ecology and Evolution Beyond the Forest Edge (Oxford Univ. Press, 2019).Bakker, E. S. et al. Combining paleo-data and modern exclosure experiments to assess the impact of megafauna extinctions on woody vegetation. Proc. Natl Acad. Sci. USA 113, 847–855 (2016).CAS 

Google Scholar 
Carpio Camargo, A. J. et al. Assessing red deer hunting management in the Iberian Peninsula: the importance of longitudinal studies. PeerJ 9, e10872 (2021).
Google Scholar 
Gordon, I. J., Manning, A. D., Navarro, L. M. & Rouet-Leduc, J. Domestic livestock and rewilding: are they mutually exclusive? Front. Sustain. Food Syst. 5, 550410 (2021).
Google Scholar 
Bond, W. J. Large parts of the world are brown or black: a different view on the ‘Green World’ hypothesis. J. Veg. Sci. 16, 261–266 (2005).
Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Google Scholar 
ESRI. ArcGIS Desktop: Release 10.3. (Environmental Systems Research Institute, Redlands, CA, 2014).Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth. BioScience 51, 933–938 (2001).
Google Scholar 
Fløjgaard, C., Pedersen, P. B. M., Sandom, C. J., Svenning, J.-C. & Ejrnæs, R. Exploring a natural baseline for large-herbivore biomass in ecological restoration. J. Appl. Ecol. 59, 18–24 (2022).
Google Scholar 
Haller, T. et al. Conflicts, security and marginalisation: institutional change of the pastoral commons in a ‘glocal’ world. Rev. Sci. Tech. Off. Int. Epiz. 35, 405–416 (2016).CAS 

Google Scholar 
Torrents-Ticó, M., Fernández-Llamazares, A., Burgas, D. & Cabeza, M. Convergences and divergences between scientific and Indigenous and Local Knowledge contribute to inform carnivore conservation. Ambio 50, 990–1002 (2021).
Google Scholar 
Griffith, E. F., Pius, L., Manzano, P. & Jost, C. C. COVID-19 in pastoral contexts in the greater Horn of Africa: implications and recommendations. Pastoralism 10, 1–12 (2020).
Google Scholar 
Schieltz, J. M. & Rubenstein, D. I. Evidence-based review: positive versus negative effects of livestock grazing on wildlife. What do we really know? Environ. Res. Lett. 11, 113003 (2016).
Google Scholar 
García Sanz, A. La ganadería española entre 1750–1865: los efectos de la reforma agraria liberal. Agricultura y Sociedad 72, 81–120 (1991).
Google Scholar 
San Miguel, A., Roig, S. & Perea, R. The pastures of Spain. Pastos 46, 6–39 (2016).
Google Scholar 
Epp, H. & Dyck, I. Early human-bison population interdependence in the Plains ecosystem. Gt. Plains Res. 12, 323–337 (2002).
Google Scholar 
Hristov, A. N. Historic, pre-European settlement, and present-day contribution of wild ruminants to enteric methane emissions in the United States. J. Anim. Sci 90, 1371–1375 (2012).CAS 

Google Scholar 
Bond, W. J. Ancient grasslands at risk. Science 351, 120–122 (2016).CAS 

Google Scholar 
Bond, W. J., Stevens, N., Midgley, G. F. & Lehmann, C. E. The trouble with trees: afforestation plans for Africa. Trends Ecol. Evol. 34, 963–965 (2019).
Google Scholar 
Ellis, E. C. et al. People have shaped most of terrestrial nature for at least 12,000 years. Proc. Natl Acad. Sci. USA 118, e2023483118 (2021).CAS 

Google Scholar 
Swette, B. & Lambin, E. F. Institutional changes drive land use transitions on rangelands: the case of grazing on public lands in the American West. Glob. Environ. Change 66, 102220 (2021).
Google Scholar 
Hayek, M. N., Harwatt, H., Ripple, W. J. & Mueller, N. D. The carbon opportunity cost of animal-sourced food production on land. Nat. Sustain. 4, 21–24 (2021).
Google Scholar 
Kristensen, J. A., Svenning, J.-C., Georgiou, K. & Mahli, Y. Can large herbivores enhance ecosystem carbon persistence? Trends Ecol. Evol. 37, 117–128 (2022).CAS 

Google Scholar 
Carmona, C. P., Azcárate, F. M., Oteros-Rozas, E., González, J. A. & Peco, B. Assessing the effects of seasonal grazing on holm oak regeneration: Implications for the conservation of Mediterranean dehesas. Biol. Cons. 159, 240–247 (2013).
Google Scholar 
García-Fernández, A. et al. Herbivore corridors sustain genetic footprint in plant populations: a case for Spanish drove roads. PeerJ 7, e7311 (2019).
Google Scholar 
Scoones, I. Living with Uncertainty: New Directions in Pastoralism Development in Africa, Ch. 1 (ITDG, 1995).Pardo, G., Casas, R., del Prado, A. & Manzano, P. Carbon footprint of transhumant sheep farms: accounting for natural baseline emissions in Mediterranean systems. Preprint at Research Square https://doi.org/10.21203/rs.3.rs-1838904/v1 (2022).Odhiambo, M. & Manzano, P. Making Way. Developing National Legal and Policy Frameworks for Pastoral Mobility (FAO, 2022).del Prado, A., Manzano, P. & Pardo, G. The role of the European small ruminant dairy sector on stabilizing global temperatures: lessons from GWP* warming-equivalent emission metrics. J. Dairy Res. 8, 8–15 (2021).
Google Scholar 
Molina-Flores, B., Manzano-Baena, P. & Coulibaly, M. A. The Role of Livestock in Food Security, Poverty Reduction and Wealth Creation in West Africa (FAO, 2020).Lasanta, T., Cortijos-López, M., Errea, M. P., Khorchani, M. & Nadal-Romero, E. An environmental management experience to control wildfires in the mid-mountain Mediterranean area: shrub clearing to generate mosaic landscapes. Land Use Policy 118, 106147 (2022).
Google Scholar 
Torres-Miralles, M. et al. Contribution of High Nature Value farming systems to sustainable livestock production: a case from Finland. Sci. Total Environ. 839, 156267 (2022).CAS 

Google Scholar 
Manzano, P. et al. Towards a holistic understanding of pastoralism. One Earth 4, 651–665 (2021).
Google Scholar 
Karlsson, J. O., Parodi, A., Van Zanten, H. H., Hansson, P. A. & Röös, E. Halting European Union soybean feed imports favours ruminants over pigs and poultry. Nat. Food 2, 38–46 (2021).
Google Scholar 
Leroy, F. et al. Transformation of animal agriculture should be evidence-driven and respectful of livestock’s benefits and contextual aspects. Animal 16, 100644 (2022).
Google Scholar 
Jackson, R. D. Grazed perennial grasslands can match current beef production while contributing to climate mitigation and adaptation. Agric. Environ. Lett. 7, e20059 (2022).
Google Scholar 
Mahli, Y. et al. The role of large wild animals in climate change mitigation and adaptation. Curr. Biol. 32, R181–R196 (2022).
Google Scholar 
O’Bryan, C. J. et al. Unrecognized threat to global soil carbon by a widespread invasive species. Glob. Change Biol. 28, 877–882 (2022).
Google Scholar 
Karp, A. T., Faith, J. T., Marlon, J. R. & Straver, A. C. Global response of fire activity to late Quaternary grazer extinctions. Science 374, 1145–1148 (2021).CAS 

Google Scholar 
Ripple, W. J. et al. World scientists’ warning of a climate emergency 2021. Bioscience 71, 894–898 (2021).
Google Scholar 
Gilbert, M. et al. Global distribution data for cattle, buffaloes, horses, sheep, goats, pigs, chickens and ducks in 2010. Sci. Data 5, 180227 (2018).
Google Scholar 
Mann, D. H., Groves, P., Kunz, M. L., Reanier, R. E. & Gaglioti, B. V. Ice-age megafauna in Arctic Alaska: extinction, invasion, survival. Quat. Sci. Rev. 70, 91–108 (2013).
Google Scholar 
Sandom, C., Faurby, S., Sandel, B. & Svenning, J. C. Global late Quaternary megafauna extinctions linked to humans, not climate change. Proc. Royal Soc. B 281, 20133254 (2014).
Google Scholar 
Fariña, R. A., Czerwonogora, A. & di Giacomo, M. Splendid oddness: revisiting the curious trophic relationships of South American Pleistocene mammals and their abundance. An. Acad. Bras. Ciênc. 86, 311–331 (2014).
Google Scholar  More

EU. Regulation (EU) No 1380/2013 of the European Parliament and of the Council. Off. J. Eur. Union 354 22–61 (2013).Sguotti, C. et al. Catastrophic dynamics limit Atlantic cod recovery. Proc. R. Soc. B Biol. Sci. 286, 20182877 (2019).Article 

Google Scholar 
ICES. in Report of the ICES Advisory Committee 2012. (2012).ICES. Advice basis. in Report of the ICES Advisory Committee 17 (2019).Emeis, K.-C. et al. The North Sea—A shelf sea in the Anthropocene. J. Mar. Syst. 141, 18–33 (2015).Article 

Google Scholar 
Beare, D. et al. An increase in the abundance of anchovies and sardines in the northwestern North Sea since 1995. Glob. Chang. Biol. 10, 1209–1213 (2004).Article 

Google Scholar 
Petitgas, P. et al. Anchovy population expansion in the North Sea. Mar. Ecol. Prog. Ser. 444, 1–13 (2012).Article 

Google Scholar 
Baudron, A. R. et al. Changing fish distributions challenge the effective management of European fisheries. Ecography (Cop.) 43, 494–505 (2020).Article 

Google Scholar 
Sguotti, C. et al. Non-linearity in stock–recruitment relationships of Atlantic cod: Insights from a multi-model approach. ICES J. Mar. Sci. 77, 1492–1502 (2020).Article 

Google Scholar 
Möllmann, C. et al. Tipping point realized in cod fishery. Sci. Rep. 11, 1–12 (2021).Article 

Google Scholar 
Sguotti, C. & Cormon, X. Regime shifts—A global challenge for the sustainable use of our marine resources. in YOUMARES 8—Oceans Across Boundaries: Learning from each other 155–166 (Springer International Publishing, 2018).Beisner, B. E., Haydon, D. T. & Cuddington, K. Alternative stable states in ecology. Front. Ecol. Environ. 1, 376–382 (2003).Article 

Google Scholar 
Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shifts in ecosystems. Nature 413, 591–596 (2001).Article 
CAS 

Google Scholar 
Sguotti, C. et al. Irreversibility of regime shifts in the North Sea. Front. Mar. Sci. 1830, (2022).Yang, Y. & Yamakawa, T. Re-examination of stock-recruitment relationships: A meta-analysis. ICES J. Mar. Sci. 79, 1380–1393 (2022).Article 

Google Scholar 
Perälä, T. A., Swain, D. P. & Kuparinen, A. Examining nonstationarity in the recruitment dynamics of fishes using Bayesian change point analysis. Can. J. Fish. Aquat. Sci. 74, 751–765 (2017).Article 

Google Scholar 
Beaugrand, G. The North Sea regime shift: Evidence, causes, mechanisms and consequences. Prog. Oceanogr. 60, 245–262 (2004).Article 

Google Scholar 
ICES. ICES Stock Assessment Database. https://standardgraphs.ices.dk/stockList.aspx (2019).Szuwalski, C. S., Vert-Pre, K. A., Punt, A. E., Branch, T. A. & Hilborn, R. Examining common assumptions about recruitment: A meta-analysis of recruitment dynamics for worldwide marine fisheries. Fish Fish. 16, 633–648 (2015).Article 

Google Scholar 
Szuwalski, C. S. et al. Global forage fish recruitment dynamics: A comparison of methods, time-variation, and reverse causality. Fish. Res. 214, 56–64 (2019).Article 

Google Scholar 
Otto, S. A., Plonus, R., Funk, S. & Keth, A. INDperform: Evaluation of Indicator Performances for Assessing Ecosystem States. R Package Version 0.1.1. (2018).ICES. Plaice (Pleuronectes platessa) in Subarea 4 (North Sea) and Subdivision 20 (Skagerrak). (2019).Støttrup, J. G., Munk, P., Kodama, M. & Stedmon, C. Changes in distributional patterns of plaice Pleuronectes platessa in the central and eastern North Sea; Do declining nutrient loadings play a role?. J. Sea Res. 127, 164–172 (2017).Article 

Google Scholar 
ICES. Hake (Merluccius merluccius) in subareas 4, 6, and 7, and divisions 3.a, 8.a–b, and 8.d, Northern stock (Greater North Sea, Celtic Seas, and the northern Bay of Biscay). (2019).ICES. Herring (Clupea harengus) in Subarea 4 and divisions 3.a and 7.d, autumn spawners (North Sea, Skagerrak and Kattegat, eastern English Channel). (2019).EU. Council Regulation (EC) No. 811/2004 of 21.4.2004 establishing measures for the recovery of the Northern hake stock. (2004).ICES. Haddock (Melanogrammus aeglefinus) in Subarea 4, Division 6.a, and Subdivision 20 (North Sea, West of Scotland, Skagerrak). (2019).ICES. Cod (Gadus morhua) in Subarea 4, Division 7.d, and Subdivision 20 (North Sea, eastern English Channel, Skagerrak). (2019).Rijnsdorp, A. D. & Vingerhoed, B. Feeding of plaice Pleuronectes platessa L. and sole Solea solea (L.) in relation to the effects of bottom trawling. J. Sea Res. 45, 219–229 (2001).Article 

Google Scholar 
Rijnsdorp, A., Van Damme, C. & Witthames, P. R. Implications of fisheries-induced changes in stock structure and reproductive potential for stock recovery of a sex-dimorphic species, North Sea plaice. ICES J. Mar. Sci. 67, 1931–1938 (2010).Article 

Google Scholar 
Engelhard, G. H., Pinnegar, J. K., Kell, L. T. & Rijnsdorp, A. D. Nine decades of North Sea sole and plaice distribution. ICES J. Mar. Sci. 68, 1090–1104 (2011).Article 

Google Scholar 
Baudron, A. R. & Fernandes, P. G. Adverse consequences of stock recovery: European hake, a new “choke” species under a discard ban?. Fish Fish. 16, 563–575 (2015).Article 

Google Scholar 
Huang, B. et al. Extended reconstructed sea surface temperature version 4 (ERSST.v4). J. Clim. 28, 911–930 (2015).Article 

Google Scholar 
Quante, M. & Colijn, F. North Sea Region Climate Change Assessment (Springer Nature, 2016).Book 

Google Scholar 
IPCC. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. (2014).Erdman, C. & Emerson, J. W. bcp: An R package for performing a Bayesian analysis of change point problems. J. Stat. Softw. 23, 1–13 (2007).Article 

Google Scholar 
Killick, R. & Eckley, I. A. changepoint: An R package for changepoint analysis. J. Stat. Softw. 58, 1–19 (2014).Article 

Google Scholar 
R Core Team. R: A language and environment for statistical computing. R Found. Stat. Comput. (2018).Zeileis, A., Leisch, F., Hornik, K. & Kleiber, C. Strucchange. An R package for testing for structural change in linear regression models. J. Stat. Softw. 7, 1–38 (2002).Article 

Google Scholar 
Ogle, D. H., Wheeler, P. & Dinno, A. FSA: Fisheries Stock Analysis. R package version 0.8.30. (2020).Beverton, R. J. H. & Holt, S. J. On the Dynamics of Exploited Fish Populations (Fisheries Investigations, 1957).
Google Scholar 
Ricker, W. E. Stock and recruitment. J. Fish. Board Can. 11, 559–623 (1954).Article 

Google Scholar 
Muggeo, V. M. R. Segmented: An R package to fit regression models with broken-line relationships. R News 8, 20–25 (2008).
Google Scholar 
James, G., Witten, D., Hastie, T. & Tibshirani, R. An Introduction to Statistical Learning (Springer, New York, 2013).Book 
MATH 

Google Scholar 
Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed effects models and extensions in ecology with R (Springer, 2009).Book 
MATH 

Google Scholar 
Wickham, H. Ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).Book 
MATH 

Google Scholar 
Pedersen, T. L. Patchwork: The Composer of Plots. R package version 1.1.1. (2020). More

Field approachOur searches uncovered 301 papers reporting field studies. After screening the titles abstracts, and full texts, we kept 130 articles for the analysis (Supplementary Fig. 1), from which we obtained 1342 observations regarding 57 Cx. mosquito species from 28 countries and 135 localities (Fig. 1A). Of these 1342 observations, 733 (54.61%) were classified as high quality, (i.e., the number of individuals tested was specified) (Supplementary Tables 1 and 2). The best represented countries were the USA (64.7%, number of observations = 869), Italy (9.3%, n = 125), and Iran (2.9%, n = 39). Based on mosquito field surveillance and individuals testing positive, we concluded that JES is distributed mainly in the Nearctic, Palearctic and Oriental regions (Fig. 1A).Figure 1(A) Weighted Minimum Infection Rates and (B) Weighted Transmission Efficiency of mosquito populations for JES. High-quality data. The size of circles represents the magnitude of the estimates. Map was generated using R software version 4.1.2 with the packages mapdata, maps and tydiverse (https://www.r-project.org) and edited with Inkscape (https://inkscape.org/es/).Full size imageWest Nile virusWNV was detected mainly in the USA (76.5%, number of observations = 826), Italy (4.9%, n = 53) and Iran (3.6%, n = 39) (Fig. 1A). We also recorded 23 species (57%, 41 species) interacting with this virus (Supplementary Tables 1 and 3). Cx. quinquefasciatus became naturally infected in North America [Infection Frequency (IF) = 2.33] (Table 1). We recorded WNV interacting with Cx. tritaeniorhynchus in Asia (IF = 1.02), with Cx. pipiens in Europe (IF = 1.74) and with Cx. antennatus, Cx. neavei, Cx. perexiguus, Cx. perfuscus, Cx. poicilipes, Cx. quinquefasciatus and Cx. tritaeniorhynchus in Africa (IF = 1) (Supplementary Table 3).Table 1 Description of the variables.Full size tableThe highest infection rates were found in North America in Cx. restuans [Standardized minimum infection rate (SMIR) = 56.01], and in Africa and Europe in Cx. pipiens (SMIR = 20.45 and 29.25, respectively). No positive SMIR values were reported in Asian mosquitoes, and Oceanic mosquitoes were not sampled for this virus (Fig. 2A and Supplementary Table 3). The highest infection risk or potential was recorded in species from the USA, such as Cx. restuans (Infection Risk (IR) = 69.50), Cx. pipiens (IR = 55) and Cx. tarsalis (IR = 52.16) (Fig. 3A and Supplementary Table 3). Finally, WNV lineage 1 was detected in Algeria, Turkey, Portugal, Mexico, Tunisia, Iran, Spain and Italy, lineage 2 in Italy, Bulgaria, Greece and the Czech Republic, and lineage 5 in India (Supplementary Table 2).Figure 2(A) Box plots for the Weighted Minimum Infection Rates and (B) Weighted Transmission Efficiency Rates for JES. Boxes indicate 2nd and 3rd quartiles, vertical lines upper and lower quartiles, and horizontal lines the median. Points indicate outliers. The Y axis was transformed to Sqrt (Square root).Full size imageFigure 3JES (A) Infection Risk and (B) Transmission Risk by mosquito species.Full size imageJapanese encephalitis virusJEV was detected mainly in Taiwan (21.6%, n = 24), Korea (18%, n = 20) and Australia (15.3%, n = 17) (Fig. 1A). We found 23 mosquito species interacting with JEV: Cx. vishnui was the one most frequently found to be positive (IF = 1.20), followed by Cx. tritaeniorhynchus (IF = 1.17), Cx. pipiens and Cx. annulus (IF = 0.98) in Asia, while the most susceptible species in Oceania were Cx. sitiens and Cx. gelidus (IF = 0.71) (Supplementary Table 3).The highest SMIR values were recorded in Asia in Cx. rubithoracis (62.38), Cx. annulus (47.68) and Cx. tritaeniorhynchus (28.16) (Fig. 2A and Supplementary Table 3), and Cx. fuscocephala had the highest estimated natural IR (Fig. 3A, Supplementary Table 3). Three genotypes were recorded: genotype I (strain VNKT/479/2007, VNKT/486/2007, and JEV Ishikawa12), genotype III (Tibet-Culex-JEV1-5), and genotype V (K12YJ1174). These were isolated in China, Vietnam, and Japan (Genotype I), Italy, China (Genotype III), and Korea (Genotype V) (Supplementary Table 2).Usutu virusField studies on USUV have been conducted in Europe, Africa, and Asia, most of them in Italy (66.6%, n = 72), Czechoslovakia (11.1%, n = 12) and Slovakia (7.4%, n = 8) (Fig. 1A). Six species were reported to be susceptible to natural infection. Cx. perexiguus had the highest IF and SMIR (1.30) (Supplementary Table 3). In Africa, Cx. antennatus (IF = 1), and in Asia Cx. pipiens (IF = 1) were the most likely to be positive, while Cx. pipiens had the highest IR (5.19) (Fig. 3A, Supplementary Table 3). The recorded strains were USU181_09/USU090-10/USU173_09 (Italy) and USU/Croatia/Zagreb-102/2018 (Italy).St. Louis encephalitis virusThe field studies on SLEV focused on North America (97.7%, n = 43) and Brazil (2.2%, n = 1). Three species were recorded interacting with this virus. Cx. erraticus had the highest IF (2.06), SMIR (2.06) and IR, followed by Cx. quinquefasciatus (North America) (IF = 0.73, SMIR = 1.97) (Fig. 2A, Supplementary Table 3).The highest estimated IR of JES was for Cx. pipiens (Europe), which can be naturally infected with WNV and USUV, followed by Cx. quinquefasciatus (North America), which can be infected with WNV and SLEV (Fig. 3A).Experimental approachExperimental studies were reported in 481 articles. After screening the titles, abstracts, and full texts, as well as opportunistic records, 95 articles remained for the analysis (Supplementary Fig. 2). From these we obtained 189 high quality observations of the TE of JES in 11 countries, 40 localities, and 21 species (Fig. 1B, Supplementary Table 1). The USA was the best represented country (54.4%, n = 103), followed by Germany (13.2%, n = 25) and Australia (12.6%, n = 24). There was, however, a notable lack of information on the vector competence of Cx. mosquitos for JES in many regions of the world, such as Central and South America, and Africa (Fig. 1B).The most common means of infection was oral (94.8%, 395 observations), while the rest were intrathoracic. Intrathoracic infection bypasses the midgut barrier so is not considered natural infection. We therefore carried out the subsequent analyses using only the data on oral infection (Supplementary Table 4).We used a generalised linear model (GLM) for the statistical analysis, which was conducted only on the WNV dataset (strain NY99), the only one with sufficient observations for the purpose (n = 63). We did not find a significant effect of viral titre, temperature, or days post infection on TE. However, more data with a wide range of values is necessary to confirm these observations. On the other hand, we found that the Extrinsic Incubation Period (as DPI) was shorter at higher temperatures (Fig. 4 and Supplementary Table 5).Figure 4Relationship between temperature and Days Post Infection for WNV strain NY99.Full size imageWest Nile virusMosquito populations from many locations on all continents have been studied for their vector competence for this virus, particularly in the USA (60.3% of observations, n = 96), Germany (15.7%, n = 25) and Australia (6.9%) (Fig. 1B). Our bibliographic research revealed 21 species of Cx. with the ability to transmit WNV under laboratory conditions (Supplementary Table 6). Cx. pipiens (North America) and Cx. tarsalis were the most frequently studied species and were the most efficient in transmitting the virus (Transmission Frequency (TF) = 2.33) (Table1). Cx. quinquefasciatus had the highest TF (1.70) in Africa, Cx. modestus in Europe (TF = 1.32), and Cx. annulirostris and Cx. quinquefasciatus in Oceania (TF = 1.48) (Supplementary Table 6).Concerning Standardized Transmission Rates (STE) estimates, Cx. quinquefasciatus had the highest values in the USA (STE = 1.63), Cx. pipiens in Europe (0.90), Cx. tritaeniorhynchus in Asia (1.8), Cx. neavei in Africa (0.17) and Cx. annulirostris in Oceania (2.45) (Fig. 2B, Supplementary Table 6). We found 20 different strains of WNV tested. The TE of the various WNV strains vary considerably, but lineage 1 was more efficient than lineage 2. There were also more studies on the lineage 1 strains (n = 11), which exhibited high variation (Fig. 5).Figure 5Box plots for WNV (A) lineages and (B) strains used to measure Weighted Transmission Rates.Full size imageJapanese encephalitis virusJEV has been studied mainly in mosquito populations from France (45%, n = 20) and Australia (34%, n = 15), but also the United Kingdom, India, Taiwan, New Zealand, and the USA (Fig. 1B). Six mosquito species are capable of transmitting JEV. Cx. pipiens (Europe) had the highest TF (1.85), while Cx. gelidus had high values of STE (1.73) (Fig. 3B and Supplementary Table 6).St. Louis encephalitis virusVector competence for SLEV has been studied in two countries: the USA (93.3%, n = 42) and Argentina (6.6%, n = 3), and 7 mosquito species have been investigated. Cx. nigripalpus was the most efficient in transmitting the virus (TF = 1.60), while Cx. pipiens had the highest STE (0.68) (Figs. 2B, 3B and Supplementary Table 6).Usutu virusStudies have also been conducted on the Usutu virus in mosquito populations in the USA (28.57%, n = 4), the United Kingdom (42.8%, n = 6) and Senegal (25%, n = 4), in particular on Cx. neavei, Cx. pipiens and Cx. quinquefasciatus (TF = 1). Cx. neavei also had the highest STE (0.79) (Fig. 3B).We found reports of JES transmission under laboratory conditions in 22 Cx. species, and natural infections in 32 species (55.1% of the total sample) in the field. Cx. pipiens complex (biotypes quinquefasciatus, pipiens, molestus and pallens) was the most common vector accounting for 36.9% (n = 660) of the experimental observations and 25.7% (n = 1342) of the field observations. With both approaches, WNV was the most common flavivirus, accounting for 80.4% of the field observations and 86.7% of the experimental data (Fig. 1A,B). Only WNV, therefore, had enough observations to make comparison between the experimental and field data possible. We were able to compare 16 mosquito species and found a high positive correlation between TF and IF (R = 0.57, p = 0.02) (Fig. 7).In summary, we found that the species with the highest infection-transmission risk (IRT) for WNV was Cx. restuans, for USUV it was Cx. pipiens (Europe), for SLEV Cx. quinquefasciatus (North America), and for JEV Cx. gelidus (Oceania) (Fig. 6 and Supplementary Tables 2 and 6).Figure 6JES infection-transmission risk by continent and flavivirus.Full size image More