Philippa Kaur

AffiliationsSanger Institute, Wellcome Trust Genome Campus, Hinxton, UKPhysilia Ying Shi ChuaLaboratory of Genomics and Molecular Medicine, Department of Biology, University of Copenhagen, Copenhagen, DenmarkJacob Agerbo RasmussenCenter for Evolutionary Hologenomics, Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, DenmarkJacob Agerbo RasmussenAuthorsPhysilia Ying Shi ChuaJacob Agerbo RasmussenCorresponding authorCorrespondence to
Physilia Ying Shi Chua. More

We selected eleven major wheat production provinces of China for the study area, which comprise the largest winter wheat-sowing fraction: Henan, Shandong, Anhui, Jiangsu, Hebei, Hubei, Shanxi, Shaanxi, Sichuan, Xinjiang, and Gansu (Fig. 1). The wheat planting area is about 22 million ha in these provinces, accounting for more than 93% of the total wheat planting area. The total wheat production in these regions contributes more than 96% of the total wheat production in China, with more than 128 million tons in 201933.We developed a methodological framework for high-resolution AGB mapping. It mainly includes three parts: (1) Data acquisition and processing. (2) The WOFOST model parameterization and calibration. (3) Data assimilation (Fig. 1). Each part is explained in more detail below.Data acquisition and processingMeteorological dataChina Meteorological Forcing Dataset34,35 is used as weather driving data for the WOFOST model. The dataset is based on the internationally existing Princeton reanalysis data, Global Land Data Assimilation System data, Global Energy and Water Cycle Experiment-Surface Radiation Budget radiation data, and Tropical Rainfall Measuring Mission precipitation data. It is made by fusing the conventional meteorological observation data of the China Meteorological Administration. It includes seven elements: near-surface air temperature, air pressure, near-surface total humidity, wind speed, ground downward shortwave radiation, ground downward longwave radiation, and ground precipitation rate. The meteorological drive elements required for WOFOST are daily radiation, minimum temperature, maximum temperature, water vapor pressure, average wind speed, and precipitation. Details of these variables that participated in the WOFOST model can be referred to in Table S1.Soil characteristics measurements and phenology observationsSoil and phenology data were collected at 177 agricultural meteorological stations (AMS) from 2007 to 2015 (Fig. 1). Soil characteristics include soil moisture content at wilting points, field capacity, and saturation. To be consistent with the corresponding units in the crop model, the original data in weight water content was converted into volume water content through the corresponding soil bulk density measurements. Winter wheat phenology observations include the date of emergence (more than 50% of the wheat seedlings in the field show the first green leaves and reached about 2 cm), anthesis (the inner and outer glumes of the middle and upper florets of more than 50% of the wheat ears in the whole field are open, and the anthers loose powder), and maturity (more than 80% of the wheat grains turn yellow, the glumes and stems turn yellow, and only the upper first and second nodes are still slightly green). In most cases, the phenological stage “anthesis” is missing. The anthesis date was calculated by adding seven days to the observed heading date (when more than 50% of the wheat in the whole field exposes the tip of the ear from the sheath of the flag leaf).County-level yield statistics dataThe county-level yield data was collected from city statistical yearbooks of the study area from 2007 to 2015. Since most statistical yearbooks do not directly record per-unit yield data, the county-level yield was obtained by dividing the total yield and planting area. It is worth noting that all yields were calculated in units of metric kilograms per cultivated hectares (kg·ha−1).The winter wheat land cover dataWe used a winter wheat land cover product from a 1 km resolution dataset named ChinaCropArea1km36. This data was derived from GLASS leaf area index products and crop phenology from 2000 to 2015. This dataset is the base map of our data production.The MODIS LAI dataWe used the improved 8-days MODIS LAI products (i.e., 1 km) generated by Yuan et al.32 to assimilate the WOFOST model. The products applied the modified temporal-spatial filter and Savitzky-Golay filter to overcome the spatial-temporal discontinuity and inconsistence of raw MODIS LAI products, which makes them more applicable for the realm of land surface and climate modeling. The products can be accessed via the Land-Atmosphere Interaction Research Group website at Sun Yat-sen University (http://globalchange.bnu.edu.cn/research/lai).The WOFOST model parameterization and calibrationThe WOFOST model introductionThe WOFOST model was initially developed as a crop growth simulation model to evaluate the yield potential of various crops in tropical countries37. In this study, we chose the WOFOST model because the model reaches a trade-off of the complexity of the crop model and is suitable for large-scale simulations3. The WOFOST model is a typical crop growth model that explains crop growth based on underlying processes such as photosynthesis and respiration and their response to changing environmental conditions38. The WOFOST model estimates phenology, LAI, aboveground biomass, and storage organ biomass (i.e., grain yield) at a daily time step39 (Fig. 2).Fig. 2Schematic overview of the major processes implemented in WOFOST. The Astronomical module calculates day length, some variables relating to solar elevation, and the fraction of diffuse radiation.Full size imageZonal parameterizationWe first divided the study area covered by AMS into seamless Thiessen polygon zones. Each Thiessen polygon contains only a single AMS. These zones represent the whole areas where any location is closer to its associated AMS point than any other AMS point. Then, we assigned parameters to the entire zone based on the AMS data, including crop calendar (date of emergence) and soil water retention parameters (soil moisture content at wilting point, field capacity, and saturation). Besides, we also optimized two main crop parameters for controlling phenological development stages, namely TSUM1 (accumulated temperature required from emergence to anthesis) and TSUM2 (accumulated temperature required from anthesis to maturity), by minimizing the cost function of the observational and simulated date corresponding to anthesis and maturity.Parameter calibration within a single zoneWe implemented the calibration of parameters within every single zone, as illustrated in Fig. 3. We calculated the average statistical yield of each county within every single zone from 2007 to 2015, then ranked the counties in descending order and divided them into three groups, namely high, medium, and low-level yield counties, by the 33% quantile and 67% quantile of the average statistical yield. The three counties corresponding to 17% quantile, 50% quantile, and 83% quantile would be used for subsequent calibration and represent the corresponding three yield level groups. We used the statistical yields (converted to dry matter mass based on the standard moisture content of 12.5%) of the three counties for multiple years and a harvest index for each province to convert the county-level yield to AGB for calibration. The harvest index of each province was mainly estimated from the AMS’s dynamic growth records on the biomass composition of the dominant winter wheat varieties of the province and a published literature40. Besides, we collected the maximum LAI observations on all agrometeorological stations in all years in the study area, according to its histogram. We found that the histogram follows a normal distribution with a mean of 6.5 and a standard deviation of 1.5. Finally, we calibrated three sets of parameters corresponding to three yield level groups in each single zone according to the three selected counties.Fig. 3Flow chart of parameter calibration within a single zone.Full size imageWe designed a three-step calibration strategy for a specific yield level group. Firstly, as winter wheat varieties did not change significantly according to information recorded by agrometeorological stations from 2007 to 2015, we assumed the crop parameters of winter wheat remain unchanged every three years to combine three years of observational data to calibrate the parameters of the WOFOST model better. We maximized a log-likelihood function based on the maximum LAI statistics and every three-year county-level yield and AGB data mentioned to optimize selected crop parameters (see Table S2 in the Supplement Materials).The log-likelihood function was constructed as follows:$$log;{{rm{L}}}_{{rm{LAI}}}=-frac{1}{2}left[dlogleft(2pi right)+logleft(left|{Sigma }_{{rm{LAI}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{LAI}}};{mu }_{{rm{LAI}}},{Sigma }_{{rm{LAI}}}right)}^{2}right]$$
(1)
$$log;{{rm{L}}}_{{rm{TWSO}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{TWSO}}};{{boldsymbol{mu }}}_{{rm{TWSO}}},{{boldsymbol{Sigma }}}_{{rm{TWSO}}}right)}^{2}right]$$
(2)
$$log;{{rm{L}}}_{{rm{AGB}}}=-frac{1}{2}left[dlog(2pi )+logleft(left|{{boldsymbol{Sigma }}}_{{rm{AGB}}}right|right)+{rm{MD}}{left({{bf{x}}}_{{rm{AGB}}};{{boldsymbol{mu }}}_{{rm{AGB}}},{{boldsymbol{Sigma }}}_{{rm{AGB}}}right)}^{2}right]$$
(3)
$$log;{rm{L}}=log;{L}_{{rm{LAI}}}+log;{L}_{{rm{TWSO}}}+log;{L}_{{rm{AGB}}}$$
(4)
Where log L is the natural logarithm of the likelihood function, d is the dimension, that is, the number of years of joint calibration, which is set to 3 in this study xLAI is the vector composed of the maximum value of the 3-year LAI simulated by WOFOST, μLAI and ΣLAI are the mean vector and error covariance matrix of maximum LAI based on observation statistics. The annual maximum LAI was assumed to be independent, and the mean and standard deviation for each year was set the same as the result of Fig. 3. Similarly, xTWSO and xAGB are the yield vector and AGB vector at maturity of 3 years simulated by WOFOST, and μTWSO, μAGB are their corresponding county-level statistic vector, ΣTWSO and ΣAGB are their corresponding error covariance matrix. In this study, we assumed that the annual yield or AGB was independent, and their corresponding standard deviation was 10% of their statistical value. |Σ| is the determinant of Σ. The expression ({rm{MD}}{({bf{x}};{boldsymbol{mu }},{boldsymbol{Sigma }})}^{2}={({bf{x}}-{boldsymbol{mu }})}^{{rm{T}}}{{boldsymbol{Sigma }}}^{-1}({bf{x}}-{boldsymbol{mu }})), where MD is the Mahalanobis distance between the point x and the mean vector μ.Secondly, we optimized the inter-annual irrigation. We optimized two parameters every year: the critical value of soil moisture (denoted as SMc) and the amount of irrigation (denoted as V). When the soil moisture simulated by WOFOST is lower than SMc, an irrigation event will be triggered, and the irrigation amount is V cm. In this study, we combined three-year data for calibration with six parameters for optimization. The optimization strategy is the same as the previous step by maximizing the log-likelihood function. Finally, we fixed the optimized irrigation parameters and repeated the first step to calibrate the selected crop parameters and obtain the final optimal parameters.Data assimilationConsidering that MODIS LAI is relatively low compared to the actual LAI of winter wheat41, we select a weak-constraint cost function based on the least square of normalized observational and simulated LAI as shown in Eq. (5), which is assimilating the trend information of MODIS LAI into the crop growth model.$$J={sum }_{{rm{t}}=1}^{{rm{n}}}{left(frac{{{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}{{{rm{LAI}}}_{{rm{MODIS}}}^{max}-{{rm{LAI}}}_{{rm{MODIS}}}^{min}}-frac{{{rm{LAI}}}_{{rm{WOFOS}}}^{{rm{t}}}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}{{{rm{LAI}}}_{{rm{WOFOS}}}^{max}-{{rm{LAI}}}_{{rm{WOFOS}}}^{min}}right)}^{2}$$
(5)
Where ({{rm{LAI}}}_{{rm{MODIS}}}^{{rm{t}}}) and .. are MODIS LAI and WOFOST simulated LAI of time t. ({{rm{LAI}}}_{{rm{MODIS}}}^{max}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{max}) are maximum of MODIS LAI and WOFOST simulated LAI. ({{rm{LAI}}}_{{rm{MODIS}}}^{min}) and ({{rm{LAI}}}_{{rm{WOFOS}}}^{min}) are minimum of MODIS LAI and WOFOST simulated LAI. J is the value of the cost function.We reinitialize the day of emergence (IDEM), the life span of leaves growing at 35 °C (SPAN), and thermal time from emergence to anthesis (TSUM1) in the WOFOST model on each 1 km winter wheat pixel according to cost function between WOFOST LAI and MODIS LAI. Besides, we applied the Subplex algorithm from the NLOPT library (https://github.com/stevengj/nlopt) for parameter optimization. More

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

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

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

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

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

Ricklefs, R. E. A comprehensive framework for global patterns in biodiversity. Ecol. Lett. 7, 1–15 (2004).Article 

Google Scholar 
Dolson, S. J. et al. Diversity and phylogenetic community structure across elevation during climate change in a family of hyperdiverse neotropical beetles (Staphylinidae). Ecography 44, 740–752 (2021).Article 

Google Scholar 
Montaño-Centellas, F. A., McCain, C. & Loiselle, B. A. Using functional and phylogenetic diversity to infer avian community assembly along elevational gradients. Glob. Ecol. Biogeogr. 29, 232–245 (2020).Article 

Google Scholar 
Wiens, J. J. et al. Niche conservatism as an emerging principle in ecology and conservation biology. Ecol. Lett. 13, 1310–1324 (2010).PubMed 
Article 

Google Scholar 
Cavender-Bares, J., Kozak, K. H., Fine, P. V. A. & Kembel, S. W. The merging of community ecology and phylogenetic biology. Ecol. Lett. 12, 693–715 (2009).PubMed 
Article 

Google Scholar 
Mayfield, M. M. & Levine, J. M. Opposing effects of competitive exclusion on the phylogenetic structure of communities. Ecol. Lett. 13, 1085–1093 (2010).PubMed 
Article 

Google Scholar 
Webb, C. O., Ackerly, D. D., McPeek, M. A. & Donoghue, M. J. Phylogenies and community ecology. Annu. Rev. Ecol. Syst. 33, 475–505 (2002).Article 

Google Scholar 
Hubbell, S. P. The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32) Vol. 32 (Princeton University Press, 2001).
Google Scholar 
Kraft, N. J. B., Cornwell, W. K., Webb, C. O. & Ackerly, D. D. Trait evolution, community assembly, and the phylogenetic structure of ecological communities. Am. Nat. 170, 271–283 (2007).PubMed 
Article 

Google Scholar 
Cadotte, M. W. & Tucker, C. M. Should environmental filtering be abandoned?. Trends Ecol. Evol. 32, 429–437 (2017).PubMed 
Article 

Google Scholar 
Mouchet, M. A. et al. Functional diversity measures: An overview of their redundancy and their ability to discriminate community assembly rules. Funct. Ecol. 24, 867–876 (2010).Article 

Google Scholar 
Graham, C. H. & Fine, P. V. A. Phylogenetic beta diversity: Linking ecological and evolutionary processes across space in time. Ecol. Lett. 11, 1265–1277 (2008).PubMed 
Article 

Google Scholar 
Qian, H., Jin, Y., Leprieur, F., Wang, X. & Deng, T. Geographic patterns and environmental correlates of taxonomic and phylogenetic beta diversity for large-scale angiosperm assemblages in China. Ecography 43, 1706–1716 (2020).Article 

Google Scholar 
Swenson, N. G. et al. Phylogenetic and functional alpha and beta diversity in temperate and tropical tree communities. Ecology 93, 112–125 (2012).Article 

Google Scholar 
Qian, H., Hao, Z. & Zhang, J. Phylogenetic structure and phylogenetic diversity of angiosperm assemblages in forests along an elevational gradient in Changbaishan, China. J. Plant Ecol. 7, 154–165 (2014).Article 

Google Scholar 
Chase, J. M. & Myers, J. A. Disentangling the importance of ecological niches from stochastic processes across scales. Philos. Trans. R. Soc. B Biol. Sci. 366, 2351–2363 (2011).Article 

Google Scholar 
Leibold, M. A., Economo, E. P. & Peres-Neto, P. Metacommunity phylogenetics: Separating the roles of environmental filters and historical biogeography. Ecol. Lett. 13, 1290–1299 (2010).PubMed 
Article 

Google Scholar 
Ricklefs, R. E. Evolutionary diversification and the origin of the diversity-environment relationship. Ecology 87, 3–13 (2006).Article 

Google Scholar 
Zhang, J. L. et al. Phylogenetic beta diversity in tropical forests: Implications for the roles of geographical and environmental distance. J. Syst. Evol. 51, 71–85 (2013).Article 

Google Scholar 
Baselga, A. The relationship between species replacement, dissimilarity derived from nestedness, and nestedness. Glob. Ecol. Biogeogr. 21, 1223–1232 (2012).Article 

Google Scholar 
Leprieur, F. et al. Quantifying phylogenetic beta diversity: Distinguishing between ‘true’ turnover of lineages and phylogenetic diversity gradients. PLoS ONE https://doi.org/10.1371/journal.pone.0042760 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Bishop, T. R., Robertson, M. P., van Rensburg, B. J. & Parr, C. L. Contrasting species and functional beta diversity in montane ant assemblages. J. Biogeogr. 42, 1776–1786 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Economo, E. P., Narula, N., Friedman, N. R., Weiser, M. D. & Guénard, B. Macroecology and macroevolution of the latitudinal diversity gradient in ants. Nat. Commun. 9, 1–8 (2018).CAS 
Article 

Google Scholar 
Lessard, J. P., Fordyce, J. A., Gotelli, N. J. & Sanders, N. J. Invasive ants alter the phylogenetic structure of ant communities. Ecology 90, 2664–2669 (2009).PubMed 
Article 

Google Scholar 
Liu, C., Dudley, K. L., Xu, Z. H. & Economo, E. P. Mountain metacommunities: climate and spatial connectivity shape ant diversity in a complex landscape. Ecography 41, 101–112 (2018).Article 

Google Scholar 
Smith, M. A., Hallwachs, W. & Janzen, D. H. Diversity and phylogenetic community structure of ants along a Costa Rican elevational gradient. Ecography 37, 720–731 (2014).Article 

Google Scholar 
Machac, A., Janda, M., Dunn, R. R. & Sanders, N. J. Elevational gradients in phylogenetic structure of ant communities reveal the interplay of biotic and abiotic constraints on diversity. Ecography 34, 364–371 (2011).Article 

Google Scholar 
Guo, Q. et al. Global variation in elevational diversity patterns. Sci. Rep. 3, 1 (2013).CAS 

Google Scholar 
Kluge, J., Kessler, M. & Dunn, R. R. What drives elevational patterns of diversity? A test of geometric constraints, climate and species pool effects for pteridophytes on an elevational gradient in Costa Rica. Glob. Ecol. Biogeogr. 15, 358–371 (2006).Article 

Google Scholar 
Sanders, N. J., Lessard, J. P., Fitzpatrick, M. C. & Dunn, R. R. Temperature, but not productivity or geometry, predicts elevational diversity gradients in ants across spatial grains. Glob. Ecol. Biogeogr. 16, 640–649 (2007).Article 

Google Scholar 
Malsch, A. K. F. et al. An analysis of declining ant species richness with increasing elevation at Mount Kinabalu, Sabah, Borneo. Asian Myrmecol. 2, 33–49 (2008).
Google Scholar 
Pérez-Toledo, G. R., Valenzuela-González, J. E., Moreno, C. E., Villalobos, F. & Silva, R. R. Patterns and drivers of leaf-litter ant diversity along a tropical elevational gradient in Mexico. J. Biogeogr. 48, 2515 (2021).Article 

Google Scholar 
Szewczyk, T. M. & McCain, C. M. A systematic review of global drivers of ant elevational diversity. PLoS ONE 11, e155040 (2016).Article 

Google Scholar 
McCain, C. M. & Grytnes, J.-A.A. Elevational gradients in species richness. In Encyclopedia of Life Sciences (ed. Wiley, J.) (Wiley, 2010). https://doi.org/10.1002/9780470015902.a0022548.Chapter 

Google Scholar 
Silva, R. R. & Brandão, C. R. F. Morphological patterns and community organization in leaf-litter ant assemblages. Ecol. Monogr. https://doi.org/10.1890/08-1298.1 (2010).Article 

Google Scholar 
Dunn, R. R. et al. Climatic drivers of hemispheric asymmetry in global patterns of ant species richness. Ecol. Lett. 12, 324–333 (2009).PubMed 
Article 

Google Scholar 
Warren, R. J. & Chick, L. Upward ant distribution shift corresponds with minimum, not maximum, temperature tolerance. Glob. Chang. Biol. 19, 2082–2088 (2013).ADS 
PubMed 
Article 

Google Scholar 
Cerdá, X. & Retana, J. Alternative strategies by thermophilic ants to cope with extreme heat: Individual versus colony level traits. Oikos 89, 155–163 (2000).Article 

Google Scholar 
Kadochová, Š & Frouz, J. Thermoregulation strategies in ants in comparison to other social insects, with a focus on red wood ants (Formica rufa group). F1000 Res. 2, 280 (2013).Article 

Google Scholar 
Moreau, C. S., Bell, C. D., Vila, R., Archibald, S. B. & Pierce, N. E. Phylogeny of the ants: diversification in the age of angiosperms. Science 312, 101–104 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rabeling, C., Brown, J. M. & Verhaagh, M. Newly discovered sister lineage sheds light on early ant evolution. Proc. Natl. Acad. Sci. 105, 14913–14917 (2008).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ward, P. S., Brady, S. G., Fisher, B. L. & Schultz, T. R. The evolution of myrmicine ants: Phylogeny and biogeography of a hyperdiverse ant clade (Hymenoptera: Formicidae). Syst. Entomol. 40, 61–81 (2015).Article 

Google Scholar 
Pie, M. R. The macroevolution of climatic niches and its role in ant diversification. Ecol. Entomol. 41, 301–307 (2016).Article 

Google Scholar 
Smith, M. R. Revision of the genus Stenamma Westwood in America north of Mexico (Hymenoptera, Formicidae). Am. Midl. Nat. 57, 133–174 (1957).Article 

Google Scholar 
Herbers, J. M. & Johnson, C. A. Social structure and winter survival in acorn ants. Oikos 116, 829–835 (2007).Article 

Google Scholar 
Kaspari, M. & Weiser, M. D. Ant activity along moisture gradients in a neotropical forest1. Biotropica 32, 703–711 (2006).Article 

Google Scholar 
Flores, O., Seoane, J., Hevia, V. & Azcárate, F. M. Spatial patterns of species richness and nestedness in ant assemblages along an elevational gradient in a Mediterranean mountain range. PLoS ONE 13, 1–16 (2018).
Google Scholar 
Almeida, R. P. S. et al. Induced drought strongly affects richness and composition of ground-dwelling ants in the eastern Amazon. BioRxiv (2020).Le Breton, J., Chazeau, J. & Jourdan, H. Immediate impacts of invasion by Wasmannia auropunctata (Hymenoptera: Formicidae) on native litter ant fauna in a New Caledonian rainforest. Austral Ecol. 28, 204–209 (2003).Article 

Google Scholar 
Vonshak, M., Dayan, T., Ionescu-Hirsh, A., Freidberg, A. & Hefetz, A. The little fire ant Wasmannia auropunctata: A new invasive species in the Middle East and its impact on the local arthropod fauna. Biol. Invasions 12, 1825–1837 (2010).Article 

Google Scholar 
Wheeler, W. M. Ants: Their Structure, Development and Behavior (Columbia University Press, 1910).
Google Scholar 
Cavender-Bares, J., Ackerly, D. D., Baum, D. A. & Bazzaz, F. A. Phylogenetic overdispersion in Floridian oak communities. Am. Nat. 163, 823–843 (2004).CAS 
PubMed 
Article 

Google Scholar 
Parr, C. L., Sinclair, B. J., Andersen, A. N., Gaston, K. J. & Chown, S. L. Constraint and competition in assemblages: A cross-continental and modeling approach for ants. Am. Nat. 165, 481–494 (2005).PubMed 
Article 

Google Scholar 
Retana, J. & Cerdá, X. Patterns of diversity and composition of Mediterranean ground ant communities tracking spatial and temporal variability in the thermal environment. Oecologia 123, 436–444 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84, 3105–3117 (2003).Article 

Google Scholar 
Graham, C. H., Parra, J. L., Rahbek, C. & McGuire, J. A. Phylogenetic structure in tropical hummingbird communities. Proc. Natl. Acad. Sci. 106, 19673–19678 (2009).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Camacho, G. P., Loss, A. C., Fisher, B. L. & Blaimer, B. B. Spatial phylogenomics of acrobat ants in Madagascar—Mountains function as cradles for recent diversity and endemism. J. Biogeogr. 1, 1706–1719. https://doi.org/10.1111/jbi.14107 (2021).Article 

Google Scholar 
Lobo, J. M. & Halffter, G. Biogeographical and ecological factors affecting the altitudinal variation of mountainous communities of coprophagous beetles (Coleoptera: Scarabaeoidea): A comparative study. Ann. Entomol. Soc. Am. 93, 115–126 (2000).Article 

Google Scholar 
Halffter, G., Favila, M. & Arellano, L. Spatial distribution of three groups of Coleoptera along an altitudinal transect in the Mexican Transition Zone and its biogeographical implications. Elytron 9, 1–10 (1995).
Google Scholar 
Blaimer, B. B. et al. Phylogenomic methods outperform traditional multi-locus approaches in resolving deep evolutionary history: a case study of formicine ants. BMC Evol. Biol. 15, 1–14 (2015).Article 

Google Scholar 
Longino, J. T., Branstetter, M. G. & Colwell, R. K. How ants drop out: ant abundance on tropical mountains. PLoS ONE 9, e104030 (2014).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Longino, J. T. & Branstetter, M. G. The truncated bell: An enigmatic but pervasive elevational diversity pattern in Middle American ants. Ecography 42, 272–283 (2019).Article 

Google Scholar 
Branstetter, M. G. Origin and diversification of the cryptic ant genus Stenamma Westwood (Hymenoptera: Formicidae), inferred from multilocus molecular data, biogeography and natural history. Syst. Entomol. 37, 478–496 (2012).Article 

Google Scholar 
Prebus, M. Insights into the evolution, biogeography and natural history of the acorn ants, genus Temnothorax Mayr (hymenoptera: Formicidae). BMC Evol. Biol. 17, 1–22 (2017).Article 

Google Scholar 
Kluge, J. & Kessler, M. Phylogenetic diversity, trait diversity and niches: Species assembly of ferns along a tropical elevational gradient. J. Biogeogr. 38, 394–405 (2011).Article 

Google Scholar 
Janzen, D. H. Why mountain passes are higher in the tropics. Am. Nat. 101, 233–249 (1967).Article 

Google Scholar 
Fernandes, G. W. et al. Cerrado to rupestrian grasslands: Patterns of species distribution and the forces shaping them along an altitudinal gradient. in Ecology and Conservation of Mountaintop Grasslands in Brazil 345–378 (2016). https://doi.org/10.1007/978-3-319-29808-5_15.Leibold, M. A. et al. The metacommunity concept: A framework for multi-scale community ecology. Ecol. Lett. 7, 601–613 (2004).Article 

Google Scholar 
Perrigo, A., Hoorn, C. & Antonelli, A. Why mountains matter for biodiversity. J. Biogeogr. 47, 315–325 (2020).Article 

Google Scholar 
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Colwell, R. K., Brehm, G., Cardelus, C. L., Gilman, A. C. & Longino, J. T. Global warming, elevational range shifts and lowland biotic attrition in the wet tropics. Science 322, 258–261 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Moreau, C. S. & Bell, C. D. Testing the museum versus cradle tropical biological diversity hypothesis: Phylogeny, diversification, and ancestral biogeographic range evolution of the ants. Evolution 67, 2240–2257 (2013).PubMed 
Article 

Google Scholar 
Borowiec, M. L. Generic revision of the ant subfamily Dorylinae (Hymenoptera, Formicidae). Zookeys 1, 280 (2016).
Google Scholar 
Lapolla, J. S., Brady, S. G. & Shattuck, S. O. Phylogeny and taxonomy of the Prenolepis genus-group of ants (Hymenoptera: Formicidae). Syst. Entomol. 35, 118–131 (2010).Article 

Google Scholar 
Schmidt, C. A. & Shattuck, S. O. The higher classification of the ant subfamily Ponerinae (Hymenoptera: Formicidae), with a review of ponerine ecology and behavior. Zootaxa 3817, 1–242 (2014).CAS 
PubMed 
Article 

Google Scholar 
Revell, L. J. phytools: An R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).Article 

Google Scholar 
Arnan, X., Arcoverde, G. B., Pie, M. R., Ribeiro-Neto, J. D. & Leal, I. R. Increased anthropogenic disturbance and aridity reduce phylogenetic and functional diversity of ant communities in Caatinga dry forest. Sci. Total Environ. 631, 429–438 (2018).ADS 
PubMed 
Article 

Google Scholar 
Divieso, R., Silva, T. S. R. & Pie, M. R. Morphological evolution in the ant reproductive caste. BioRxiv https://doi.org/10.1101/2020.07.18.210302 (2020).Article 

Google Scholar 
Paradis, E. et al. Package ‘ape’. Anal. Phylogenet. Evol. 2, 1–10 (2019).
Google Scholar 
Faith, D. P. Conservation evaluation and phylogenetic diversity. Biol. Conserv. 61, 1–10 (1992).Article 

Google Scholar 
Tucker, C. M. et al. A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biol. Rev. 92, 698–715 (2017).PubMed 
Article 

Google Scholar 
Webb, C. O. Exploring the phylogenetic structure of ecological communities: An example for rain forest trees. Am. Nat. 156, 145–155 (2000).PubMed 
Article 

Google Scholar 
Tucker, C. M. et al. Assessing the utility of conserving evolutionary history. Biol. Rev. 94, 1740–1760 (2019).PubMed 
Article 

Google Scholar 
Kembel, S. W. et al. Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26, 1463–1464 (2010).CAS 
PubMed 
Article 

Google Scholar 
R Core Team. A language and environment for statistical computing. R Found. Stat. Comput. 2, https://www.R-project.org (2021).Baselga, A. & Orme, C. D. L. Betapart: An R package for the study of beta diversity. Methods Ecol. Evol. 3, 808–812 (2012).Article 

Google Scholar 
Dobrovolski, R., Melo, A. S., Cassemiro, F. A. S. & Diniz-Filho, J. A. F. Climatic history and dispersal ability explain the relative importance of turnover and nestedness components of beta diversity. Glob. Ecol. Biogeogr. 21, 191–197 (2012).Article 

Google Scholar 
Peixoto, F. P. et al. Geographical patterns of phylogenetic beta-diversity components in terrestrial mammals. Glob. Ecol. Biogeogr. 26, 573–583 (2017).Article 

Google Scholar 
Körner, C. The use of ‘altitude’ in ecological research. Trends Ecol. Evol. 22, 569–574 (2007).PubMed 
Article 

Google Scholar 
Sundqvist, M. K., Sanders, N. J. & Wardle, D. A. Community and ecosystem responses to elevational gradients: Processes, mechanisms, and insights for global change. Annu. Rev. Ecol. Evol. Syst. 44, 261–280 (2013).Article 

Google Scholar 
Cuervo-Robayo, A. P. et al. An update of high-resolution monthly climate surfaces for Mexico. Int. J. Climatol. 34, 2427–2437 (2014).Article 

Google Scholar 
Hijmans, R. J., Phillips, S., Leathwick, J., Elith, J. & Hijmans, M. R. J. Package ‘dismo’. Circles 9, 1–68 (2017).
Google Scholar 
Schwarz, G. Estimating the dimension of a model. Ann. Stat. 6, 461–464 (1978).MathSciNet 
MATH 
Article 

Google Scholar 
Guthery, F. S., Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach Vol. 67 (Springer, 2003).
Google Scholar 
Mazerolle, M. J. Improving data analysis in herpetology: Using Akaike’s information criterion (AIC) to assess the strength of biological hypotheses. Amphib. Reptil. 27, 169–180 (2006).Article 

Google Scholar 
Ferrier, S., Manion, G., Elith, J. & Richardson, K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Divers. Distrib. 13, 252–264 (2007).Article 

Google Scholar 
Fitzpatrick, M. C. et al. Environmental and historical imprints on beta diversity: Insights from variation in rates of species turnover along gradients. Proc. R. Soc. B Biol. Sci. 280, 20131201 (2013).Article 

Google Scholar 
Manion, G. et al. gdm: Generalized dissimilarity modeling. R Packag. version (2018).Wickham, H. Ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 3, 180–185 (2011).Article 

Google Scholar  More