More stories

  • in

    Urban blue–green space landscape ecological health assessment based on the integration of pattern, process, function and sustainability

    Study areaHarbin is located in the centre of Northeast Asia, between 44°04’–46° 40′ N and 125° 42′–130° 10′ E24,26. The site has a mid-temperate continental monsoon climate, with an average annual temperature of 3.6° C and an average annual precipitation is 569.1 mm. The main precipitation months being from June to September, accounting for about 60% of the annual precipitation, the main snow months are from November to January24,25. The overall topography is high in the east and low in the west, with mountains and hills predominating in the east and plains predominating in the west27. In this study, we identified the central district of Harbin, where urban construction activities are frequent and the population is dense, as the study area. According to the “Harbin City Urban Master Plan (2011–2020)” (revised draft in 2017), the specific scope includes Daoli District, Daowai District, Nangang District, Xiangfang District, Pingfang District, Songbei District’s administrative district, Hulan District, and Acheng District part of the area, with a total area of 4187 km2 (Fig. 2). The blue–green space in this study included woodland, grassland, cultivated land, wetland and water that permeate inside and outside the construction sites. They all have integrated functions such as ecology, supply, beautification, culture, and disaster prevention and avoidance, and have a decisive influence on the urban ecological environment.Figure 2Schematic of study area. The Figure is created using ArcGIS ver.10.2 (https://www.esri.com/).Full size imageData sourcesThe data used in this research included the following: land-cover date (30 m × 30 m) of two periods (2011, 2020) spported by the China Geographic National Conditions Data Cloud Platform (http://www.dsac.cn/), Meteorological datasets (1 km × 1 km) were obtained from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http:∥www.resdc.cn/), including air temperature, precipitation, and surface runoff. ASTER GDFM elevation data (30 m × 30 m) came from the Geospatial Data Cloud (http:∥www.gscloud.cn), from which the slope was extracted. Soil data (1 km × 1 km) were from the World Soil Database (HWSD) China Soil Data Set (v1.1). The normalized difference vegetation index (NDVI) and modified normalized difference water index (MNDWI) data (30 m × 30 m) came from the National Comprehensive Earth Observation Data Sharing Platform (http://www.chinageoss.org/), ET datasets (30 m × 30 m) were drawn from the NASA-USGS (https://lpdaac.usgs.gov/). Social and economic data were mainly obtained through the Harbin statistical yearbook and the Harbin social and economic bulletin.Framework of urban blue–green space LEH assessmentUrban blue–green space is a politically defined man-land coupling region composed of ecological, economic, and social systems, which is greatly disturbed by human activities11. The essence of urban blue–green space LEH is that the landscape ecological function sustainably meets human needs28,29. The landscape ecological function reflects the value orientation of human beings to blue–green space, and to a large extent affects the blue–green landscape ecological pattern and process. The interaction between the blue–green landscape ecological pattern and process drives the overall dynamics of blue–green space. Meanwhile, presenting certain landscape ecological function characteristics, which provide ecological support for various human activities30,31,32. While the pattern and process of blue–green space both profoundly influence and are influenced by human activities33,34. This influence is long-term, the standard of LEH should not be fixed in real-time health, but should fully consider the sustainability of the health state.In summary, the landscape ecological pattern, process, function, and sustainability are not separate, but a complex of mutual integration, and organic unity. In this study, we constructed an integrated assessment framework of blue–green space LEH that included four units: pattern, process, service, and sustainability (Fig. 3). In the assessment framework, the LEH of urban blue–green space involves two dimensions: the first is the health status of the urban blue–green space itself, emphasizing the maintenance of the ecological conditions, thereby potentially satisfying a series of diversity goals. The other is that urban blue–green space, as a part of social and economic development, could sustainably provide the ability to meet (subject) needs and goals.Figure 3Key units, interactions of urban blue–green space LEH.Full size imageLandscape ecological patternThe landscape ecological pattern of urban blue–green space is a spatial mosaic combination of landscape elements at different levels or the same level. Affected by human activities interference31, the landscape ecological pattern shows the changing trend of landscape structure complexity, landscape type diversification, and landscape fragmentation. The assessment of urban landscape ecological pattern should be a comprehensive reflection of this changing trend1. Landscape pattern indexes are the most frequently applied which could reflect the structural composition and spatial configuration characteristics of the landscape4,35. This study took landscape ecology as the entry point and selected the landscape pattern indexes that can quantitatively reflect the change characteristics of landscape structural composition and spatial configuration under the disturbance. In this way, the landscape disturbance index (U), landscape connectivity index (CON), and landscape adaptability index (LAI) were used as the indexes for the assessment of landscape ecological pattern health.

    (1)

    Landscape disturbance index (U)

    There are two kinds of relationships between the landscape ecological pattern and the external disturbance: compatibility and conflict. As the landscape ecological pattern has accommodating characteristics, the disturbance beyond the accommodating capacity will degrade the landscape ecological pattern36,37. The landscape disturbance index (U) could characterize the degree of fragmentation, dispersion, and morphological changes in landscape pattern38. The index is a comprehensive index that can reflect the health of the landscape pattern by quantifying the ability of ecosystems to accommodate external disturbances. It consists of the landscape fragmentation index, the inverse of the fractional dimension, and the dominance index. They measure the response of the landscape pattern to external disturbance from the perspective of different landscape types, the same landscape type, and landscape diversity, respectively36,38, and their weights were determined by the entropy weight method. The formula is as follows:$$ U = alpha N_{{{Fi}}} + bD_{{{Fi}}} + cD_{{{Oi}}} $$
    (1)
    where NFi is the landscape fragmentation index, DFi is the inverse of the fractional dimension, DOi is the dominance index, and a, b, and c are the corresponding weights, which were 0.20, 0.5, and 0.3 in this study, respectively.

    (2)

    Landscape connectivity index (CON)

    The most direct result of landscape ecological pattern degradation caused by external disturbance is that the flow of energy, material, and information among ecological patches is reduced or even blocked, ultimately the stability of the landscape pattern is decreased. The connectivity could characterize the ability of landscape ecological pattern to mitigate risk transmission, which is significant for the dynamic stability of landscape ecological pattern39,40. The landscape connectivity index (CON) could measure the connectivity between ecosystem components through the aggregation or dispersion trend of patches41. The better the connectivity, the stronger the stability of landscape ecological pattern. The formula is as follows:$$ CON = frac{{100sumlimits_{s = 1}^{q} {sumlimits_{h ne l}^{p} {C_{{{shl}}} } } }}{{sumlimits_{s = 1}^{s} {left[ {q_{{s}} (q_{{s}} – 1)/2} right]} }} $$
    (2)
    where qs is the number of plaques of patch type s, Cshl is the link between patch h and patch l in s within the delimited distance.

    (3)

    Landscape Restorability Index (LRI)

    The ability to recover to its original structure when subjected to disturbances is an important criterion for the landscape ecological pattern42. Research confirmed that the restorability of the landscape ecological pattern is closely related to the structure, function, diversity, and uniformity of distribution. The landscape restorability index (LRI) combines the above landscape information and could indicate the restorability of the landscape ecological pattern in response to disturbance43. The index consists of the patch density, Shannon diversity index, and the landscape evenness, the patch density is the number of patches per square kilometer. The Shannon diversity index reflects the change in the proportion of landscape types. The landscape evenness index shows the distribution evenness of patches in terms of area. The larger the LRI index, the more complex and evenly distributed the structure is, and the more recovery ability of the landscape pattern against disturbance is. The formula is as follows:$$ LRI = PD times SHDI times SHEI $$
    (3)
    where PD is the patch density, SHDI is the Shannon diversity index, and SHEI is the landscape evenness index.Landscape ecological processThe landscape ecological process of urban blue–green space is extremely complex for it involves multiple factors such as natural ecology, economy, and culture. Landscape ecological process assessment is the measure of the self-organized capacity and the efficiency of ecological processes within and among patches44. A blue–green space with a healthy landscape ecological process should have the ability to adapt to conventional land use under human management and maintain physiological integrity while maintaining the balance of ecological components. Specifically, the landscape ecological process could quickly restore its balance after severe disturbances, with strong organization, suitability, recoverability, and low sensitivity45,46. A single model hardly to gets good research on landscape ecological process under the urban scale. The comprehensive application of multidisciplinary methods is effective means to solve the problem. Regarding this, we selected ecological indexes and models from four aspects: organization, suitability, restoration, and sensitivity to assess the landscape ecological process of urban blue–green space.

    (1)

    Organization index (O)

    The organization of the landscape ecological process is the maintenance ability of stable and orderly material cycling and energy flow within and between landscapes47. The normalized vegetation index (NDVI) and the modified normalized difference water index (MNDWI) could reflect the efficiency and order of ecological processes. Such as accumulation of organic matter, fixation of solar energy, nutrient cycling, regeneration, and metabolism13. The indexes are the external performance of the internal dynamics and organizational capabilities of the ecological process. In recent years, it has been widely used in the assessment of related to landscape ecological process. The formulas are as follows:$$ NDVI = frac{NIR – R}{{NIR + R}} $$$$ MNDWI = frac{p(green) – p(MIR)}{{p(green) + p(MIR)}} $$
    (4)
    where (NDVI) is the normalized vegetation index, (MNDWI) is the modified water body index, (NIR) is the reflectance value in the near-infrared band, (R) is the reflectance value in the visible channel, (p(green)) and (p(MIR)) are the normalized values in the green and mid-infrared bands.

    (2)

    Suitability index (Q)

    The suitability of the landscape ecological process is a measurement of the self-regulating ability of the landscape ecosystem. That is, to effectively maintain the ecological process in a state of being protected from disturbance during the occasional changes caused by the external environment2. The water conservation amount index (Q) can measure the operating capacity of ecosystems to maintain ecological balance, water conservation, climate regulation, and other ecological processes by integrating the water balance of rainfall, surface runoff, and evaporation41. It could reflect the suitability of landscape ecological process to regional environment and developmental conditions. The formula is as follows:$$ Q = R – J – ET $$
    (5)
    where Q is the water conservation amount, R is the annual rainfall, J is the surface runoff, ET is the evapotranspiration.

    (3)

    Recoverability index (ECO)

    The recoverability of the landscape ecological process refers to the ability of an ecosystem to return to its original operating state after being subjected to external impacts. Land-use types play an essential role in landscape ecological recoverability48. The ecological recoverability index (ECO) uses the resilience coefficients of land-use types to reflect the level of ecosystem resilience38. Based on previous studies, the resilience coefficient of land-use types was assigned (Table 1).

    (4)

    Sensitivity index(A)

    Table 1 Resilience coefficients of different land use types.Full size tableThe sensitivity index (A) could be used to indicate landscape ecological process formation, change, and vulnerability to disturbance31. We started from the physical effects of blue–green space on sand production, water confluence, and sediment transport, introduced the Soil Erosion Modulus to characterize the sensitivity of landscape ecological processes to disturbance. The index effectively combines landscape ecology, erosion mechanics, soil science, and sediment dynamics49. The formula is as follows:$$ begin{gathered} A = R_{{i}} cdot K cdot LS cdot C cdot P hfill \ L = (l/22.1)^{m} hfill \ S = left{ begin{gathered} 10.8sin theta + 0.03,theta < 5^{ circ } hfill \ 16.8sin theta - 0.50,5^{ circ } le theta < 10^{ circ } hfill \ 21.9sin theta - 0.96,theta ge 10^{ circ } hfill \ end{gathered} right. hfill \ C = left{ begin{gathered} 1,c = 0 hfill \ 0.6508 - 0.3436lg c,0 < c le 78.3% hfill \ 0,c > 78.3% hfill \ end{gathered} right. hfill \ end{gathered} $$
    (6)
    where A is the soil erosion modulus. Ri is the rainfall erosion factor, K is the soil erosion factor, L and S are slope the length factor and the slope factor respectively, C is the vegetation coverage and management factor, P is the soil and water conservation factor, l is the slope length value, m is the slope length index, and θ the is slope value.Landscape ecological functionThe landscape ecological function determines the ability of ecological service50,51,52, the ecological service of urban blue–green space depends on the human value orientation48. It includes four categories: supply, support, regulation, and culture. Based on Maslow’s Hierarchy of Needs and Alderfer’s ERG theory, scholars have summarized the three major needs of human beings for urban blue–green space. Namely, securing the living environment to meet the survival needs, improving social relationships to meet the interaction needs, and cultivating cultural cultivation to meet the development needs53. Specifically corresponding to the landscape ecological function of urban blue–green space, supply is not the main function, only plays a subsidiary role, support is the basic guarantee, regulation is the basic need for urban environmental construction, and culture is an important element of high-quality social life. Ecosystem service value (ESV) can realize the measurement of ecological service function by calculating the specific value of life support products and services produced by the ecosystem54,55,56. Considering the human value orientation of the urban blue–green space landscape ecological function, the weights were given by consulting 16 experts, with supply, regulation, support, and culture weights of 0.2, 0.3, 0.3, 0.2, respectively. The formula is as follows:$$ ESV = sumlimits_{k = 1}^{n} {S_{k} times V_{k}^{{}} } $$
    (7)
    where Sk is the area of landscape type k, Vk is the value coefficient of the ecosystem service function of landscape type k .Landscape ecological sustainabilityWu (2013) proposed a research framework for landscape sustainability based on a summary of related studies, stating that landscape ecological sustainability is the ability to provide ecosystem services in a long-term and stable manner34. The framework emphasized that landscape sustainability should focus on the analysis of ecosystem service trade-offs effect34,57. In the process of dynamic change of urban blue–green space ecosystem, there are complex trade-offs among various ecosystem services. This is important for promoting the optimal overall benefits of various ecosystem services and achieving sustainable development of urban ecology58. In addition, as a special type of human-centered ecosystem developed by humans based on nature, human well-being is also very important for the landscape ecological sustainability of urban blue–green space. For this reason, we introduced ecosystem service trade-offs (EST) and ecological construction input (IEC) as assessment indexes of landscape ecological sustainability.

    (1)

    Ecosystem service trade-offs (EST)

    This study applied the root mean square deviation of ecological services to quantify ecosystem service trade-offs (EST). The index could effectively measure the average difference in standard deviation between individual ecosystem services and the average ecosystem services. It is a simple and effective way to evaluate the trade-offs among ecosystem services. The formula is as follows:$$ EST = sqrt {frac{1}{n – 1}sumnolimits_{i = 1}^{n} {(ES_{std} – overline{ES}_{std} } } )^{2} $$
    (8)
    where ESstd is the normalized ecosystem services, n is the number of ecosystem services , and (overline{ES}_{std}) is the mean value of normalized ecosystem services.

    (2)

    Ecological construction input (ECI)

    Human well-being is a premise for the landscape ecological sustainability of urban blue–green spaces, it is closely related to government investment in ecological construction planning34. From the perspective of economics, this study assessed the human well-being obtained by urban blue–green space with the ratio of urban ecological construction investment to GDP, that is, the ecological construction input (ECI). The formula is as follows:$$ ECI = EI/G $$
    (9)
    where EI is the amount of ecological construction investment, and G is the gross regional product.Evaluation methodThe index weight determines its relative importance in the index system, and the selection of the weight calculation method in the decision-making of multi-attribute problems has an important impact on the assessment results21. Traditional weighting methods can be divided into two categories, subjective weighting method and objective weighting method21,38. The subjective weighting method is represented by the analytic hierarchy process (AHP), Delphi method, and so on. It has the advantage of simplicity, but the disadvantage is too subjective and randomness because it was completely dependent on the knowledge and experience of decision makers. The objective weighting method is represented by the entropy weighting method (EWM), principal component analysis, variation coefficient method, and so on. And it has been widely recognized for reflecting the variability of assessment results18, but the values of indexes have significant influence and the calculation results are not stable. Considering the limitations of the single weighting method, the weights of each assessment index in this study were determined by the combination of subjective weight and objective weight. Among them, the subjective weighting selected the AHP, and the objective weighting selected the EWM (Table 2). The formula is as follows:$$ w_{{j}} = alpha w_{{j}}^{{{AHP}}} + (1 – alpha )w_{{j}}^{{{EWM}}} $$
    (10)
    $$ w_{{j}}^{{{EWM}}} = d_{{j}} /sumlimits_{i = 1}^{m} {d_{{j}} } $$
    (11)
    $$ d_{{j}} = 1 – e_{{j}} $$
    (12)
    $$ e_{{j}} = – ksumlimits_{i = 1}^{n} {f_{{{ij}}} ln (f_{{{ij}}} )} ,;k = 1/ln (n) $$
    (13)
    $$ f_{{{ij}}} = X^{prime}_{{{ij}}} /sumlimits_{i = 1}^{n} {X^{prime}_{{{ij}}} } $$
    (14)
    where (W_{{j}}^{{}}) is the combined weight. (W_{{j}}^{{_{AHP} }}) is the weight of the j-th index of the AHP, (W_{{j}}^{{{EWM}}}) is the weight of the j-th index of the EWM, dj is the information entropy of the j-th index, ej is the entropy value of the j-th index, (f_{{{ij}}}) is the proportion of the index value of the j-th sample under the i-th indexm, (X^{prime}_{{{ij}}}) is the standardized value of the i-th sample of the j-th index, m is the number of index, n is the number of samples, and (alpha) was taken as 0.5.Table 2 Weight of assessment index.Full size tableSince the dimensions of indexes are different, it is necessary to unify the dimensions of the index to avoid the errors caused by direct calculation to make the evaluation results inaccurate. The range standardization was used to normalize the index data and bound its value in the interval [0, 1], the range standardization can be expressed as follows15,23:$$ {text{Positive indicator}}left( + right):A_{{{ij}}} = (X_{{{ij}}} – X_{{{jmin}}} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
    (15)
    $$ {text{Negative indicator}}left( – right):A_{{{ij}}} = (X_{{{jmax}}} – X_{ij} )/(X_{{{jmax}}} – X_{{{jmin}}} ) $$
    (16)
    Additionally, we divided the LEH index into five levels from high to low using an equal-interval approach as follows40: [1–0.8) healthy, [0.8–0.6) sub-healthy, [0.6–0.4) moderately healthy, [0.4–0.2) unhealthy, [0.2–0] pathological, corresponding level I–V. And the level transfer of LEH in different periods was divided into three types: improvement type, degradation type, and stabilization type. For example, III-II means that the transfer from level III to level II is the improvement type.Spatial autocorrelation analysisSpatial autocorrelation analysis is one of the basic methods in theoretical geography. It could deeply investigate the spatial correlation characteristics of data, including global spatial autocorrelation and local spatial autocorrelation23. The global spatial autocorrelation uses global Moran’s I to evaluate the degree of their spatial agglomeration or differentiation of an attribute value in the study area. The local spatial autocorrelation is a decomposed form of the global spatial autocorrelation18,21, including four types: HH(High-High), LL(Low-Low), HL(High-Low), LH(Low–High). In this study, spatial autocorrelation analysis was applied to study the spatial correlation characteristics of blue–green space LEH. The calculation formulas are as follows:$$ I = frac{{Nsumlimits_{i} {sumlimits_{v} {W_{iv} (Y_{i} – overline{Y} )(Y_{v} – overline{Y} )} } }}{{(sumlimits_{i} {sumlimits_{v} {W_{iv} } } )sumlimits_{i} {(Y_{i} – overline{Y} )} }} $$
    (17)
    $$ I_{i} = frac{{Y_{i} – overline{Y} }}{{S_{x}^{2} }}sumlimits_{v} {left[ {W_{iv} (Y_{i} – overline{Y} )} right]} $$
    (18)
    where N is the number of space units, (W_{iv}) is the spatial weight, (Y_{i} ,Y_{v}) are the variable attribute values of the area (i,v), (overline{Y}) is the variable mean, (S_{x}^{2}) is the variance, (I) is the global Moran’s I index, and (I_{i}) is the local Moran’s I index. More

  • in

    Bacterial matrix metalloproteases and serine proteases contribute to the extra-host inactivation of enteroviruses in lake water

    Virus propagation and enumerationEchovirus-11 (E11, Gregory strain, ATCC VR737) and Coxsackievirus-A9 (CVA9, environmental strain from sewage, kindly provided by the Finnish National Institute for Health and Welfare) stocks were produced by infecting sub-confluent monolayers of BGMK cells as described previously [7]. Viruses were released from infected cells by freezing and thawing the culture flasks three times. To eliminate cell debris, the suspensions were centrifuged at 3000 × g for 5 min. Each stock solution was stored at −20 °C until use. Infectious virus concentrations were enumerated by a most probable number (MPN) infectivity assay as described in the Supplementary Information. The assay limit of detection (LoD), defined as the concentration corresponding to one positive cytopathic effect in the lowest dilution of the MPN assay under the experimental conditions used, corresponding to 2 MPN/mL.Inactivation of enteroviruses by bacterial consortia from lake waterTo study the inactivation of CVA9 and E11 by a bacterial consortium from lake water, four surface water samples were collected from Lake Geneva (Ecublens, Switzerland) during the summer 2021. Each sampling event was conducted on warm and sunny days, to minimize biological variation. Immediately after sampling, large particles of the sample were removed by filtering 500 mL of water on a 8 μm nitrocellulose filter membrane (Merck Millipore, Cork, Ireland). The sample was then filtered through a 0.8 μm nitrocellulose filter membrane (Merck Millipore) to remove large microorganisms such as protists. The resulting water sample corresponds to the bacterial fraction used to study virus inactivation.For inactivation experiments, each virus was spiked into individual 1 mL aliquots of fractionated lake water to a final concentration of 106 MPN/mL, and samples were incubated for 48 h at 30 °C without shaking. Duplicate experiments were conducted for each virus and each lake water sample. Experiments to control for thermal inactivation were conducted using the same procedure but by replacing the fractionated lake water with sterile milliQ water. Viral infectivity at times 0 h and 48 h was determined by MPN as described above. Virus decay was calculated as log10 (C/C0), where C is the residual titer after 48 h of incubation, and C0 is the initial titer. The experimental LoD was approximately 5-log10.These same experiments were conducted for three new water samples in the presence of four protease inhibitors with the following final concentrations: E64—10 μM (E3132, Sigma–Aldrich, Saint-Louis, MO, USA), GM6001—4 μM (CC1010, Sigma–Aldrich), Chymostatin—100 μM (C7268, Sigma–Aldrich), and PMSF—100 μM (P7626, Sigma–Aldrich). Each inhibitor was added to 1 mL of fractionated lake water, vortexed for 30 seconds, and incubated at room temperature for 15 min, before adding the two viral strains under the same conditions as described above.Bacterial isolation, cultivation, and storageBacteria were isolated from two water samples from Lake Geneva’s Ecublens beach, taken in November 2019 (Fall, 89 isolates) and May 2020 (Spring, 47 isolates). Bacteria recovery was performed on R2A agar plate (BD Difco, Franklin Lakes, NJ, USA) as described previously [15]. Briefly, successive dilutions from 10−1 to 10−5 were carried out in sterile water for each sample. For each dilution, a volume of 1 mL was deposited on three separate R2A plates, before being incubated at 22, 30, and 37 °C. After 5 days of incubation, each colony was picked and enriched on a new R2A plate. To ensure purity, each isolate was successively plated five times on R2A plate and incubated at the same temperature as the initial isolation. Each purified isolate was cryopreserved in R2A / 20% glycerol at −80 °C. The isolates were named based on the water body (Lake (L)), isolation temperature, and the isolation order (L-T°C-number).Bacterial identificationThe identification of each isolate was performed by 16 S rRNA gene sequencing using the pair of primers 27 F (5’- AGA GTT TGA TCM TGG CTC AG- 3’, Microsynth AG, Balgach, Switzerland) / 786 R (5’- CTA CCA GGG TAT CTA ATC – 3’, Microsynth AG), following a methodology previously described [15]. The thermocycling conditions and the purification of PCR products are described in the Supplementary Information. The complete list of isolated bacteria and associated accession numbers is given in Supplementary Table 1.Phylogenetic inference and metadata visualizationThe consensus from 16 S rRNA gene sequences of the 136 isolates was aligned using the MUSCLE algorithm [16]. The phylogenetic analysis of 566 bp aligned sequences from the V2-V4 16 S rRNA gene regions (Positions: 152–717) was performed using Molecular Evolutionary Genetics Analysis X software [17]. Phylogeny was inferred by maximum likelihood, with 1000 bootstrap iterations to test the robustness of the nodes. The resulting tree was uploaded and formatted using iTOL [18].Virus incubation with bacterial isolatesFor the preparation of the bacteria before co-incubation, each one was first cultured on R2A agar for 48 h at their initial isolation temperature. Overnight suspensions of each bacterial isolate were grown in R2A broth at room temperature under constant agitation (180 rpm). For co-incubation experiments, 200 μL of each bacterial suspension were mixed with 100 μL of a 105 MPN/mL stock of E11 or CVA9. Then, each condition was supplemented with 600 μL of R2A broth. Incubation was carried out for 96 h at room temperature, without shaking. At the end of the co-incubation, each tube was centrifuged for 15 min at 9000 × g (4 °C) to eliminate bacteria, and the residual infectious viral titer was enumerated by MPN assay as described above [7]. Each co-incubation experiment was carried out in triplicate. Control experiments were performed under the same conditions but using sterile R2A. Virus decay was quantified as log10 (Cexp/Cctrl), where Cexp is the residual titer after a co-incubation for 96 h, and Cctrl is the titer after incubation of the virus in sterile R2A for 96 h. The experimental LoD was 3-log10.Protease activity measurement using casein and gelatin agar platesCasein agar was prepared as follows: 20 g of skim milk (BD Difco), supplemented with 1 g glucose were reconstituted with 200 mL of distilled water. Likewise, a 10% bacteriological agar solution was prepared in a final volume of 200 mL. Finally, a solution consisting of 0.8% NaCl, 0.02% KCl, 0.144% Na2HPO4, and 0.024% KH2PO4 was reconstituted in 600 mL of water. All solutions were autoclaved for 15 min at 110 °C. The solutions were mixed, and 25 mL were poured into each Petri dish. Gelatin agar was composed of 0.4% peptone, 0.1% yeast extract, 1.5% gelatin and 1.5% bacteriological agar. The mixture was autoclaved 15 min at 120 °C, and 25 mL of medium was poured into each Petri dish.For each isolate, an overnight suspension was performed in R2A broth at room temperature, before spotting 15 μL of each suspension at the center of both gelatin and casein agar plates. Each plate was incubated at 22, 30, or 37 °C for 72 h, depending on the initial isolation temperature of the bacteria. Casein-degrading activity (cas), which is exerted by many different protease classes, and gelatin-degrading activity (gel), which is mostly caused by MMPs, were revealed by a hydrolysis halo around the producing bacteria. Hydrolysis diameters were measured in millimeters (mm) to report the extent of the proteolytic effect of each strain on both substrates.Protease activity quantification in cell-free supernatantUsing the same bacterial suspensions as for bacterial/virus co-incubation, 200 μL of each suspension was inoculated into 600 μL of R2A broth and incubated without shaking for 96 h at room temperature. Each culture was centrifuged for 15 min at 9000 × g at 4 °C. The resulting cell-free supernatants (CFS) were stored at −20 °C until use. For each CFS, protease activity was measured using the Protease Activity Assay Kit (ab112152, Abcam, Cambridge, UK), which measures general protease activity (pgen) except MMPs, and the MMP Activity Assay Kit (ab112146, Abcam), which selectively measures MMP activity (mmp). Briefly, for the Protease Activity Assay kit, 50 μL of the substrate was added into each well of a dark-bottom plate containing 50 μL of each CFS. Standard trypsin provided by the kit was used as a positive control. For the MMP Activity Assay kit, 50 μL of each CFS was incubated with 50 μL of 2 mM APMA for 3 h at 37 °C, prior to the activity test. Collagenase I (C0130, Sigma–Aldrich) was used as a positive control. R2A broth was used as a negative control for each assay. Protease activity was measured at time 0 and after 60 min, using a Synergy MX fluorescence reader (BioTek). The excitation and emission wavelengths were set to 485 and 530 nm, respectively. The emitted fluorescence, generated by proteolytic cleavage of the substrate of each kit, was calculated as follows: ∆RFU = RFU (60 min) − RFU (0 min). Proteolytic activity was calculated in mmol/min/μL based on the emitted fluorescence measured for trypsin and collagenase I at known proteolytic activities.Data analysisStatistical analyses to compare inactivation data were performed by one-way t-test or one-way ANOVA with Dunnett’s post-hoc test in GraphPad Prism v.9. An alpha value of 0.05 was used as a threshold for statistical significance. For each dataset we confirmed that data were normally distributed.To analyze a potential correlation between protease activity and viral decay, the decay values for each virus strain was related to the four protease activity tests of this study using a scatterplot combined with a Kernel density estimation. The analyses were performed with R v.3.6.1 using the SmoothScatter function of the R Base package.A Left-Censored Tobit model (CTM) with mixed effects was chosen to investigate interactions between protease activity and the decay measured for each virus strain. Briefly, the CTM with mixed effect was chosen for three reasons: (1) The protocol used to measure viral decay had a limit of quantification of −3-log10, and 152 measurement points reached the detection limit, requiring the use of this value as the left-censored value of the model; (2) The two virus strains used in the study showed distinct responses after exposure to environmental bacteria, preventing the use of a multiple linear regression model; (3) Among biological replicates of co-incubation experiments, inactivation variability was observed, suggesting the concomitant action of random biological effects (e.g., production of other compounds than proteases by bacteria, or differences in protease production rate between replicates for each bacterial isolate). The resulting statistical model was then formulated as follows:$$log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = ; beta _0 + beta _1;{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2} + beta _2sqrt {left[ {pgen} right]_i} + beta _3sqrt {left[ {mmp} right]_i} + beta _4sqrt {left[ {cas} right]_i} \ + beta _5sqrt {left[ {gel} right]_i} + beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {pgen} right]_i} + beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {mmp} right]_i} \ + beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {cas} right]_i} + beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}sqrt {left[ {gel} right]_i} + alpha _{{{{{{{{mathrm{id}}}}}}}}_i} + varepsilon _i$$$${{{mbox{where}}}}; log left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) = left{ {begin{array}{*{20}{c}} { – 3} & {{{{{{{{mathrm{if}}}}}}}};{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right) le – 3} \ {{{{{{{{mathrm{log}}}}}}}}left( {frac{{C_{{{{{{mathrm{exp}}}}}}}}}{{C_{{{{{{mathrm{ctrl}}}}}}}}}} right)} & {{{{{{{{mathrm{otherwise}}}}}}}}} end{array}} right.$$$$alpha _{{{{{{{{mathrm{id}}}}}}}}_i}sim {{{{{{{mathrm{i}}}}}}}}.{{{{{{{mathrm{i}}}}}}}}.;{{{{{{{mathrm{d}}}}}}}}.;{rm N}left( {0,;sigma _{{{{{{{{mathrm{id}}}}}}}}}^2} right)$$$${{{{{{{mathrm{for}}}}}}}};i in left{ {1,2, ldots } right}$$for which β0 defines the model intercept, (beta _1{rm I}_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the main effect of the virus factor on the viral decay, (beta _2,;beta _3,;beta _4,;{{{{{{{mathrm{and}}}}}}}};beta _5) corresponds to the main effects of the different protease activity measurements on viral decay, (beta _6I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _7I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},;beta _8I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2},{{{{{{{mathrm{and}}}}}}}};beta _9I_{{{{{{{{mathrm{virus}}}}}}}}_i = 2}) corresponds to the interaction effects between each of these variables and the viral decay, (alpha _{{{{{{{{mathrm{id}}}}}}}}_i}) corresponds to the mixed effect of the model and (varepsilon _i) corresponds to the error term of the model. The selection of the model is further detailed in the Supplementary Information (Supplementary Material and Figs. S1 and S2).The full dataset included in the correlation analysis and the CTM is provided in Supplementary Table 2. A description of the variables used is given in the Supplementary Information. The dataset was analyzed using the censReg package in R [19]. The R code is given in the Supplementary Information. More

  • in

    Modern aridity in the Altai-Sayan mountain range derived from multiple millennial proxies

    1500-year stable carbon and oxygen isotopes in larch tree-ring celluloseThe δ13Ccell (Fig. 1a, Fig. S2) and δ18Ocell (Fig. 1b, Fig. S3) records span 516–2016 CE, at annual resolution. The δ13Ccell timeseries shows mostly increasing trends during the first millennium of the Common Era (516–1120 CE), and similarly at the end of the last millennium (1720–2016 CE). The maximum δ13Ccell value occurs in 2016 CE (−19.6‰; + 3.2σ), while the minimum occurs in 686 CE (−24.7‰, −3.6σ) relative to the average for the period 516–2016 CE (−22.04‰) (Table S2, Fig. S2). The standard error (SE) for the whole analysed period is 0.02.Figure 1Annually resolved δ13Ccell (a) and δ18O cell (b) in Siberian larch tree-ring cellulose chronologies for the period from 516 to 2016 CE. Chronologies are smoothed by a 101-year Hamming window to highlight a centennial scale. The dotted and dashed lines indicate the number of trees analysed.Full size imageThe δ18Ocell timeseries (Fig. 1b, Fig. S3) showed two positive and one negative extreme over the past 1500 years, with the minimum value (19.9‰; −6.3σ), occurring in 536 CE, and maximum values (31.9‰; + 3.8σ and 32.2‰; + 4.4σ), occurring in 1266 and 2008 CE, respectively (Table S2, Fig. S3). The SE for the whole analysed period is 0.03. The δ18Ocell data has higher standard deviation (SD) (1.15) than δ13Ccell (0.75).Less than 1% of values in the δ18Ocell record are classified as extreme, with the standard deviation ≥  ± 3σ. The δ13Ccell and δ18Ocell records are significantly correlated (r = 0.1, p = 0.0001, n = 1500).Local climate signals preserved in δ13Ccell and δ18Ocell recordsWe used weather observations from the local Mugur-Aksy weather station (50°N, 90°E, 1850 m asl) (Table S1) to derive quantitative paleoclimatic reconstructions from our δ13Ccell and δ18Ocell timeseries. A multiple linear regression analysis revealed significant correlations between δ13Ccell and July precipitation (r = −0.58; p  More

  • in

    Enhanced spring warming in a Mediterranean mountain by atmospheric circulation

    Foster, G. & Rahmstorf, S. Global temperature evolution 1979–2010. Environ. Res. Lett. 6, 044022 (2011).ADS 
    Article 

    Google Scholar 
    Cahill, N., Rahmstorf, S. & Parnell, A. C. Change points of global temperature. Environ. Res. Lett. 10, 084002 (2015).ADS 
    Article 

    Google Scholar 
    Yan, X. H. et al. The global warming hiatus: Slowdown or redistribution?. Earth’s Future 4, 472–482 (2016).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Karl, T. R. et al. Possible artifacts of data biases in the recent global surface warming hiatus. Science 348, 5632 (2015).Article 

    Google Scholar 
    Cohen, J. L., Furtado, J. C., Barlow, M., Alexeev, V. A. & Cherry, J. E. Asymmetric seasonal temperature trends. Geophys. Res. Lett. 39, 04705. https://doi.org/10.1029/2011GL050582 (2012).ADS 
    Article 

    Google Scholar 
    Pepin, N. C. & Lundquist, J. D. Temperature trends at high elevations: patterns across the globe. Geophys. Res. Lett. 35, 14 (2008).Article 

    Google Scholar 
    Rangwala, I. & Miller, J. R. Climate change in mountains: a review of elevation-dependent warming and its possible causes. Clim. Change 114, 527–547 (2012).ADS 
    Article 

    Google Scholar 
    Wang, Q., Fan, X. & Wang, M. Recent warming amplification over high elevation regions across the globe. Clim. Dyn. 43, 87–101 (2014).CAS 
    Article 

    Google Scholar 
    Fan, X., Wang, Q., Wang, M. & Jiménez, C. V. Warming amplification of minimum and maximum temperatures over high-elevation regions across the globe. PLoS ONE 10, e0140213 (2015).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Pepin, N. et al. Elevation-dependent warming in mountain regions of the world. Nat. Clim. Chang. 5, 424 (2015).ADS 
    Article 

    Google Scholar 
    Piccarreta, M., Lazzari, M. & Pasini, A. Trends in daily temperature extremes over the Basilicata region (southern Italy) from 1951 to 2010 in a Mediterranean climatic context. Int. J. Climatol. 35, 1964–1975 (2015).Article 

    Google Scholar 
    Gonzalez-Hidalgo, J. C., Peña-Angulo, D., Brunetti, M. & Cortesi, N. Recent trend in temperature evolution in Spanish mainland (1951–2010): from warming to hiatus. Int. J. Climatol. 36, 2405–2416 (2016).Article 

    Google Scholar 
    McCullough, I. M. et al. High and dry: high elevations disproportionately exposed to regional climate change in Mediterranean-climate landscapes. Landsc. Ecol. 31, 1063–1075 (2016).Article 

    Google Scholar 
    Sanz-Elorza, M., Dana, E. D., González, A. & Sobrino, E. Changes in the high-mountain vegetation of the central Iberian Peninsula as a probable sign of global warming. Ann. Bot. 92, 273–280 (2003).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Peñuelas, J. & Boada, M. A global change induced biome shift in the Montseny mountains (NE Spain). Glob. Change Biol. 9, 131–140 (2003).ADS 
    Article 

    Google Scholar 
    Linares, J. C. & Tíscar, P. A. Climate change impacts and vulnerability of the southern populations of Pinus nigra subsp. salzmannii. Tree Physiol. 30, 795–806 (2010).PubMed 
    Article 

    Google Scholar 
    Giorgi, F., Hurrell, J. W., Marinucci, M. R. & Beniston, M. Elevation dependency of the surface climate change signal: a model study. J. Clim. 10, 288–296 (1997).ADS 
    Article 

    Google Scholar 
    Palazzi, E., Mortarini, L., Terzago, S. & Von Hardenberg, J. Elevation-dependent warming in global climate model simulations at high spatial resolution. Clim. Dyn. 52, 2685–2702 (2019).Article 

    Google Scholar 
    Poyatos, R., Latron, J. & Llorens, P. Land use and land cover change after agricultural abandonment. Mt. Res. Dev. 23, 362–368 (2003).Article 

    Google Scholar 
    Mouillot, F., Ratte, J. P., Joffre, R., Mouillot, D. & Rambal, S. Long-term forest dynamic after land abandonment in a fire prone Mediterranean landscape (central Corsica, France). Landsc. Ecol. 20, 101–112 (2005).Article 

    Google Scholar 
    Zellweger, F. et al. Forest microclimate dynamics drive plant responses to warming. Science 368, 772–775 (2020).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Ríos-Cornejo, D., Penas, Á., Álvarez-Esteban, R. & Del Río, S. Links between teleconnection patterns and mean temperature in Spain. Theor. Appl. Climatol. 122, 1–18 (2015).ADS 
    Article 

    Google Scholar 
    Nogués-Bravo, D., Araújo, M. B., Errea, M. P. & Martinez-Rica, J. P. Exposure of global mountain systems to climate warming during the 21st Century. Glob. Environ. Chang. 17, 420–428 (2007).Article 

    Google Scholar 
    Vicente-Serrano, S. M., Beguería, S., López-Moreno, J. I., El Kenawy, A. M. & Angulo-Martínez, M. Daily atmospheric circulation events and extreme precipitation risk in northeast Spain: Role of the North Atlantic Oscillation, the Western Mediterranean Oscillation, and the Mediterranean Oscillation. J. Geophys. Res. Atmos. 114, D8 (2009).Article 

    Google Scholar 
    Guzman-Morales, J. & Gershunov, A. Climate change suppresses Santa Ana winds of southern California and sharpens their seasonality. Geophys. Res. Lett. 46, 2772–2780. https://doi.org/10.1029/2018GL080261 (2019).ADS 
    Article 

    Google Scholar 
    Yu, M. & Ruggieri, E. Change point analysis of global temperature records. Int. J. Climatol. 39, 3679–3688 (2019).Article 

    Google Scholar 
    Giorgi, F. Climate change hot-spots. Geophys. Res. Lett. 33, 08707. https://doi.org/10.1029/2006GL025734 (2006).ADS 
    Article 

    Google Scholar 
    García, M. J. L. Recent warming in the Balearic Sea and Spanish Mediterranean coast: Towards an earlier and longer summer. Atmósfera 28, 149–160 (2015).Article 

    Google Scholar 
    Toreti, A., Desiato, F., Fioravanti, G. & Perconti, W. Seasonal temperatures over Italy and their relationship with low-frequency atmospheric circulation patterns. Clim. Change 99, 211–227 (2010).ADS 
    Article 

    Google Scholar 
    Scorzini, A. R. & Leopardi, M. Precipitation and temperature trends over central Italy (Abruzzo Region): 1951–2012. Theor. Appl. Climatol. 135, 959–977 (2019).ADS 
    Article 

    Google Scholar 
    Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. Nature 479, 384–387 (2011).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Juang, J.-Y., Katul, G., Siqueira, M., Stoy, P. & Novick, K. Separating the effects of albedo from eco-physiological changes on surface temperature along a successional chronosequence in the southeastern United States. Geophys. Res. Lett. 34, 21408. https://doi.org/10.1029/2007.GL03129 (2007).ADS 
    Article 

    Google Scholar 
    Boulant, N., Kunstler, G., Rambal, S. & Lepart, J. Seed supply, drought, and grazing determine spatio-temporal patterns of recruitment for native and introduced invasive pines in grasslands. Divers. Distrib. 14, 862–874 (2008).Article 

    Google Scholar 
    Améztegui, A. Land-use changes as major drivers of mountain pine (Pinus uncinata Ram.) expansion in the Pyrenees. Glob. Ecol. Biogeogr. 19, 632–641 (2010).
    Google Scholar 
    Rambal, S. Relations entre couverts végétaux des parcours et cycle de l’eau. In L’eau des troupeaux en alpages et sur parcours: une ressource à gérer, aménager, partager (ed. Lepart, J.) 25–37 (Association Française de Pastoralisme et Cardère éditeur, 2015).
    Google Scholar 
    Fonderflick, J., Lepart, J., Caplat, P., Debussche, M. & Marty, P. Managing agricultural change for biodiversity conservation in a Mediterranean upland. Biol. Conserv. 143, 737–746 (2010).Article 

    Google Scholar 
    Abadie, J. et al. Forest recovery since 1860 in a Mediterranean region: Drivers and implications for land use and land cover spatial distribution. Landsc. Ecol. 33, 289–305 (2018).Article 

    Google Scholar 
    Cervera, T., Pino, J., Marull, J., Padró, R. & Tello, E. Understanding the long-term dynamics of forest transition: From deforestation to afforestation in a Mediterranean landscape (Catalonia, 1868–2005). Land Use Policy 80, 318–331 (2019).Article 

    Google Scholar 
    Wolpert, F., Quintas-Soriano, C. & Plieninger, T. Exploring land-use histories of tree-crop landscapes: a cross-site comparison in the Mediterranean Basin. Sustain. Sci. 15, 1267–1283 (2020).Article 

    Google Scholar 
    Lasanta-Martínez, T., Vicente-Serrano, S. M. & Cuadrat-Prats, J. M. Mountain Mediterranean landscape evolution caused by the abandonment of traditional primary activities: A study of the Spanish Central Pyrenees. Appl. Geogr. 25, 47–65 (2005).Article 

    Google Scholar 
    Malandra, F., Vitali, A., Urbinati, C., Weisberg, P. J. & Garbarino, M. Patterns and drivers of forest landscape change in the Apennines range, Italy. Reg. Environ. Change 19, 1973–1985 (2019).Article 

    Google Scholar 
    Zhang, Q. et al. Reforestation and surface cooling in temperate zones: Mechanisms and implications. Glob. Change Biol. 26, 3384–3401 (2020).ADS 
    Article 

    Google Scholar 
    Rambal, S., Lacaze, B. & Winkel, T. Testing an area-weighted model for albedo or surface temperature of mixed pixels in Mediterranean woodlands. Int. J. Remote Sens. 11, 1495–1499 (1990).Article 

    Google Scholar 
    Luyssaert, S. et al. Land management and land-cover change have impacts of similar magnitude on surface temperature. Nat. Clim. Change 4, 389–393. https://doi.org/10.1038/nclimate2196 (2014).ADS 
    Article 

    Google Scholar 
    Novick, K. A. & Katul, G. G. The duality of reforestation impacts on surface and air temperature. J. Geophys. Res. Biogeosci. 125, e05543 (2020).Article 

    Google Scholar 
    Davy, R. & Esau, I. Differences in the efficacy of climate forcings explained by variations in atmospheric boundary layer depth. Nat. Commun. 7, 1–8 (2016).Article 

    Google Scholar 
    Serafin, S. et al. Exchange processes in the atmospheric boundary layer over mountainous terrain. Atmosphere 9, 102. https://doi.org/10.3390/atmos9030102 (2018).ADS 
    CAS 
    Article 

    Google Scholar 
    Perugini, L. et al. Biophysical effects on temperature and precipitation due to land cover change. Environ. Res. Lett. 12, 053002 (2017).ADS 
    Article 

    Google Scholar 
    Visbeck, M. H., Hurrell, J. W., Polvani, L. & Cullen, H. M. The North Atlantic oscillation: Past, present, and future. Proc. Natl. Acad. Sci. 98, 12876–12877 (2001).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Hurrell, J. W. Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation. Science 269, 676–679 (1995).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Martín, P., Sabatés, A., Lloret, J. & Martin-Vide, J. Climate modulation of fish populations: the role of the Western Mediterranean Oscillation (WeMO) in sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) production in the north-western Mediterranean. Clim. Change 110, 925–939 (2012).ADS 
    Article 

    Google Scholar 
    Schwingshackl, C., Hirschi, M. & Seneviratne, S. I. Global contributions of incoming radiation and land surface conditions to maximum near surface air temperature variability and trend. Geophys. Res. Lett. 45, 5034–5044 (2018).ADS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Philipona, R., Behrens, K. & Ruckstuhl, C. How declining aerosols and rising greenhouse gases forced rapid warming in Europe since the 1980s. Geophys. Res. Lett. 36, L02806. https://doi.org/10.1029/2008GL036350 (2009).ADS 
    Article 

    Google Scholar 
    Schwarz, M., Folini, D., Yang, S., Allan, R. P. & Wild, M. Changes in atmospheric shortwave absorption as important driver of dimming and brightening. Nat. Geosci. 13, 110–115 (2020).ADS 
    CAS 
    Article 

    Google Scholar 
    Norris, J. R. & Wild, M. Trends in aerosol radiative effects over Europe inferred from observed cloud cover, solar “dimming”, and solar “brightening”. J. Geophys. Res. Atmos. 112, D08214. https://doi.org/10.1029/2006JD007794 (2007).ADS 
    Article 

    Google Scholar 
    Mateos, D. et al. Quantifying the respective roles of aerosols and clouds in the strong brightening since the early 2000s over the Iberian Peninsula. J. Geophys. Res. Atmos. 119, 10–382 (2014).Article 

    Google Scholar 
    Sanchez-Lorenzo, A. et al. Reassessment and update of long-term trends in downward surface shortwave radiation over Europe (1939–2012). J. Geophys. Res. Atmos. 120, 9555–9569 (2015).ADS 
    Article 

    Google Scholar 
    Kambezidis, H. D., Kaskaoutis, D. G., Kalliampakos, G. K., Rashki, A. & Wild, M. The solar dimming/brightening effect over the Mediterranean Basin in the period 1979–2012. J. Atmos. Solar Terr. Phys. 150, 31–46 (2016).ADS 
    Article 

    Google Scholar 
    Chiacchio, M. & Wild, M. Influence of NAO and clouds on long-term seasonal variations of surface solar radiation in Europe. J. Geophys. Res. Atmos. 115, 0022. https://doi.org/10.1029/2009JD012182 (2010).Article 

    Google Scholar 
    Wild, M. Decadal changes in radiative fluxes at land and ocean surfaces and their relevance for global warming. Wiley Interdiscipl. Rev. Clim. Change 7, 91–107 (2016).Article 

    Google Scholar 
    Held, I. M. & Soden, B. J. Water vapor feedback and global warming. Annu. Rev. Energy Environ. 25, 441–475 (2000).Article 

    Google Scholar 
    Dessler, A. E. & Sherwood, S. C. A matter of humidity. Science 323, 1020–1021 (2009).CAS 
    PubMed 
    Article 

    Google Scholar 
    Ruckstuhl, C., Philipona, R., Morland, J. & Ohmura, A. Observed relationship between surface specific humidity, integrated water vapor, and longwave downward radiation at different altitudes. J. Geophys. Res. Atmos. 112(D03302), 2007. https://doi.org/10.1029/2006JD007850 (2007).Article 

    Google Scholar 
    Parras-Berrocal, I. M. et al. The climate change signal in the Mediterranean Sea in a regionally coupled atmosphere–ocean model. Ocean Sci. 16, 743–765. https://doi.org/10.5194/os-16-743-2020 (2020).ADS 
    Article 

    Google Scholar 
    Reale, M. et al. The regional earth system model RegCM-ES: Evaluation of the Mediterranean climate and marine biogeochemistry. J. Adv. Model. Earth Syst. 12, e001812 (2020).Article 

    Google Scholar 
    Sen, P. K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 63, 1379–1389 (1968).MathSciNet 
    MATH 
    Article 

    Google Scholar 
    Kelliher, F. M., Leuning, R. & Schulze, E. D. Evaporation and canopy characteristics of coniferous forests and grasslands. Oecologia 95, 153–163 (1993).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Linacre, E. T. Simpler empirical expression for actual evapotranspiration rates-a discussion. Agric. Meteorol. 11, 451–452 (1973).Article 

    Google Scholar 
    Jones, P. D., Jónsson, T. & Wheeler, D. Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland. Int. J. Climatol. 17, 1433–1450 (1997).Article 

    Google Scholar 
    Palutikof, J. P. Analysis of Mediterranean climate data: measured and modelled. In Mediterranean Climate: Variability and Trends (ed. Bolle, H. J.) (Springer, 2003).
    Google Scholar 
    Martin-Vide, J. & Lopez-Bustins, J. A. The western Mediterranean oscillation and rainfall in the Iberian Peninsula. Int. J. Climatol. 26, 1455–1475 (2006).Article 

    Google Scholar  More

  • in

    Unravelling seasonal trends in coastal marine heatwave metrics across global biogeographical realms

    Smale, D. A. et al. Marine heatwaves threaten global biodiversity and the provision of ecosystem services. Nat. Clim. Change https://doi.org/10.1038/s41558-019-0412-1 (2019).Article 

    Google Scholar 
    Smith, K. E. et al. Socioeconomic impacts of marine heatwaves: Global issues and opportunities. Science 374, eabj3593 (2021).CAS 
    Article 

    Google Scholar 
    Frolicher, T. L., Fischer, E. M. & Gruber, N. Marine heatwaves under global warming. Nature 560, 360–364. https://doi.org/10.1038/s41586-018-0383-9 (2018).ADS 
    CAS 
    Article 
    PubMed 

    Google Scholar 
    Oliver, et al. Longer and more frequent marine heatwaves over the past century. Nat. Commun. https://doi.org/10.1038/s41467-018-03732-9 (2018).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Oliver, et al. Projected marine heatwaves in the 21st century and the potential for ecological impact. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00734 (2019).Article 

    Google Scholar 
    Hobday, A. J. et al. A hierarchical approach to defining marine heatwaves. Prog. Oceanogr. 141, 227–238 (2016).ADS 
    Article 

    Google Scholar 
    Banzon, V., Smith, T. M., Chin, T. M., Liu, C. & Hankins, W. A long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies. Earth Syst. Sci. Data 8, 165–176 (2016).ADS 
    Article 

    Google Scholar 
    Wernberg, T. et al. An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat. Clim. Change 3, 78–82. https://doi.org/10.1038/nclimate1627 (2013).ADS 
    Article 

    Google Scholar 
    Arias-Ortiz, A. et al. A marine heatwave drives massive losses from the world’s largest seagrass carbon stocks. Nat. Clim. Change https://doi.org/10.1038/s41558-018-0096-y (2018).Article 

    Google Scholar 
    Smale, D. A. Impacts of ocean warming on kelp forest ecosystems. New Phytol. 225, 1447–1454 (2020).Article 

    Google Scholar 
    Couch, C. S. et al. Mass coral bleaching due to unprecedented marine heatwave in Papahānaumokuākea Marine National Monument (Northwestern Hawaiian Islands). PLoS ONE 12, e0185121. https://doi.org/10.1371/journal.pone.0185121 (2017).CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Oliver, E. C. J. et al. The unprecedented 2015/16 Tasman Sea marine heatwave. Nat. Commun. https://doi.org/10.1038/ncomms16101 (2017).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Montie, S., Thomsen, M. S., Rack, W. & Broady, P. A. Extreme summer marine heatwaves increase chlorophyll a in the Southern Ocean. Antarct. Sci. 32, 508–509 (2020).ADS 
    Article 

    Google Scholar 
    Gupta, A. S. et al. Drivers and impacts of the most extreme marine heatwaves events. Sci. Rep. 10, 1–15 (2020).ADS 
    Article 

    Google Scholar 
    Holbrook, N. J. et al. A global assessment of marine heatwaves and their drivers. Nat. Commun. 10, 1–13 (2019).CAS 
    Article 

    Google Scholar 
    La Sorte, F. A., Johnston, A. & Ault, T. R. Global trends in the frequency and duration of temperature extremes. Clim. Change 166, 1. https://doi.org/10.1007/s10584-021-03094-0 (2021).ADS 
    Article 

    Google Scholar 
    Thomsen, et al. Local extinction of bull kelp (Durvillaea spp.) due to a marine heatwave. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00084 (2019).Article 

    Google Scholar 
    Strydom, S. et al. Too hot to handle: Unprecedented seagrass death driven by marine heatwave in a World Heritage Area. Glob. Change Biol. 26, 3525–3538. https://doi.org/10.1111/gcb.15065 (2020).ADS 
    Article 

    Google Scholar 
    Leggat, W. P. et al. Rapid coral decay is associated with marine heatwave mortality events on reefs. Curr. Biol. 29, 2723-2730.e2724. https://doi.org/10.1016/j.cub.2019.06.077 (2019).CAS 
    Article 
    PubMed 

    Google Scholar 
    Wernberg, T. et al. Climate-driven regime shift of a temperate marine ecosystem. Science 353, 169–172. https://doi.org/10.1126/science.aad8745 (2016).ADS 
    CAS 
    Article 
    PubMed 

    Google Scholar 
    Thomsen, M. S. & McGlathery, K. Facilitation of macroalgae by the sedimentary tube forming polychaete Diopatra cuprea. Estuar. Coast. Shelf Sci. 62, 63–73. https://doi.org/10.1016/j.ecss.2004.08.007 (2005).ADS 
    Article 

    Google Scholar 
    Spalding, M. D. et al. Marine ecoregions of the world: A bioregionalization of coastal and shelf areas. Bioscience 57, 573–583 (2007).Article 

    Google Scholar 
    Costello, M. J. & Chaudhary, C. Marine biodiversity, biogeography, deep-sea gradients, and conservation. Curr. Biol. 27, R511–R527. https://doi.org/10.1016/j.cub.2017.04.060 (2017).CAS 
    Article 
    PubMed 

    Google Scholar 
    Tait, L. W., Thoral, F., Pinkerton, M. H., Thomsen, M. S. & Schiel, D. R. Loss of giant kelp, Macrocystis pyrifera, driven by marine heatwaves and exacerbated by poor water clarity in New Zealand. Front. Mar. Sci. https://doi.org/10.3389/fmars.2021.721087 (2021).Article 

    Google Scholar 
    Marin, M., Feng, M., Phillips, H. E. & Bindoff, N. L. A global, multiproduct analysis of coastal marine heatwaves: Distribution, characteristics, and long-term trends. J. Geophys. Res. Oceans 126, e2020JC016708. https://doi.org/10.1029/2020JC016708 (2021).ADS 
    Article 

    Google Scholar 
    Kain, J. M. The seasons in the subtidal. Brit. Phycol. J. 24, 203–215 (1989).Article 

    Google Scholar 
    Atkinson, J., King, N. G., Wilmes, S. B. & Moore, P. J. Summer and winter marine heatwaves favor an invasive over native seaweeds. J. Phycol. 56, 1591–1600. https://doi.org/10.1111/jpy.13051 (2020).CAS 
    Article 
    PubMed 

    Google Scholar 
    Salinger, M. J. et al. The unprecedented coupled ocean-atmosphere summer heatwave in the New Zealand region 2017/18: Drivers, mechanisms and impacts. Environ. Res. Lett. 14, 044023 (2019).ADS 
    Article 

    Google Scholar 
    Amaya, D. J., Miller, A. J., Xie, S.-P. & Kosaka, Y. Physical drivers of the summer 2019 North Pacific marine heatwave. Nat. Commun. 11, 1–9 (2020).Article 

    Google Scholar 
    Di Lorenzo, E. & Mantua, N. Multi-year persistence of the 2014/15 North Pacific marine heatwave. Nat. Clim. Change 6, 1042–1047. https://doi.org/10.1038/nclimate3082 (2016).ADS 
    Article 

    Google Scholar 
    Cayan, D. R. Large-scale relationships between sea surface temperature and surface air temperature. Mon. Weather Rev. 108, 1293–1301 (1980).ADS 
    Article 

    Google Scholar 
    Hipel, K. W. & McLeod, A. I. Time Series Modelling of Water Resources and Environmental Systems (Elsevier, 1994).
    Google Scholar 
    trend: non-parametric trend tests and changepoint detection.–R package ver. 1.1. 2 (2020).Costanza, R. et al. The value of the world’s ecosystem services and natural capital. Nature 387, 253–260 (1997).ADS 
    CAS 
    Article 

    Google Scholar 
    Halpern, B. S. et al. A global map of human impact on marine ecosystems. Science 319, 948–952 (2008).ADS 
    CAS 
    Article 

    Google Scholar 
    Harley, C. D. et al. The impacts of climate change in coastal marine systems. Ecol. Lett. 9, 228–241. https://doi.org/10.1111/j.1461-0248.2005.00871.x (2006).ADS 
    Article 
    PubMed 

    Google Scholar 
    Thomsen, M. S. & South, P. M. Communities and attachment networks associated with primary, secondary and alternative foundation species; a case study of stressed and disturbed stands of southern bull kelp. Diversity 11, 56. https://doi.org/10.3390/d11040056 (2019).Article 

    Google Scholar 
    Smale, D. A. & Wernberg, T. Extreme climatic event drives range contraction of a habitat-forming species. Proc. R. Soc. B Biol. Sci. 280, 20122829 (2013).Article 

    Google Scholar 
    Thomsen, M. S. et al. Cascading impacts of earthquakes and extreme heatwaves have destroyed populations of an iconic marine foundation species. Divers. Distrib. (2021).Rogers-Bennett, L. & Catton, C. Marine heat wave and multiple stressors tip bull kelp forest to sea urchin barrens. Sci. Rep. 9, 1–9 (2019).CAS 
    Article 

    Google Scholar 
    Filbee-Dexter, K. et al. Marine heatwaves and the collapse of marginal North Atlantic kelp forests. Sci. Rep. 10, 1–11 (2020).Article 

    Google Scholar 
    Thomson, J. A. et al. Extreme temperatures, foundation species, and abrupt ecosystem change: An example from an iconic seagrass ecosystem. Glob. Change Biol. 21, 1463–1474. https://doi.org/10.1111/gcb.12694 (2015).ADS 
    Article 

    Google Scholar 
    Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377. https://doi.org/10.1038/nature21707 (2017).ADS 
    CAS 
    Article 
    PubMed 

    Google Scholar 
    Le Nohaïc, M. et al. Marine heatwave causes unprecedented regional mass bleaching of thermally resistant corals in northwestern Australia. Sci. Rep. 7, 14999. https://doi.org/10.1038/s41598-017-14794-y (2017).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Smale, D. A., Wernberg, T. & Vanderklift, M. A. Regional-scale variability in the response of benthic macroinvertebrate assemblages to a marine heatwave. Mar. Ecol. Prog. Ser. 568, 17–30. https://doi.org/10.3354/meps12080 (2017).ADS 
    Article 

    Google Scholar 
    Cavole, L. et al. Biological impacts of the 2013–2015 warm-water anomaly in the Northeast Pacific: Winners, losers, and the future. Oceanography (Washington D.C.) https://doi.org/10.5670/oceanog.2016.32 (2016).Article 

    Google Scholar 
    Coleman, M. A., Minne, A. J. P., Vranken, S. & Wernberg, T. Genetic tropicalisation following a marine heatwave. Sci. Rep. 10, 12726. https://doi.org/10.1038/s41598-020-69665-w (2020).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Collette, B. B. in Reproduction and sexuality in marine fishes 21–64 (University of California Press, 2010).Yatsu, A. & Shimada, H. Distributions of Epipelagic Fishes, Squids, Marine Mammals. Bulletin 53 North Pacific Commission, 111–146.Hirst, A., Roff, J. & Lampitt, R. A synthesis of growth rates in marine epipelagic invertebrate zooplankton. Adv. Mar. Biol. 44, 1–142 (2003).CAS 
    Article 

    Google Scholar 
    Smale, D. A. & Wernberg, T. Satellite-derived SST data as a proxy for water temperature in nearshore benthic ecology. Mar. Ecol. Prog. Ser. 387, 27–37 (2009).ADS 
    Article 

    Google Scholar 
    Bernardello, R., Serrano, E., Coma, R., Ribes, M. & Bahamon, N. A comparison of remote-sensing SST and in situ seawater temperature in near-shore habitats in the western Mediterranean Sea. Mar. Ecol. Prog. Ser. 559, 21–34 (2016).ADS 
    Article 

    Google Scholar 
    Brewin, R. J. et al. Evaluating operational AVHRR sea surface temperature data at the coastline using benthic temperature loggers. Remote Sens. 10, 925 (2018).ADS 
    Article 

    Google Scholar 
    Smit, A. J. et al. A coastal seawater temperature dataset for biogeographical studies: Large biases between in situ and remotely-sensed data sets around the Coast of South Africa. PLoS ONE 8, e81944. https://doi.org/10.1371/journal.pone.0081944 (2013).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Marin, M., Bindoff, N. L., Feng, M. & Phillips, H. E. Slower long-term coastal warming drives dampened trends in coastal marine heatwave exposure. J. Geophys. Res. Oceans https://doi.org/10.1029/2021jc017930 (2021).Article 

    Google Scholar 
    Lourenço, C. R. et al. Upwelling areas as climate change refugia for the distribution and genetic diversity of a marine macroalga. J. Biogeogr. 43, 1595–1607. https://doi.org/10.1111/jbi.12744 (2016).Article 

    Google Scholar 
    Riegl, B. & Piller, W. E. Possible refugia for reefs in times of environmental stress. Int. J. Earth Sci. 92, 520–531. https://doi.org/10.1007/s00531-003-0328-9 (2003).Article 

    Google Scholar 
    El Glynn, P. W. Nino-southern oscillation 1982–1983: Nearshore population, community, and ecosystem responses. Annu. Rev. Ecol. Syst. 19, 309–346. https://doi.org/10.1146/annurev.es.19.110188.001521 (1988).Article 

    Google Scholar 
    Glynn, P. W. & D’Croz, L. Experimental evidence for high temperature stress as the cause of El Niño-coincident coral mortality. Coral Reefs 8, 181–191. https://doi.org/10.1007/bf00265009 (1990).ADS 
    Article 

    Google Scholar 
    Glynn, P. W., Maté, J. L., Baker, A. C. & Calderón, M. O. Coral bleaching and mortality in panama and ecuador during the 1997–1998 El Niño-Southern Oscillation Event: Spatial/temporal patterns and comparisons with the 1982–1983 event. Bull. Mar. Sci. 69, 79–109 (2001).
    Google Scholar 
    Podestá, G. P. & Glynn, P. W. The 1997–98 El Niño event in Panama and Galápagos: An update of thermal stress indices relative to coral bleaching. Bull. Mar. Sci. 69, 43–59 (2001).
    Google Scholar 
    Ladah, L. B. & Zertuche-Gonzalez, J. A. Giant kelp (Macrocystis pyrifera) survival in deep water (25–40 m) during El Nino of 1997–1998 in Baja California, Mexico. Bot. Marina 47, 367–372. https://doi.org/10.1515/bot.2004.054 (2004).Article 

    Google Scholar 
    Kayanne, H. Validation of degree heating weeks as a coral bleaching index in the northwestern Pacific. Coral Reefs 36, 63–70 (2017).ADS 
    Article 

    Google Scholar 
    Le Nohaïc, M. et al. Marine heatwave causes unprecedented regional mass bleaching of thermally resistant corals in northwestern Australia. Sci. Rep. 7, 1–11 (2017).ADS 
    Article 

    Google Scholar 
    Marba, N. & Duarte, C. M. Mediterranean warming triggers seagrass (Posidonia oceanica) shoot mortality. Glob. Change Biol. 16, 2366–2375. https://doi.org/10.1111/j.1365-2486.2009.02130.x (2010).ADS 
    Article 

    Google Scholar 
    Bennett, S., Wernberg, T., Arackal Joy, B., de Bettignies, T. & Campbell, A. H. Central and rear-edge populations can be equally vulnerable to warming. Nat. Commun. 6, 10280. https://doi.org/10.1038/ncomms10280 (2015).ADS 
    CAS 
    Article 
    PubMed 

    Google Scholar 
    Filbee-Dexter, K. et al. Marine heatwaves and the collapse of marginal North Atlantic kelp forests. Sci. Rep. 10, 13388. https://doi.org/10.1038/s41598-020-70273-x (2020).ADS 
    CAS 
    Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Snover, M. L. Ontogenetic habitat shifts in marine organisms: Influencing factors and the impact of climate variability. Bull. Mar. Sci. 83, 53–67 (2008).
    Google Scholar 
    Harley, C. D. Climate change, keystone predation, and biodiversity loss. Science 334, 1124–1127 (2011).ADS 
    CAS 
    Article 

    Google Scholar 
    Kelaher, B. P., Coleman, M. A. & Bishop, M. J. Ocean warming, but not acidification, accelerates seagrass decomposition under near-future climate scenarios. Mar. Ecol. Prog. Ser. 605, 103–110 (2018).ADS 
    CAS 
    Article 

    Google Scholar 
    De Senerpont Domis, L. N. et al. Plankton dynamics under different climatic conditions in space and time. Freshw. Biol. 58, 463–482 (2013).Article 

    Google Scholar 
    Fossheim, M. et al. Recent warming leads to a rapid borealization of fish communities in the Arctic. Nat. Clim. Change 5, 673–677. https://doi.org/10.1038/nclimate2647 (2015).ADS 
    Article 

    Google Scholar 
    Morales-Nin, B. & Panfili, J. Seasonality in the deep sea and tropics revisited: What can otoliths tell us?. Mar. Freshw. Res. 56, 585–598 (2005).Article 

    Google Scholar 
    Alongi, D. M. Ecology of tropical soft-bottom benthos: A review with emphasis on emerging concepts. Rev. Biol. Trop. 37, 85–100 (1989).
    Google Scholar 
    Hobday, A. J., Spillman, C. M., Paige Eveson, J. & Hartog, J. R. Seasonal forecasting for decision support in marine fisheries and aquaculture. Fish. Oceanogr. 25, 45–56. https://doi.org/10.1111/fog.12083 (2016).Article 

    Google Scholar 
    Spillman, C. M., Smith, G. A., Hobday, A. J. & Hartog, J. R. Onset and decline rates of marine heatwaves: Global trends, seasonal forecasts and marine management. Front. Clim. https://doi.org/10.3389/fclim.2021.801217 (2021).Article 

    Google Scholar 
    Schlegel, R. W., Oliver, E. C. J., Wernberg, T. & Smit, A. J. Nearshore and offshore co-occurrence of marine heatwaves and cold-spells. Prog. Oceanogr. 151, 189–205. https://doi.org/10.1016/j.pocean.2017.01.004 (2017).ADS 
    Article 

    Google Scholar 
    Huang, B. et al. Improvements of the daily optimum interpolation sea surface temperature (DOISST) Version 2.1. J. Clim. 34, 2923–2939. https://doi.org/10.1175/jcli-d-20-0166.1 (2021).ADS 
    Article 

    Google Scholar 
    OBPG, N. & Stumpf, R. P. Distance to Nearest Coastline: 0.01-Degree Grid. Distributed by the Pacific Islands Ocean Observing System (PacIOOS). http://pacioos.org/metadata/dist2coast_1deg.html and https://data.noaa.gov/dataset/dataset/distance-to-nearest-coastline-0-01-degree-grid2http://www.pacioos.hawaii.edu/metadata/dist2coast_1deg.html (2012).Schlegel, R. W. & Smit, A. J. heatwaveR: A central algorithm for the detection of heatwaves and cold-spells. J. Open Source Softw. 3(27), 821 (2018).ADS 
    Article 

    Google Scholar 
    Sakurai, T., Yukio, K. & Kuragano, T. in Proceedings. 2005 IEEE International Geoscience and Remote Sensing Symposium, 2005. IGARSS’05. 2606–2608 (IEEE). More

  • in

    Temporal patterns in the soundscape of a Norwegian gateway to the Arctic

    Ellison, W. T., Southall, B. L., Clark, C. W. & Frankel, A. S. A new context-based approach to assess marine mammal behavioral responses to anthropogenic sounds. Conserv. Biol. 26, 21–28 (2012).CAS 
    PubMed 
    Article 

    Google Scholar 
    Williams, R., Clark, C. W., Ponirakis, D. & Ashe, E. Acoustic quality of critical habitats for three threatened whale populations. Anim. Conserv. 17, 174–185 (2014).Article 

    Google Scholar 
    Halliday, W. D., Pine, M. K. & Insley, S. J. Underwater noise and arctic marine mammals: Review and policy recommendations. Environ. Rev. 28, 438–448 (2020).Article 

    Google Scholar 
    Kvadsheim, P. H. et al. Impact of Anthropogenic Noise on the Marine Environment: Status of Knowledge and Management (Springer, 2020).
    Google Scholar 
    Weilgart, L. S. & Whitehead, H. Distinctive vocalizations from mature male sperm whales (Physeter macrocephalus). Can. J. Zool. 66, 1931–1937 (1988).Article 

    Google Scholar 
    Simon, M., Stafford, K. M., Beedholm, K., Lee, C. M. & Madsen, P. T. Singing behavior of fin whales in the Davis Strait with implications for mating, migration and foraging. J. Acoust. Soc. Am. 128, 3200 (2010).ADS 
    PubMed 
    Article 

    Google Scholar 
    Alves, D., Amorim, M. C. P. & Fonseca, P. J. Assessing acoustic communication active space in the Lusitanian toadfish. J. Exp. Biol. 219, 1122–1129 (2016).PubMed 

    Google Scholar 
    Linnenschmidt, M., Teilmann, J., Akamatsu, T., Dietz, R. & Miller, L. A. Biosonar, dive, and foraging activity of satellite tracked harbor porpoises (Phocoena phocoena). Mar. Mamm. Sci. 29, E77–E97 (2013).Article 

    Google Scholar 
    Giorli, G. & Goetz, K. T. Foraging activity of sperm whales (Physeter macrocephalus) off the east coast of New Zealand. Sci. Rep. 9, 12182 (2019).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Baumgartner, M. F. & Fratantoni, D. M. Diel periodicity in both sei whale vocalization rates and the vertical migration of their copepod prey observed from ocean gliders. Limnol. Oceanogr. 53, 2197–2209 (2008).ADS 
    Article 

    Google Scholar 
    Urazghildiiev, I. R. & Van Parijs, S. M. Automatic grunt detector and recognizer for Atlantic cod (Gadus morhua). J. Acoust. Soc. Am. 139, 2532–2540 (2016).ADS 
    PubMed 
    Article 

    Google Scholar 
    Ladich, F. Ecology of sound communication in fishes. Fish Fish. 20, 552–563 (2019).Article 

    Google Scholar 
    Radford, C. A., Stanley, J. A., Simpson, S. D. & Jeffs, A. G. Juvenile coral reef fish use sound to locate habitats. Coral Reefs 30, 295–305 (2011).ADS 
    Article 

    Google Scholar 
    Pierpoint, C. Harbour porpoise (Phocoena phocoena) foraging strategy at a high energy, near-shore site in south-west Wales, UK. J. Mar. Biol. Assoc. UK 88, 1167–1173 (2008).Article 

    Google Scholar 
    Pijanowski, B. C. et al. Soundscape ecology: The science of sound in the landscape. Bioscience 61, 203–216 (2011).Article 

    Google Scholar 
    Stanley, J. A., Radford, C. A. & Jeffs, A. G. Location, location, location: Finding a suitable home among the noise. Proc. R. Soc. B Biol. Sci. 279, 3622–3631 (2012).Article 

    Google Scholar 
    Buscaino, G. et al. Temporal patterns in the soundscape of the shallow waters of a Mediterranean marine protected area. Sci. Rep. 6, 1–13 (2016).Article 
    CAS 

    Google Scholar 
    Gasc, A., Francomano, D., Dunning, J. B. & Pijanowski, B. C. Future directions for soundscape ecology: The importance of ornithological contributions. Auk 134, 215–228 (2017).Article 

    Google Scholar 
    Putland, R. L., Constantine, R. & Radford, C. A. Exploring spatial and temporal trends in the soundscape of an ecologically significant embayment. Sci. Rep. 7, 5713 (2017).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Pieretti, N., LoMartire, M., Farina, A. & Danovaro, R. Marine soundscape as an additional biodiversity monitoring tool: A case study from the Adriatic Sea (Mediterranean Sea). Ecol. Indic. 83, 13–20 (2017).Article 

    Google Scholar 
    Gillespie, D., Palmer, L., Macaulay, J., Sparling, C. & Hastie, G. Passive acoustic methods for tracking the 3D movements of small cetaceans around marine structures. PLoS ONE 15, e0229058 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Van Parijs, S. et al. Management and research applications of real-time and archival passive acoustic sensors over varying temporal and spatial scales. Mar. Ecol. Prog. Ser. 395, 21–36 (2009).ADS 
    Article 

    Google Scholar 
    Warren, V. E., McPherson, C., Giorli, G., Goetz, K. T. & Radford, C. A. Marine soundscape variation reveals insights into baleen whales and their environment: a case study in central New Zealand. R. Soc. Open Sci. 8, 201503 (2021).ADS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Ahonen, H. et al. The underwater soundscape in western Fram Strait: Breeding ground of Spitsbergen’s endangered bowhead whales. Mar. Pollut. Bull. 123, 97–112 (2017).CAS 
    PubMed 
    Article 

    Google Scholar 
    Hildebrand, J. A. Anthropogenic and natural sources of ambient noise in the ocean. Mar. Ecol. Prog. Ser. 395, 5–20 (2009).ADS 
    Article 

    Google Scholar 
    Ross, D. Ship sources of ambient noise. IEEE J. Ocean. Eng. 30, 257–261 (2005).ADS 
    Article 

    Google Scholar 
    Popper, A. N. & Hawkins, A. The Effects of Noise on Aquatic Life Vol. 730 (Springer, 2012).Book 

    Google Scholar 
    Hubert, J. et al. Effects of broadband sound exposure on the interaction between foraging crab and shrimp: A field study. Environ. Pollut. 243, 1923–1929 (2018).CAS 
    PubMed 
    Article 

    Google Scholar 
    Weilgart, L. The Impact of Ocean Noise Pollution on Fish and Invertebrates (Springer, 2018).
    Google Scholar 
    Kvadsheim, H., Sivle, L. D., Hansen, R. R. & Karlsen, H. E. Effekter av Menneskeskapt støy på Havmiljø Rapport til Miljødirektoratet om Kunnskapsstatus FFI-RAPPORT. (2017).Parks, S. E., Johnson, M., Nowacek, D. & Tyack, P. L. Individual right whales call louder in increased environmental noise. Biol. Lett. 7, 33–35 (2011).PubMed 
    Article 

    Google Scholar 
    Meh, F. et al. Humpback whales Megaptera novaeangliae alter calling behavior in response to natural sounds and vessel noise. Mar. Ecol. Prog. Ser. 607, 251–268 (2018).Article 

    Google Scholar 
    Leroy, E. C., Royer, J.-Y., Bonnel, J. & Samaran, F. Long-term and seasonal changes of large whale call frequency in the Southern Indian ocean. J. Geophys. Res. Ocean. 123, 8568–8580 (2018).ADS 
    Article 

    Google Scholar 
    Parks, S. E., Clark, C. W. & Tyack, P. L. Short- and long-term changes in right whale calling behavior: The potential effects of noise on acoustic communication. J. Acoust. Soc. Am. 122, 3725–3731 (2007).ADS 
    PubMed 
    Article 

    Google Scholar 
    Clark, C. et al. Acoustic masking in marine ecosystems: intuitions, analysis, and implication. Mar. Ecol. Prog. Ser. 395, 201–222 (2009).ADS 
    Article 

    Google Scholar 
    PAME. Underwater Noise in the Arctic: A State of Knowledge Report (PAME, 2019).
    Google Scholar 
    Beszczynska-Möller, A., Woodgate, R., Lee, C., Melling, H. & Karcher, M. A synthesis of exchanges through the main oceanic gateways to the Arctic Ocean. Oceanography 24, 82–99 (2011).Article 

    Google Scholar 
    Ramm, T. Hungry During Migration? Humpback Whale Movement from the Barents Sea to a Feeding Stopover in Northern Norway Revealed by Photo-ID Analysis. (MSc thesis. UiT The Arctic University of Norway, 2020).Broms, F. et al. Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway. in Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway (ed. Wenzel, F. W.) (European Cetacean Society Special Publication Series, 2015).Aniceto, A. S. et al. Arctic marine data collection using oceanic gliders: Providing ecological context to cetacean vocalizations. Front. Mar. Sci. 7, 547 (2020).Article 

    Google Scholar 
    Jourdain, E. & Vongraven, D. Humpback whale (Megaptera novaeangliae) and killer whale (Orcinus orca) feeding aggregations for foraging on herring (Clupea harengus) in Northern Norway. Mamm. Biol. 86, 27–32 (2017).Article 

    Google Scholar 
    Christiansen, J. S., Mecklenburg, C. W. & Karamushko, O. V. Arctic marine fishes and their fisheries in light of global change. Glob. Chang. Biol. 20, 352–359 (2014).ADS 
    PubMed 
    Article 

    Google Scholar 
    Rødland, E. S. & Bjørge, A. Residency and abundance of sperm whales (Physeter macrocephalus) in the Bleik Canyon, Norway. Mar. Biol. Res. 11, 974–982 (2015).Article 

    Google Scholar 
    Nøttestad, L. et al. Prey selection of offshore killer whales Orcinus orca in the Northeast Atlantic in late summer: spatial associations with mackerel. Mar. Ecol. Prog. Ser. 499, 275–283 (2014).ADS 
    Article 

    Google Scholar 
    Bjørge, A., Aarefjord, H., Kaarstad, S., Kleivane, L. & Øien, N. Harbour porpoise (Phocoena phocoena) in Norwegian waters (Springer, 1991).
    Google Scholar 
    Gjøseter, H. et al. Fisken og Havet. https://doi.org/10.1111/maec.12075 (2010).Article 

    Google Scholar 
    ICES. ICES Report on Ocean Climate 2009 No.304. (2010).ICES. Report of the Working Group on Widely Distributed Stocks (WGWIDE). (2010).Rey, F. Phytoplankton: The grass of the sea. In The Norwegian Sea Ecosystem (ed. Skjoldal, H. R.) 97–136 (Academic Press, 2004).
    Google Scholar 
    Huse, G. et al. Effects of interactions between fish populations on ecosystem dynamics in the Norwegian Sea : Results of the INFERNO project. Mar. Biol. Res. 8, 415–419 (2012).Article 

    Google Scholar 
    Godø, O. R., Johnsen, S. & Torkelsen, T. The LoVe ocean observatory is in operation. Mar. Technol. Soc. J. 48, 24–30 (2014).Article 

    Google Scholar 
    Cooke, J. G. Balaenoptera physalus. The IUCN Red List of Threatened Species: e.T2478A50349982 (2018).Leonard, D. & Øien, N. Estimated abundances of Cetacean species in the Northeast Atlantic from Norwegian Shipboard Surveys Conducted in 2014–2018. NAMMCO Sci. Publ. 11, 4694 (2020).
    Google Scholar 
    Øygard, S. H. Simulations of Acoustic Transmission Loss of Fin Whale Calls Reaching the LoVe Ocean Observatory. (MSc thesis. University of Bergen, 2018).Steiner, L. et al. A link between male sperm whales, Physeter macrocephalus, of the Azores and Norway. J. Mar. Biol. Assoc. UK 92, 1751–1756 (2012).Article 

    Google Scholar 
    Olafsen, T., Winther, U., Olsen, Y. & Skjermo, J. Verdiskaping Basert på Produktive hav i 2050 1–76 (Springer, 2012).
    Google Scholar 
    Wenz, G. M. Acoustic ambient noise in the ocean: Spectra and sources. J. Acoust. Soc. Am. 34, 1936–1956 (1962).ADS 
    Article 

    Google Scholar 
    Klinck, H. et al. Seasonal presence of cetaceans and ambient noise levels in polar waters of the North Atlantic. J. Acoust. Soc. Am. 132, 176–181 (2012).Article 

    Google Scholar 
    Burnham, R. E., Duffus, D. A. & Mouy, X. The presence of large whale species in Clayoquot Sound and its offshore waters. Cont. Shelf Res. 177, 15–23 (2019).ADS 
    Article 

    Google Scholar 
    Romagosa, M. et al. Baleen whale acoustic presence and behaviour at a Mid-Atlantic migratory habitat, the Azores Archipelago. Sci. Rep. 10, 61489 (2020).
    Google Scholar 
    Tervo, O. Acoustic Behaviour of Bowhead Whales Balaena mysticetus in Disko Bay, Western Greenland. PhD thesis. (2011).Magnúsdóttir, E. E. & Lim, R. Subarctic singers: Humpback whale (Megaptera novaeangliae) song structure and progression from an Icelandic feeding ground during winter. PLoS ONE 14, e0210057 (2019).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Samaran, F. et al. Seasonal and geographic variation of southern blue whale subspecies in the Indian Ocean. PLoS ONE 8, e70 (2013).Article 

    Google Scholar 
    Norris, T. F., Dunleavy, K. J., Yack, T. M. & Ferguson, E. L. Estimation of minke whale abundance from an acoustic line transect survey of the Mariana Islands. Mar. Mammal Sci. 33, 574 (2017).Article 

    Google Scholar 
    Marques, T. A. et al. Estimating animal population density using passive acoustics. Biol. Rev. Camb. Philos. Soc. 88, 287–309 (2013).PubMed 
    Article 

    Google Scholar 
    Dunlop, R. A. The effects of vessel noise on the communication network of humpback whales. R. Soc. Open Sci. 6, 190967 (2019).ADS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Christensen, I., Haug, T. & Øien, N. A review of feeding and reproduction in large baleen whales (Mysticeti) and sperm whales Physeter macrocephalus in Norwegian and adjacent waters. ICES J. Mar. Sci. 49, 341–355 (1992).Article 

    Google Scholar 
    Aniceto, A. S. et al. Monitoring marine mammals using unmanned aerial vehicles: quantifying detection certainty. Ecosphere 9, e02122 (2018).Article 

    Google Scholar 
    Pedersen, G., Storheim, E., Sivle, L. D., Godø, O. R. & Ødegaard, L. A. Concurrent passive and active acoustic observations of high-latitude shallow foraging sperm whales (Physeter macrocephalus) and mesopelagic prey layer. J. Acoust. Soc. Am. 141, 1–10 (2017).Article 

    Google Scholar 
    Vogel, E. F. The influence of herring (Clupea harengus) biomass and distribution on killer whale (Orcinus orca) movements on the Norwegian shelf (UiT The Arctic University of Norway, 2020).
    Google Scholar 
    Williams, R. et al. Impacts of anthropogenic noise on marine life: Publication patterns, new discoveries, and future directions in research and management. Ocean Coast. Manag. 115, 17–24 (2015).Article 

    Google Scholar 
    Garibbo, S. et al. Low-frequency ocean acoustics: Measurements from the Lofoten-Vesterålen Ocean Observatory, Norway. (2020). https://doi.org/10.1121/2.0001324.Dekeling, R. P. A. et al. Monitoring Guidance for Underwater Noise in European Seas, Part I: Executive Summary (Springer, 2014).
    Google Scholar 
    Erbe, C. International regulation of underwater noise. Acoust. Aust. 41, 1–10 (2013).
    Google Scholar 
    Halliday, W. D., Insley, S. J., Hilliard, R. C., de Jong, T. & Pine, M. K. Potential impacts of shipping noise on marine mammals in the western Canadian Arctic. Mar. Pollut. Bull. 123, 73–82 (2017).CAS 
    PubMed 
    Article 

    Google Scholar 
    Halliday, W. D. et al. Underwater sound levels in the Canadian Arctic, 2014–2019. Mar. Pollut. Bull. 168, 112437 (2021).CAS 
    PubMed 
    Article 

    Google Scholar 
    Ødegaard, L., Pedersen, G. & Johnsen, E. Underwater Noise From Wind At the High North Love Ocean Observatory. UACE 2019 Conf. Proc. 359–366 (2019).Zhang, G., Forland, T. N., Johnsen, E., Pedersen, G. & Dong, H. Measurements of underwater noise radiated by commercial ships at a cabled ocean observatory. Mar. Pollut. Bull. 153, 110948 (2020).CAS 
    PubMed 
    Article 

    Google Scholar 
    Maystrenko, Y. P., Olesen, O., Gernigon, L. & Gradmann, S. Deep structure of the Lofoten–Vesterålen segment of the Mid-Norwegian continental margin and adjacent areas derived from 3-D density modeling. J. Geophys. Res. Solid Earth 122, 1402–1433 (2017).ADS 
    Article 

    Google Scholar 
    Gillespie, D. et al. PAMGUARD: Semiautomated, open source software for real-time acoustic detection and localization of cetaceans. J. Acoust. Soc. Am. 125, 2547–2547 (2009).ADS 
    Article 

    Google Scholar 
    Hollander, M. & Wolfe, D. A. Nonparametric Statistical Methods (Wiley, 1973).MATH 

    Google Scholar 
    Vogel, E. F. et al. Killer whale movements on the Norwegian shelf are associated with herring density. Mar. Ecol. Prog. Ser. 665, 217–231 (2021).ADS 
    Article 

    Google Scholar 
    Garcia, H. A. et al. Temporal-spatial, spectral, and source level distributions of fin whale vocalizations in the Norwegian Sea observed with a coherent hydrophone array. ICES J. Mar. Sci. 76, 268–283 (2019).Article 

    Google Scholar 
    Davis, G. E. et al. Exploring movement patterns and changing distributions of baleen whales in the western North Atlantic using a decade of passive acoustic data. Glob. Chang. Biol. 26, 4812–4840 (2020).ADS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Risch, D. et al. Minke whale acoustic behavior and multi-year seasonal and diel vocalization patterns in Massachusetts Bay, USA. Mar. Ecol. Prog. Ser. 489, 279–295 (2013).ADS 
    Article 

    Google Scholar 
    Le Tixerant, M., Le Guyader, D., Gourmelon, F. & Queffelec, B. How can Automatic Identification System (AIS) data be used for maritime spatial planning?. Ocean Coast. Manag. 166, 18–30 (2018).Article 

    Google Scholar 
    Team, R. C. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).Sumner, M. D. The Tag Location Problem. 133 (2011).Sumner, M. D., Wotherspoon, S. J. & Hindell, M. A. Bayesian estimation of animal movement from archival and satellite tags. PLoS ONE 4, e7324 (2009).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Halliday, W. D. et al. The coastal Arctic marine soundscape near Ulukhaktok, Northwest Territories, Canada. Polar Biol. 43, 623–636 (2020).Article 

    Google Scholar 
    Ezzet, F. & Pinheiro, J. Linear, generalized linear, and nonlinear mixed effects models. Pharm. Sci. Quant. Pharmacol. 1, 103–135. https://doi.org/10.1002/9780470087978.ch4 (2006).Article 

    Google Scholar 
    Mazerolle, M. J. AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c). (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 
    Book 

    Google Scholar 
    Pante, E. & Simon-Bouhet, B. marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS ONE 8, e73051 (2013).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Wessel, P. & Walter, H. F. S. A global self-consistent, hierarchical, high-resolution shoreline database. J. Geophys. Res. 101, 8741–8743 (1996).ADS 
    Article 

    Google Scholar 
    Sueur, J., Aubin, T. & Simonis, C. Equipment review: Seewave, a free modular tool for sound analysis and synthesis. Bioacoustics 18, 213–226 (2008).Article 

    Google Scholar 
    The Mathworks Inc. MATLAB (R2019a). (MathWorks Inc., 2019).Merchant, N. D. et al. Measuring acoustic habitats. Methods Ecol. Evol. 6, 257–265 (2015).PubMed 
    PubMed Central 
    Article 

    Google Scholar  More

  • in

    The dynamical complexity of seasonal soundscapes is governed by fish chorusing

    Data collectionThe acoustic recordings were collected during 2017 off the Changhua coast (24° 4.283 N/120° 19.102 E) (Fig. 5) by deploying a passive acoustic monitoring (PAM) device from Wildlife Acoustics, which was an SM3M recorder moored at a depth of 18–20 m. The hydrophone recorded continuously with a sampling frequency of 48 kHz and a sensitivity of −164.2 dB re:1 v/µPa. The acoustic files were recorded in the.WAV format with a duration of 60 minutes. The hydrophone setup prior to deployment is shown in Fig. 6. Table 2 contains the details for the monitoring period with the corresponding season and the number of hours of recordings each season used in this study. Studies have shown that the presence of seasonal chorusing at this monitoring site in the frequency range of 500–2500 Hz is caused by two types of chorusing15,38, with chorusing starting in early spring, reaching a peak in summer, and starting to decline late autumn, and silencing in winter6. Previous studies6,15,38 at this monitoring site have derived the details of two types of fish calls responsible for chorusing (Type 1 and Type 2); Supplementary Fig. 1 shows the spectrogram, waveform, and power spectrum density of the individual calls. Supplementary Table 1 tabulated are the acoustic features of the two call types. The monitoring region, Changhua, lies in the Eastern Taiwan Strait (ETS). The ETS is ~350 km in length and ~180 km wide64. The ETS experiences diverse oceanographic and climatic variations influenced by monsoons in summer and winter65 and extreme events caused by tropical storms, typhoons in summer, and wind/cold bursts occurring in winter66,67,68.Fig. 5: Study area located off the Taiwan Strait.Map of the Changhua coast located in Taiwan Strait, Taiwan depicting the deployed passive acoustic monitoring recorder at site A1. The map was produced in Matlab 9.11 (The Mathworks, Natick, MA; http://www.mathworks.com/) using mapping toolbox function geobasemap(). Full global basemap composed of high-resolution satellite imagery hosted by Esri (https://www.esri.com/).Full size imageFig. 6: Setup of the SM3M submersible recorder.SM3M recorder fastened to the steel frame (length and breadth = 1.22 m, height = 0.52 m) with plastic cable zip ties prior to deployment.Full size imageTable 2 Passive acoustic monitoring device specifications and monitoring duration during different seasons.Full size tableAcoustic data analysisThe acoustic data were analyzed using the PAMGuide toolbox in Matlab60. The seasonal spectrograms were computed with an FFT size of 1024 points and a 1 s time segment averaged to a 60 s resolution. The sound pressure levels (SPL) were computed in the frequency band of 500–3500 Hz and programmed to provide a single value every hour, thus resulting in 984, 1344, and 1440 data points in spring, summer, and winter, respectively (Supplementary Data 1).Determining the regularity and complexity with the complexity-entropy planeThe complexity-entropy plane was utilized in this study to quantify the structural statistical complexity and the regularity in the hourly acoustical and seasonal SPL time series data. The C-H plane is a 2D plane representation of the permutation entropy on the horizontal axis that quantifies the regularity, and the vertical axis is represented by the statistical complexity quantifying the correlation structure in the temporal series.For a given time series ({{x(t)}}_{t=1}^{N}), the N’ ≡ N − (m − 1) the values of the vectors for the length m  > 1 are ranked as$${X}_{s}=left({x}_{s-(m-1)},{x}_{s-(m-2)},ldots ,{x}_{s}right),s=1,ldots ,,{N}^{{prime} }$$
    (1)
    Within each vector, the values are reordered in the ascending order of their amplitude, yielding the set of ordering symbols (patterns) ({r}_{0},{r}_{1},ldots ,{r}_{m-1}) such that$${x}_{s-{r}_{0}}le {x}_{s-{r}_{1}}le ..,..le {x}_{s-{r}_{(m-1)}}$$
    (2)
    This symbolization scheme was introduced by Bandt and Pompe69. The scheme performs the local ordering of a time series to construct a probability mass function (PMF) of the ordinal patterns of the vector. The corresponding vectors (pi ={r}_{0},{r}_{1},ldots ,{r}_{(m-1)}) may presume any of the m! possible permutations of the set ({{{{{mathrm{0,1}}}}},ldots ,m-1}) and symbolically represent the original vector. For instance, for a given time series {9, 4, 5, 6, 1,…} with length m = 3, provides 3! different order patterns with six mutually exclusive permutation symbols are considered. The first three-dimensional vector is (9, 4, 5), following the Eq. (1), this vector corresponds to ((,{x}_{s-2},{x}_{s-1},{x}_{s})). According to Eq. (2), it yields ({x}_{s-1}le {x}_{s}le {x}_{s-2}). Then, the ordinal pattern satisfying the Eq. (2) will be (1, 0, 2). The second 3-dimensional vector is (4, 5, 6), and (2, 1, 0) will be its associated permutation, and so on.The permutation entropy (H) of order m ≥ 2 is defined as the Shannon entropy of the Brandt-Pompe probability distribution p(π)69$$Hleft(mright)=,-{mathop{sum}limits _{{pi }}}pleft(pi right){{{log }}}_{2}p(pi )$$
    (3)
    where ({pi }) represents the summation over all possible m! permutations of order m, (p(pi )) is the relative frequency of each permutation (pi), and the binary logarithm (base of 2) is evaluated to quantify the entropy in bits. H(m) attains the maximum ({{log }}m!) for (p(pi )=1/m!). Then the normalized Shannon entropy is given by$$0le H(m)/{{{{{rm{ln}}}}}},m!le 1$$
    (4)
    where the lower bound H = 0 corresponds to more predictable signals with fewer fluctuations, an either strictly increasing or decreasing series (with a single permutation), and the upper bound H = 1 corresponds to an unpredictable random series for which all the m! possible permutations are equiprobable. Thus, H quantifies the degree of disorder inherent in the time series. The choice of the pattern length m is essential for calculating the appropriate probability distribution, particularly for m, which determines the number of accessible states given by m!70,71. As a rule of thumb, the length T of the time series must satisfy the condition T (gg) m! in order to obtain reliable statistics, and for practical purposes, Bandt and Pompe suggested choosing m = 3,…,7 69.The statistical complexity measure is defined with the product form as a function of the Bandt and Pompe probability distribution P associated with the time series. (Cleft[Pright]) is represented as33$$Cleft[Pright]=frac{J[P,U]}{{J}_{{max }}}{H}_{s}[P]$$
    (5)
    where ({H}_{s}left[Pright]=Hleft[Pright]/{{log }}m!) is the normalized permutation entropy. (J[P,U]) is the Jensen divergence$$Jleft[P,Uright]=left{Hleft[frac{P+U}{2}right]-frac{H[P]}{2}-frac{H[U]}{2}right}$$
    (6)
    which quantifies the difference between the uniform distributions U and P, and ({J}_{{max }})is the maximum possible value of (Jleft[P,Uright]) that is obtained from one of the components of P = 1, with all the other components being zero:$$Jleft[P,Uright]=-frac{1}{2}left[frac{m!+1}{m!}{{log }}left(m!+1right)-2{{log }}left(2m!right)+{{log }}(m!)right]$$
    (7)
    For each value of the normalized permutation entropy (0le {H}_{s}[P]le 1) there is a corresponding range of possible statistical complexity (Cleft[Pright]) values. Thus, the upper (({C}_{{max }})) and lower ((C_{{min }})) complexity bounds in the C-H plane are formed. The periodic sequences such as sine and series with increasing and decreasing (with ({H}_{s}[P]=0)) and completely random series such as white noise (for which (Jleft[P,Uright]=0) and ({H}_{s}[P]=1)) will have zero complexity. Furthermore, for each given value of the (0le {H}_{s}[P]le 1), there exists a range of possible values of the statistical complexity, ({C}_{{min }}le C[P]le {C}_{{max }}). The procedure for evaluating the complexity bounds ({C}_{{min }}) and ({C}_{{max }}) is given in Martin et al.72. We utilized the R package ‘statcomp’73 to evaluate the statistical complexity (C) and the permutation entropy (H) using the command global-complexity() for the order m = 6, and the command limit_curves(m, fun = ‘min/max’) was utilized to evaluate the complexity boundaries ({C}_{{min }}) and ({C}_{{max }}). In this study, we constructed two C-H planes: (1) C and H was computed for each hourly acoustic file during different seasons. The resulting lengths of C and H during spring, summer, and autumn-winter are similar to the number of hours in the particular season (Table 2). (2) C and H was computed every 4–5 days for the seasonal SPL. The resulting length of C and H obtained during spring was 9 points (each value of C and H for every 109 h), and in summer and autumn-winter was 12 points (each value of C and H for every 112 and 120 h).Determining predictability and dynamics (linear/nonlinear) using EDMIn this study, we utilized EDM to quantify the predictability (forecasting) and distinguish between the linear stochastic and nonlinear dynamics in the seasonal soundscape SPL. EDM involves phase-space reconstruction via delay coordinate embeddings to make forecasts and to determine the ‘predictability portrait’ of time series data40. The first step in EDM is to determine the optimal embedding dimension (E), and this is obtained using a method based on simplex projection41. The simplex projection is carried out by dividing the dataset into two equal parts, of which the first part is called the library and the other part is called the target. The library set is used to build a series of non-parametric models (known as predictors) for the one step ahead predictions for the E varying between 1 and 10. Then the model’s accuracies are determined when the model is applied to the target dataset and the prediction skill (⍴) for the actual and predicted datasets is measured. The embedding dimension with the highest predictive skill is the optimal E.For the appropriate optimal E chosen, the predictability profile of the time series data is obtained by evaluating ⍴ at Tp = 1, 2, 3, … steps ahead. The flat prediction profile of the ⍴–Tp curve indicates that the time series is purely random (low ⍴) or regularly oscillating (high ⍴). In contrast, a decreasing ⍴ as Tp increases may reject the possibility of an underlying uncorrelated stochastic process and indicate the presence of low-dimensional deterministic dynamics. However, the concern with the predictability profile is that it may exhibit predictability even if time series are purely stochastic (such as autocorrelated red noise). Hence, a nonlinear test can be performed by using S-maps (sequential locally weighted global linear maps) to distinguish between linear stochastic and nonlinear dynamics in the time series dataset by fitting a local linear map. S-maps similar to simplex projects provide the forecasts in phase-space by quantifying the degree to which points are weighted when fitting the local linear map, which is given by the nonlinear localization parameter θ. When θ = 0, the entire library set will exhibit equal weights regardless of the target prediction, which mathematically resembles the model of a linear autoregressive process. In contrast, if θ  > 0, the forecasts of the library provided by the S-map depend on the local state of the target prediction, thus producing large weights, and the unique local fittings can vary in phase-space to incorporate nonlinear behavior. Consequently, if the (⍴–θ) dynamics profile shows the highest ⍴ at θ = 0 that is reduced as θ increases, it represents linear stochastic dynamics. If the ⍴ achieves the highest value at θ  > 0, then the dynamics are represented by nonlinear dynamics.In this study, the R package “rEDM”74 was used to evaluate the optimal E, prediction profile (⍴–Tp), and dynamics profile (⍴–θ) for the seasonal SPL dataset. While evaluating these entities, the data points are equally into two parts, and sequentially the first half is chosen as the library set and the other as the target set. The length of the library and the target set for spring, summer, and autumn-winter are 492, 672, and 720. The command EmbedDimension() was used to determine the forecast skill for the E ranging from 1 to 10 and the optimal E with the highest forecast skill (Supplementary Fig. 2) was chosen. In this study, we found that for all seasons, the optimal E was 2. The (⍴–Tp) was evaluated for Tp varying between 1 and 100 using the command PredictInterval() and the (⍴–θ) was evaluated using the command PredictNonlinear() for θ = 0, 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5,0.75, 1.0, 1.5, 2, and 3 to 20.StatisticsThe nonparametric Kruskal–Wallis test, followed by post hoc Bonferroni’s multiple comparisons, was used to test differences in the seasonal H and C that were obtained directly from the hourly acoustic data during chorusing hours, as well as the H and C obtained for the seasonal SPL and the seasonal forecast skill. The statistical calculations were performed using the R package “agricolae” 75. More

  • in

    Global forest management data for 2015 at a 100 m resolution

    Reference data collectionIn February 2019, we involved forest experts from different regions around the world and organized a workshop to (1) discuss the variety of forest management practices that take place in various parts of the world; (2) explore what types of forest management information could be collected by visual interpretation of very high-resolution images from Google Maps and Microsoft Bing Maps, in combination with Sentinel time series and Normalized Difference Vegetation Index (NDVI) profiles derived from Google Earth Engine (GEE); (3) generalize and harmonize the definitions at global scale; (4) finalize the Geo-Wiki interface for the crowdsourcing campaigns; and (5) build a data set of control points (or the expert data set), which we used later to monitor the quality of the crowdsourced contributions by the participants. Based on the results of this analysis, we launched the crowdsourcing campaigns by involving a broader group of participants, which included people recruited from remote sensing, geography and forest research institutes and universities. After the crowdsourcing campaigns, we collected additional data with the help of experts. Hence, the final reference data consists of two parts: (1) a randomly stratified sample collected by crowdsourcing (49,982 locations); (2) a targeted sample collected by experts (176,340 locations, at those locations where the information collected from the crowdsourcing campaign was not large enough to ensure a robust classification).DefinitionsTable 1 contains the initial classification used for visual interpretation of the reference samples and the aggregated classes presented in the final reference data set. For the Geo-Wiki campaigns, we attempted to collect information (1) related to forest management practices and (2) recognizable from very high-resolution satellite imagery or time series of vegetation indices. The final reference data set and the final map contain an aggregation of classes, i.e., only those that were reliably distinguishable from visual interpretation of satellite imagery.Table 1 Forest management classes and definitions.Full size tableSampling design for the crowdsourcing campaignsInitially, we generated a random stratified sample of 110,000 sites globally. The total number of sample sites was chosen based on experiences from past Geo-Wiki campaigns12, a practical estimation of the potential number of volunteer participants that we could engage in the campaign, and the expected spatial variation in forest management. We used two spatial data sets for the stratification of the sample: World Wildlife Fund (WWF) Terrestrial Ecoregions13 and Global Forest Change14. The samples were stratified into three biomes, based on WWF Terrestrial Ecoregions (Fig. 2): boreal (25 000 sample sites), temperate (35,000 sample sites) and tropical (50,000 sample sites). Within each biome, we used Hansen’s14 Global Forest Change maps to derive areas with “forest remaining forest” 2000–2015, “forest loss or gain”, and “permanent non-forest” areas.Fig. 2Biomes for sampling stratification (1 – boreal, 2 – temperate, 3 – sub-tropical and tropical).Full size imageThe sample size was determined from previous experiences, taking into account the expected spatial variation in forest management within each biome. Tropical forests had the largest sample size because of increasing commodity-driven deforestation15, the wide spatial extent of plantations, and slash and burn agriculture. Temperate forests had a larger sample compared to boreal forests due to their higher fragmentation. Each sample site was classified by at least three different participants, thus accounting for human error and varying expertise16,17,18. At a later stage, following a preliminary analysis of the data collected, we increased the number of sample sites to meet certain accuracy thresholds for every mapped class (aiming to exceed 75% accuracy).The Geo‐Wiki applicationGeo‐Wiki.org is an online application for crowdsourcing and expert visual interpretation of satellite imagery, e.g., to classify land cover and land use. This application has been used in several data collection campaigns over the last decade16,19,20,21,22,23. Here, we implemented a new custom branch of Geo‐Wiki (‘Human impact on Forest’), which is devoted to the collection of forest management data (Fig. 3). Various map overlays (including satellite images from Google Maps, Microsoft Bing Maps and Sentinel 2), campaign statistics and tools to aid interpretation, such as time series profiles of NDVI, were provided as part of this Geo‐Wiki branch, giving users a range of options and choices to facilitate image classification and general data collection. Google Maps and Microsoft Bing Maps include mosaics of very high-resolution satellite and aerial imagery from different time periods and multiple image providers, including the Landsat satellites operated by NASA and USGS as base imagery to commercial image providers such as Digital Globe. More information on the spatial and temporal distribution of very high-resolution satellite imagery can be found in Lesiv et al.24. This collection of images was supplied as guidance for visual interpretation16,20. Participants could analyze time series profiles of NDVI from Landsat, Sentinel 2 and MODIS images, which were derived from Google Earth Engine (GEE). More information on tools can be found in Supplementary file 1.Fig. 3Screenshot of the Geo‐Wiki interface showing a very high-resolution image from Google Maps and a sample site as a 100 mx100 m blue square, which the participants classified based on the forest management classes on the right.Full size imageThe blue box in Fig. 3 corresponds to 100 m × 100 m pixels aligned with the Sentinel grid in UTM projection. It is the same geometry required for the classification workflow that is used to produce the Copernicus Land Cover product for 201511.Before starting the campaign, the participants were shown a series of slides designed to help them gain familiarity with the interface and to train them in how to visually determine and select the most appropriate type of land use and forest management classes at each given location, thereby increasing both consistency and accuracy of the labelling tasks among experts. Once completed, the participants were shown random locations (from the random stratified sample) on the Geo‐Wiki interface and were then asked to select one of the forest management classes outlined in the Definition section (see Table 1 above).Alternatively, if there was either insufficient quality in the available imagery, or if a participant was unable to determine the forest management type, they could skip such a site (Fig. 3). If a participant skipped a sample site because it was too difficult, other participants would then receive this sample site for classification, whereas in the case of the absence of high-resolution satellite imagery, i.e., Google Maps and Microsoft Bing Maps, this sample site was then removed from the pool of available sample sites. The skipped locations were less than 1% of the total amount of locations assigned for labeling. Table 2 shows the distribution of the skipped locations by countries, based on the subset of the crowdsourced data where all the participants agreed.Table 2 Distribution of the skipped locations by countries.Full size tableQuality assurance and data aggregation of the crowdsourced dataBased on the experience gained from previous crowdsourcing campaigns12,19, we invested in the training of the participants (130 persons in total) and overall quality assurance. Specifically, we provided initial guidelines for the participants in the form of a video and a presentation that were shown before the participants could start classifying in the forest management branch (Supplementary file 1). Additionally, the participants were asked to classify 20 training samples before contributing to the campaign. For each of these training samples, they received text‐based feedback regarding how each location should be classified. Summary information about the participants who filled in the survey at the end of the campaign (i.e., gender, age, level of education, and their country of residence) is provided in the Supplementary file 2. We would like to note that 130 participants is a high number, especially taking the complexity of the task into consideration.Furthermore, during the campaign, sample sites that were part of the “control” data set were randomly shown to the participants. The participants received text-based feedback regarding whether the classification had been made correctly or not, with additional information and guidance. By providing immediate feedback, our intention was that participants would learn from their mistakes, increasing the quality and classification accuracy over time. If the text‐based feedback was not sufficient to provide an understanding of the correct classification, the participants were able to submit a request (“Ask the expert”) for a more detailed explanation by email.The control set was independent of the main sample, and it was created using the same random stratified sampling procedure within each biome and the stratification by Global Forest Change maps14 (see “Sample design” section). To determine the size of the control sample, we considered two aspects: (a) the maximum number of sample sites that one person could classify during the entire campaign; (b) the frequency at which control sites would appear among the task sites (defined at 15%, which is a compromise between the classification of as many unknown locations as possible and a sufficient level of quality control, based on previous experience). Our control sample consisted of 5,000 sites. Each control sample site was classified twice by two different experts. When the two experts agreed, these sample sites were added to the final control sample. Where disagreement occurred (in 25% of cases), these sample sites were checked again by the experts and revised accordingly. During the campaign, participants had the option to disagree with the classification of the control site and submit a request with their opinion and arguments. They received an additional quality score in the situation when they were correct, but the experts were not. This procedure also ensured an increase in the quality of the control data set.To incentivize participation and high-quality classifications, we offered prizes as part of the campaign design. The ranking system for the prize competition considered both the quality of the classifications and the number of classifications provided by a participant. The quality measure was based on the control sample discussed above. The participants randomly received a control point, which was classified in advance by the experts. For every control point, a participant could receive a maximum of +30 points (fully correct classification) to a minimum of −30 points (incorrect classification). In the case where the answer was partly correct (e.g., the participant correctly classified that the forest is managed, but misclassified the regeneration type), they received points ranging from 5 to 25.The relative quality score for each participant was then calculated as the total sum of gained points divided by the maximum sum of points that this participant could have earned. For any subsequent data analysis, we excluded classifications from those participants whose relative quality score was less than 70%. This threshold corresponds to an average score of 10 points at each location (out of a maximum of 30 points), i.e., where participants were good at defining the aggregated forest management type but may have been less good at providing the more detailed classification.Unfortunately, we observed some imbalance in the proportion of participants coming from different countries, e.g. there were not so many participants from the tropics. This could have resulted in interpretation errors, even when all the participants agreed on a classification. To address this, we did an additional quality check. We selected only those sample sites where all the participants agreed and then randomly checked 100 sample sites from each class. Table 3 summarizes the results of this check and explains the selection of the final classes presented in Table 1.Table 3 Qualitative analysis of the reference sample sites with full agreement.Full size tableAs a result of the actions outlined in Table 3, we compiled the final reference data set, which consisted of 49,982 consistent sample sites.Additional expert data collectionWe used the reference data set to produce a test map of forest management (the classification algorithm used is described in the next section). By checking visually and comparing against the control data set, we found that the map was of insufficient quality for many locations, especially in the case of heterogeneous landscapes. While several reasons for such an unsatisfactory result are possible, the experts agreed that a larger sample size would likely increase the accuracy of the final map, especially in areas of high heterogeneity and for forest management classes that only cover a small spatial extent. To increase the amount of high-quality training data and hence to improve the map, we collected additional data using a targeted approach. In practice, the map was uploaded to Geo-Wiki, and using the embedded drawing tools, the experts randomly checked locations on the map, focusing on their region of expertise and added classified polygons in locations where the forest management was misclassified. To limit model overfitting and oversampling of certain classes, the experts also added points for correctly mapped classes to keep the density of the points the same. This process involved a few iterations of collecting additional points and training the classification algorithm until the map accuracy reached 75%. In total, we collected an additional 176,340 training points. With the 49,982 consistent training points from the Geo-Wiki campaigns, this resulted in 226,322 (Fig. 4). This two-pronged approach would not have been possible without the exhaustive knowledge obtained from running the initial Geo-Wiki campaigns, including numerous questions raised by the campaign participants. Figure 4 also highlights in yellow the areas of very high sampling density, I.e., those collected by the experts. The sampling intensity of these areas is much higher in comparison with the randomly distributed crowdsourced locations, and these are mainly areas with very mixed forest classes or small patches, in most cases, including plantations.Fig. 4Distribution of reference locations.Full size imageClassification algorithmTo produce the forest management map for the year 2015, we applied a workflow that was developed as part of the production of the Copernicus Global Land Services land cover at 100 m resolution (CGLS-LC100) collection 2 product11. A brief description of the workflow (Fig. 5), focusing on the implemented changes, is given below. A more thorough explanation, including detailed technical descriptions of the algorithms, the ancillary data used, and the intermediate products generated, can be found in the Algorithm Theoretical Basis Document (ATBD) of the CGLS-LC100 collection 2 product25.Fig. 5Workflow overview for the generation of the Copernicus Global Land Cover Layers. Adapted from the Algorithm Theoretical Basis Document25.Full size imageThe CGLS-LC100 collection 2 processing workflow can be applied to any satellite data, as it is unspecific to different sensors or resolutions. While the CGLS-LC100 Collection 2 product is based on PROBA-V sensor data, the workflow has already been tested with Sentinel 2 and Landsat data, thereby using it for regional/continental land cover (LC) mapping applications11,26. For generating the forest management layer, the main Earth Observation (EO) input was the PROBA-V UTM Analysis Ready Data (ARD) archive based on the complete PROBA-V L1C archive from 2014 to 2016. The ARD pre-processing included geometric transformation into a UTM coordinate system, which reduced distortions in high northern latitudes, as well as improved atmospheric correction, which converted the Top-of-Atmosphere reflectance to surface reflectance (Top-of-Canopy). In a further processing step, gaps in the 5-daily PROBA-V UTM multi-spectral image data with a Ground Sampling Distance (GSD) of ~0.001 degrees (~100 m) were filled using the PROBA-V UTM daily multi-spectral image data with a GSD of ~0.003 degrees (~300 m). This data fusion is based on a Kalman filtering approach, as in Sedano et al.27, but was further adapted to heterogonous surfaces25. Outputs from the EO pre-processing were temporally cleaned by using the internal quality flags of the PROBA-V UTM L3 data, a temporal cloud and outlier filter built on a Fourier transformation. This was done to produce consistent and dense 5-daily image stacks for all global land masses at 100 m resolution and a quality indicator, called the Data Density Indicator (DDI), used in the supervised learning process of the algorithm.Since the total time series stack for the epoch 2015 (a three-year period including the reference year 2015 +/− 1 year) would be composed of too many proxies for supervised learning, the time and spectral dimension of the data stack had to be condensed. The spectral domain was condensed by using Vegetation Indices (VIs) instead of the original reflectance values. Overall, ten VIs based on the four PROBA-V reflectance bands were generated, which included: Normalized Difference Vegetation Index (NDVI); Enhanced Vegetation Index (EVI); Structure Intensive Pigment Index (SIPI); Normalized Difference Moisture Index (NDMI); Near-Infrared reflectance of vegetation (NIRv); Angle at NIR; HUE and VALUE of the Hue Saturation Value (HSV) color system transformation. The temporal domain of the time series VI stacks was then condensed by extracting metrics, which are used as general descriptors to enable distinguishing between the different LC classes. Overall, we extracted 266 temporal, descriptive, and textual metrics from the VI times series stacks. The temporal descriptors were derived through a harmonic model, fitted through the time series of each of the VIs based on a Fourier transformation28,29. In addition to the seven parameters of the harmonic model that describe the overall level and seasonality of the VI time series, 11 descriptive statistics (mean, standard deviation, minimum, maximum, sum, median, 10th percentile, 90th percentile, 10th – 90th percentile range, time step of the first minimum appearance, and time step of the first maximum appearance) and one textural metric (median variation of the center pixel to median of the neighbours) were generated for each VI. Additionally, the elevation, slope, aspect, and purity derived at 100 m from a Digital Elevation Model (DEM) were added. Overall, 270 metrics were extracted from the PROBA-V UTM 2015 epoch.The main difference to the original CGLS-LC100 collection 2 algorithms is the use of forest management training data instead of the global LC reference data set, as well as only using the discrete classification branch of the algorithm. The dedicated regressor branch of the CGLS-LC100 collection 2 algorithm, i.e., outputting cover fraction maps for all LC classes, was not needed for generating the forest management layer.In order to adapt the classification algorithm to sub-continental and continental patterns, the classification of the data was carried out per biome cluster, with the 73 biome clusters defined by the combination of several global ecological layers, which include the ecoregions 2017 dataset30, the Geiger-Koeppen dataset31, the global FAO eco-regions dataset32, a global tree-line layer33, the Sentinel-2 tiling grid and the PROBA-V imaging extent;30,31 this, effectively, resulted in the creation of 73 classification models, each with its non-overlapping geographic extent and its own training dataset. Next, in preparation for the classification procedure, the metrics of all training points were analyzed for outliers, as well as screened via an all-relevant feature selection approach for the best metric combinations (i.e., best band selection) for each biome cluster in order to reduce redundancy between parameters used in the classification. The best metrics are defined as those that have the highest separability compared to other metrics. For each metric, the separability is calculated by comparing the metric values of one class to the metric values of another class; more details can be found in the ATBD25. The optimized training data set, together with the quality indicator of the input data (DDI data set) as a weight factor, were used in the training of the Random Forest classifier. Moreover, a 5-fold cross-validation was used to optimize the classifier parameters for each generated model (one per biome).Finally, the Random Forest classification was used to produce a hard classification, showing the discrete class for each pixel, as well as the predicted class probability. In the last step, the discrete classification results (now called the forest management map) are modified by the CGLS-LC100 collection 2 tree cover fraction layer29. Therefore, the tree cover fraction layer, showing the relative distribution of trees within one pixel, was used to remove areas with less than 10% tree cover fraction in the forest management layer, following the FAO definition of forest. Figure 6 shows the class probability layer that illustrates the model behavior, highlighting the areas of class confusion. This layer shows that there is high confusion between forest management classes in heterogeneous landscapes, e.g., in Europe and the Tropics while homogenous landscapes, such as Boreal forests, are mapped with high confidence. It is important to note that a low probability does not mean that the classification is wrong.Fig. 6The predicted class probability by the Random Forest classification.Full size image More