Ecology

Two databases including yield, management, and weather data for maize (n = 17,013) and soybean (n = 24,848) involving US crop performance trials conducted in 28 states between 2016 to 2018 for maize and between 2014 to 2018 for soybean, were developed (Fig. 1). Crop yield and management data were obtained from publicly available variety performance trials which are typically performed yearly in several locations across each state (see methods for more information). Final databases were separated in training (80% of database) and testing (20% of database) datasets using stratified sampling by year, use of irrigation, and soil type. For each crop, an extreme gradient boosting (XGBoost, see methods for more information) algorithm to estimate yield based on soil type and weather conditions (E), seed traits (G) and management practices (M) was developed (see variables listed in Tables S1 and S2 for maize and soybean, respectively, and data science workflow in Fig. S1).Figure 1Locations where maize and soybean trials were performed during the examined period. The map was developed in ArcGIS Pro 2.8.0 (https://www.esri.com).Full size imageThe developed algorithms exhibited a high degree of accuracy when estimating yield in independent datasets (test dataset not used for model calibration) (Fig. 2). For maize, the root mean square error (RMSE) and mean absolute error (MAE) was a respective 4.7 and 3.6% of the dataset average yield (13,340 kg/ha). For soybean, the respective RMSE and MAE was 6.4 and 4.9% of the dataset average yield (4153 kg/ha). As is evident in the graphs (Fig. 2), estimated yields exhibited a high degree of correlation with actual yields for both algorithms in the independent datasets. For maize and soybean, 72.3 and 60% of cases in the test dataset deviated less than 5% from actual yields, respectively. Maximum deviation for maize and soybean reached 43 and 70%, respectively. Data points with deviations greater than 15% from actual yield were 1.5% in maize and 3.6% in soybean databases. These results suggest that the developed algorithms can accurately estimate maize and soybean yields utilizing database-generated information involving reported environmental, seed genetic, and crop management variables.Figure 2Actual versus algorithm-derived maize (left) and soybean (right) yield in test datasets. Black solid line indicates y = x, red short-dashed lines, black dashed lines, and red long-dashed lines indicate ± 5, 10, and 15% deviation from the y = x line. RMSE, root mean square error; MAE, mean absolute error; r2, coefficient of determination; n = number of observations. Each observation corresponds to a yield of an individual cropping system in a specific environment (location-year).Full size imageIn contrast to statistical models, ML algorithms can be complex, and the effect of single independent variables may not obvious. However, accumulated local effects (ALE) plots14 can aid the understanding and visualization of important and possibly correlated features in ML algorithms. For both crops, indicatively important variables included sowing date, seeding rate, nitrogen fertilizer (for maize), row spacing (for soybean) and June to September cumulative precipitation (Fig. 3). Across the entire region and for both crops, the algorithm-derived trends suggest that above average yields occur in late April to early May sowing dates, but sharply decrease thereafter. Similar responses have been observed in many regional studies across the US for both, maize15,16,17,18 and soybean19. Similarly, simulated yield curves due to increasing seeding rate are in close agreement with previous maize20,21 and soybean22 studies. The maize algorithm has captured the increasing yield due to increasing N fertilizer rate. The soybean algorithm suggests that narrower row spacing resulted in above average yield compared to wider spacing. Such response has been observed in many regions across the US23. Season cumulative precipitation between 400 and 700 mm resulted in above average yields for both crops.Figure 3Accumulated local effect plots for maize sowing date (A), seeding rate (B), Nitrogen fertilizer rate (C), and cumulative precipitation between June and September (mm) (D), and soybean sowing date (E), seeding rate (F), row spacing (G), and cumulative precipitation between June and September (mm) (H).Full size imageThe responses in the ALE plots (Fig. 3) suggest that these algorithms have captured the general expected average responses for important single features. Nevertheless, our databases include hundreds of locations with diverse environments across the US and site-specific crop responses which may vary due to components of the G × E × M interaction. We argue that, instead of examining a single or low-order management interactions, site-specific evaluation of complex high order interactions (a.k.a. cropping systems) can reveal yield differences that current research approaches cannot fully explore and quantify. For example, sowing date exerts a well-known impact on maize and soybean yield. For each crop separately, by creating a hypothetical cropping system (a single combination of all management and traits in Tables S1 and S2) in a randomly chosen field in south central Wisconsin (latitude = 43.34, longitude = -89.38), and by applying the developed algorithms, we can generate estimates of maize and soybean yield. For that specific field and cropping system (out of the vast number of management combinations a farmer can choose from), maize yield with May 1st sowing was 711 kg/ha greater (6% increase) than June sowing (Fig. 4A). By creating scenarios with 256 background cropping system choices (Table S3), the resultant algorithm-derived yield estimate difference for the same sowing date contrast (averaged across varying cropping systems) was smaller but still positive (3% increase), although the range of possible yield differences was wider (Fig. 4B). However, when comparing, instead of averaging, the estimated yield potential among the simulated cropping systems, a 2903 kg/ha yield difference (25% difference) was observed (Fig. 4C). Interestingly, when focusing on the early sown fields that were expected to exhibit the greatest yield, the same yield difference was observed (Fig. 4D). This result shows that sub-optimal background management can mitigate the beneficial effect of early sowing (Table S4).Figure 4Maize yield difference (in kg/ha and percentage) due to sowing date (May 1st vs. June 1st) for a single identical background cropping system (A), maize yield difference due to sowing date when averaged across 256 (3 years × 256 cropping systems = 768 year-specific yields) (B), maize yield variability in each of the 256 cropping systems (C), and maize yield variability in each of the 128 cropping systems with early sowing (D). Soybean yield difference due to sowing date (May 1st vs June 1st) for a single identical background cropping system (E), soybean yield difference due to sowing date when averaged across 128 (5 years × 128 cropping systems = 640 year-specific yields) (F), soybean yield in each of the 128 cropping systems (G) and soybean yield variability due in each of the 64 cropping systems with early sowing (H). Within each panel, the horizontal red and grey lines indicate the boxplot with maximum and minimum yield, respectively. In the left four panels, boxes delimit first and third quartiles; solid lines inside boxes indicate median and green triangles indicate means. Upper and lower whiskers extend to maximum and minimum yields. Each maize and soybean cropping system is a respective 8-way and a 7-way interaction of management practices in a randomly chosen field in Wisconsin, USA (Table S3 and S5, respectively).Full size imageIn the case of soybean, a May 1st sowing resulted in greater yield (588 kg/ha; a 14% increase) than a June 1st in the single background cropping system (Fig. 4E). The result was consistent when yield differences due to sowing date were averaged across 128 background cropping system choices (Table S5) (Fig. 4F). Similar to what was observed in maize, among all cropping systems, yield varied by 1704 kg/ha (44% difference) (Fig. 4G). When focusing only on the early sown fields, a 1181 kg/ha yield difference (27% yield increase) was observed (Fig. 4H). In agreement with maize, this result highlights the importance of accounting for sub-optimal background management which can mitigate the beneficial effect of early sowing (Table S6).We note here the ability of farmers to change management practices can be limited due to an equipment constraint (e.g., change planter unit row width) or simply impossible (e.g., change the previous year’s crop). Thus, recommended management practices that were evaluated in studies that used specific background management may not be applicable in some instances. The benefits of the foregoing approach, which involves extensive up-to-date agronomic datasets and high-level computational programing, can have important and immediate implications in future agricultural trials. Our approach allows for more precise examination of complex management interactions in specific environments (soil type and growing season weather) across the US (region covered in Fig. 1). The ability to extract single management practice information (even across cropping systems) is also possible by utilizing ALE plots, or by calculation of the frequency at which a given level/rate of a management practice appeared among the highest yielding cropping systems (Tables S4 and S6).Among all available 30-d weather variables, many were strongly correlated in both crop databases (Figs. S2 and S3 for maize and soybean, respectively). Models using all 30-d interval variables with r  More

Sampling and genotypingWe collected hair and cheek swab samples from 77 men from geographically and linguistically distinct groups of the Caucasus: Kartvelian speakers from Georgia and Turkey, Northeast Caucasian speakers and Turkic speakers from the Russian Federation and Armenian speakers from Georgia’s southern province of Javakheti, descendants of the families displaced from Mush and Erzurum provinces of eastern Turkey in the early nineteenth century (Table 1, Fig. 1). To maximize the representativeness of the genetic signature of each population, the samples were collected from locals with no ancestors from outside of the respective ethnic/geographic population over the last three generations. DNA was extracted from follicles of 10–12 male chest hairs and cheek swab samples. Extraction was performed using Qiagen DNeasy Blood and Tissue kit, following the manufacturer’s recommendations (Qiagen, Valencia, CA, USA). The DNA samples were genotyped for 693,719 autosomal and 17,678 X-chromosomal SNPs by Family Tree DNA (FTDNA—Gene By Gene, Ltd, Houston, TX, www.familytreedna.com).Table 1 Modern study populations of the Caucasus. Latitude and longitude georeference population hubs.Full size tableFigure 1The distribution of the study populations: averaged centroids of ancient populations (uniquely colored points in the main map, see Table 2 for details) and hubs of the modern Caucasian populations (identified in the inset map, see Table 1 for details). Glacial human refugia extracted from Gavashelishvili and Tarkhnishvili5 are shaded in purple. The map is generated using QGIS Desktop 3.10.6-A Coruña (https://qgis.org).Full size imageOur dataset of modern Caucasian genotypes was supplemented with published 10 modern Mbuti (Supplementary Table S1) and 122 Upper Paleolithic-Mesolithic human genotypes, retrieved as a part of 1240 K dataset from David Reich’s Lab website, Harvard University (https://reich.hms.harvard.edu/downloadable-genotypes-present-day-and-ancient-dna-data-compiled-published-papers; see Supplementary Table S2 for details). The ancient genotypes were selected such that they either dated from the LGM or fell within the glacial refugia identified by Gavashelishvili and Tarkhnishvili5. We did so in order to maximize the genetic signature of potential refugial populations in our analysis. We divided the ancient genotypes into 2000-year-long intervals, and then grouped each of these intervals into geographic units (hereafter ancient populations, Table 2, Fig. 1). The modern and ancient genotypes were merged using PLINK 1.9 (PLINK 1.9: www.cog-genomics.org/plink/1.9/27.Table 2 Ancient study populations. The ancient genotypes are divided into 2000-year-long intervals, and then each of these intervals is grouped into geographic units (i.e. ancient populations). Age, latitude and longitude are averaged across each ancient population (see Supplementary Table S2 for details).Full size tableEthics statementThe research team members, through their contacts in the studied communities, inquired whether locals would voluntarily participate in genetic research that would help clarify the genetic makeup of the Caucasus. A verbal agreement was made with volunteer donors of DNA samples, according to which the results would be communicated, electronically or in hard copy, with participants individually. Participants were informed that, upon the completion of the lab work, the research would be published without mentioning the names of sample donors. Those who agreed provided us with the envelopes containing their chest hairs or cheek swab samples, with the birthplace of their ancestors (last three generations) written on the envelope or a piece of paper. In accordance with the preferences of the sample donors, the agreement was verbal and not written. The envelopes and papers are stored as evidence of voluntary provision of the samples and the related information. Analysis of data was done anonymously, using only location and ethnic information; only the first and third authors of the manuscript had access to names associated with the samples. Therefore, this study was based on noninvasive and nonintrusive sampling (volunteers provided hair and swab samples they collected themselves), and the information destined for open publication does not contain any personal information. The study methodology and the procedure of verbal consent was discussed in detail with and approved by the members of the Ilia State University Commission for Ethical Issues before the field survey started, and the commission decided that formal ethical approval was not needed for conducting this study. This is confirmed in a letter from the commission chairman, a copy of which has been provided to the journal editor as part of the submission process.Genetic affinity and geographyFirst, we measured genetic affinity between the modern Caucasian populations, and between the modern populations and the ancient populations of hunter-gatherers, and then tested whether the genetic affinity between these populations was determined by geographic features. Data were mapped using QGIS Desktop 3.10.6-A Coruña, whereas graphs were created using the “ggplot2” package28 in R version 3.5.229.To evaluate genetic affinities and structure of the modern populations, we used Wright’s fixation index (Fst), inbreeding coefficient, admixture analysis and the principal component analysis (PCA). For these procedures we filtered the raw SNP genotypes in PLINK 1.9, first removing all SNPs with the minor allele frequency  0.3, calculated in windows of 50 bp size and 10 bp steps (–maf 0.05 –indep-pairwise 50 10 0.3). Since all individuals in our dataset possess a single copy of the X-chromosome, we did not expect any differential ploidy bias, and SNPs on the X were treated similarly to those on the autosomes. Fst pairwise values were calculated using the smartpca program of EIGENSOFT30 with default parameters, inbreed: YES, and fstonly: YES. The relationship between the modern populations based on Fst values was visualized by constructing a neighbor-joining tree using the “ape” package31 in R version 3.5.2. The average and standard deviation of the inbreeding coefficient for each population was calculated using “fhat2” estimate of PLINK 1.9. The LD pruned genotypes were used in ADMIXTURE 1.3.032, performed in unsupervised mode in order to infer the population structure from the modern individuals. The number of clusters (k) was varied from 2 to 7 and the fivefold cross-validation error was calculated for each k33. We conducted principal components analysis in the smartpca program of EIGENSOFT30, using default parameters and the lsqproject: YES and numoutlieriter: 0 options. Eigenvectors of principal components were inferred with the modern populations from the Caucasus, while the ancient populations were then projected onto the PCA plots. We also assessed the relatedness between sampled individuals using kinship coefficients estimated by KING34.To quantify genetic affinities between the modern and ancient populations, we used the programs qp3Pop and qpDstat in the ADMIXTOOLS suite (https://github.com/DReichLab35 for f3- and f4-statistics, respectively. f3-statistics of the form f3(X,Y,Outgroup) measure the amount of shared genetic drift of populations X and Y after their divergence from an outgroup. We used an ancient population and a modern Caucasian population for X, Y and Mbuti as an outgroup. f4-statistics of the form f4(Outgroup,Test;X,Y) show if population Test is equally related to X and Y or shares an excess of alleles with either of the two. In the f4-statistic calculation we used Mbuti for Outgroup, a modern population of the Caucasus for Test, and X and Y for contemporaneous ancient populations. This meant that f4  0 indicated higher genetic affinity between the test population and Y.To quantify geographic features, we derived least-cost paths and measured least-cost distances (LCD) between the modern and ancient populations using the Least Cost Path Plugin for QGIS. The computation of LCD considers a friction grid that is a raster map where each cell indicates the relative difficulty (or cost) of moving through that cell. A least-cost path minimizes the sum of frictions of all cells along the path, and this sum is the least-cost distance (LCD). For impedance to human movement and expansion, we used 15 geographic features (Table 3). All gridded geographic features (i.e. raster layers) were resampled to a resolution of 1 km using the nearest-neighbor assignment technique. All possible subsets of the 15 geographic features, that did not cancel out each other, were used to calculate different variables of LCD. We assumed that most human movements occurred during climate warming events when the earth’s surface was not dramatically different from that of today, and hence used the current data of the geographic features.Table 3 Geographic features used in combinations to calculate least-cost distances (LCD) between ancient populations and modern Caucasians.Full size tableLinking genetic affinity and geographyGeneralized additive models (GAMs) were used to fit the outgroup f3-statistic to time and variously calculated LCD between the modern and ancient populations using the “mgcv” package36 in R version 3.5.2. Time between the modern and ancient populations was measured in BP (years before present, defined by convention as years before 1950 CE). We used GAMs because without any assumptions they are able to find nonlinear and non-monotonic relationships. GAMs were fitted using a Gamma family with a log link function. Penalized thin plate regression splines were used to represent all the smooth terms. The restricted maximum likelihood (REML) estimation method was implemented to estimate the smoothing parameter because it is the most robust of the available GAM methods36.Model and variable selection were performed by exploring LCD, time BP and the interaction term. The predictive power of the models was evaluated through a tenfold cross-validation. The cross-validation of many models was handled through R’s parallelization capabilities37,38. The best model was selected by the mean squared error of the cross-validation. Akaike’s Information Criterion (AIC) is generally used as a means for model selection. However, we preferred cross-validation for model selection because AIC a priori assumes that simpler models with the high goodness of fit are more likely to have the higher predictive power, while cross-validation without any a priori assumptions measures the predictive performance of a model by efficiently running model training and testing on the available data.We additionally validated the effect of different subsets of geographic features by assessing the relationship between statistically significant values of f4-statistic (i.e. |Z| > 3) and each subset. The relationship between f4-statistic of the form of f4(Outgroup,Test;X,Y) and geographic features was determined by measuring the agreement between the negative/positive signs of f4-statistic and the difference in LCD (LCD.D) for each pair of contemporaneous ancient populations X and Y. LCD.D was calculated as (LCD1–LCD2), where LCD1 was least-cost distance between the test population and X, and LCD2 was least-cost distance between the test population and Y. LCD.D  0 indicated less least-cost distance between Test and Y. So, the same sign of f4 and LCD.D values indicated agreement between geographic proximity and genetic affinity. We used Cohen’s kappa39 to measure the agreement.In order to test if geographic features (Table 3) accounted for present-day genetic differentiation in the Caucasus, we measured the relationship between Fst and LCD across the modern populations using the Mantel test in the “vegan” package40 in R version 3.5.2. In addition, we checked whether contribution from ancient samples was related to today’s genetic differentiation. To do so, we calculated median of f3-statistic of ancient populations of each geographic grouping (e.g. the following 6 populations made up one group: Balkans 39,950–41,950 BP, Balkans 37,950–39,950 BP, Balkans 31,950–33,950 BP, Balkans 9950–11,950 BP, Balkans 7950–9950 BP, Balkans 5950–7950 BP). Then we measured the manhattan distance of f3 median values of all combinations of the geographic groupings between the modern populations and compared the results to Fst and LCD using the Mantel test. More

Regular structures and isothermal theoremFor networks where each node represents a single individual, the isothermal theorem of evolutionary graph theory shows that the fixation probability is the same as the fixation probability of a well-mixed population if the temperature distribution is homogeneous across the whole population1. The temperature of a node defined as the sum over all the weights leads to that node. This theorem extends to structured meta-populations for any migration probability (lambda ): If the underlying structure of the meta-population that connects the patches is a regular network and the local population size is identical in each patch, the temperature of all individuals is identical, regardless of the value of the migration probability. Therefore, the fixation probability in a population with such a structure is the same as the fixation probability in a well-mixed population of the same total population size (N=sum _{j=1}^M N_j), given by ( phi _{mathrm{wm}}^N(r)).Small migration regimeIf the migration probability is small enough such that the time between two subsequent migration events (( sim frac{1}{lambda } )) is much longer than the absorption time within any patch, then at the time of each migration event we may suppose that the meta-population is in a homogeneous configuration22,28. In other words, the low migration regime is an approximation in which we neglect the probability that the meta-population is not in a homogeneous configuration at the time of migration events. We define a homogeneous configuration of the meta-population as a configuration in which in all patches either all individuals are mutants, or all are wild-types.Therefore, instead of having (2^N) states, where N is the population size, the system has only (2^M) states, where M is the number of patches. Thus, we can calculate the fixation probability exactly as in the case of a standard evolutionary graph model where each node represents a single individual but with a modified transition probabilities.In a network with homogeneous patches, in order to increase the number of homogeneous mutant-patches one individual mutant needs to migrate to one of its neighbouring homogeneous wild-type-patches and reaches fixation there. For example if node j is occupied by mutants and one of its neighbouring patches, node k, is occupied by wild-types, the probability that one mutant individual from patch j migrates to patch k and reaches fixation there is (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(r) ), where (mathrm{deg} (j) ) is the degree of node j to take into account that the mutant can move to different patches. This is analogous to the probability that one mutant in node j replaces one wild-type in node k ,(T^{jrightarrow k}), in the network of individuals.Similarly, if node j is occupied by wild-types and one of its neighbouring patches, node j, is occupied by mutants the probability that one wild-type individual from patch j migrates to patch k and reaches fixation there equals to (frac{lambda }{mathrm{deg} (j)}phi _{mathrm{wm}}^{N_{k}}(1/r) ) where (mathrm{deg} (j) ). Overall, we can move from network of individuals to the network of homogeneous patches by replacing the transition probabilities with the product of migration and fixation probabilities.Two-patch meta-populationThe simplest non-trivial case is the fixation probability in a two-patch meta-population with different local size for small migration probability (lambda ). If the migration probability (lambda ) is very small, a new mutant first needs to take over its own patch and only then the first migrant arrives in the second patch. To be more precise, the time between two migration events has to be much higher than the typical time that it takes for the migrant to take over the patch or go extinct again38. In this case, we can divide the dynamics into two phases: A first phase in which a mutant invades one patch and a second phase in which a homogeneous patch of mutants invades the whole meta-population. Assume a new mutation arises in patch 1. Only if this mutant reaches fixation in patch 1, it also has a chance to reach fixation in patch 2. When patch 1 consists of only mutants and patch 2 consists of only wild-types, there are two possibilities for the ultimate fate of the mutant:

(i)

Eventually, the offspring of one mutant selected from patch 1 for reproduction will migrate to patch 2 and reach fixation there. The wild-type goes extinct. This happens with probability ( frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)).

(ii)

Eventually, the offspring of one wild-type selected from patch 2 for reproduction will migrate to patch 1 and the mutant goes extinct. This occurs with probability ( frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})).

