Several studies have documented a global pattern of phenological advancement that is consistent with ongoing climate change. However, the magnitude of these phenological shifts is highly variable across taxa and locations. This variability of phenological responses has been difficult to explain mechanistically. To examine how the evolution of multi-trait cueing strategies could produce variable responses to climate change, we constructed a model in which organisms evolve strategies that integrate multiple environmental cues to inform anticipatory phenological decisions. We simulated the evolution of phenological cueing strategies in multiple environments, using historic climate data from 78 locations in North America and Hawaii to capture features of climatic correlation structures in the real world. Organisms in our model evolved diverse strategies that were spatially autocorrelated across locations on a continental scale, showing that similar strategies tend to evolve in similar climates. Within locations, organisms often evolved a wide range of strategies that showed similar response phenotypes and fitness outcomes under historical conditions. However, these strategies responded differently to novel climatic conditions, with variable fitness consequences. Our model shows how the evolution of phenological cueing strategies can explain observed variation in phenological shifts and unexpected responses to climate change.

Please see the methods description in our associated publication. 

This following code guide will be published in the Appendix of the associated paper:

Code guide

The simulation code for this project is designed to work with the following file structure, which is necessary for the code to run. Files are italicized while folders are not. Exact file names are contained in quotes. A ZIP file with the scripts and data files in the appropriate file structure can be found on Dryad. (Update 2/25/2021: the file "Phenology_simulations_directory.zip" contains the complete data; "phenology_simulations_dir_slim.zip" is a smaller version that only has one location in data-years (ABERDEEN.MS) for previewing purposes).

• Project directory
  • Climate data
    • “location summary with distance to coast.csv”
  • data-years
    • all data files to be used
    • “location summary with photoperiod.Rdata”
  • data-years-lisvski
    • all data files to be used
  • fitcurve
    • “skewgauss.R”
  • parameters
    • Any parameter files to run. The following are the parameter files used to generate the results presented here:
    • “cu-ccshift-mut005-parameters.R”
    • “cu-cctemp-mut005-parameters.R”
    • “cu-cuphotoperiod-mut005-parameters.R”
    • “sensitivity-win50-mut005-parameters.R”
  • results
  • scripts
    • “master-script.R”
    • “rate-setup.R”
    • “simulation.R”
    • “windows-subs.R”
  • yearinds

Code is initiated by running master-script.R. The “set_wrkdir()” function at the beginning of this file should be updated to point to the home folder you are using, and “prefix” and “suffix” should be updated to match the name of the parameter file to run (by default, set “prefix” to everything before the “.R” of the parameter file name, and set “suffix” to be an empty string). The “sum.merge.only” and “test.only” parameters are to analyze existing runs and test the code using the first two locations, respectively. These should be set to FALSE for a complete run. Outside of the parameter files, these are the only lines of code that may require modification to implement our model.

After reading in parameter (below) and function files, “master-script.R” calls “simulation.R” script for every location. “simulation.R” reads in the imputed climate data files, and calculates derived climate measures (cumulative cues, fitness, etc). It then calls “rate-setups.R” to initialize mutation rates and starting trait ranges. “simulation.R” then uses the function year_var_analyze() to calculate climate metrics, which are saved in “…-stats.Rdata” in the data-years folder. The script then saves the current scripts (into the results folder), parameters (into the results folder), and the random order of years to be used (into the yearinds folder). These are useful for tracking changes to the model as well as making comparisons between model runs. Climate change scenarios are then generated, creating a changed version of the locations climate data, which is saved in the results folder. “simulation.R” then uses parallel processing (using between 1 and 7 clusters depending on the computer’s number of cores) to carry out the actual simulation, which is run through the run_sim() function. This and most other functions are defined in the “windows-subs.R” file. Results of each simulation run (along with intermediate objects and functions) are saved in the corresponding results folder in a file labeled “dat.Rdata”. The realized relative cue use is then calculated and saved in “finalpop-yrtest-acteff-dat.Rdata” (see Methods for an explanation of this metric). As a reminder, this is calculated for each individual of the final population exposed to each year of the historic climate. The climate change scenario is then evaluated using the climate_master() function. This function takes the final population of the simulation and evaluates its emergence and performance on each year of the changed climate, saving the results in “temp.Rdata” in the climate folder within the corresponding results folder. Finally, fitness parameters used for this location are saved in “…_fitparms.Rdata” in the corresponding results folder.

