Philippa Kaur

LookingGlass design and optimizationDataset generationThe taxonomic organization of representative Bacterial and Archaeal genomes was determined from the Genome Taxonomy Database, GTDB51 (release 89.0). The complete genome sequences were downloaded via the NCBI Genbank ftp52. This resulted in 24,706 genomes, comprising 23,458 Bacterial and 1248 Archaeal genomes.Each genome was split into read-length chunks. To determine the distribution of realistic read lengths produced by next-generation short-read sequencing machines, we obtained the BioSample IDs52 for each genome, where they existed, and downloaded their sequencing metadata from the MetaSeek53 database using the MetaSeek API. We excluded samples with average read lengths less than 60 or greater than 300 base pairs. This procedure resulted in 7909 BioSample IDs. The average read lengths for these sequencing samples produced the read-length distribution (Supplementary Fig. 1) with a mean read length of 136 bp. Each genome was split into read-length chunks (with zero overlap in order to maximize information density and reduce data redundancy in the dataset): a sequence length was randomly selected with replacement from the read-length distribution and a sequence fragment of that length was subset from the genome, with a 50% chance that the reverse complement was used. The next sequence fragment was chosen from the genome starting at the end point of the previous read-length chunk, using a new randomly selected read length, and so on. These data were partitioned into a training set used for optimization of the model; a validation set used to evaluate model performance during parameter tuning and as a benchmark to avoid overfitting during training; and a test set used for final evaluation of model performance. To ensure that genomes in the training, validation, and test sets had low sequence similarity, the sets were split along taxonomic branches such that genomes from the Actinomycetales, Rhodobacterales, Thermoplasmata, and Bathyarchaeia were partitioned into the validation set; genomes from the Bacteroidales, Rhizobiales, Methanosarcinales, and Nitrososphaerales were partitioned into the test set; and the remaining genomes remained in the training set. This resulted in 529,578,444 sequences in the training set, 57,977,217 sequences in the validation set, and 66,185,518 sequences in the test set. We term this set of reads the GTDB representative set (Table 1).Table 1 Summary table of datasets used.Full size tableThe amount of data needed for training was also evaluated (Supplementary Fig. 2). Progressively larger amounts of data were tested by selecting at random 1, 10, 100, or 500 read-length chunks from each of the GTDB representative genomes in the GTDB representative training set. Additionally, the performance of smaller but more carefully selected datasets, representing the diversity of the microbial tree of life, were tested by selecting for training one genome at random from each taxonomic class or order in the GTDB taxonomy tree. In general, better accuracy was achieved in fewer epochs with a greater amount of sequencing data (Supplementary Fig. 2); however, a much smaller amount of data performed better if a representative genome was selected from each GTDB taxonomy class.The final LookingGlass model was trained on this class-level partition of the microbial tree of life. We term this dataset the GTDB class set (Table 1). The training, validation, and test sets were split such that no classes overlapped across sets: the validation set included 8 genomes from each of the classes Actinobacteria, Alphaproteobacteria, Thermoplasmata, and Bathyarchaeia (32 total genomes); the test set included 8 genomes from each of the classes Bacteroidia, Clostridia, Methanosarcinia, and Nitrososphaeria (32 total genomes); and the training set included 1 genome from each of the remaining classes (32 archaeal genomes and 298 bacterial genomes for a total of 330 genomes). This resulted in a total of 6,641,723 read-length sequences in the training set, 949,511 in the validation set, and 632,388 in the test set (Supplementary Data 1).Architecture design and trainingRecurrent neural networks (RNNs) are a type of neural network designed to take advantage of the context dependence of sequential data (such as text, video, audio, or biological sequences), by passing information from previous items in a sequence to the current item in a sequence54. Long short-term memory networks (LSTMs)55 are an extension of RNNs, which better learn long-term dependencies by handling the RNN tendency to “forget” information farther away in a sequence56. LSTMs maintain a cell state which contains the “memory” of the information in the previous items in the sequence. LSTMs learn additional parameters which decide at each step in the sequence which information in the cell state to “forget” or “update”.LookingGlass uses a three-layer LSTM encoder model with 1152 units in each hidden layer and an embedding size of 104 based on the results of hyperparameter tuning (see below). It divides the sequence into characters using a kmer size of 1 and a stride of 1, i.e., is a character-level language model. LookingGlass is trained in a self-supervised manner to predict a masked nucleotide, given the context of the preceding nucleotides in the sequence. For each read in the training sequence, multiple training inputs are considered, shifting the nucleotide that is masked along the length of the sequence from the second position to the final position in the sequence. Because it is a character-level model, a linear decoder predicts the next nucleotide in the sequence from the possible vocabulary items “A”, “C”, “G”, and “T”, with special tokens for “beginning of read”, “unknown nucleotide” (for the case of ambiguous sequences), “end of read” (only “beginning of read” was tokenized during LookingGlass training), and a “padding” token (used for classification only).Regularization and optimization of LSTMs require special approaches to dropout and gradient descent for best performance57. The fastai library58 offers default implementations of these approaches for natural language text, and so we adopt the fastai library for all training presented in this paper. We provide the open source fastBio python package59 which extends the fastai library for use with biological sequences.LookingGlass was trained on a Pascal P100 GPU with 16GB memory on Microsoft Azure, using a batch size of 512, a back propagation through time (bptt) window of 100 base pairs, the Adam optimizer60, and utilizing a Cross Entropy loss function (Supplementary Table 1). Dropout was applied at variable rates across the model (Supplementary Table 1). LookingGlass was trained for a total of 12 days for 75 epochs, with progressively decreasing learning rates based on the results of hyperparameter optimization (see below): for 15 epochs at a learning rate of 1e−2, for 15 epochs at a learning rate of 2e−3, and for 45 epochs at a learning rate of 1e−3.Hyperparameter optimizationHyperparameters used for the final training of LookingGlass were tuned using a randomized search of hyperparameter settings. The tuned hyperparameters included kmer size, stride, number of LSTM layers, number of hidden nodes per layer, dropout rate, weight decay, momentum, embedding size, bptt size, learning rate, and batch size. An abbreviated dataset consisting of ten randomly selected read-length chunks from the GTDB representative set was created for testing many parameter settings rapidly. A language model was trained for two epochs for each randomly selected hyperparameter combination, and those conditions with the maximum performance were accepted. The hyperparameter combinations tested and the selected settings are described in the associated Github repository61.LookingGlass validation and analysis of embeddingsFunctional relevanceDataset generation
In order to assess the ability of the LookingGlass embeddings to inform the molecular function of sequences, metagenomic sequences from a diverse set of environments were downloaded from the Sequence Read Archive (SRA)62. We used MetaSeek53 to choose ten metagenomes at random from each of the environmental packages defined by the MIxS metadata standards63: built environment, host-associated, human gut, microbial mat/biofilm, miscellaneous, plant-associated, sediment, soil, wastewater/sludge, and water, for a total of 100 metagenomes. The SRA IDs used are available in (Supplementary Table 2). The raw DNA reads for these 100 metagenomes were downloaded from the SRA with the NCBI e-utilities. These 100 metagenomes were annotated with the mi-faser tool27 with the read-map option to generate predicted functional annotation labels (to the fourth digit of the Enzyme Commission (EC) number), out of 1247 possible EC labels, for each annotatable read in each metagenome. These reads were then split 80%/20% into training/validation candidate sets of reads. To ensure that there was minimal overlap in sequence similarity between the training and validation set, we compared the validation candidate sets of each EC annotation to the training set for that EC number with CD-HIT64, and filtered out any reads with >80% DNA sequence similarity to the reads of that EC number in the training set (the minimum CD-HIT DNA sequence similarity cutoff). In order to balance EC classes in the training set, overrepresented ECs in the training set were downsampled to the mean count of read annotations (52,353 reads) before filtering with CD-HIT. After CD-HIT processing, any underrepresented EC numbers in the training set were oversampled to the mean count of read annotations (52,353 reads). The validation set was left unbalanced to retain a distribution more realistic to environmental settings. The final training set contained 61,378,672 reads, while the validation set contained 2,706,869 reads. We term this set of reads and their annotations the mi-faser functional set (Table 1).
As an external test set, we used a smaller number of DNA sequences from genes with experimentally validated molecular functions. We linked the manually curated entries of Bacterial or Archaeal proteins from the Swiss-Prot database65 corresponding to the 1247 EC labels in the mi-faser functional set with their corresponding genes in the EMBL database66. We downloaded the DNA sequences, and selected ten read-length chunks at random per CDS. This resulted in 1,414,342 read-length sequences in the test set. We term this set of reads and their annotations the Swiss-Prot functional set (Table 1).