Therefore, the probability that a single mutant arising in patch 1 reaches fixation in the entire population is $$begin{aligned} phi _{mathrm{wm}}^{N_1}(r) frac{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)}{frac{N_1 r}{N_1 r+N_2} phi _{mathrm{wm}}^{N_2}(r)+frac{N_2}{N_1r+N_2} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) }=phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) frac{1 }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1} left( tfrac{1}{r}right) }. end{aligned}$$
(3a)
Similarly the probability that a mutant arising in patch 2 takes over the whole population equals$$begin{aligned} phi _{mathrm{wm}}^{N_2}(r) phi _{mathrm{wm}}^{N_1}(r) frac{1 }{phi _{mathrm{wm}}^{N_1}(r)+frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }. end{aligned}$$
(3b)
If we assume that the mutant arises in a patch with a probability proportional to the patch size, the average fixation probability (phi _{bullet !!-!!bullet }) in a two patch population for small migration probability is the weighted sum of Eqs. (3a) and (3b),$$begin{aligned} phi _{bullet !!-!!bullet }&= phi _{mathrm{wm}}^{N_1}(r) phi _{mathrm{wm}}^{N_2}(r) nonumber \&quad times left( frac{frac{N_1}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_2}(r) +frac{N_2}{N_1} frac{1}{r}phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } +frac{frac{N_2}{N_1+N_2} }{ phi _{mathrm{wm}}^{N_1}(r) +frac{N_1}{N_2} frac{1}{r} phi _{mathrm{wm}}^{N_2}left( tfrac{1}{r}right) }right) . end{aligned}$$
(4)
In the case of neutrality, (r=1), we recover (phi _{bullet !!-!!bullet } = frac{1}{N_1+N_2})—the fixation probability in a population of the total size of the two patches. For identical patch sizes, ( N_1=N_2 ), Eq. (4) simplifies to$$begin{aligned} phi _{bullet !!-!!bullet } = left( phi _{mathrm{wm}}^{N_1}(r)right) ^2 frac{1}{phi _{mathrm{wm}}^{N_1}(r)+frac{1}{r} phi _{mathrm{wm}}^{N_1}left( tfrac{1}{r}right) } = phi _{mathrm{wm}}^{2 N_1}(r), end{aligned}$$
(5)
where the simplification to the fixation probability within a single population of size (2N_1) reflects the validity of the isothermal theorem.For (N_1 ne N_2), we approximate Eq. (4) for weak and strong selection. Let us first consider highly advantageous mutants, (r gg 1). In this case, we have (phi _{mathrm{wm}}^{N_1}(r) gg phi _{mathrm{wm}}^{N_1}(tfrac{1}{r})) and thus we can neglect the possibility that a wild-type takes over a mutant patch if patch sizes are sufficiently large. The probability (phi _{bullet !!-!!bullet } ) then becomes a weighted average reflecting patch sizes. For identical patch size (N_1=N_2 = N/2), it reduces to (phi _{bullet !!-!!bullet } approx phi _{mathrm{wm}}^{N_1}(r)=phi _{mathrm{wm}}^{N/2}(r)). In other words, taking over the first patch is sufficient to make fixation in the entire population certain. For patches of very different size, (N_1 gg N_2), we have (N approx N_1) and find (phi _{bullet !!-!! bullet } approx phi _{mathrm{wm}}^{N}(r), ) which implies that fixation is driven by the fixation process in the larger patch, regardless of where the mutant arises. Note that there is a difference between the case of identical patch size and very different patch size . The case of highly disadvantageous mutants, (r ll 1), can be handled in a very similar way.Next, we consider weak selection, (r approx 1). We can approximate the fixation probability as (phi _{mathrm{wm}}^{N}(r^{pm 1}) approx frac{1}{N} pm frac{N-1}{2N} (r-1)). With this, we find$$begin{aligned} phi _{bullet !!-!!bullet } approx frac{1}{N_1+N_2} +frac{1}{2} left( 1 – frac{1}{N_1+N_2} -frac{(N_1-N_2)^2}{(N_1^2+N_2^2)^2} N_1 N_2right) (r-1). end{aligned}$$
(6)
For identical patch size (N_1=N_2 = N/2), this reduces to$$begin{aligned} phi _{bullet !!-!!bullet } approx tfrac{1}{N} +tfrac{N-1}{2N} (r-1), end{aligned}$$
(7)
which is the known result for a single population of size (N=N_1+N_2). When patches have very different size, (N_1 gg N_2) such that (N approx N_1), we recover the same result. Thus, the difference between the fixation probability of a two-patch meta-population with identical patch size and the fixation probability of a two-patch meta-population with very different patch size that we found for highly advantageous mutants is no longer observed for weak selection.When migration probabilities become larger, our approximation is no longer valid and we need to rely on numerical approaches. Figure 2 illustrates the difference between the fixation probability of a two-patch structure meta-population and the equivalent well-mixed population of size (N_1+N_2 ) when migration is low using Eq. (4) and comparing with the numerical approach in Ref.39.While the fixation probability of the two-patch meta-population is very close to the fixation probability of the well-mixed population40, a close inspection reveals an interesting property: For low migration probabilities and (N_1 ne N_2), the two patch structure is a suppressor of selection in the original sense of Lieberman et al.1: For advantageous mutations, (r >1), it decreases the fixation probability, whereas for disadvantageous mutations, (r1) and negative for (r1 ) the minimum fixation probability occurs when the two patch sizes are identical, ( N_1=N_2=N/2 ). Similarly, for fitness values ( r More

1.Liu, M. et al. Evaluation of high-resolution satellite rainfall products using rain gauge data over complex terrain in southwest China. Theoret. Appl. Climatol. 119, 203–219. https://doi.org/10.1007/s00704-014-1092-4 (2014).ADS 
Article 

Google Scholar 
2.Piao, S. et al. Leaf onset in the northern hemisphere triggered by daytime temperature. Nat. Commun. 6, 6911. https://doi.org/10.1038/ncomms7911 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
3.Zheng, Z. et al. Continuous but diverse advancement of spring-summer phenology in response to climate warming across the Qinghai-Tibetan Plateau. Agric. For. Meteorol. 223, 194–202. https://doi.org/10.1016/j.agrformet.2016.04.012 (2016).ADS 
Article 

Google Scholar 
4.Shen, M. et al. Can changes in autumn phenology facilitate earlier green-up date of northern vegetation?. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2020.108077 (2020).Article 

Google Scholar 
5.Zhang, G., Zhang, Y., Dong, J. & Xiao, X. Green-up dates in the Tibetan Plateau have continuously advanced from 1982 to 2011. Proc. Natl. Acad. Sci. U. S. A. 110, 4309–4314. https://doi.org/10.1073/pnas.1210423110 (2013).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
6.Zhu, Z. C. et al. Greening of the Earth and its drivers. Nat. Clim. Chang. 6, 791–795. https://doi.org/10.1038/nclimate3004 (2016).ADS 
CAS 
Article 

Google Scholar 
7.Sun, Q., Li, B., Zhou, G., Jiang, Y. & Yuan, Y. Delayed autumn leaf senescence date prolongs the growing season length of herbaceous plants on the Qinghai-Tibetan Plateau. Agric. For. Meteorol. 284, 1. https://doi.org/10.1016/j.agrformet.2019.107896 (2020).Article 

Google Scholar 
8.Gao, Q. et al. Climatic change controls productivity variation in global grasslands. Sci. Rep. 6, 26958. https://doi.org/10.1038/srep26958 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
9.Zhang, K. et al. Vegetation greening and climate change promote multidecadal rises of global land evapotranspiration. Sci. Rep. 5, 15956. https://doi.org/10.1038/srep15956 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
10.Li, Z., Chen, Y., Wang, Y. & Fang, G. Dynamic changes in terrestrial net primary production and their effects on evapotranspiration. Hydrol. Earth Syst. Sci. 20, 2169–2178. https://doi.org/10.5194/hess-20-2169-2016 (2016).ADS 
Article 

Google Scholar 
11.Jeong, S.-J., Ho, C.-H., Gim, H.-J. & Brown, M. E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008. Global Change Biol. 17, 2385–2399. https://doi.org/10.1111/j.1365-2486.2011.02397.x (2011).ADS 
Article 

Google Scholar 
12.Wang, Y., Gao, Q., Liu, T., Tian, Y. & Yu, M. The greenness of major shrublands in china increased from 2001 to 2013. Remote Sens. https://doi.org/10.3390/rs8020121 (2016).Article 

Google Scholar 
13.Xu, X. et al. Plant community structure regulates responses of prairie soil respiration to decadal experimental warming. Glob. Change Biol. 21, 3846–3853. https://doi.org/10.1111/gcb.12940 (2015).ADS 
Article 

Google Scholar 
14.Gang, C. et al. Drought-induced dynamics of carbon and water use efficiency of global grasslands from 2000 to 2011. Ecol. Ind. 67, 788–797. https://doi.org/10.1016/j.ecolind.2016.03.049 (2016).CAS 
Article 

Google Scholar 
15.Yao, J., Yang, Q., Mao, W., Zhao, Y. & Xu, X. Precipitation trend–Elevation relationship in arid regions of the China. Glob. Planet. Change 143, 1–9. https://doi.org/10.1016/j.gloplacha.2016.05.007 (2016).ADS 
Article 

Google Scholar 
16.Yuan, X. et al. Vegetation changes and land surface feedbacks drive shifts in local temperatures over Central Asia. Sci. Rep. 7, 3287. https://doi.org/10.1038/s41598-017-03432-2 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Liu, Y., Kumar, M., Katul, G. G., Feng, X. & Konings, A. G. Plant hydraulics accentuates the effect of atmospheric moisture stress on transpiration. Nat. Clim. Chang. 10, 691–695. https://doi.org/10.1038/s41558-020-0781-5 (2020).ADS 
CAS 
Article 

Google Scholar 
18.Li, Y. et al. Evaluation and projection of snowfall changes in High Mountain Asia based on NASA’s NEX-GDDP high-resolution daily downscaled dataset. Environ. Res. Lett. https://doi.org/10.1088/1748-9326/aba926 (2020).Article 

Google Scholar 
19.Yao, T. et al. Chained impacts on modern environment of interaction between Westerlies and Indian Monsoon on Tibetan Plateau. Bull. Chin. Acad. Sci. 32, 976–984. https://doi.org/10.16418/j.issn.1000-3045.2017.09.007 (2017).Article 

Google Scholar 
20.Pritchard, H. D. Asia’s shrinking glaciers protect large populations from drought stress. Nature 569, 649–654. https://doi.org/10.1038/s41586-019-1240-1 (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
21.He, Z. et al. Assessing temperature sensitivity of subalpine shrub phenology in semi-arid mountain regions of China. Agric. For. Meteorol. 213, 42–52. https://doi.org/10.1016/j.agrformet.2015.06.013 (2015).ADS 
Article 

Google Scholar 
22.Zhou, J. et al. Alpine vegetation phenology dynamic over 16years and its covariation with climate in a semi-arid region of China. Sci. Total Environ. 572, 119–128. https://doi.org/10.1016/j.scitotenv.2016.07.206 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
23.Zhao, J. et al. Increased precipitation offsets the negative effect of warming on plant biomass and ecosystem respiration in a Tibetan alpine steppe. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107761 (2019).Article 
PubMed 

Google Scholar 
24.Deng, H., Pepin, N. C. & Chen, Y. Changes of snowfall under warming in the Tibetan Plateau. J. Geophys. Res. Atmos. 122, 7323–7341. https://doi.org/10.1002/2017jd026524 (2017).ADS 
Article 

Google Scholar 
25.Yao, T. Tackling on environmental changes in Tibetan Plateau with focus on water, ecosystem and adaptation. Sci. Bull. 64, 1. https://doi.org/10.1016/j.scib.2019.03.033 (2019).Article 

Google Scholar 
26.Shen, M. et al. Strong impacts of daily minimum temperature on the green-up date and summer greenness of the Tibetan Plateau. Glob. Change Biol. 22, 3057–3066. https://doi.org/10.1111/gcb.13301 (2016).ADS 
Article 

Google Scholar 
27.Piao, S. et al. Plant phenology and global climate change: Current progresses and challenges. Glob. Change Biol. 25, 1922–1940. https://doi.org/10.1111/gcb.14619 (2019).ADS 
Article 

Google Scholar 
28.Xu, M. & Xue, X. A research on summer vegetation characteristics & short-time responses to experimental warming of alpine meadow in the Qinghai-Tibetan Plateau. Acta Ecol. Sin. 33, 2071–2083. https://doi.org/10.5846/stxb201112201935 (2013).Article 

Google Scholar 
29.Huang, N., He, J. S., Chen, L. & Wang, L. No upward shift of alpine grassland distribution on the Qinghai-Tibetan Plateau despite rapid climate warming from 2000 to 2014. Sci. Total Environ. 625, 1361–1368. https://doi.org/10.1016/j.scitotenv.2018.01.034 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
30.Piao, S. et al. Evidence for a weakening relationship between interannual temperature variability and northern vegetation activity. Nat. Commun. 5, 5018. https://doi.org/10.1038/ncomms6018 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
31.Huang, M. et al. Air temperature optima of vegetation productivity across global biomes. Nat. Ecol. Evol. 3, 1. https://doi.org/10.1038/s41559-019-0838-x (2019).CAS 
Article 

Google Scholar 
32.Liu, H., Zhang, M., Lin, Z. & Xu, X. Spatial heterogeneity of the relationship between vegetation dynamics and climate change and their driving forces at multiple time scales in Southwest China. Agric. For. Meteorol. 256–257, 10–21. https://doi.org/10.1016/j.agrformet.2018.02.015 (2018).ADS 
Article 

Google Scholar 
33.Chen, Z., Wang, W. & Fu, J. Vegetation response to precipitation anomalies under different climatic and biogeographical conditions in China. Sci. Rep. 10, 830. https://doi.org/10.1038/s41598-020-57910-1 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
34.Guo, H. et al. Space-time characterization of drought events and their impacts on vegetation in Central Asia. J. Hydrol. 564, 1165–1178. https://doi.org/10.1016/j.jhydrol.2018.07.081 (2018).ADS 
Article 

Google Scholar 
35.Li, P., Hu, Z. & Liu, Y. Shift in the trend of browning in Southwestern Tibetan Plateau in the past two decades. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2020.107950 (2020).Article 

Google Scholar 
36.Liu, Z., Li, C., Zhou, P. & Chen, X. A probabilistic assessment of the likelihood of vegetation drought under varying climate conditions across China. Sci. Rep. 6, 35105. https://doi.org/10.1038/srep35105 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
37.Gao, Q.-Z., Li, Y., Xu, H.-M., Wan, Y.-F. & Jiangcun, W.-Z. Adaptation strategies of climate variability impacts on alpine grassland ecosystems in Tibetan Plateau. Mitig. Adapt. Strat. Glob. Change 19, 199–209. https://doi.org/10.1007/s11027-012-9434-y (2012).CAS 
Article 

Google Scholar 
38.Guo, Y. & Wang, C. Trends in precipitation recycling over the Qinghai-Xizang Plateau in last decades. J. Hydrol. 517, 826–835. https://doi.org/10.1016/j.jhydrol.2014.06.006 (2014).ADS 
Article 

Google Scholar 
39.Schlaepfer, D. R. et al. Climate change reduces extent of temperate drylands and intensifies drought in deep soils. Nat. Commun. 8, 14196. https://doi.org/10.1038/ncomms14196 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
40.Yao, J. et al. Climatic and associated atmospheric water cycle changes over the Xinjiang, China. J. Hydrol. 585, 1. https://doi.org/10.1016/j.jhydrol.2020.124823 (2020).Article 

Google Scholar 
41.Sun, A. et al. Quantified hydrological responses to permafrost degradation in the headwaters of the Yellow River (HWYR) in High Asia. Sci. Total Environ. 712, 135632. https://doi.org/10.1016/j.scitotenv.2019.135632 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
42.Brun, F., Berthier, E., Wagnon, P., Kaab, A. & Treichler, D. A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000–2016. Nat. Geosci. 10, 668–673. https://doi.org/10.1038/NGEO2999 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
43.Luo, D., Liu, L., Jin, H., Wang, X. & Chen, F. Characteristics of ground surface temperature at Chalaping in the Source Area of the Yellow River, northeastern Tibetan Plateau. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107819 (2020).Article 

Google Scholar 
44.Che, M. et al. Spatial and temporal variations in the end date of the vegetation growing season throughout the Qinghai-Tibetan Plateau from 1982 to 2011. Agric. For. Meteorol. 189–190, 81–90. https://doi.org/10.1016/j.agrformet.2014.01.004 (2014).ADS 
Article 

Google Scholar 
45.Ji, Z. et al. Investigation of mineral aerosols radiative effects over High Mountain Asia in 1990–2009 using a regional climate model. Atmos. Res. 178–179, 484–496. https://doi.org/10.1016/j.atmosres.2016.05.003 (2016).CAS 
Article 

Google Scholar 
46.Wang, X. et al. Spring temperature change and its implication in the change of vegetation growth in North America from 1982 to 2006. Proc. Natl. Acad. Sci. U. S. A. 108, 1240–1245. https://doi.org/10.1073/pnas.1014425108 (2011).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Piao, S. et al. Responses and feedback of the Tibetan Plateau’s alpine ecosystem to climate change. Chin. Sci. Bull. 64, 2842–2855. https://doi.org/10.1360/TB-2019-0074 (2019).Article 

Google Scholar 
48.Zeng, Z. et al. Climate mitigation from vegetation biophysical feedbacks during the past three decades. Nat. Clim. Change 7, 432–436. https://doi.org/10.1038/nclimate3299 (2017).ADS 
Article 

Google Scholar 
49.Xu, H. J., Wang, X. P. & Yang, T. B. Trend shifts in satellite-derived vegetation growth in Central Eurasia, 1982–2013. Sci. Total Environ. 579, 1658–1674. https://doi.org/10.1016/j.scitotenv.2016.11.182 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
50.Zhang, Y. et al. Satellite-observed global terrestrial vegetation production in response to water availability. Remote Sens. 13, 1. https://doi.org/10.3390/rs13071289 (2021).Article 

Google Scholar 
51.Curio, J. & Scherer, D. Seasonality and spatial variability of dynamic precipitation controls on the Tibetan Plateau. Earth Syst. Dyn. Discus. https://doi.org/10.5194/esd-2016-1,10.5194/esd-2016-1 (2016).Article 

Google Scholar 
52.Li, J., Sun, C. & Jin, F. F. NAO implicated as a predictor of Northern Hemisphere mean temperature multidecadal variability. Geophys. Res. Lett. 40, 5497–5502. https://doi.org/10.1002/2013gl057877 (2013).ADS 
Article 

Google Scholar 
53.Turner, A. G. & Annamalai, H. Climate change and the South Asian summer monsoon. Nat. Clim. Chang. 2, 587–595. https://doi.org/10.1038/nclimate1495 (2012).ADS 
Article 

Google Scholar 
54.Crimmins, T. M., Crimmins, M. A. & DavidBertelsen, C. Complex responses to climate drivers in onset of spring flowering across a semi-arid elevation gradient. J. Ecol. 98, 1042–1051. https://doi.org/10.1111/j.1365-2745.2010.01696.x (2010).Article 

Google Scholar 
55.Du, J. et al. Interacting effects of temperature and precipitation on climatic sensitivity of spring vegetation green-up in arid mountains of China. Agric. For. Meteorol. 269–270, 71–77. https://doi.org/10.1016/j.agrformet.2019.02.008 (2019).ADS 
Article 

Google Scholar 
56.Huang, J. et al. Global semi-arid climate change over last 60 years. Clim. Dyn. 46, 1131–1150. https://doi.org/10.1007/s00382-015-2636-8 (2015).Article 

Google Scholar 
57.Sun, J., Qin, X. & Yang, J. The response of vegetation dynamics of the different alpine grassland types to temperature and precipitation on the Tibetan Plateau. Environ. Monit. Assess. 188, 20. https://doi.org/10.1007/s10661-015-5014-4 (2016).Article 
PubMed 

Google Scholar 
58.Ganjurjav, H. et al. Differential response of alpine steppe and alpine meadow to climate warming in the central Qinghai-Tibetan Plateau. Agric. For. Meteorol. 223, 233–240. https://doi.org/10.1016/j.agrformet.2016.03.017 (2016).ADS 
Article 

Google Scholar 
59.Xu, M. et al. Year-round warming and autumnal clipping lead to downward transport of root biomass, carbon and total nitrogen in soil of an alpine meadow. Environ. Exp. Bot. 109, 54–62. https://doi.org/10.1016/j.envexpbot.2014.07.012 (2015).CAS 
Article 

Google Scholar 
60.Xie, J. et al. Land surface phenology and greenness in Alpine grasslands driven by seasonal snow and meteorological factors. Sci. Total Environ. 725, 138380. https://doi.org/10.1016/j.scitotenv.2020.138380 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
61.Zhang, Y. et al. Vegetation dynamics and its driving forces from climate change and human activities in the Three-River Source Region, China from 1982 to 2012. Sci. Total Environ. 563–564, 210–220. https://doi.org/10.1016/j.scitotenv.2016.03.223 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
62.Liu, L. et al. Elevation-dependent decline in vegetation greening rate driven by increasing dryness based on three satellite NDVI datasets on the Tibetan Plateau. Ecol. Indic. https://doi.org/10.1016/j.ecolind.2019.105569 (2019).Article 

Google Scholar 
63.Piao, S. et al. Altitude and temperature dependence of change in the spring vegetation green-up date from 1982 to 2006 in the Qinghai-Xizang Plateau. Agric. For. Meteorol. 151, 1599–1608. https://doi.org/10.1016/j.agrformet.2011.06.016 (2011).ADS 
Article 

Google Scholar 
64.Zhang, X., Tarpley, D. & Sullivan, J. T. Diverse responses of vegetation phenology to a warming climate. Geophys. Res. Lett. https://doi.org/10.1029/2007gl031447 (2007).Article 

Google Scholar 
65.Gao, Y. et al. Vegetation net primary productivity and its response to climate change during 2001–2008 in the Tibetan Plateau. Sci. Total Environ. 444, 356–362. https://doi.org/10.1016/j.scitotenv.2012.12.014 (2013).ADS 
CAS 
Article 
PubMed 

Google Scholar 
66.Shen, M. et al. Influences of temperature and precipitation before the growing season on spring phenology in grasslands of the central and eastern Qinghai-Tibetan Plateau. Agric. For. Meteorol. 151, 1711–1722. https://doi.org/10.1016/j.agrformet.2011.07.003 (2011).ADS 
Article 

Google Scholar 
67.Chen, N. et al. The compensation effects of post-drought regrowth on earlier drought loss across the tibetan plateau grasslands. Agric. For. Meteorol. https://doi.org/10.1016/j.agrformet.2019.107822 (2020).Article 

Google Scholar 
68.Zhao, W. et al. Contributions of climatic factors to interannual variability of the vegetation index in Northern China Grasslands. J. Clim. 33, 175–183. https://doi.org/10.1175/jcli-d-18-0587.1 (2020).ADS 
Article 

Google Scholar 
69.Liang, J. et al. Where will threatened migratory birds go under climate change? Implications for China’s national nature reserves. Sci. Total Environ. 645, 1040–1047. https://doi.org/10.1016/j.scitotenv.2018.07.196 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
70.Qu, S. et al. What drives the vegetation restoration in Yangtze River basin, China: Climate change or anthropogenic factors?. Ecol. Ind. 90, 438–450. https://doi.org/10.1016/j.ecolind.2018.03.029 (2018).Article 

Google Scholar 
71.Yin, L. et al. What drives the vegetation dynamics in the Hengduan Mountain region, southwest China: Climate change or human activity?. Ecol. Ind. 112, 106013. https://doi.org/10.1016/j.ecolind.2019.106013 (2020).Article 

Google Scholar 
72.Zhou, X., Yamaguchi, Y. & Arjasakusuma, S. Distinguishing the vegetation dynamics induced by anthropogenic factors using vegetation optical depth and AVHRR NDVI: A cross-border study on the Mongolian Plateau. Sci. Total Environ. 616–617, 730–743. https://doi.org/10.1016/j.scitotenv.2017.10.253 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
73.Li, Y. et al. The effects of fencing on carbon stocks in the degraded alpine grasslands of the Qinghai-Tibetan Plateau. J. Environ. Manag. 128, 393–399. https://doi.org/10.1016/j.jenvman.2013.05.058 (2013).CAS 
Article 

Google Scholar 
74.Liu, X. et al. How does grazing exclusion influence plant productivity and community structure in alpine grasslands of the Qinghai-Tibetan Plateau?. Glob. Ecol. Conserv. https://doi.org/10.1016/j.gecco.2020.e01066 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
75.Li, W. et al. Effects of grazing regime on vegetation structure, productivity, soil quality, carbon and nitrogen storage of alpine meadow on the Qinghai-Tibetan Plateau. Ecol. Eng. 98, 123–133. https://doi.org/10.1016/j.ecoleng.2016.10.026 (2017).Article 

Google Scholar 
76.Deng, L. et al. Effects of grazing exclusion on carbon sequestration in China’s grassland. Earth Sci. Rev. 173, 84–95. https://doi.org/10.1016/j.earscirev.2017.08.008 (2017).ADS 
CAS 
Article 

Google Scholar 
77.Yu, L. et al. Effects of grazing exclusion on soil carbon dynamics in alpine grasslands of the Tibetan Plateau. Geoderma 353, 133–143. https://doi.org/10.1016/j.geoderma.2019.06.036 (2019).ADS 
CAS 
Article 

Google Scholar 
78.Shao, Q. et al. Effects of an ecological conservation and restoration project in the Three-River Source Region, China. J. Geograph. Sci. 27, 183–204. https://doi.org/10.1007/s11442-017-1371-y (2016).Article 

Google Scholar 
79.Sun, Q. et al. A systematic review of research studies on the estimation of net primary productivity in the Three-River Headwater Region, China. J. Geograph. Sci. 27, 161–182. https://doi.org/10.1007/s11442-017-1370-z (2016).Article 

Google Scholar 
80.Shen, X. et al. Marshland loss warms local land surface temperature in China. Geophys. Res. Lett. https://doi.org/10.1029/2020GL087648 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
81.Shen, X. et al. Aboveground biomass and its spatial distribution pattern of herbaceous marsh vegetation in China. Sci. China Earth Sci. 64, 1115–1125. https://doi.org/10.1007/s11430-020-9778-7 (2021).ADS 
Article 

Google Scholar 
82.Wang, Y. et al. Spatiotemporal change of aboveground biomass and its response to climate change in marshes of the Tibetan Plateau. Int. J. Appl. Earth Observ. Geoinf. https://doi.org/10.1016/j.jag.2021.102385 (2021).Article 

Google Scholar 
83.Jeong, S.-J., Ho, C.-H. & Jeong, J.-H. Increase in vegetation greenness and decrease in springtime warming over east Asia. Geophys. Res. Lett. https://doi.org/10.1029/2008gl036583 (2009).Article 

Google Scholar 
84.Shen, M. et al. Evaporative cooling over the Tibetan Plateau induced by vegetation growth. Proc. Natl. Acad. Sci. U. S. A. 112, 9299–9304. https://doi.org/10.1073/pnas.1504418112 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
85.Shen, X. et al. Asymmetric effects of daytime and nighttime warming on spring phenology in the temperate grasslands of China. Agric. For. Meteorol. 259, 240–249. https://doi.org/10.1016/j.agrformet.2018.05.006 (2018).ADS 
Article 

Google Scholar 
86.Shen, X. et al. Spatiotemporal variation in vegetation spring phenology and its response to climate change in freshwater marshes of Northeast China. Sci. Total Environ. 666, 1169–1177. https://doi.org/10.1016/j.scitotenv.2019.02.265 (2019).ADS 
CAS 
Article 
PubMed 

Google Scholar 
87.Niittynen, P. et al. Fine-scale tundra vegetation patterns are strongly related to winter thermal conditions. Nat. Clim. Chang. https://doi.org/10.1038/s41558-020-00916-4 (2020).Article 

Google Scholar 
88.Wu, D. et al. Evaluation of spatiotemporal variations of global fractional vegetation cover based on GIMMS NDVI data from 1982 to 2011. Remote Sens. 6, 4217–4239. https://doi.org/10.3390/rs6054217 (2014).ADS 
Article 

Google Scholar 
89.Zhang, H. et al. Calculation of evapotranspiration in different climatic zones combining the long-term monitoring data with bootstrap method. Environ. Res. 191, 110200. https://doi.org/10.1016/j.envres.2020.110200 (2020).CAS 
Article 
PubMed 

Google Scholar 
90.Kalisa, W. et al. Assessment of climate impact on vegetation dynamics over East Africa from 1982 to 2015. Sci. Rep. 9, 16865. https://doi.org/10.1038/s41598-019-53150-0 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
91.Chen, Y. Geographical data analysis with Matlab 202–220 (Chen, 2012). More

Technological advancements improve our ability to manage natural resources. This is particularly relevant for small scale fisheries, where there is a need for low-cost data sources to improve our understanding of fishing effort, catch, and the associated sustainability of fish resources required for global food security. GPS trackers have now been widely used to study the behaviour of small-scale fisheries25,40. We found focusing on patterns of vessel movement to be a low-cost, reliable approach to identify fishing grounds, as well as to understand both the spatial and temporal usage of AFADs, and ultimately predicting the resulting catch.We acknowledge that the number of actual AFADs used by our tracked vessels is likely much higher than the number estimated in this study. This is in part due to our requirement for a potential AFAD to have been visited at least two times before we considered it a confirmed AFAD. These criteria significantly reduced the number of AFADs reported (from 139 to 72 AFADs). However, we erred on the cautious side as we were unable to distinguish between AFAD fishing and other non-AFAD fishing behaviours, such as bait fishing, that might involve vessels being stationary. Furthermore, given that the length of trip for a vessel is 5 to 20 days, the one-month period over which a SPOT Trace tracker is deployed means there is a maximum of two fishing trips possible during our observation period. This leads to the potential that even the tracked vessels may have additional AFADs they use outside of the fishing trips observed in this study period.Another source of underestimation in AFAD numbers may come from the distance parameter we employed in our analysis. During the ground-truthing, only two out of three visited AFADs were detected by DBSCAN. This is because the radius of movement between two of the FADs (Fig. 2) was overlapping. This is possible, as currents and winds displace AFADs synchronously, and thus tangling is reduced, allowing AFADs to be deployed closer together than the sum of their surface radii. Therefore, the distance among vessel positions clustered two AFADs, identifying them as a single AFAD, given the criteria we applied. The implication of this potential for multiple AFADs within a DBSCAN cluster is that the locations we detected could actually represent a much larger number of AFADs that are deployed close together. Future extensions of this work could include estimating the number of AFADs within clusters using the geometric pattern of the boundary of the cluster. For instance, a figure-eight shaped boundary would indicate there are two FADs in a cluster rather than one. However, SPOT Trace deployments would need to be longer to provide adequate data to distinguish this subtlety.Since our study did not include records from the first time each of the AFADs were deployed, we were unable to determine the absolute lifetime of AFADs in the region. However, based on the vessel tracking data, only a few AFADs were visited for nearly one year implying that AFADs might be failing in less than one year. Because the record of AFAD usage is from the vessel perspective, when the tracker on a vessel is removed at the end of its month long deployment, the record stops while the AFADs may still exist and remain in use. If other vessels in the study use the same AFAD, the record for that AFAD will continue, but if not, it ends with the removal of the tracker from the vessel using it. Hence, the lifespan of AFADs we report is an estimate that should be treated as a minimum lifespan. Moreover, since fishers tend to deploy AFADs in a particular fishing location, it is also possible that the fisher has deployed a new AFAD in the same location. However, given the deployment precision required this may not be as big of a source of error as underestimation.Conversely, from long periods of inactivity at individual AFADs (as shown in Fig. 4), we suspect that some AFADs may have been lost and replaced over the course of the longer use patterns we observed. These inactivity periods take place during the wet season, which typically has rougher weather and poorer fishing conditions, particularly for small vessels. Hence, we might anticipate fewer vessel days at sea or the loss of AFADs due to failure of their moorings during periods of high swell. The asynchrony in the time at which inactivity patterns begin and end, however, suggests that a lack of fishing activity is unlikely to be the sole source of the observed inactivity periods and that there is likely a contribution from AFAD loss and replacement. With additional tracking data on individual vessels, it might be possible to disentangle these differences by looking for subtle shifts in the centres of the spatial clusters, indicating a new deployment. However, the current observations are inadequate to provide this level of resolution.The AFAD sharing practices identified in our study reveal a management opportunity to reduce the number of AFADs deployed. The use of AFADs can be maximized by extending the users beyond the owners of an individual AFAD, or by considering AFADs a community resource. While perhaps not suitable in all areas, given that sharing AFAD is relatively widespread, this presents a viable option. Developing a management system that allows limits on the total number of AFADs but provides for a system of rotating access may allow for the establishment of a biologically sustainable system of AFADs whilst minimizing social and economic disruption to the fishers. Moreover, it may also reduce the incentives for fishers to keep AFAD locations private.The catch data obtained from the port sampling allowed us to identify the factors that influence the total catch. The number of AFADs visited is the main factor that significantly affects the weight of catch by a vessel on a fishing trip, given the average catch of a vessel. Trip success increased as more AFADs were visited, but then declined sharply beyond 3 AFADs. Similarly, for a given vessel, as trip lengths increased, catches were lower.This pattern might be expected if fishers are considered as central place foragers in the context of the optimal foraging theory41. Vessels typically leave and return to the same port. Presumably while at sea, they attempt to either maximize their catch or at least satisfy a minimum required catch to meet their fixed costs. In either event, one would expect fishers to extend their trip length if catch rates are low to try to meet their objective, subject to other constraints such as fuel supply or adverse weather. In this context, if they visit an AFAD and have a low catch rate, one would expect fishers to move to another AFAD. Thus together, the number of AFADs visited and the length of the trip provide a reliable predictor of the quality of a fishing trip, in terms of variation around the average for a given vessel. This information is very useful, as it suggests that the SPOT Tracking data, or other vessel tracking information, can be used as a proxy for port sampling. Thus, remote monitoring of the vessels can be used to get some measure of stock status, via catch rates, or as a check against port sampling or logbooks to check their veracity. Given the rapidly falling cost of technologies, such as the SPOT trackers used in this study, proxies for catch rates such as the one we developed here could facilitate fleet-wide monitoring. In Indonesia, with a quarter-million small scale vessels spread across thousands of islands this scalability is critical, and given Indonesia has the third highest marine catch in the world42, the resulting management improvements have global ramifications.The case of Indonesian FAD management challenges reflects current global FAD management challenges, especially in artisanal coastal fisheries in Pacific island countries where AFADs are commonly used43. We found that AFAD deployments in Indonesia are very dense, and frequently well inside the minimum ten nautical miles spacing required by law. Based on our study, it is also clear that vessels are using more than the three AFADs limit allowed in current regulations. These high densities and usage rates could be reducing the effectiveness of AFADs to aggregate the fish by dividing the fish concentration among close AFADs and thus decreasing catch rates. Moreover, the current concentrated use of AFADs could also be leading to large numbers of lost and abandoned AFAD structures, with significant impacts on the ecosystem and local habitats44,45. Fishers could deploy fewer AFADs, thus decreasing their potential impacts. The regulation of AFADs in Indonesia, which has been in place since 2014, is still not effectively enforced due to technical issues. Moreover, the users of this type of FAD are dominated by small-scale fishers whose livelihoods and food supplies likely depend on the additional efficiency, making management more problematic.Expansion of the current study from a monthly sampling approach to continuous monitoring of vessels would greatly improve our ability to discern AFAD use patterns, infer catch dynamics, and ultimately investigate the potential for management strategies that could balance maximizing the benefits from AFAD deployments and controlling their environmental and social impacts. Ultimately minor technological improvements which extend tracking device lifetimes, along with links to other electronic monitoring approaches such as low-cost onboard cameras or electronic logbooks and landing records could allow cost effective monitoring of the vast small scale fleet in Indonesia, leading to better fishery outcomes at a significantly reduced cost. Expanding these approaches, particularly in the case of rapidly falling technology costs, has significant promise for improving management across the many fisheries and sectors in Indonesia, and elsewhere.Most of the global FADs are managed by the regional fisheries management organizations (RFMOs), and not all member countries have implemented regulations regarding FAD use46,47 (IOTC, 2018). Given the large proportion of world tuna production which is dominated by floating object fishing, compared to fishing on free schooling tuna48, more investment in FAD management will likely yield an overall improvement in fisheries management and catch sustainability. Paired with addressing management of Indonesia’s very large small scale tuna sector, which lands half the national catch, these regulations could significantly improve sustainability in the Indo-Pacific region. More

1.Walther, G. R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
2.Mooij, W. M. et al. The impact of climate change on lakes in the Netherlands: A review. Aquat. Ecol. 39, 381–400 (2005).CAS 
Article 

Google Scholar 
3.Walter, B., Peters, J. & van Beusekom, J. E. E. The effect of constant darkness and short light periods on the survival and physiological fitness of two phytoplankton species and their growth potential after re-illumination. Aquat. Ecol. 51, 591–603 (2017).CAS 
Article 

Google Scholar 
4.Woodward, G., Perkins, D. M. & Brown, L. E. Climate change and freshwater ecosystems: Impacts across multiple levels of organization. Philos. Trans. R. Soc. B Biol. Sci. 365, 2093–2106 (2010).Article 

Google Scholar 
5.Wagner, H., Fanesi, A. & Wilhelm, C. Title: Freshwater phytoplankton responses to global warming. J. Plant Physiol. 203, 127–134 (2016).CAS 
PubMed 
Article 

Google Scholar 
6.Gilbert, J. A. Some phytoplankton like it hot. Nat. Clim. Change 3, 954–955 (2013).ADS 
Article 

Google Scholar 
7.Hense, I., Meier, H. E. M. & Sonntag, S. Projected climate change impact on Baltic Sea cyanobacteria: Climate change impact on cyanobacteria. Clim. Change 119, 391–406 (2013).CAS 
Article 

Google Scholar 
8.Trombetta, T. et al. Water temperature drives phytoplankton blooms in coastal waters. PLoS One 14, e0214933 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
9.Jin, P. & Agustí, S. Fast adaptation of tropical diatoms to increased warming with trade-offs. Sci. Rep. 8, 17771 (2018).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Pinceel, T., Buschke, F., Weckx, M., Brendonck, L. & Vanschoenwinkel, B. Climate change jeopardizes the persistence of freshwater zooplankton by reducing both habitat suitability and demographic resilience. BMC Ecol. 18, 1–9 (2018).Article 

Google Scholar 
11.Shin, H. R. & Kneitel, J. M. Warming interacts with inundation timing to influence the species composition of California vernal pool communities. Hydrobiologia 843, 93–105 (2019).Article 

Google Scholar 
12.Montrone, A. et al. Climate change impacts on vernal pool hydrology and vegetation in northern California. J. Hydrol. 574, 1003–1013 (2019).ADS 
Article 

Google Scholar 
13.Williams, D. D. The biology of temporary waters. Biol. Tempor. Waters https://doi.org/10.1093/acprof:oso/9780198528128.001.0001 (2007).Article 

Google Scholar 
14.Waterkeyn, A., Grillas, P., Vanschoenwinkel, B. & Brendonck, L. Invertebrate community patterns in Mediterranean temporary wetlands along hydroperiod and salinity gradients. Freshw. Biol. 53, 1808–1822 (2008).CAS 
Article 

Google Scholar 
15.Lemmens, P. et al. How to maximally support local and regional biodiversity in applied conservation? Insights from pond management. PLoS One 8, e72538 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Lischeid, G. et al. Natural ponds in an agricultural landscape: External drivers, internal processes, and the role of the terrestrial-aquatic interface. Limnologica 68, 5–16 (2018).CAS 
Article 

Google Scholar 
17.Mancinelli, G., Mali, S. & Belmonte, G. Species richness and taxonomic distinctness of zooplankton in ponds and small lakes from Albania and North Macedonia: The role of bioclimatic factors. Water (Switzerland) 11, 2384 (2019).
Google Scholar 
18.Gołdyn, B., Kowalczewska-Madura, K. & Celewicz-Gołdyn, S. Drought and deluge: Influence of environmental factors on water quality of kettle holes in two subsequent years with different precipitation. Limnologica 54, 14–22 (2015).Article 
CAS 

Google Scholar 
19.Salmaso, N. & Tolotti, M. Phytoplankton and anthropogenic changes in pelagic environments. Hydrobiologia https://doi.org/10.1007/s10750-020-04323-w (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
20.Celewicz, S., Czyż, M. J. & Gołdy, B. Feeding patterns in Eubranchipus grubii (Dybowski 1860) (Branchiopoda: Anostraca) and its potential influence on the phytoplankton communities of vernal pools. J. Limnol. 77, 276–284 (2018).Article 

Google Scholar 
21.Rasconi, S., Winter, K. & Kainz, M. J. Temperature increase and fluctuation induce phytoplankton biodiversity loss—Evidence from a multi-seasonal mesocosm experiment. Ecol. Evol. 7, 2936–2946 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
22.Celewicz-Goldyn, S. & Kuczynska-Kippen, N. Ecological value of macrophyte cover in creating habitat for microalgae (diatoms) and zooplankton (rotifers and crustaceans) in small field and forest water bodies. PLoS One 12, e0177317 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
23.Kozak, A., Celewicz-Gołdyn, S. & Kuczyńska-Kippen, N. Cyanobacteria in small water bodies: The effect of habitat and catchment area conditions. Sci. Total Environ. 646, 1578–1587 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Iacarella, J. C., Barrow, J. L., Giani, A., Beisner, B. E. & Gregory-Eaves, I. Shifts in algal dominance in freshwater experimental ponds across differing levels of macrophytes and nutrients. Ecosphere 9, e02086 (2018).Article 

Google Scholar 
25.Toseland, A. et al. The impact of temperature on marine phytoplankton resource allocation and metabolism. Nat. Clim. Change 3, 979–984 (2013).ADS 
CAS 
Article 

Google Scholar 
26.Richardson, J. et al. Response of cyanobacteria and phytoplankton abundance to warming, extreme rainfall events and nutrient enrichment. Glob. Change Biol. 25, 3365–3380 (2019).ADS 
Article 

Google Scholar 
27.De Senerpont Domis, L. N., Mooij, W. M. & Huisman, J. Climate-induced shifts in an experimental phytoplankton community: A mechanistic approach. Hydrobiologia 584, 403–413 (2007).Article 

Google Scholar 
28.Boyce, D. G., Lewis, M. R. & Worm, B. Global phytoplankton decline over the past century. Nature 466, 591–596 (2010).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Hinder, S. L. et al. Changes in marine dinoflagellate and diatom abundance under climate change. Nat. Clim. Change 2, 271–275 (2012).ADS 
Article 

Google Scholar 
30.Winder, M. & Sommer, U. Phytoplankton response to a changing climate. Hydrobiologia 698, 5–16 (2012).Article 

Google Scholar 
31.Machado, K. B., Vieira, L. C. G. & Nabout, J. C. Predicting the dynamics of taxonomic and functional phytoplankton compositions in different global warming scenarios. Hydrobiologia 830, 115–134 (2019).CAS 
Article 

Google Scholar 
32.O’Neil, J. M., Davis, T. W., Burford, M. A. & Gobler, C. J. The rise of harmful cyanobacteria blooms: The potential roles of eutrophication and climate change. Harmful Algae 14, 313–334 (2012).Article 
CAS 

Google Scholar 
33.Rasconi, S., Gall, A., Winter, K. & Kainz, M. J. Increasing water temperature triggers dominance of small freshwater plankton. PLoS One 10, e0140449 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
34.Wirth, C., Limberger, R. & Weisse, T. Temperature × light interaction and tolerance of high water temperature in the planktonic freshwater flagellates Cryptomonas (Cryptophyceae) and Dinobryon (Chrysophyceae). J. Phycol. 55, 404–414 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Wang, H. et al. High antioxidant capability interacts with respiration to mediate two Alexandrium species growth exploitation of photoperiods and light intensities. Harmful Algae 82, 26–34 (2019).CAS 
PubMed 
Article 

Google Scholar 
36.Fakhri, M., Arifin, N. B., Budianto, B., Yuniarti, A. & Hariati, A. M. Effect of salinity and photoperiod on growth of microalgae Nannochloropsis sp. and Tetraselmis sp. Nat. Environ. Pollut. Technol. 14, 563–566 (2015).
Google Scholar 
37.Torzillo, G., Sacchi, A. & Materassi, R. Temperature as an important factor affecting productivity and night biomass loss in Spirulina platensis grown outdoors in tubular photobioreactors. Bioresour. Technol. 38, 95–100 (1991).Article 

Google Scholar 
38.Shatwell, T., Köhler, J. & Nicklisch, A. Temperature and photoperiod interactions with phosphorus-limited growth and competition of two diatoms. PLoS One 9, e102367 (2014).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Li, G., Talmy, D. & Campbell, D. A. Diatom growth responses to photoperiod and light are predictable from diel reductant generation. J. Phycol. 53, 95–107 (2017).CAS 
PubMed 
Article 

Google Scholar 
40.Reynolds, C. S. Vegetation Processes in the Pelagic: A Model for Ecosystem Theory. Excellence in Ecology Vol. 77 (Ecology Institute, 1997).
Google Scholar 
41.Elliott, J. A., Jones, I. D. & Thackeray, S. J. Testing the sensitivity of phytoplankton communities to changes in water temperature and nutrient load, in a temperate lake. Hydrobiologia 559, 401–411 (2006).CAS 
Article 

Google Scholar 
42.Jöhnk, K. D. et al. Summer heatwaves promote blooms of harmful cyanobacteria. Glob. Change Biol. 14, 495–512 (2008).ADS 
Article 

Google Scholar 
43.Elliott, J. A. Is the future blue–green? A review of the current model predictions of how climate change could affect pelagic freshwater cyanobacteria. Water Res. 46, 1364–1371 (2012).CAS 
PubMed 
Article 

Google Scholar 
44.Ullah, H., Nagelkerken, I., Goldenberg, S. U. & Fordham, D. A. Climate change could drive marine food web collapse through altered trophic flows and cyanobacterial proliferation. PLoS Biol. 16, e2003446 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Hansson, L. A. et al. Food-chain length alters community responses to global change in aquatic systems. Nat. Clim. Change 3, 228–233 (2013).ADS 
Article 

Google Scholar 
46.Burgmer, T. & Hillebrand, H. Temperature mean and variance alter phytoplankton biomass and biodiversity in a long-term microcosm experiment. Oikos 120, 922–933 (2011).Article 

Google Scholar 
47.Hillebrand, H., Burgmer, T. & Biermann, E. Running to stand still: Temperature effects on species richness, species turnover, and functional community dynamics. Mar. Biol. 159, 2415–2422 (2012).Article 

Google Scholar 
48.Lewandowska, A. M. et al. Responses of primary productivity to increased temperature and phytoplankton diversity. J. Sea Res. 72, 87–93 (2012).ADS 
Article 

Google Scholar 
49.Lewandowska, A. M., Hillebrand, H., Lengfellner, K. & Sommer, U. Temperature effects on phytoplankton diversity—The zooplankton link. J. Sea Res. 85, 359–364 (2014).ADS 
Article 

Google Scholar 
50.Bergkemper, V., Stadler, P. & Weisse, T. Moderate weather extremes alter phytoplankton diversity—A microcosm study. Freshw. Biol. 63, 1211–1224 (2018).CAS 
Article 

Google Scholar 
51.McMinn, A. & Martin, A. Dark survival in a warming world. Proc. R. Soc. B Biol. Sci. 280, 20122909 (2013).CAS 
Article 

Google Scholar 
52.Waibel, A., Peter, H. & Sommaruga, R. Importance of mixotrophic flagellates during the ice-free season in lakes located along an elevational gradient. Aquat. Sci. 81, 1–10 (2019).CAS 
Article 

Google Scholar 
53.Chen, B. Patterns of thermal limits of phytoplankton. J. Plankton Res. 37, 285–292 (2015).Article 

Google Scholar 
54.Reeves, S., McMinn, A. & Martin, A. The effect of prolonged darkness on the growth, recovery and survival of Antarctic sea ice diatoms. Polar Biol. 34, 1019–1032 (2011).Article 

Google Scholar 
55.van de Poll, W. H., Abdullah, E., Visser, R. J. W., Fischer, P. & Buma, A. G. J. Taxon-specific dark survival of diatoms and flagellates affects Arctic phytoplankton composition during the polar night and early spring. Limnol. Oceanogr. 65, 903–914 (2020).ADS 
Article 

Google Scholar 
56.Poniewozik, M. & Juráň, J. Extremely high diversity of euglenophytes in a small pond in eastern Poland. Plant Ecol. Evol. 151, 18–34 (2018).Article 

Google Scholar 
57.Shafik, H. M., Herodek, S., Présing, M. & Vörös, L. Factors effecting growth and cell composition of cyanoprokaryote Cylindrospermopsis raciborskii (Wołoszyńska) Seenayya et Subba Raju. Algol. Stud. Hydrobiol. Suppl. 103, 75–93 (2001).
Google Scholar 
58.Tang, E. P. Y. & Vincent, W. F. Effects of daylength and temperature on the growth and photosynthesis of an Arctic cyanobacterium, Schizothrix calcicola (Oscillatoriaceae). Eur. J. Phycol. 35, 263–272 (2000).Article 

Google Scholar 
59.Agasild, H., Zingel, P., Tõnno, I., Haberman, J. & Nõges, T. Contribution of different zooplankton groups in grazing on phytoplankton in shallow eutrophic Lake Võrtsjärv (Estonia). Hydrobiologia 584, 167–177 (2007).Article 

Google Scholar 
60.Gołdyn, R. & Kowalczewska-Madura, K. Interactions between phytoplankton and zooplankton in the hypertrophic Swarzȩdzkie Lake in western Poland. J. Plankton Res. 30, 33–42 (2008).Article 
CAS 

Google Scholar 
61.Tovar-Sanchez, A., Duarte, C. M., Hernández-León, S. & Sañudo-Wilhelmy, S. A. Krill as a central node for iron cycling in the Southern Ocean. Geophys. Res. Lett. 34, L11601 (2007).ADS 
Article 
CAS 

Google Scholar 
62.Hunt, R. J. & Matveev, V. F. The effects of nutrients and zooplankton community structure on phytoplankton growth in a subtropical Australian reservoir: An enclosure study. Limnologica 35, 90–101 (2005).Article 

Google Scholar 
63.Yvon-Durocher, G. et al. Five years of experimental warming increases the biodiversity and productivity of phytoplankton. PLoS Biol. 13, e1002324 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
64.Gołdyn, B., Chudzińska, M., Barałkiewicz, D. & Celewicz-Gołdyn, S. Heavy metal contents in the sediments of astatic ponds: Influence of geomorphology, hydroperiod, water chemistry and vegetation. Ecotoxicol. Environ. Saf. 118, 103–111 (2015).PubMed 
Article 
CAS 

Google Scholar 
65.IPCC. Climate Change 2007: The Physical Science Basis (Cambridge University Press, 2007).
Google Scholar 
66.Christensen, J. H. & Christensen, O. B. A summary of the PRUDENCE model projections of changes in European climate by the end of this century. Clim. Change 81, 7–30 (2007).ADS 
Article 

Google Scholar 
67.Beniston, M. et al. Future extreme events in European climate: An exploration of regional climate model projections. Clim. Change 81, 71–95 (2007).Article 

Google Scholar 
68.Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2001).
Google Scholar 
69.Arbizu, P. M. pairwiseAdonis: Pairwise Multilevel Comparison Using Adonis. R Packag. version 0.0.1. (2017).70.Rink, B. & Raak, C. J. F. Principal response curves: Analysis of time-dependent multivariate responses of biological community to stress. Environ. Toxicol. Chem. 18, 138–148 (1999).Article 