After each simulations for each location have been run, “master-script.R” uses trait_eff_summarize_small()to aggregate the realized relative cue use information to the individual (“acteff-agg-byindiv.Rdata”) and the simulation level (“acteff-agg-bysim.Rdata”) level in the corresponding results folder, and aggregates all simulation information to the simulation level in “all-finalpop-yrtest-simlevel-….Rdata” in the appropriate results folder. After all simulations and aggregations have been carried out, “master-script.R” calls dat_sum_merge()  to summarize climate metrics, and then merge these metrics with location data, simulation results, and realized relative cue use in historic and changed climate regimes, which is saved in into a single data frame with simulation-level summaries which is saved in “climateVsPops-….Rdata”. This file contains the majority of the information presented in this paper. Below is a description of each column’s contents.

In these descriptions and in the parameter guide that follows, “simulation” refers to one of several identically parameterized instances of the model in a location; these are best imagined as alternate realities. The results we present had 30 simulations for each location for each set of parameters. “Runs” refers to different executions of the code (typically an execution of the code applies the parameters to every location, producing multiple simulations for each of 78 locations). Presumably these runs are carried out with the purpose of comparing consequences of differences in parameters; thus each run is likely to have a different parameter file.

The results file has a row for each simulation. Each simulation in a location had a different random sequence of years (1000 by default) drawn from the same original collection of years of recorded daily climate data with missing values imputed (see Environmental Data in Methods and Supplements). Unless stated otherwise, climate metrics were calculated based on the sequence of years actually experienced, producing slightly different values between simulations.

Results file column contents