Fine-tuning procedure
We fine-tuned the LookingGlass language model to predict the functional annotation of DNA reads, to demonstrate the speed with which an accurate model can be trained using our pretrained LookingGlass language model. The architecture of the model retained the 3-layer LSTM encoder and the weights of the LookingGlass language model encoder, but replaced the language model decoder with a new multiclass classification layer with pooling (with randomly initialized weights). This pooling classification layer is a sequential model consisting of the following layers: a layer concatenating the output of the LookingGlass encoder with min, max, and average pooling of the outputs (for a total dimension of 104*3 = 312), a batch normalization67 layer with dropout, a linear layer taking the 312-dimensional output of the batch norm layer and producing a 50-dimensional output, another batch normalization layer with dropout, and finally a linear classification layer that is passed through the log(Softmax(x)) function to output the predicted functional annotation of a read as a probability distribution of the 1247 possible mi-faser EC annotation labels. We then trained the functional classifier on the mi-faser functional set described above. Because the >61 million reads in the training set were too many to fit into memory, training was done in 13 chunks of ~5-million reads each until one total epoch was completed. Hyperparameter settings for the functional classifier training are seen in Supplementary Table 1.

Encoder embeddings and MANOVA test
To test whether the LookingGlass language model embeddings (before fine-tuning, above) are distinct across functional annotations, a random subset of ten reads per functional annotation was selected from each of the 100 SRA metagenomes (or the maximum number of reads present in that metagenome for that annotation, whichever was greater). This also ensured that reads were evenly distributed across environments. The corresponding fixed-length embedding vectors for each read was produced by saving the output from the LookingGlass encoder (before the embedding vector is passed to the language model decoder) for the final nucleotide in the sequence. This vector represents a contextually relevant embedding for the overall sequence. The statistical significance of the difference between embedding vectors across all functional annotation groups was tested with a MANOVA test using the R stats package68.
Evolutionary relevance
Dataset generation
The OrthoDB database69 provides orthologous groups (OGs) of proteins at various levels of taxonomic distance. For instance, the OrthoDB group “77at2284” corresponds to proteins belonging to “Glucan 1,3-alpha-glucosidase at the Sulfolobus level”, where “2284” is the NCBI taxonomy ID for the genus Sulfolobus.
We tested whether embedding similarity of homologous sequences (sequences within the same OG) is higher than that of nonhomologous sequences (sequences from different OGs). We tested this in OGs at multiple levels of taxonomic distance—genus, family, order, class, and phylum. At each taxonomic level, ten individual taxa at that level were chosen from across the prokaryotic tree of life (e.g., for the genus level, Acinetobacter, Enterococcus, Methanosarcina, Pseudomonas, Sulfolobus, Bacillus, Lactobacillus, Mycobacterium, Streptomyces, and Thermococcus were chosen). For each taxon, 1000 randomly selected OGs corresponding to that taxon were chosen; for each of these OGs, five randomly chosen genes within this OG were chosen.
OrthoDB cross-references OGs to UniProt65 IDs of the corresponding proteins. We mapped these to the corresponding EMBL CDS IDs66 via the UniProt database API65; DNA sequences of these EMBL CDSs were downloaded via the EMBL database API. For each of these sequences, we generated LookingGlass embedding vectors.