Google Scholar 
71.Lepš, J. & Šmilauer, P. Multivariate Analysis of Ecological Data using CANOCO. Bulletin of the Ecological Society of America Vol. 87 (Cambridge University Press, 2003).MATH 
Book 

Google Scholar 
72.Jongman, R. H. G., Ter Braak, C. J. F. & van Tongeren, O. F. R. Data Analysis in Community and Landscape Ecology. Data Analysis in Community and Landscape Ecology (Cambridge University Press, 1995). https://doi.org/10.1017/cbo9780511525575.Book 

Google Scholar 
73.ter Braak, J. F. C. & Šmilauer, P. Canoco Reference Manual and CanoDraw for Windows User’s Guide (Microcomputer Power, 2002).
Google Scholar 
74.R Development Core Team. R: A Language and Environment for Statistical Computing (2020).75.Oksanen, J. et al. vegan: Community Ecology Package. R Packag. version 2.5-7 (2020). More

Insect preparationWe collected 200–300 larvae and pupae of the diamondback moth from cabbage and cauliflower fields in Wuhan, Beijing, and Shenyang in late September from 2008 to 2012. We mixed all collections into one stock colony because of no geographic differentiation in this species from these sites59. We reared all individuals on cabbage leaves spread evenly across five screen cages (35 × 35 × 15 cm) in growth chambers at constant temperature (25 ± 1 °C) with 15-h light:9-h dark photoperiod, and relative humidity set at 60 ± 10%. We moved any new pupae to new screen cages for adult emergence. Emerged adults were fed with 10% honey solution in cotton balls. To collect eggs, we dipped four small pieces of laboratory film (7 × 5 cm) in fresh cabbage juice for 3–5 s and hung the treated film pieces on the top of each screen cage. To further enlarge the population for our experiments, we reared these insects in the artificial diet in plastic boxes at 25 ± 1 °C. We transferred 200 eggs to the surface of 120 g artificial diet (Southland Products Incorporated, USA) in each plastic box (10 × 10 × 9 cm). The hatched larvae dropped to the surface of the artificial diet and fed on it. Once individuals developed into 3rd or 4th instar larvae, pupae or adults, they were exposed to 10 °C for 24 h (to simulate gradually reduced temperatures in late autumn and allow a thermal acclimation) just before they were placed to the low-temperature regimes for overwintering tests. Overall, we obtained >7000 larvae and >8000 pupae for the laboratory experiment, and >8000 larvae, >8000 pupae, >4000 adults for the field experiment. We have compared the life history traits of insects reared on the artificial diet with natural host plants (cabbage leaves), they performed similarly (Peng and Li, unpublished data).Laboratory experiment of winter survivalSite selectionTo identify what factor determines winter survival under different winter thermal conditions, we conducted a laboratory experiment that simulated temperature regimes of 10 selected sites across a latitudinal gradient in China (Fig. 1, Supplementary Table 1) at which this species is known to occur and damage cruciferous crops during the growing season.Temperature treatmentTo simulate the winter temperatures in the 10 geographically distinct sites (Fig. 1a, Supplementary Table 1), we collected daily mean temperatures during winter (November to next April of 1966–2010) at each site from China Meteorological Data Service Centre (http://data.cma.cn/en). Then, we fitted a polynomial model to the temporal changes of winter daily mean temperatures for each site (Fig. 1b). To simplify the logistics of temperature control procedures, we set all temperature regimes in combinations of linear decline, horizontal maintenance and linear increase to mimic the polynomial changes of winter temperatures in the 10 sites, and adjusted temperature every 10 days as needed (Fig. 1c). We controlled the winter temperature changes of the 10 sites with climate chambers (RXZ-280B, Jiangnan Ltd., Ningbo, China) and refrigerators (Royalstar BCD-246GER) according to curves in Fig. 1c.Experimental protocolsWe conducted a winter survival experiment with 10 low-temperature regimes (Fig. 1c). We exposed 6050 larvae and 7150 pupae to 10 low-temperature regimes according to the experiment design (Fig. 1c). Then we sampled 55 larvae and 65 pupae every 10 days from each temperature regime resulting in 11 sampling points. Sampled larvae were placed at 25 °C for 1.5 days to observe the survival based on if their body kept fresh green60 and appendage moved after touching with a brush34. The pupae were placed at 25 °C, RH 70–80% and photoperiod of 16 L:8D for emergence to determine the survival (emergence rate). These samples were not returned to the temperature treatments. Thus, no individual was measured more than once and each sample interval represents an independent observation.Field survival experiments across 12 geographic sitesTo verify the cold survivals from the laboratory simulation and identify the best predictor under natural conditions, we conducted field experiments to explore the winter survival for multiple years at various geographic sites in China (Fig. 1a, Supplementary Table 1). The diamondback moth overwinters either in remaining cabbage plants or in fallen leaves (post-harvest conditions) in regions without standing cabbage crops in the winter. We tested the winter survival of larvae, pupae and adults in the caged cabbage plants or post-harvest conditions in fallen leaves on the soil surface at each site for 3–4 months. We transferred 30 larvae, 30 pupae or 30 adults from our stock rearing to a cabbage plant in the field. Then each plant was covered with a screen cage to avoid disturbance and contain focal individuals (see Supplementary Fig. 2). We set 6–8 cages for larvae, pupae and adults, respectively, in a field in November or early December. After an exposure of 1, 2, 3 and 4 months, we collected 2 cages of larvae, 2 cages of pupae and 2 cages of adults from the field at each sampling point and kept individuals in the laboratory (25 ± 1 °C, RH 65–75%, L:D = 16:8 h) for two days. We checked the survival status of the larvae based on the change in body coloration (i.e. if the larval body kept fresh green colour)60, pupal survival based on whether adults could emerge from the pupae, and adult survival based on if their appendage moved after touching with a brush.To simulate the field microenvironment of post-harvest conditions in winter, we filled half of a glass jar (diameter = 5.5 cm, height = 14 cm) with moist soil. Then, we transferred 30 larvae or 30 pupae to the soil surface, covered the insects with leaves, and then covered the glass jar with a nylon net (see Supplementary Fig. 2). We buried 6–8 jars for larvae and pupae, respectively, and kept the top of the jar at ground surface level at each site in November and early December. Because almost all adults died in few days within the jar, we did not test the adult survival in post-harvest conditions. After an exposure of 1, 2, 3 and 4 months, we took 2 jars of larvae and 2 jars of pupae per sampling period from the field and placed them in the laboratory with 25 ± 1 °C, RH 65–75%, L:D = 16:8 h for 2 days. The survival status of the larvae and pupae was checked with the same procedures as the overwintering tests on caged cabbage plants. Note that as in the standing plant experiment, no individual was tested more than once assuming that each observation is independent at the replicate level.Modelling and predicting winter survivalModel developmentOur goal was to identify key metrics that best predict the winter survival of the diamondback moth across a climatic gradient. To achieve this goal, we took several steps. First, we fit a set of predictive models to the laboratory experiments to identify which metric and model best describes survival under controlled conditions. We focused on three alternative predictors: the lowest daily mean temperature (MinDTmean), mean temperature (DTmean) combined with exposure days, and low-temperature degree-days (LTDD). The MinDTmean model assumes that survival can simply be predicted as a function of the lowest temperature an individual experienced during its exposure time. The DTmean model assumes that survival depends on both the average temperature individuals experience below the cold threshold for survival (11 °C)32 and exposure duration (note that exposure time varied systematically in 10-day increments). Finally, the LTDD model predicts survival depending on coldness below the cold threshold. We calculated LTDD by summing up negative deviations of daily mean temperatures from the cold threshold (11.0 °C) during each exposure period for each simulated geographic site (Fig. 1c). To detect potential relationships, we fit each model using three different functions, i.e. linear, exponential and sigmoid models to describe the survival probability (Supplementary Table 2, Fig. 1). We estimated parameters of models in SigmaStat 3.5 and compared model fit using R2 and AIC values (see detailed models in Supplementary Table 2).Field validation of survival modelsTo validate winter survival models derived from the laboratory (see models in Supplementary Table 2) for complex and variable field conditions, we compared model predictions to observed survivals in field experiments across 12 different geographic sites over 5 years (Fig. 1a, Supplementary Table 1). To make the connection, we first collected daily mean temperatures recorded at the nearest weather stations to our field sites from China Meteorological Data Service Centre. We then calculated MinDTmean, DTmean and exposure days, and LTDD for each site for each treated period and input these values into these laboratory models to predict winter survival. Note, that because the coefficients were calculated from the laboratory experiment, predictions are completely independent of survival observed under field conditions. During the model validation, we excluded the field data of south China, e.g. Guangzhou, Changsha and Wuhan where the warmer temperatures allowed moths to continue their regular life cycle during the whole winter, resulting in unrealistic winter survival. We also excluded replicates in which glass jars were filled with water and destroyed the tested insects. We used linear regression to compare predicted survival with field observations. The validity of each model was evaluated based on the variance explained, slopes of linear regressions and prediction bias (i.e. deviation from unity slope). Finally, we selected the exponential model driven by LTDD as the model to predict the global distribution of winter survival due to its lowest AIC value (Supplementary Table 2) and the least bias (Supplementary Table 3, Fig. 2) among all models.Global prediction of overwintering range shiftTo extrapolate our winter survival predictions to a global scale under present and future climate conditions, we downloaded global historical daily mean temperature data for 50 years (1967–2016) from Berkeley Earth (1° × 1° grid, http://berkeleyearth.org/data/). We added 1, 2, 3, 4, 5 and 6 °C to mean temperatures of 2012–2016, respectively, to represent the different future warming scenarios37. Then, we calculated the annual LTDD in the northern hemisphere with Eq. (1) and in the southern hemisphere with Eq. (2). For xi,j  x0, we excluded the xi,j for the calculation LTDD. For Eq. (1), we started the calculation of LTDD from July 1st (Julian date 182), ended on June 30th of next year (Julian date 181) to cover the whole low-temperature season in the northern hemisphere cross the calendar year. We used LTDD for every year during past conditions to our validated survival model (LTDD-dependent exponential model) and further calculated the expected corresponding yearly winter survival and 5-year mean survival. Since the diamondback moth only feeds on Brassicaceae plants61, we incorporated host availability to refine the pest distributions. We retrieved Brassicaceae occurrence data during 1967–2016 (3,720,971 records) from the Global Biodiversity Information Facility (GBIF) database (www.gbif.org), and excluded unknown and duplicate records; 919,808 records were retained to model the global distribution of host plants. We used a dataset of eight selected bioclimatic variables as described in a previous Brassicaceae biogeographic study62, including isothermality (bio3), temperature seasonality (bio4), min temperature of coldest month (bio6), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation seasonality (bio15), precipitation of warmest quarter (bio18), precipitation of coldest quarter (bio19) from Worldclim dataset63 (http://worldclim.org). We ran the species distribution model using the Maxent algorithm in R package dismo64. Model outputs were presented in grid ranks of host plant presence probability from 0 (unsuitable) to 1 (most suitable). Based on the known distribution of Brassicaceae, we only included grid cells with Brassicaceae presence probability ≥0.3 for our final survival and distribution analysis to ensure the presence of the host plant and mapped them with Arcmap 10.2 (Environmental Systems Research Institute) (see Fig. 3a, e). To show spatial-temporal changes in the geographic distribution of winter survival, we quantified the historical change (expansions or contractions) in the overwintering range based on the total numbers of grids for each year between 1967 and 2016 relative to the baseline area in 1967 (see Fig. 3f) and further calculated average changes of every 5 years (see Fig. 3b). We selected sites in the overwintering marginal belt (with winter survival between 1 and 5%) in the baseline year (1967), calculated the annual LTDD of these sites from 1967 to 2016, and built the linear trend of annual LTDD for years 1967–2016 (see Fig. 3g). We predicted distribution changes for future scenarios (added 1, 2, 3, 4, 5 and 6 °C to the current mean temperatures of 2012–2016) relative to the baseline area of 1967–1971 (see Fig. 3c, d).Meta-analysis linking pesticide resistance to overwintering typeData preparation: literature search and selection criteriaWe performed a comprehensive literature survey to collect data on pesticide resistance of the diamondback moth worldwide. We searched for publications in databases of ISI Web of Science, Scopus and China National Knowledge Infrastructure (CNKI) using keywords “pesticide resistance” in combination with “diamondback moth” or “Plutella xylostella” and expanded references in the selected papers. We reviewed titles, abstracts and in many cases the full articles for relevance and agreement with our inclusion criteria. Studies were included if they (1) monitored the pesticide resistance of field populations, (2) used the leaf dip bioassay method to test pesticide resistance which is the most commonly used method recommended by Insecticide Resistance Action Committee (IRAC, http://www.irac-online.org); (3) provided resistance ratio of field populations. Resistance ratio (abbreviated as RR) is the magnitude of pesticide resistance and is commonly calculated by dividing the median lethal concentration (LC50) of a tested field population by LC50 of the susceptible population (without exposure to pesticide). The LC50 is commonly estimated from a concentration-mortality curve of a given pesticide. The preliminary literature search resulted in 2151 studies out of which 62 matched these criteria. A PRISMA diagram describing details of our literature search is available in Supplementary Fig. 4.Data preparation: data extractionWe extracted data from each selected publication, the names of pesticides, sampling locations and years of field populations, number of tested individuals in a bioassay, resistance ratio of field populations (RR), LC50 of field populations (LC50field) and susceptible populations (LC50susceptible), and 95% confidence intervals (CIs) of LC50field and LC50susceptible. Some studies generated results from multiple types of pesticides with the same field population, each of which was considered as a different entry. Finally, we gathered 1806 entries for pesticide resistance of field populations of the diamondback moth.Data preparation: calculation of the weighted effect sizeWe conducted a meta-analysis to test if pesticide resistance levels vary across different types of overwintering sites. To account for differences in sample sizes and variances in resistance ratios across studies, we calculated the corrected (weighted) resistance ratio for each study following the method in Hedges et al.65. We calculated the logarithm of resistance ratio (logRR) to present the effect size for each entry and further calculated the weighted effect size (wlogRR) by$${{{{{rm{wlogRR}}}}}}={{{{{rm{logRR}}}}}}times w$$
(3)
where w is the weighting factor of each entry, with w = 1/sqrt(VlogRR)66. To consider the contribution from both field and susceptible population, the pooled variance VlogRR was calculated as follows65:$${{{{{{rm{V}}}}}}}_{{{{{{rm{logRR}}}}}}}=frac{{{{{{{{rm{SE}}}}}}}_{{{{{{rm{field}}}}}}}}^{2}}{{n}_{{{{{{rm{field}}}}}}}times {{{{{{{rm{LC50}}}}}}}_{{{{{{rm{field}}}}}}}}^{2}}+frac{{{{{{{{rm{SE}}}}}}}_{{{{{{rm{susceptible}}}}}}}}^{2}}{{n}_{{{{{{rm{susceptible}}}}}}}times {{{{{{{rm{LC50}}}}}}}_{{{{{{rm{susceptible}}}}}}}}^{2}}$$
(4)
where LC50field and LC50susceptible, SEfield and SEsusceptible, nfield and nsusceptible, are LC50, the standard error of LC50 and sample size for field population and susceptible population, respectively. SEfield and SEsusceptible can be calculated from their own confidence intervals (95% CI)67:$${{{{{rm{SE}}}}}}=frac{{{{{{{rm{CI}}}}}}}_{{{{{{rm{upper}}}}}}{{{{{rm{limit}}}}}}}-{{{{{{rm{CI}}}}}}}_{{{{{{rm{lower}}}}}}{{{{{rm{limit}}}}}}}}{2times 1.96}$$
(5)
where CIupper limit is the upper limit and CIlower limit is the lower limit of the 95% CI for LC50.We used the prognostic method68 to estimate VlogRR for entries that miss either 95% CI or LC50 based on the average VlogRR of the other complete entries.Data preparation: potential moderator variablesSeveral factors could influence pesticide resistance besides overwintering temperatures. The effective temperature degree-days (ETDD) may change the annual number of generations, the intensity of pesticide application, and thus the selection stress47, e.g. between 7.4 and 33 °C for the diamondback moth32. In addition, the variety of pesticides used in a study may also affect the resistance levels through their mode of actions (the lethal mechanism) and cross-resistance69,70. To account for these potentially confounding factors, we collected the mode of action for each variety of pesticides from IRAC, and calculated LTDD, ETDD and overwintering type for each of the 1806 original records. We collected data for daily mean temperatures for each site from Berkeley Earth. For each location, we calculated the mean annual LTDD, ETDD and winter survival average across the 5 years before the sample. We split the sampling sites into three types based on predicted winter survival of the diamondback moth: (1) the permanent (overwintering) sites: locations with the mean winter survivals ≥5%, (2) marginal sites: locations with the mean winter survivals 1–5%, (3) transient (non-overwintering) sites: locations with the mean winter survivals More

1.Elmqvist, T. et al. Benefits of restoring ecosystem services in urban areas. Curr. Opin. Environ. Sustain. 14, 101–108 (2015).Article 

Google Scholar 
2.Jones, H. P. et al. Restoration and repair of Earth’s damaged ecosystems. Proc. R. Soc. B Biol. Sci. https://doi.org/10.1098/rspb.2017.2577 (2018).Article 

Google Scholar 
3.Rey Benayas, J. M., Newton, A. C., Diaz, A. & Bullock, J. M. Enhancement of biodiversity and ecosystem services by ecological restoration: A meta-analysis. Science 325, 1121–1124. https://doi.org/10.1126/science.1172460 (2009).ADS 
CAS 
Article 
PubMed 

Google Scholar 
4.Brudvig, L. A. et al. Interpreting variation to advance predictive restoration science. J. Appl. Ecol. 54, 1018–1027. https://doi.org/10.1111/1365-2664.12938 (2017).Article 

Google Scholar 
5.Suding, K. N. Toward an era of restoration in ecology: Successes, failures, and opportunities ahead. Annu. Rev. Ecol. Evol. Syst. 42, 465–487. https://doi.org/10.1146/annurev-ecolsys-102710-145115 (2011).Article 

Google Scholar 
6.Reynolds, H. L., Packer, A., Bever, J. D. & Clay, K. Grassroots ecology: Plant-microbe-soil interactions as drivers of plant community structure and dynamics. Ecology 84, 2281–2291 (2003).Article 

Google Scholar 
7.Van Der Heijden, M. G. A., Bardgett, R. D. & Van Straalen, N. M. The unseen majority: Soil microbes as drivers of plant diversity and productivity in terrestrial ecosystems. Ecol. Lett. 11, 296–310 (2008).Article 

Google Scholar 
8.Hoeksema, J. D. et al. A meta-analysis of context-dependency in plant response to inoculation with mycorrhizal fungi. Ecol. Lett. 13, 394–407. https://doi.org/10.1111/j.1461-0248.2009.01430.x (2010).Article 
PubMed 

Google Scholar 
9.Schultz, P. A. et al. Evidence of a mycorrhizal mechanism for the adaptation of Andropogon gerardii (Poaceae) to high- and low-nutrient prairies. Am. J. Bot. 88, 1650–1656. https://doi.org/10.2307/3558410 (2001).CAS 
Article 
PubMed 

Google Scholar 
10.Koske, R. E., & Gemma, J. N. Mycorrhizae and succession in plantings of beachgrass in sand dunes. Am. J. Bot. 84(1), 118–130 (1997).Article 

Google Scholar 
11.Smith, M. E., Facelli, J. M. & Cavagnaro, T. R. Interactions between soil properties, soil microbes and plants in remnant-grassland and old-field areas: a reciprocal transplant approach. Plant Soil 433, 127–145. https://doi.org/10.1007/s11104-018-3823-2 (2018).CAS 
Article 

Google Scholar 
12.Tipton, A. G., Middleton, E. L., Spollen, W. G. & Galen, C. Anthropogenic and soil environmental drivers of arbuscular mycorrhizal community composition differ between grassland ecosystems. Botany 97, 85–99. https://doi.org/10.1139/cjb-2018-0072 (2019).Article 

Google Scholar 
13.Hamman, S. T. & Hawkes, C. V. Biogeochemical and microbial legacies of non-native grasses can affect restoration success. Restor. Ecol. 21, 58–66. https://doi.org/10.1111/j.1526-100X.2011.00856.x (2013).Article 

Google Scholar 
14.Emery, S. M. & Rudgers, J. A. Beach restoration efforts influenced by plant variety, soil inoculum, and site effects. J. Coast. Res. 27, 636. https://doi.org/10.2112/jcoastres-d-10-00120.1 (2010).Article 

Google Scholar 
15.Sylvia, D. M., Jarstfer, A. G. & Vosátka, M. Comparisons of vesicular-arbuscular mycorrhizal species and inocula formulations in a commercial nursery and on diverse Florida beaches. Biol. Fertil. Soils 16, 139–144. https://doi.org/10.1007/BF00369416 (1993).Article 

Google Scholar 
16.Sylvia, D. M. & Will, M. E. Establishment of vesicular-arbuscular mycorrhizal fungi and other microorganisms on a beach replenishment site in Florida. Appl. Environ. Microbiol. 54, 348–352 (1988).ADS 
CAS 
Article 

Google Scholar 
17.Wubs, E. R. J., van der Putten, W. H., Bosch, M. & Bezemer, T. M. Soil inoculation steers restoration of terrestrial ecosystems. Nat. Plants 2, 16107. https://doi.org/10.1038/nplants.2016.107 (2016).Article 
PubMed 

Google Scholar 
18.Bothe, H., Turnau, K. & Regvar, M. The potential role of arbuscular mycorrhizal fungi in protecting endangered plants and habitats. Mycorrhiza 20, 445–457. https://doi.org/10.5586/asbp.2008.019 (2010).Article 
PubMed 

Google Scholar 
19.Middleton, E. L. & Bever, J. D. Inoculation with a native soil community advances succession in a grassland restoration. Restor. Ecol. 20, 218–226. https://doi.org/10.1111/j.1526-100X.2010.00752.x (2012).Article 

Google Scholar 
20.Crawford, K. M., Busch, M. H., Locke, H. & Luecke, N. C. Native soil microbial amendments generate trade-offs in plant productivity, diversity, and soil stability in coastal dune restorations. Restor. Ecol. https://doi.org/10.1111/rec.13073 (2019).Article 

Google Scholar 
21.Eom, A. H., Hartnett, D. C. & Wilson, G. W. T. Host plant species effects on arbuscular mycorrhizal fungal communities in tallgrass prairie. Oecologia 122, 435–444. https://doi.org/10.1007/s004420050050 (2000).ADS 
Article 
PubMed 

Google Scholar 
22.Brundrett, M. C. & Tedersoo, L. Evolutionary history of mycorrhizal symbioses and global host plant diversity. New Phytol. 220, 1108–1115 (2018).Article 

Google Scholar 
23.Bever, J. D., Mangan, S. A. & Alexander, H. M. Maintenance of plant species diversity by pathogens. Annu. Rev. Ecol. Evol. Syst. 46, 305–325. https://doi.org/10.1146/annurev-ecolsys-112414-054306 (2015).Article 

Google Scholar 
24.Crawford, K. M. et al. When and where plant-soil feedback may promote plant coexistence: a meta-analysis. Ecol. Lett. 22, 13278. https://doi.org/10.1111/ele.13278 (2019).Article 

Google Scholar 
25.Mills, K. E. & Bever, J. D. Maintenance of diversity within plant communities: Soil pathogens as agents of negative feedback. Ecology 79, 1595–1601. https://doi.org/10.1890/0012-9658(1998)079[1595:MODWPC]2.0.CO;2 (1998).Article 

Google Scholar 
26.Koziol, L. et al. The plant microbiome and native plant restoration: The example of native mycorrhizal fungi. Bioscience 68, 996–1006 (2018).Article 

Google Scholar 
27.Maltz, M. R. & Treseder, K. K. Sources of inocula influence mycorrhizal colonization of plants in restoration projects: A meta-analysis. Restor. Ecol. 23, 625–634. https://doi.org/10.1111/rec.12231 (2015).Article 

Google Scholar 
28.Koziol, L. & Bever, J. D. AMF, phylogeny, and succession: Specificity of response to mycorrhizal fungi increases for late-successional plants. Ecosphere https://doi.org/10.1002/ecs2.1555 (2016).Article 

Google Scholar 
29.Middleton, E. L. et al. Locally adapted arbuscular mycorrhizal fungi improve vigor and resistance to herbivory of native prairie plant species. Ecosphere 6, 276. https://doi.org/10.1890/ES15-00152.1 (2015).Article 

Google Scholar 
30.Solís-Domínguez, F. A., Valentín-Vargas, A., Chorover, J. & Maier, R. M. Effect of arbuscular mycorrhizal fungi on plant biomass and the rhizosphere microbial community structure of mesquite grown in acidic lead/zinc mine tailings. Sci. Total Environ. 409, 1009–1016. https://doi.org/10.1016/j.scitotenv.2010.11.020 (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
31.Vogelsang, K. M., Reynolds, H. L. & Bever, J. D. Mycorrhizal fungal identity and richness determine the diversity and productivity of a tallgrass prairie system. New Phytol. 172, 554–562. https://doi.org/10.1111/j.1469-8137.2006.01854.x (2006).Article 
PubMed 

Google Scholar 
32.Larimer, A. L., Bever, J. D. & Clay, K. Consequences of simultaneous interactions of fungal endophytes and arbuscular mycorrhizal fungi with a shared host grass. Oikos 121, 2090–2096. https://doi.org/10.1111/j.1600-0706.2012.20153.x (2012).Article 

Google Scholar 
33.Sikes, B. A., Cottenie, K. & Klironomos, J. N. Plant and fungal identity determines pathogen protection of plant roots by arbuscular mycorrhizas. J. Ecol. 97, 1274–1280. https://doi.org/10.1111/j.1365-2745.2009.01557.x (2009).Article 

Google Scholar 
34.Defeo, O. et al. Threats to sandy beach ecosystems: A review. Estuar. Coast. Shelf Sci. 81, 1–12 (2009).ADS 
Article 

Google Scholar 
35.Feagin, R. A. et al. Going with the flow or against the grain? The promise of vegetation for protecting beaches, dunes, and barrier islands from erosion. Front. Ecol. Environ. 13, 203–210 (2015).Article 

Google Scholar 
36.Feagin, R. A. et al. The role of beach and sand dune vegetation in mediating wave run up erosion. Estuar Coast Shelf Sci. 219, 97–106. https://doi.org/10.1016/j.ecss.2019.01.018 (2019).ADS 
Article 

Google Scholar 
37.Sigren, J. M., Figlus, J. & Armitage, A. R. Coastal sand dunes and dune vegetation: Restoration, erosion, and storm protection. Shore Beach 82, 5–12 (2014).
Google Scholar 
38.Sigren, J. M. et al. The effects of coastal dune volume and vegetation on storm-induced property damage: Analysis from Hurricane Ike. J. Coast Res. 341, 164–173. https://doi.org/10.2112/jcoastres-d-16-00169.1 (2018).Article 

Google Scholar 
39.Silva, R. et al. Response of vegetated dune-beach systems to storm conditions. Coast. Eng. 109, 53–62. https://doi.org/10.1016/j.coastaleng.2015.12.007 (2016).Article 

Google Scholar 
40.Lane, C., Wright, S. J., Roncal, J. & Maschinski, J. Characterizing environmental gradients and their influence on vegetation zonation in a subtropical coastal sand dune system. J. Coast. Res. 4, 213–224. https://doi.org/10.2112/07-0853.1 (2008).CAS 
Article 

Google Scholar 
41.Miller, T. E., Gornish, E. S. & Buckley, H. L. Climate and coastal dune vegetation: Disturbance, recovery, and succession. Plant Ecol. 206, 97–104. https://doi.org/10.1007/s11258-009-9626-z (2010).Article 

Google Scholar 
42.Hewitt, E. J. & Eden, A. Sand and water culture methods used in the study of plant nutrition. Analyst 78, 329–330 (1953).
Google Scholar 
43.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria). https://www.R-project.org/ (2020).
Google Scholar 
44.Farrer, E. C. & Goldberg, D. E. Litter drives ecosystem and plant community changes in cattail invasion. Ecol. Appl. 19, 398–412. https://doi.org/10.1890/08-0485.1 (2009).Article 
PubMed 