• name: Location ID (see supplements table 1)
• sim: simulation number
• eff.day: simulation day effect, averaged across all individuals and all years using the acomp function in the compositions package (see realized relative cue use in Methods)
• eff.cutemp: as day, for cumulative temperature
• eff.cuprecip: as day, for cumulative precipitation
• emerge: emergence day of final population averaged across all years and all individuals
• emerge.cc: emergence day of final population under climate change scenario averaged across all years and all individuals
• geofit: raw units of “fitness” obtained by final generation (used to determine relative). Averaged across all individuals in all years
• geofit.cc: as geofit, but using climate change climate data
• lis.tpred: climate metric that represents the predictability of temperature⁵. This was calculated for the original sequence of years (for each location, the years were ordered chronologically and treated as sequential). As such, the same value was used for all simulations with the same location.
• lis.ppred: as lis.tpred, but for precipitation
• lis.mpred: as lis.tpred, but for moisture
• lis.tseason: climate metric that represents the seasonality of temperature⁵
• lis.pseason: as lis.tseason, but for precipitation
• lis.mseason: as listseason, but for moisture
• pau.tvar: climate metric that represents variability in temperature⁶
• pau.mvar: as pau.tvar, but for moisture
• pau.pvar: as pau.tvar, but for precipitation
• pau.tseason: climate metric associated with each simulation. This metric represents seasonality of temperature⁶
• pau.mseason: as pau.tseason, but for moisture
• pau.pseason: as pau.tseason, but for precipitation
• b.day.mn: day trait of final population, averaged across all individuals
• b.cutemp.mn: as b.day.mn, but cumulative temperature
• b.cuprecip.mn: as b.day.mn, but cumulative precipitation
• med.cv: coefficient of variation in the “median day”. For each year of the climate data, we determined the first day the cumulative temperature reached or exceeded half the maximum cumulative temperature of the mean year for that location. Coefficient of variation was calculated from these median days
• imput.temp.mn: fraction of the daily temperature values that originated form imputation
• imput.prec.mn: as imput.temp.mn, but for precipitation
• med.var: as med.cv, but variance instead of coefficieint of variation.
• within.mn: The “within-year unpredictability” was calculated for each year by comparing the daily temperature to the daily temperatures of the mean year. Each year was allowed to advance or retreat by discrete days, and the temperatures could increase or decrease by a constant; using sum of squared errors, we found the daily shift and temperature change constant to make the current year best match the average year. The “within-year unpredictability” was then the mean of the square of the lag 1 difference (“diff”) of the residuals of this fit. This represents day-to-day inconsistencies that weren’t present in the average year. This metric was calculated for each year of the climate, and then averaged across all years.
• within.nodiff.mn: as within.mn, but using the mean of the square of differences (skipping the “diff” step). This represents divergences between the current year and the average year.
• Var.dailyfit.mn: For each year, we calculated the variance in the relative fitness available on each day. This was averaged across all years.
• lr.var: The “left-right shift” for a year was the number of days advancement needed to best fit the current year’s daily temperatures to the average year of that location. For each location, the lr.var was the variance of these measures across all years for that location.
• Updown.var: as lr.var, but using the temperature shift needed to best fit the current year’s daily temperatures to the average year for that location.
• Temp.mn: mean of all yearly mean temps for the simulation (Note: this was calculated with the temperature values that had been shifted so that the minimum temperature in each location was 0)
• Temp.orig.mn: mean of all yearly mean temps for the simulation using the original, unshifted temperatures
• Temp.predtemp.mn: mean correlation of daily temperature on day n to sum of temperatures on days (n+1):(n+duration). As with all other X.predY.mn, this was a metric for how well X acted as a cue to predict the entire span of Y experienced by an individual that chose to emerge on day n.
• Temp.predfit.mn: as temp.predtemp.mn, but how well temperature predicted sum of fitness
• Temp.predmoist.mn: as temp.predtemp.mn, but how well temperature predicted sum of moisture
• Moist.mn: mean moisture across all years
• Moist.predmoist.mn: as temp.predtemp.mn, but how well moisture predicted sum of moisture
• Moist.predfit.mn: as temp.predtemp.mn, but how well moisture predicted sum of fitness Moist.predtemp.mn: as temp.predtemp.mn, but how well moisture predicted sum of temperature
• Moist.cv.mn: average across all years of the within-year coefficient of variation of daily moisture
• Precip.mn: average across all years of within-year mean precipitation
• Precip.var.btwn: between-year variance in the yearly mean precipitations
• Precip.var.mn: average across all years of within-year variance in precipitation
• Precip.cv.mn: as precip.var.mn, but using coefficient of variation
• Temp.varbtwn: between-year variance in the yearly mean temperatures
• Temp.cvbtwn: between-year coefficient of variation in the yearly mean temperatures
• Precip.cvbtwn: between-year coefficient of variation in the yearly mean precipitation
• Moist.varbtwn: between-year variance in the yearly mean moisture
• Moist.cvbtwn: between-year coefficient of variation in the yearly mean moisture
• Temp.lag1temp.mn: mean lag 1 autocorrelation of daily temperature. As with all X.lag1Y.mn metrics, this was a simple metric for how well X acted as cue to predict Y
• Temp.lag1moist.mn: mean correlation between temperature of day n and moisture of day n+1
• Precip.lag1precip.mn: mean lag 1 autocorrelation of precipitation
• Moist.lag1temp.mn: mean correlation between moisture on day n and temperature on day n+1
• Moist.lag1moist.mn: mean lag1 autocorrelation of moisture
• Moist.mn.var: same as moist.varbtwn
• Temp.cv.mn: mean across years of the within-year coefficient of variation of daily temperatures
• Temp.varwin.mn: as temp.cv.mean, but variance instead of coefficient of variation
• Moist.varwin.mn: as temp.varwin.mn, but moisture instead of temperature
• Lat: latitude of location
• Lon: longitude of location
• Elev: elevation of location
• Coast.dist: distance from location to nearest coast
• Count.good: number of years of climate data available.

Parameter file

The specifics of any simulation run are controlled by the parameters defined in the parameter file. Below are descriptions of each parameter