Homologous and nonhomologous sequence pairs
To create a balanced dataset of homologous and nonhomologous sequence pairs, we compared all homologous pairs of the five sequences in an OG (total of ten homologous pairs) to an equal number of randomly selected out-of-OG comparisons for the same sequences; i.e., each of the five OG sequences was compared to 2 other randomly selected sequences from any other randomly selected OG (total of ten nonhomologous pairs). We term this set of sequences, and their corresponding LookingGlass embeddings, the OG homolog set (Table 1).

Embedding and sequence similarity
For each sequence pair, the sequence and embedding similarity were determined. The embedding similarity was calculated as the cosine similarity between embedding vectors. The sequence similarity was calculated as the Smith-Waterman alignment score using the BioPython70 pairwise2 package, with a gap open penalty of −10 and a gap extension penalty of −1. The IDs of chosen OGs, the cosine similarities of the embedding vectors, and sequence similarities of the DNA sequences are available in the associated Github repository61.

Comparison to HMM-based domain searches for distant homology detection
Distantly related homologous sequences that share, e.g., Pfam71, domains can be identified using HMM-based search methods. We used hmmscan25 (e-val threshold = 1e−10) to compare homologous (at the phylum level) sequences in the OG homolog set, for which the alignment score was less than 50 bits and the embedding similarity was greater than 0.62 (total: 21,376 gene pairs). Specifically, we identified Pfam domains in each sequence and compared whether the most significant (lowest e-value) domain for each sequence was identified in common for each homologous pair.
Environmental relevance
Encoder embeddings and MANOVA test
The LookingGlass embeddings and the environment of origin for each read in the mi-faser functional set were used to test the significance of the difference between the embedding vectors across environmental contexts. The statistical significance of this difference was evaluated with a MANOVA test using the R stats package68.
Oxidoreductase classifier
Dataset generation
The manually curated, reviewed entries of the Swiss-Prot database65 were downloaded (June 2, 2020). Of these, 23,653 entries were oxidoreductases (EC number 1.-.-.-) of Archaeal or Bacterial origin (988 unique ECs). We mapped their UniProt IDs to both their EMBL CDS IDs and their UniRef50 IDs via the UniProt database mapper API. Uniref50 IDs identify clusters of sequences with >50% amino acid identity. This cross-reference identified 28,149 EMBL CDS IDs corresponding to prokaryotic oxidoreductases, belonging to 5451 unique UniRef50 clusters. We split this data into training, validation, and test sets such that each UniRef50 cluster was contained in only one of the sets, i.e., there was no overlap in EMBL CDS IDs corresponding to the same UniRef50 cluster across sets. This ensures that the oxidoreductase sequences in the validation and test sets are dissimilar to those seen during training. The DNA sequences for each EMBL CDS ID were downloaded via the EMBL database API. These data generation process were repeated for a random selection of non-oxidoreductase UniRef50 clusters, which resulted in 28,149 non-oxidoreductase EMBL CDS IDs from 13,248 unique UniRef50 clusters.
Approximately 50 nucleotide read-length chunks (selected from the representative read-length distribution, as above) were selected from each EMBL CDS DNA sequence, with randomly selected start positions on the gene and a 50% chance of selecting the reverse complement, such that an even number of read-length sequences with “oxidoreductase” and “not oxidoreductase” labels were generated for the final dataset. This procedure produced a balanced dataset with 2,372,200 read-length sequences in the training set, 279,200 sequences in the validation set, and 141,801 sequences in the test set. We term this set of reads and their annotations the oxidoreductase model set (Table 1). In order to compare the oxidoreductase classifier performance to an HMM-based method, reads with “oxidoreductase” labels in the oxidoreductase model test set (71,451 reads) were 6-frame translated and searched against the Swiss-Prot protein database using phmmer25 (reporting e-val threshold = 0.05, using all other defaults).