Google Scholar 
45.Bauer, J. T., Koziol, L. & Bever, J. D. Local adaptation of mycorrhizae communities changes plant community composition and increases aboveground productivity. Oecologia https://doi.org/10.1007/s00442-020-04598-9 (2020).Article 
PubMed 

Google Scholar 
46.Ohsowski, B. M., Klironomos, J. N., Dunfield, K. E. & Hart, M. M. The potential of soil amendments for restoring severely disturbed grasslands. Appl. Soil. Ecol. 60, 77–83. https://doi.org/10.1016/j.apsoil.2012.02.006 (2012).Article 

Google Scholar 
47.Koziol, L. & Bever, J. D. The missing link in grassland restoration: arbuscular mycorrhizal fungi inoculation increases plant diversity and accelerates succession. J. Appl. Ecol. 54, 1301–1309. https://doi.org/10.1111/1365-2664.12843 (2017).Article 

Google Scholar 
48.Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193. https://doi.org/10.1016/0169-5347(94)90088-4 (1994).CAS 
Article 
PubMed 

Google Scholar 
49.Heneghan, L. et al. Integrating soil ecological knowledge into restoration management. Restor. Ecol. 16, 608–617. https://doi.org/10.1111/j.1526-100X.2008.00477.x (2008).Article 

Google Scholar 
50.Wubs, E. R. J. et al. Single introductions of soil biota and plants generate long-term legacies in soil and plant community assembly. Ecol. Lett. 22, 1145–1151 (2019).Article 

Google Scholar 
51.Hestrin, R., Hammer, E. C., Mueller, C. W. & Lehmann, J. Synergies between mycorrhizal fungi and soil microbial communities increase plant nitrogen acquisition. Commun. Biol. 2, 233–242. https://doi.org/10.1038/s42003-019-0481-8 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More