• runs.types: the name of each location file to use, in vector of characters
• traits: used for sensitivity analysis. Changing from “day”, “cutemp”, and “cuprecip” may cause the current code to break
• num.sims (default: 30): the number of simulations to run for each location. Each simulation has identical parameterizations, but the randomly generated starting genotypes, randomly generated sequence of years, and randomly generated mutations of each generation will differ.
• N (default: 500): number of individuals per simulation
• num.years (default 1000): number of years to simulate.
• duration (default: 10): number of days organism gathers fitness
• lag (default: 1): number of days after organism decides to emerge before it begins gathering fitness
• mut.dist (default: .005): fraction of trait range to use as standard deviation of mutation rate
• fattail (default: FALSE): indicator to determine whether to use a normal distribution or Cauchy distribution to generate mutation values. Early testing showed that the two distributions produced indistinguishable dynamics in the model, and subsequent simulations used the normal distribution. Setting this to TRUE may cause the current code to break.
• base.temp (default: 0): threshold below which temperature isn’t added to cumulative temperature
• decay (default: .8): decay parameter alpha for use when calculating moisture (see Methods)
• fit.shape (default: “skewgauss”): defines which fitness shape script to use. The “skewgauss.R” is the only fitness shape script designed to work with the current code.
• Other.name (default: “moist”): determines what daily environmental measure combines with temperature to define fitness.
• shape.temp (default: 15): initial guess at shape parameter for the temperature fitness curve. The final shape parameter is determined by minimizing sum of squared errors.
• shape.other (default: 15): as shape.temp, but for the other fitness-determining environmental measure (by default, moisture)
• best.temp.quant (default: .9): Used to fit fitness curves. At  this quantile of temperatures for each location, the skew normal has its maximum value (see Methods)
• best.other.quant (default: 0.9): as best.temp.quant, but for the other measure (by default, moisture)
• min.quant (default: 0.1): Used to fit fitness curves. At this quantile of temperature and moisture, skew normal function has a value that is (ratio.min) of the peak. (see Methods)
• ratio.min (default: 0.1): determines ratio of values at the min.quant value and best.*.quant value for each skew normal curve.
• plot.extra (default: FALSE): If true, plot emergence of a subset of the generations of each simulation
• plot.pheno (default: FALSE): If true, plot phenotype of a subset of the generations of each simulation
• burnin (default: 100): Used for emergence and phenotype plots, defines how many of the initial generations to skip in the plots
• year.label (default: “A”): When the simulation script is generating a random sequence of years for all simulations in a given location, it checks to see if there exists a file storing random sequences for this location with the same num.sims, num.years, year.set (below) and year.label; if it finds such a file, it uses those randomly generated year sequences. That is, no two simulations within a run will have the same sequence of years, but using the same year.label, num.sims, num.years, and year.set values when running two different parameter files will re-use the same num.sim (default: 30) randomly generated sequences of years for each of the separate parameter file runs. This allows for fair comparisons between different model structures or parameterizations (since ABERDEEN.MS run 1 of each parameterization will have the same sequence of years). Changing the year.label between parameter files causes each to have a separate sequence of years.
• year.set (default: “all.years”): as year.label – shows up in a different part of the naming scheme but plays the same role.
• temp.increase (default: 0 for “shift” climate change, 3 for “warming” climate change): how much should each daily temperature be increased in the climate change scenario?
• precip.shift (default: 5 for “shift” climate change, 0 for “warming climate change”): how many days should precipitation be advanced in the climate change scenario?
• Lockstep (default: TRUE): should temperature advance with precipitation in the climate change scenario?
• tempsd.increase (default: 1.2): for use in deprecated climate change scenario
• phensd.increase (default: 1.2): for use in deprecated climate change scenario
• intra.increase (default: 1.2): for use in deprecated climate change scenario
• scen.temp (default: TRUE): with change in climate change functions, this should always be TRUE
• scen.tempsd (default: FALSE): with change in climate change functions, this should always be FALSE
• scen.phensd (default: FALSE): with change in climate change functions, this should always be FALSE
• scen.intra (default: FALSE): with change in climate change functions, this should always be FALSE
• save.small (default: TRUE): If FALSE, saves information on every individual of every generation in each simulation. If TRUE, only saves information on every individual of the first and final generations.
• photo.flag (default: FALSE): if TRUE, use cumulative photoperiod instead of day of year for the “day” cue.
• do.lisovski (default: FALSE): if TRUE, calculates metrics based on Lisovski et al.⁵ , and saves them in separate data files in data-years-lisovski. As this is a lengthy calculation, it’s best to re-use existing calculations of this (and the data files of these calculations can be found in X).
• stats.lisovski: if TRUE, uses the metrics associated with do.lisovski. This should be TRUE unless you (a) are missing the lisovski metric data files, (b) don’t want to wait the additional hours for the scripts to calculate them, and (c) don’t want to use those metrics.