Fine-tuning procedure
Since our functional annotation classifier addresses a closer classification task to the oxidoreductase classifier than LookingGlass itself, the architecture of the oxidoreductase classifier was fine-tuned starting from the functional annotation classifier, replacing the decoder with a new pooling classification layer (as described above for the functional annotation classifier) and with a final output size of 2 to predict “oxidoreductase” or “not oxidoreductase”. Fine tuning of the oxidoreductase classifier layers was done successively, training later layers in isolation and then progressively including earlier layers into training, using discriminative learning rates ranging from 1e−2 to 5e−4, as previously described72. The fine-tuned model was trained for 30 epochs, over 18 h, on a single P100 GPU node with 16GB memory.

Model performance in metagenomes
Sixteen marine metagenomes from the surface (SRF, ~5 meters) and mesopelagic (MES, 175–800 meters) from eight stations sampled as part of the TARA expedition37 were downloaded from the SRA62 (Supplementary Table 3, SRA accession numbers ERR598981, ERR599063, ERR599115, ERR599052, ERR599020, ERR599039, ERR599076, ERR598989, ERR599048, ERR599105, ERR598964, ERR598963, ERR599125, ERR599176, ERR3589593, and ERR3589586). Metagenomes were chosen from a latitudinal gradient spanning polar, temperate, and tropical regions and ranging from −62 to 76 degrees latitude. Mesopelagic depths from four out of the eight stations were sampled from oxygen minimum zones (OMZs, where oxygen More

Listen to the latest from the world of science, with Shamini Bundell and Benjamin Thompson.

Your browser does not support the audio element.

Download MP3

In this episode:00:57 Reviving retinas to understand eyesResearch efforts to learn more about diseases of the human eye have been hampered as these organs degrade rapidly after death, and animal eyes are quite different to those from humans. To address this, a team have developed a new method to revive retinas taken from donors shortly after their death. They hope this will provide tissue for new studies looking into the workings of the human eye and nervous system.Research article: Abbas et al.08:05 Research HighlightsA technique that simplifies chocolate making yields fragrant flavours, and 3D imaging reveals some of the largest-known Native American cave art.Research Highlight: How to make a fruitier, more floral chocolateResearch Highlight: Cramped chamber hides some of North America’s biggest cave art10:54 Did life emerge in an ‘RNA world’?How did the earliest biochemical process evolve from Earth’s primordial soup? One popular theory is that life began in an ‘RNA world’ from which proteins and DNA evolved. However, this week a new paper suggests that a world composed of RNA alone is unlikely, and that life is more likely to have begun with molecules that were part RNA and part protein.Research article: Müller et al.News and Views: A possible path towards encoded protein synthesis on ancient Earth17:52 Briefing ChatWe discuss some highlights from the Nature Briefing. This time, the ‘polarised sunglasses’ that helped astronomers identify an ultra-bright pulsar, and how a chemical in sunscreen becomes toxic to coral.Nature: A ‘galaxy’ is unmasked as a pulsar — the brightest outside the Milky WayNature: A common sunscreen ingredient turns toxic in the sea — anemones suggest whySubscribe to Nature Briefing, an unmissable daily round-up of science news, opinion and analysis free in your inbox every weekday.Never miss an episode: Subscribe to the Nature Podcast on Apple Podcasts, Google Podcasts, Spotify or your favourite podcast app. Head here for the Nature Podcast RSS feed. More

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

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