DOI: 10.25338/b8tg95
title: Data from: Evolved phenological cueing strategies show variable responses to climate change

description: 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.
description: Please see the methods description in our associated
 publication. 
description: This following code guide
 will be published in the Appendix of the associated paper: <i>Code
 guide</i> 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
 
 <i>“location summary with distance to
 coast.csv”</i> 
 data-years 
 <i>all data files to be used</i>
 <i>“location summary with
 photoperiod.Rdata”</i> 
 data-years-lisvski <i>all data
 files to be used</i> 
 fitcurve 
 <i>“skewgauss.R”</i> 
 parameters Any parameter
 files to run. The following are the parameter files used to generate the
 results presented here:
 <i>“cu-ccshift-mut005-parameters.R”</i>
 <i>“cu-cctemp-mut005-parameters.R”</i>
 <i>“cu-cuphotoperiod-mut005-parameters.R”</i> <i>“sensitivity-win50-mut005-parameters.R”</i> results scripts <i>“master-script.R”</i> <i>“rate-setup.R”</i> <i>“simulation.R”</i> <i>“windows-subs.R”</i> yearinds Code is initiated by running <i>master-script.R</i>. 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, <i>“master-script.R”</i> calls “<i>simulation.R”</i> script for every location. “<i>simulation.R” </i>reads in the imputed climate data files, and calculates derived climate measures (cumulative cues, fitness, etc). It then calls <i>“rate-setups.R”</i> to initialize mutation rates and starting trait ranges. <i>“simulation.R” </i>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. <i>“simulation.R” </i>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 “<i>windows-subs.R”</i> file. Results of each simulation run (along with intermediate objects and functions) are saved in the corresponding results folder in a file labeled <i>“dat.Rdata”</i>. The realized relative cue use is then calculated and saved in <i>“finalpop-yrtest-acteff-dat.Rdata”</i> (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 “<i>temp.Rdata</i>” in the climate folder within the corresponding results folder. Finally, fitness parameters used for this location are saved in “<i>…_fitparms.Rdata</i>” in the corresponding results folder. After each simulations for each location have been run, “<i>master-script.R” </i>uses trait_eff_summarize_small()to aggregate the realized relative cue use information to the individual (“<i>acteff-agg-byindiv.Rdata”</i>) and the simulation level <i>(“acteff-agg-bysim.Rdata”</i>) level in the corresponding results folder, and aggregates all simulation information to the simulation level in <i>“all-finalpop-yrtest-simlevel-….Rdata”</i> in the appropriate results folder. After all simulations and aggregations have been carried out, “<i>master-script.R”</i> calls dat_sum_merge()<i> </i> 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 <i>“climateVsPops-….Rdata”</i>. 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. <i>Results file column contents</i> name: Location ID (see supplements table 1) sim: simulation number eff.day: simulation day effect, averaged across all individuals and all years using the <i>acomp</i> function in the <i>compositions</i> 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<sup>5</sup>. 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<sup>5</sup> lis.pseason: as lis.tseason, but for precipitation lis.mseason: as listseason, but for moisture pau.tvar: climate metric that represents variability in temperature<sup>6</sup> 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<sup>6</sup> 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. <i>Parameter file</i> 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 “<i>skewgauss.R</i>” 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.<sup>5</sup><sup> </sup>, 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 <b>X</b>). 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.