library(tidyverse)
library(lubridate)

library(lme4)
library(lmerTest)
library(MuMIn)
library(glmmTMB)
library(emmeans)
library(broom)
library(broom.mixed)
library(vegan)

library(RColorBrewer)
library(viridis)
library(gridExtra)
library(patchwork)
library(GGally)
library(ggnewscale)

library(knitr)
knitr::opts_chunk$set(comment="", cache=T, warning = F, message = F, 
                      fig.path = "figures/", 
                      dev="svglite", dev.args=list(fix_text_size=FALSE))#dev="cairo_pdf")
options(digits=4) # for kables

load("data/maxfield_data.rda")

options(contrasts = c("contr.sum", "contr.poly")) #needed for Type III SS

anova_table <- function(df, mod) { #for lmer models and lmerTest::anova
  df %>% select(trait, term, p.value) %>% 
    pivot_wider(names_from=term, values_from=p.value) %>% 
    mutate(R2m = map_dbl(mod, ~ MuMIn::r.squaredGLMM(.x)[[1,"R2m"]]))
}
Anova_table <- function(df) { #for glmmTMB models and car::Anova
  df %>% select(Trait=trait, term, p.value) %>% 
    pivot_wider(names_from=term, values_from=p.value,names_repair="unique")
}
select <- dplyr::select
lptraits <- c(leaftraits, phystraits)

Metadata

alldata <- list("treatments"=treatments %>% select(year, plotid, water, water4, snow, sun_date,precip_est_mm),
                "leaftraits"=lt.plantyr %>% 
                  select(year, round, day, plotid, plant, VWC.sm, any_of(c(leaftraits, fitnesstraits))) %>% 
                  arrange(year, round, plotid) %>% drop_na(plant) %>% filter(plant!="NA"),
                "phystraits"=pt.plantyr %>% 
                  select(year, round, plotid, plant, VWC.plant, VWC.sm, any_of(c(phystraits, fitnesstraits))) %>% 
                  arrange(year, round, plotid) %>% drop_na(plant) %>% filter(plant!="NA"))

purrr::walk(names(alldata), ~write_tsv(alldata[[.]], paste0("data/",., ".tsv"), na=""))

The following datasets are available in tab-delimited (TSV) format under the ‘data’ folder.

Treatments

The precipitation and snowmelt treatments in the split-plot design, with snowmelt dates for each plot and precipitation totals for each subplot.
Column	Description
year	year
plotid	plot number and subplot letter
water	precipitation treatment, with mock rainouts combined with controls
water4	precipitation treatment, with mock rainouts separate from controls
snow	snowmelt treatment (in 2019, plot 2 was covered but coded as normal)
sun_date	snowmelt date in each plot, determined by sunlight threshold
precip_est_mm	estimated summer precipitation in each subplot (see Methods S1)

Morphological traits

Column	Description
year	year
round	the first or second round of measurements (early/late summer)
day	the day of the year measurement was taken
plotid	plot number and subplot letter
plant	unique plant identifier
VWC.sm	the soil moisture averaged across the subplot on the date nearest the measurement
sla	specific leaf area (cm² g⁻¹)
trichome_density	trichome density (cm⁻²)
water_content	leaf water content (g H₂O g⁻¹)
survived	whether the plant survived to the next year
RGR	relative growth rate from this year to the next (year⁻¹)
flowered	whether the plant flowered the next year

Physiological traits

Column	Description
year	year
round	the date the measurement was taken
plotid	plot number and subplot letter
plant	unique plant identifier
VWC.plant	the soil moisture next to the plant on the date of measurement
VWC.sm	the soil moisture averaged across the subplot on the date nearest the measurement
photosynthesis	photosynthetic rate (µmol CO₂ m⁻² s⁻¹)
conductance	stomatal conductance (mol H₂O m⁻² s⁻¹)
WUE	intrinsic water-use efficiency (µmol CO₂ mol⁻¹ H₂O)
survived	whether the plant survived to the next year
RGR	relative growth rate from this year to the next (year⁻¹)
flowered	whether the plant flowered the next year

Sample sizes

count_measurements <- function(df, tr) summarize(df, across(all_of(tr), ~ sum(!is.na(.x)))) %>% deframe()

data.frame(dataset = c("Morphological traits", "Physiological traits"),
           leaves = c(count_measurements(lt, "leaf_area_cm2"), count_measurements(pt, "photosynthesis")),
           plants = c(lt %>% group_by(plantid) %>% summarize(across(all_of("leaf_area_cm2"), mean, na.rm=T)) %>%
                        count_measurements("leaf_area_cm2"), 
                      pt %>% group_by(plantid) %>% summarize(across(all_of("photosynthesis"), mean, na.rm=T)) %>%
                        count_measurements("photosynthesis"))) %>% kable()

dataset	leaves	plants
Morphological traits	600	264
Physiological traits	327	275

Trait correlations

full_join(pt %>% select(year, date, plantid, conductance, WUE, photosynthesis),
          lt %>% select(year, year.round, plantid, sla, water_content, trichome_density)) %>% 
  ggcorr(hjust=0.85, layout.exp=2, label=T, label_round=2)

Replicated split-plot models

Table 1

lt.subplot.r2 <- lt.subplotround %>% filter(round=="2") 
lt.plantyr.r2 <- lt.plantyr %>% filter(round=="2")

options(contrasts = c("contr.sum", "contr.poly"))
mod.split <- c(map(set_names(leaftraits), 
                   ~ lmer(lt.plantyr.r2[[.x]] ~ year * snow * water + (1|snow:plot) + (1|snow:plotid),
                          data=lt.plantyr.r2)),
               map(set_names(phystraits), 
                   ~ lmer(pt.plantyr[[.x]]    ~ year * snow * water + (1|snow:plot) + (1|snow:plotid) + (1|round),
                          data=pt.plantyr)))

mod.split.coefs <- map_dfr(mod.split, tidy, .id="trait")
mod.split.tests <- map(mod.split, anova, ddf="Ken") %>% map_dfr(tidy, .id="trait")
mod.split.emm <- map_dfr(mod.split, ~ summary(emmeans(ref_grid(.x), ~ year * snow * water)), .id = "trait")
  anova_table(mod.split.tests, mod.split) %>% kable()

trait	year	snow	water	year:snow	year:water	snow:water	year:snow:water	R2m
sla	0.0000	0.1519	0.6661	0.0120	0.7884	0.2773	0.1535	0.6027
trichome_density	0.0000	0.4894	0.5712	0.5986	0.1153	0.9208	0.5156	0.2094
water_content	0.1868	0.6817	0.0037	0.0756	0.0068	0.3675	0.4283	0.2171
photosynthesis	0.0620	0.6319	0.3075	0.4575	0.4112	0.9876	0.9137	0.1086
conductance	0.0009	0.4670	0.0023	0.4766	0.0024	0.6613	0.6894	0.3456
WUE	0.0000	0.2108	0.0002	0.4673	0.0026	0.3262	0.8511	0.4422

mod.split.tests.water <- map(mod.split, multcomp::glht, 
                             linfct = multcomp::mcp(water = c("Addition - Control = 0", "Reduction - Control = 0"))) %>% 
  map_dfr(tidy, .id="trait")

anova_table(mod.split.tests, mod.split) %>% 
  left_join(mod.split.tests.water %>% select(trait, contrast, adj.p.value) %>% 
              pivot_wider(names_from=contrast, values_from=adj.p.value)) %>% kable()

trait	year	snow	water	year:snow	year:water	snow:water	year:snow:water	R2m	Addition - Control	Reduction - Control
sla	0.0000	0.1519	0.6661	0.0120	0.7884	0.2773	0.1535	0.6027	0.5794	0.9137
trichome_density	0.0000	0.4894	0.5712	0.5986	0.1153	0.9208	0.5156	0.2094	0.5794	0.9603
water_content	0.1868	0.6817	0.0037	0.0756	0.0068	0.3675	0.4283	0.2171	0.0034	0.2324
photosynthesis	0.0620	0.6319	0.3075	0.4575	0.4112	0.9876	0.9137	0.1086	0.5205	0.5977
conductance	0.0009	0.4670	0.0023	0.4766	0.0024	0.6613	0.6894	0.3456	0.1368	0.0040
WUE	0.0000	0.2108	0.0002	0.4673	0.0026	0.3262	0.8511	0.4422	0.0051	0.0011

Figure 1

plot_split_plot <- function(emm, data, traits.plot, geom = "point", plot.emm = TRUE, facets="year") {
  
  data.long <- data %>% select(c(all_of(traits.plot), all_of(facets), water, snow)) %>% 
    pivot_longer(all_of(traits.plot), names_to="trait") %>% 
    mutate(trait = fct_relevel(trait, traits.plot), snow=fct_relevel(snow, "Normal")) %>% 
    drop_na(value, water) %>% filter(trait %in% traits.plot)
  
  traitnames.multiline <- ifelse(nchar(traitnames.units) > 30, str_replace(traitnames.units, fixed("("),"\n("), traitnames.units)
  
  split_plot <- emm %>% filter(trait %in% traits.plot) %>% 
    mutate(trait = fct_relevel(trait, traits.plot), snow=fct_relevel(snow, "Normal")) %>%  
  ggplot(aes(x=water, color=snow)) +
    labs(x="Precipitation", y="Standardized trait", color="Snowmelt") + 
    scale_y_continuous(position="right") + scale_x_discrete(guide=guide_axis(angle=90)) + 
    scale_color_manual(values=snow_pal_grey, guide=guide_legend(override.aes = list(shape=15, size=ifelse(plot.emm,5,1)))) + 
    theme_minimal() + theme(text=element_text(color="black", size=14), axis.text = element_text(color="black"),
                            axis.title.y= element_blank(),
                            panel.grid.major.x = element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y = element_blank(),
                            panel.border = element_rect(fill=NA, colour = "black"), 
                            panel.spacing=unit(0, "pt"), plot.margin = margin(0,0,0,0, "pt"), legend.position = "top") +
    switch(facets, 
           year = facet_grid(trait ~ year, scales="free_y", switch="y", 
                             labeller = as_labeller(c(traitnames.multiline, set_names(2018:2020)))),
           year.round = facet_grid(trait ~ year.round, scales="free_y", switch="y", 
                                   labeller = as_labeller(c(traitnames.multiline, set_names(levels(lt$year.round)))))) + 
    switch(geom,
           boxplot = geom_boxplot(data=data.long, aes(y=value), position=position_dodge(width=0.8), show.legend=!plot.emm,
                                  fatten=ifelse(plot.emm, NULL, 1), outlier.size=0.5),
           violin =   geom_violin(data=data.long, aes(y=value), position=position_dodge(width=0.8), show.legend=F),
           point =     geom_point(data=data.long, aes(y=value), position=position_dodge(width=0.8), shape=3))
  
  if(plot.emm) { return(split_plot + 
                          geom_linerange(aes(ymin=emmean-SE, ymax=emmean+SE,), position=position_dodge(width=0.8), size=1, show.legend=F) +
                          geom_point(aes(y=emmean), position=position_dodge(width=0.8), shape="-", size=16))
  } else return(split_plot)
}

#Vegetative and physiological traits - plot of EMMs and subplot means
grid.arrange(
  plot_split_plot(mod.split.emm, lt.subplotround %>% filter(round==2), leaftraits), #second round only 
  plot_split_plot(mod.split.emm, pt.subplot, phystraits)+ 
    guides(color = guide_legend(override.aes = list(color="white", shape=15, size=5))) + 
    theme(legend.title = element_text(color = "white"), legend.text = element_text(color = "white")), nrow=1)

Repeated measures

Table S1

# get plants that were sampled more than once for the first three rounds
lt.long.3 <- lt %>% pivot_longer(all_of(leaftraits), names_to = "trait") %>% arrange(date) %>% drop_na(value) %>%
  filter(year.round %in% c("2018.1","2018.2","2019.1"))
lt.long.repeated.3 <- lt.long.3  %>% filter(paste(trait, plantid) %in% 
                     (lt.long.3 %>% group_by(plantid, trait) %>% summarize(n = length(unique(year.round))) %>% 
                        filter(n>1) %>% mutate(trait_plantid = paste(trait, plantid)) %>% pull(trait_plantid))) %>% 
  group_by(trait, plot, plotid, plantid, water, snow, year.round) %>% summarize(value=mean(value)) 
lt.repeated.3 <- lt.long.repeated.3 %>% pivot_wider(names_from="trait")
lt.repeated.3.subplotround <- lt.repeated.3 %>% group_by(year.round, water, snow, plotid) %>% 
  summarize(across(all_of(leaftraits), mean, na.rm=T))

lt.mod.split.rm.3 <- map(set_names(leaftraits), 
                  ~ lmer(lt.repeated.3[[.x]] ~  year.round * water * snow + (1|snow:plot) + (1|snow:plotid) + (1|plantid), data=lt.repeated.3))
lt.mod.split.rm.3.tests <- map(lt.mod.split.rm.3, anova, ddf="Ken") %>% map_dfr(tidy, .id="trait")
anova_table(lt.mod.split.rm.3.tests, lt.mod.split.rm.3) %>% kable()

trait	water	snow	year.round:water	year.round:snow	water:snow	year.round:water:snow	R2m
sla	0.8730	0.3975	0.7454	0.0636	0.7876	0.8369	0.5824
trichome_density	0.9137	0.1228	0.0874	0.3010	0.7593	0.7196	0.1529
water_content	0.0271	0.6476	0.0000	0.1614	0.7026	0.5422	0.6841

Figure S2

lt.mod.split.rm.3.emm <- map_dfr(lt.mod.split.rm.3, ~ summary(emmeans(.x, ~  year.round * water * snow)), .id="trait")
plot_split_plot(lt.mod.split.rm.3.emm, lt.repeated.3.subplotround, leaftraits, facets="year.round")

Models with snowmelt date and total precipitation

Table 2

P-values

lt.subplot.r2 <- lt.subplotround %>% filter(round=="2")
lt.plantyr.r2 <- lt.plantyr %>% filter(round=="2")

options(contrasts = c("contr.sum", "contr.poly"))
mod.abs <- c(map(set_names(leaftraits), 
                 ~ lmer(lt.plantyr.r2[[.x]] ~ sun_date * precip_est_mm + (1|plot) + (1|plotid), data=lt.plantyr.r2)),
             map(set_names(phystraits), 
                 ~ lmer(pt.plantyr[[.x]]    ~ sun_date * precip_est_mm + (1|plot) + (1|plotid) + (1|round),
                        data=pt.plantyr)))

mod.abs.coefs <- map_dfr(mod.abs, tidy, .id="trait")
mod.abs.tests <- map(mod.abs, anova, ddf="Ken") %>% map_dfr(tidy, .id="trait")
anova_table(mod.abs.tests, mod.abs) %>% kable()

trait	sun_date	precip_est_mm	sun_date:precip_est_mm	R2m
sla	0.0000	0.4675	0.4396	0.5560
trichome_density	0.0015	0.0066	0.0814	0.1893
water_content	0.9537	0.2542	0.3505	0.1781
photosynthesis	0.0052	0.0973	0.1531	0.0820
conductance	0.0003	0.0247	0.0821	0.2300
WUE	0.0022	0.0099	0.0554	0.2054

Regression slopes

mod.abs.coefs %>% filter(effect=="fixed", term!="(Intercept)") %>% pivot_wider(id_cols=trait, names_from=term, values_from=c("estimate","std.error")) %>% kable()

trait	estimate_sun_date	estimate_precip_est_mm	estimate_sun_date:precip_est_mm	std.error_sun_date	std.error_precip_est_mm	std.error_sun_date:precip_est_mm
sla	6.6467	-0.7503	0.0067	0.9560	1.0007	0.0084
trichome_density	-1.3730	-1.7834	0.0090	0.4267	0.6503	0.0051
water_content	-0.0001	0.0024	0.0000	0.0018	0.0020	0.0000
photosynthesis	0.2752	0.1949	-0.0013	0.0927	0.1165	0.0009
conductance	0.0072	0.0040	0.0000	0.0017	0.0017	0.0000
WUE	-1.7948	-1.1076	0.0063	0.5112	0.4252	0.0033

Figure 2

plot_abs <- function(mod, data, traits.plot, tests=NULL, facets = "trait") {
  data.long <- data %>% select(all_of(traits.plot), water, snow, year, sun_date, precip_est_mm) %>% 
    pivot_longer(all_of(traits.plot), names_to="trait") %>% 
    mutate(trait = fct_relevel(trait, traits.plot), snow=fct_relevel(snow, "Normal"),
                      value = (value - alltraits.mean[as.character(trait)])/alltraits.sd[as.character(trait)]) %>% 
    drop_na(value, water) %>% filter(trait %in% traits.plot)
  
  precip.breaks <- seq(25,175, by=25)
  grid_emmeans <- function(mods, use.rounds = F) {
    grid.points <- list(sun_date=range(treatments$sun_date), precip_est_mm=precip.breaks)
    if(use.rounds)  map_dfr(mods, ~ summary(emmeans(.x, ~ precip_est_mm * sun_date * round, at=grid.points)), .id="trait")
    else            map_dfr(mods, ~ summary(emmeans(.x, ~ precip_est_mm * sun_date,         at=grid.points)), .id="trait")
  }
  emm.grid <- mod[traits.plot] %>% grid_emmeans(use.rounds= facets=="trait.round") %>% 
    mutate(trait = fct_relevel(trait, traits.plot),
           emmean = (emmean - alltraits.mean[as.character(trait)])/alltraits.sd[as.character(trait)])
  
  traitnames.multiline <- ifelse(nchar(traitnames) > 30, str_replace(traitnames, fixed("("),"\n("), traitnames)
  if(!is.null(tests)) {
    tests.sig <- tests %>% filter(p.value < 0.05) %>% 
          mutate(term.abbr = str_replace_all(term, c(sun_date="Sn", precip_est_mm="Pr",`:`="\U00D7"))) %>% 
          group_by(trait) %>% summarize(terms.sig = paste0("(",paste(term.abbr, collapse=", "),")")) %>% deframe
    traitnames.multiline[] <- paste(traitnames.multiline, replace_na(tests.sig[alltraits], ""))
  }
  
  abs_plot <- ggplot(emm.grid, aes(x=sun_date, y=emmean, color=precip_est_mm, group=precip_est_mm)) + 
    geom_point(data = data.long, aes(y=value, shape=year), size=2) + geom_line() + 
    scale_color_gradientn(colors=rev(water_pal), 
                          values=rev(c(1, (treatments %>% group_by(water) %>% summarize(across(precip_est_mm, mean)) %>%
                                             pull(precip_est_mm) %>% scales::rescale(from=range(precip.breaks)))[2],0)),
                          breaks=precip.breaks) + 
    labs(x="Snowmelt date (day of year)", y="Standardized trait", color="Summer\nprecipitation\n(mm)", shape="Year") + 
    scale_y_continuous(position="right")+
    theme_minimal() + theme(text=element_text(color="black", size=14), axis.text = element_text(color="black"),
                            panel.border = element_rect(fill=NA, colour = "black")) +
    switch(facets, 
           trait =     facet_wrap(vars(trait), scales="fixed", dir = "v", ncol=2, labeller = as_labeller(traitnames.multiline)),
           trait.round = facet_grid(trait ~ round, scales="free_y", switch="y", 
                                   labeller = as_labeller(c(traitnames.multiline, 
                                                            set_names(paste("Round",1:2), 1:2)))))
    
  return(abs_plot)
}

plot_abs(mod.abs, bind_rows(lt.subplot.r2, pt.subplot), c(leaftraits[2],phystraits[c(3,1,2)])) +
  scale_y_continuous(limits=c(-1.3, 1.8), breaks=c(-2:1))#subplot means

Models with soil moisture

Table 3

lt.subplot.r2 <- lt.subplotround %>% filter(round=="2")
lt.plantyr.r2 <- lt.plantyr %>% filter(round=="2")

pt.plantyr <- pt.plantyr %>% mutate(VWC = ifelse(is.na(VWC.plant), VWC.sm, VWC.plant)) #phys traits had the VWC measured next to the plant in 2019 and 2020 only - substitute VWC in subplot on closest date
pt <- pt %>% mutate(VWC = ifelse(is.na(VWC.plant), VWC.sm, VWC.plant)) 

options(contrasts = c("contr.sum", "contr.poly"))
mod.sm <- c(map(set_names(leaftraits), 
                  ~ lmer(lt.plantyr.r2[[.x]] ~ year * VWC + (1|plotid), data=lt.plantyr.r2)),
            map(set_names(phystraits[1:2]), 
                  ~ lmer(pt.plantyr[[.x]]    ~ year * VWC + I(VWC^2) + (1|plotid) + (1|round), data=pt.plantyr)),
            map(set_names(phystraits[3]), 
                  ~ lmer(pt.plantyr[[.x]]    ~ year * VWC            + (1|plotid) + (1|round), data=pt.plantyr)))
#significant effect of quadratic term for sepal_width, photosynthesis, conductance, but not year * VWC^2

mod.sm.coefs <- map_dfr(mod.sm, tidy, .id="trait")
mod.sm.tests <- map(mod.sm, anova, ddf="Ken") %>% map_dfr(tidy, .id="trait")
anova_table(mod.sm.tests, mod.sm) %>% kable()

trait	year	VWC	year:VWC	I(VWC^2)	R2m
sla	0.0116	0.1658	0.3953	NA	0.5847
trichome_density	0.0001	0.9058	0.0330	NA	0.1950
water_content	0.1778	0.0001	0.0829	NA	0.1835
photosynthesis	0.4488	0.0002	0.6687	0.0055	0.1497
conductance	0.0009	0.0001	0.0219	0.0389	0.3858
WUE	0.0000	0.0000	0.0030	NA	0.4346

Figure 3

plot_sm <- function(mod, data, traits.plot, tests=NULL) {
  data.long <- data %>% select(all_of(traits.plot), water, snow, year, sun_date, precip_est_mm, VWC) %>% 
    pivot_longer(all_of(traits.plot), names_to="trait") %>% 
    mutate(trait = fct_relevel(trait, traits.plot), snow=fct_relevel(snow, "Normal"),
                      value = (value - alltraits.mean[as.character(trait)])/alltraits.sd[as.character(trait)]) %>% 
    drop_na(value, water) %>% filter(trait %in% traits.plot)
  
  VWC.range <- range(data$VWC, na.rm=T)
  emm.grid <- mod[traits.plot] %>% 
    map_dfr(~ emmeans(.x, ~ year*VWC, at=list(VWC=seq(from=floor(VWC.range[1]*10)/10, to=ceiling(VWC.range[2]*10)/10, by=0.1))) %>% as_tibble(), .id="trait") %>% 
    mutate(trait = fct_relevel(trait, traits.plot),
           emmean = (emmean - alltraits.mean[as.character(trait)])/alltraits.sd[as.character(trait)]) %>% 
    left_join(data %>% drop_na(VWC) %>% group_by(year) %>% summarize(VWC.min = min(VWC), VWC.max=max(VWC))) %>%
    group_by(year) %>% filter(VWC >= VWC.min-0.1, VWC <= VWC.max+0.1)
  
  traitnames.multiline <- ifelse(nchar(traitnames) > 30, str_replace(traitnames, fixed("("),"\n("), traitnames)
  if(!is.null(tests)) {
    tests.sig <- tests %>% filter(p.value < 0.05) %>% 
          mutate(term.abbr = str_replace_all(term, c(year="Yr", VWC="SM",`:`="\U00D7"))) %>% 
          group_by(trait) %>% summarize(terms.sig = paste0("(",paste(term.abbr, collapse=", "),")")) %>% deframe
    print(tests.sig[traits.plot])
    traitnames.multiline[] <- paste(traitnames.multiline, replace_na(tests.sig[alltraits], ""))
  }
  
  sm_plot <- ggplot(emm.grid, aes(x=VWC, y=emmean, color=year, group=year)) + 
    facet_wrap(vars(trait), scales="fixed", dir="v", ncol=ifelse(length(traits.plot)>3,2,4), labeller = as_labeller(traitnames.multiline))+
    geom_point(data = data.long, aes(y=value, color=year), size=1) +  
    geom_line(size=2) + 
    scale_color_manual(values=year_pal) + ylab("Standardized trait") + labs(color="Year") + 
    theme_minimal() + theme(text=element_text(color="black", size=14), axis.text = element_text(color="black"),
                            panel.border = element_rect(fill=NA, colour = "black"))
  return(sm_plot)
}

grid.arrange(
  plot_sm(mod.sm, lt.plantyr.r2, leaftraits) + 
    xlab("Summer mean of soil moisture in subplot (%VWC)"),
  plot_sm(mod.sm, pt.plantyr, phystraits) + 
    xlab("Soil moisture next to plant or in subplot on closest date (%VWC)") +
    guides(color = guide_legend(override.aes = list(color="white", shape=15, size=5))) + 
    theme(legend.title = element_text(color = "white"), legend.text = element_text(color = "white")), ncol=1)

Fitness vs. treatments

Table S2

#Model effect of treatments on survival and flowering of vegetative plants - both years combined
cen.status.fitness <- cen.status.long %>% filter(status=="vegetative", census %in% c("2018","2019","2020")) %>%
  full_join(cen.RGR.long %>% drop_na(RGR) %>% filter(census!="18r2") %>% select(-census)) %>% left_join(treatments) %>% left_join(sm.subplotyear)

options(contrasts = c("contr.sum", "contr.poly"))
fitness.mod <- map(set_names(fitnesstraits), 
                    ~ glmmTMB(cen.status.fitness[[.x]] ~ year * snow * water + (1|snow:plot) + (1|snow:plotid), data=cen.status.fitness, 
                              family=ifelse(.x =="RGR", "gaussian", "binomial"))) 

fitness.mod.coefs <- map_dfr(fitness.mod, tidy, .id="trait")
fitness.mod.tests <- map(fitness.mod, car::Anova, type=3) %>% map_dfr(tidy, .id="trait")
fitness.mod.emm <- map_dfr(fitness.mod, ~ summary(emmeans(ref_grid(.x), ~ year * snow * water)), .id = "trait") %>% 
  mutate(uSE = ifelse(trait=="RGR", emmean+SE, plogis(emmean+SE)), 
         lSE = ifelse(trait=="RGR", emmean-SE, plogis(emmean-SE)), 
         emmean = ifelse(trait=="RGR", emmean, plogis(emmean)),
         trait = fct_relevel(trait, fitnesstraits)) 
Anova_table(fitness.mod.tests) %>% kable()

Trait	year	snow	water	year:snow	year:water	snow:water	year:snow:water
survived	0.0654	0.7972	0.4211	0.5528	0.3943	0.5450	0.0168
RGR	0.0000	0.0897	0.0785	0.0485	0.0000	0.5893	0.0082
flowered	0.0000	0.8817	0.0885	0.3232	0.4297	0.6055	0.0903

Figure S3

ggplot(fitness.mod.emm %>% mutate(snow=fct_relevel(snow, "Normal", "Early")), aes(x=water, color=snow, y=emmean)) + 
  facet_grid(trait ~ year, scales="free_y", labeller = as_labeller(c(fitnessnames, set_names(2018:2020))), switch="y", as.table=F)+
  geom_hline(data=data.frame(trait = fitnesstraits, plot=c(NA,NA,"yes")), aes(linetype=plot, yintercept=0), color="grey")+
  guides(linetype="none")+
  geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge(width=0.8), size=1, show.legend=F) +
  geom_point(position=position_dodge(width=0.8), shape="-", size=16) + 
  labs(x="Precipitation", color="Snowmelt") + 
  scale_y_continuous(position="right", expand=expansion(mult=c(0,0.05)))  + scale_x_discrete(guide=guide_axis(angle=90))+ 
  scale_color_manual(values=snow_pal_grey, guide=guide_legend(override.aes = list(shape=15, size=5))) + 
  theme_minimal() + theme(text=element_text(color="black", size=14), axis.text = element_text(color="black"),
                          axis.title.y= element_blank(),
                          panel.border = element_rect(fill=NA, colour = "black"),
                          panel.spacing.x=unit(0, "lines"), panel.spacing.y=unit(1, "lines"),
                          plot.margin = margin(0,0,0,0, "pt"), legend.position = "top")

Selection

lt.plantyr.r2 <- lt.plantyr %>% filter(round=="2")

step_selection_mod <- function(trait, response, treatment, set, drop=TRUE) { 
  family <- ifelse(response=="RGR", "gaussian", "binomial") 
  re <- ifelse(treatment =="snow", "snow:plot", "plot")
  t.plantyr <- list(lt = lt.plantyr.r2, pt = pt.plantyr)[[set]] %>% rename("trait"=trait)
  options(contrasts = c("contr.sum", "contr.poly"))
  mod <- glmmTMB(formula(paste(response, "~ trait * year * ", treatment," + (1|",re,")")), 
                 data=t.plantyr, family=family)
  if(drop) {
    droporder <- c(paste0("trait:year:", treatment)) 
    for(term_test in droporder) {
      pval <- car::Anova(mod) %>% tidy %>% filter(term==term_test) %>% pull(p.value)
      if(pval < 0.05) return(mod)
      else mod <- update(mod, formula(paste0(". ~ . -", term_test)))
    }
  }
  return(mod)
}

selection.combos <- expand_grid(response = fitnesstraits, treatment=c("water","snow"), set = c("lt", "pt")) 
mod.selection <- selection.combos  %>% 
  full_join(tibble(trait=c(leaftraits, phystraits), 
                   set = c(rep("lt", length(leaftraits)), rep("pt", length(phystraits))))) %>% 
  mutate(full_model = pmap(., step_selection_mod, drop=F), 
         model =      pmap(., step_selection_mod, drop=T),
         coefs = map(model, tidy),
         tests = map(model, ~ tidy(car::Anova(., type=3))))

emt.selection <- mod.selection %>% 
  mutate(yrtrend = map2(full_model, treatment, 
                        ~tidy(emtrends(.x, specs=c("year",as.character(.y)), var="trait")))) %>%
  select(response, treatment, trait, yrtrend) %>% unnest(yrtrend) %>% 
  mutate(sd.trait=alltraits.sd[trait], 
         trait.trend.sd = trait.trend*sd.trait, std.error.sd=std.error*sd.trait) #multiply by SDto get change in fitness per SD

Table S3

mod.selection %>% group_by(treatment, response) %>%
  group_walk(function(.x,.y) {cat(paste("Effect of trait on whether a plant",
                                        ifelse(.y$response=="RGR", "grew faster", as.character(.y$response)),
                                        "under different",.y$treatment,"conditions")); 
    Anova_table(.x %>% unnest(tests)) %>% kable() %>% print})

Effect of trait on whether a plant flowered under different snow conditions

Trait	(Intercept)	trait	year	snow	trait:year	trait:snow	year:snow	trait:year:snow
sla	0.3317	0.1159	0.4337	0.7854	0.4741	0.9760	0.0315	0.0292
trichome_density	0.0001	0.5036	0.1254	0.2700	0.9762	0.5908	0.8727	NA
water_content	0.9481	0.6631	0.1364	0.4622	0.1781	0.5110	0.3869	NA
photosynthesis	0.0049	0.8196	0.0241	0.6761	0.7605	0.8938	0.9809	NA
conductance	0.0000	0.1004	0.0802	0.6564	0.5563	0.2719	0.7909	NA
WUE	0.2426	0.2001	0.0113	0.0846	0.1994	0.0825	0.3988	NA

Effect of trait on whether a plant grew faster under different snow conditions

Trait	(Intercept)	trait	year	snow	trait:year	trait:snow	year:snow
sla	0.7052	0.3816	0.8329	0.2508	0.5087	0.3092	0.5174
trichome_density	0.2207	0.5082	0.0726	0.5553	0.5807	0.4760	0.6075
water_content	0.0360	0.0185	0.3471	0.7810	0.2428	0.8452	0.6853
photosynthesis	0.9038	0.3558	0.5644	0.3127	0.1033	0.9416	0.2709
conductance	0.4514	0.8889	0.8290	0.3335	0.0692	0.5832	0.2541
WUE	0.1721	0.5663	0.0367	0.6402	0.6375	0.7271	0.2155

Effect of trait on whether a plant survived under different snow conditions

Trait	(Intercept)	trait	year	snow	trait:year	trait:snow	year:snow
sla	0.0815	0.4025	0.8758	0.2835	0.9105	0.3375	0.5697
trichome_density	0.0068	0.2875	0.6033	0.5454	0.9035	0.3011	0.6245
water_content	0.4344	0.7356	0.7422	0.3973	0.8272	0.4326	0.9464
photosynthesis	0.0000	0.2788	0.3959	0.3410	0.4093	0.4677	0.5303
conductance	0.0000	0.5287	0.9881	0.9740	0.9567	0.4986	0.9033
WUE	0.0001	0.2857	0.2494	0.0283	0.2041	0.0354	0.6202

Effect of trait on whether a plant flowered under different water conditions

Trait	(Intercept)	trait	year	water	trait:year	trait:water	year:water
sla	0.5329	0.1689	0.9884	0.2451	0.9449	0.3231	0.0555
trichome_density	0.0007	0.5407	0.1706	0.7915	0.9078	0.9320	0.3337
water_content	0.8205	0.5025	0.1750	0.6309	0.2268	0.7004	0.3672
photosynthesis	0.0495	0.6821	0.0300	0.1912	0.8794	0.3623	0.8151
conductance	0.0004	0.5845	0.0654	0.2286	0.4345	0.1494	0.7986
WUE	0.3530	0.2789	0.0076	0.5942	0.1261	0.8670	0.8587

Effect of trait on whether a plant grew faster under different water conditions

Trait	(Intercept)	trait	year	water	trait:year	trait:water	year:water
sla	0.6134	0.9749	0.7442	0.7135	0.9586	0.7040	0.4100
trichome_density	0.3216	0.4079	0.1476	0.5695	0.2681	0.9229	0.0053
water_content	0.8444	0.6500	0.4799	0.8091	0.6138	0.9147	0.1149
photosynthesis	0.6730	0.5055	0.1163	0.0069	0.6106	0.3151	0.0400
conductance	0.0320	0.3380	0.3055	0.0144	0.6306	0.4573	0.0343
WUE	0.7578	0.2953	0.5070	0.0731	0.9072	0.2128	0.0087

Effect of trait on whether a plant survived under different water conditions

Trait	(Intercept)	trait	year	water	trait:year	trait:water	year:water
sla	0.0344	0.2066	0.9756	0.6896	0.8533	0.7193	0.5622
trichome_density	0.0085	0.3653	0.8689	0.6346	0.9392	0.8250	0.2619
water_content	0.2525	0.4137	0.8880	0.9807	0.9246	0.9974	0.2139
photosynthesis	0.0000	0.1342	0.4205	0.5193	0.4314	0.3896	0.9250
conductance	0.0000	0.6037	0.9599	0.0738	0.9287	0.0133	0.4460
WUE	0.0019	0.3059	0.1869	0.6154	0.1268	0.5005	0.8476

plot_selection_facets <- function(response, treatment) {
  emt.selection.filtered <- emt.selection %>% 
    select(!!!treatment, response, year, trait, p.value) %>% 
    drop_na(!!!treatment) %>% filter(response==!!response) %>% 
    mutate(sig = ifelse(p.value < 0.05,"signif","ns"))
  
  p <- bind_rows(lt.plantyr.r2, pt.plantyr) %>% 
    pivot_longer(all_of(c(leaftraits, phystraits)), names_to="trait") %>% 
    left_join(emt.selection.filtered) %>% 
    mutate(trait = fct_relevel(trait, c(leaftraits, phystraits))) %>% 
    ggplot(aes_string(y=response, x="value", color=treatment, linetype="year")) + 
    facet_wrap(vars(trait), scales="free_x", labeller = as_labeller(traitnames.units))+
    scale_color_manual(values=list(water=water_pal, snow=snow_pal)[[treatment]]) + 
    scale_linetype_manual(values=c("22","solid","42")) +
    labs(linetype="Year", shape="Year", x="", y=fitnessnames[response],
         color=c(water="Precipitation",snow="Snowmelt")[treatment]) + 
    theme_minimal() + 
    theme(text = element_text(color="black"), axis.text = element_text(color="black"),
          panel.border = element_rect(fill=NA, colour = "black"), 
          panel.spacing=unit(0, "pt"), plot.margin = margin(0,0,0,0, "pt"), 
          legend.position = "top", legend.key.width = unit(2, "lines"))
  
  if(response =="RGR") {
    p <- p + geom_hline(yintercept=0, color="grey") +
      geom_point(aes(shape=year)) + 
      geom_line(stat="smooth", size=1, se=F, method="lm")
  }  else {
    p <- p + geom_point(shape="|") + 
      geom_line(stat="smooth", size=1, se=F, 
                  method="glm", method.args = list(family="binomial")) +
      scale_y_continuous(labels=scales::label_percent(accuracy=1))
  }
  #get rightmost edge of geom_smooths to add signif labels
  smooth.ends <- layer_data(p, ifelse(response=="RGR", 3, 2)) %>% #hline adds a layer
    group_by(group, PANEL) %>% filter(x==max(x)) %>% ungroup() %>% select(-PANEL) %>% 
    bind_cols(emt.selection.filtered) %>% 
    filter(sig=="signif") %>% 
    mutate(trait = factor(trait, levels=c(leaftraits, phystraits))) 
  p + geom_point(data=smooth.ends, aes(x=x,y=y), color="black", shape=8, size=2, inherit.aes = F)
}

selection_plots <- expand_grid(treatment=c("water","snow"), response = fitnesstraits) %>% 
  mutate(plot = pmap(., plot_selection_facets))

Figure 4

(selection_plots$plot[[1]] + labs(tag="(a)")) + 
(selection_plots$plot[[3]] + labs(tag="(b)")) + 
  plot_layout(guides = "collect", ncol=1) & 
  theme(legend.position="top", plot.tag=element_text(face="bold"))

Figure 5

selection_plots$plot[[2]]

Figure S4

(selection_plots$plot[[4]] + guides(linetype="none", color="none") + labs(tag="(a)")) +
(selection_plots$plot[[5]]  + labs(tag="(b)"))+
(selection_plots$plot[[6]] + guides(linetype="none", color="none") + labs(tag="(c)")) +
  plot_layout(guides = "collect", ncol=1) & 
  theme(legend.position="top", plot.tag=element_text(face="bold"))

Timing

Figure S1

  timings_labels <- c(snowcloth="Snow cloths applied to early plots", 
                     meltdates="Mean early and normal snowmelt timing", 
                     waterdates="Summer precipitation treatments applied",
                     sm="Soil moisture recorded",         ph="Inflorescence height recorded", 
                     mt="Floral morphology recorded",     nt="Nectar traits recorded", 
                     pt= "Physiology traits recorded",    lt="Vegetative traits recorded", 
                     sds="Flower collection period",      cen="Rosette size and survival recorded")

plot_timings <- function(data) {
  data %>% mutate(variable=fct_reorder(variable, begin, na.rm=T, .desc=T),
                  drawline = ifelse(variable %in% c("lt","meltdates"), NA, TRUE)) %>% 
ggplot(aes(y=variable, color=year))+ 
  geom_linerange(aes(xmin=begin, xmax=end, linetype=drawline), 
                 position = position_dodge2(width=0.7, reverse = T), size=1.2, show.legend=FALSE) + 
  geom_text(data=data %>% drop_na(plots) %>% filter(year=="2019"),
            aes(x = end+7, label=paste(ifelse(nchar(plots)>1, "Plots","Plot"), str_replace(plots,","," & "))), 
            position = position_dodge2(width=0.8, reverse = T), size=4, hjust=0, show.legend=FALSE)+
  geom_point(aes(x=begin), position = position_dodge2(width=0.7, reverse = T), shape=15, size=2)+
  geom_point(aes(x=end), position = position_dodge2(width=0.7, reverse = T), shape=15, size=2) +
  scale_color_manual("Year", values=c(year_pal, `2021`=brewer.pal(8, name="Set2")[[4]]), guide=guide_legend(override.aes = list(size=5))) + 
  scale_y_discrete("", labels=timings_labels)+
  scale_x_continuous("Day of year", breaks=seq(80,240, by=20)) +
  theme_minimal() + theme(legend.position = "top", text=element_text(size=14, color="black"), axis.text = element_text(color="black"),
                            panel.border = element_rect(fill=NA, colour = "black"))
}

plot_timings(filter(timings, !variable %in% c("sds","mt","nt","ph")))

Experiment map

Methods S1

subplot_offset <- 1/6
size_meter <- 6
ggplot(treatments_map) + coord_fixed(clip="off")+
  geom_point(aes(x=plot_x, y=plot_y, fill=snow), size=7*size_meter, shape=22, stroke=1)+
  scale_fill_manual("Snowmelt", values=snow_pal, guide=guide_legend(override.aes = list(size = 2*size_meter)))+
  new_scale_fill()+
  geom_point(aes(x=plot_x+subplot_x*subplot_offset, y=plot_y+subplot_y*subplot_offset, fill=water4), size=2*size_meter, shape=22, stroke=1) + 
    geom_text(aes(x=plot_x+subplot_x*subplot_offset, y=plot_y+subplot_y*subplot_offset, label=subplot)) + 
  scale_fill_manual("Precipitation", values=water4_pal) +
  geom_text(aes(x=plot_x, y=plot_y, label=plot))+ theme_minimal() + 
  theme(legend.position="right", axis.title=element_blank(), axis.text = element_blank(), 
        axis.ticks = element_blank(), panel.background = element_blank(), panel.grid = element_blank())

Code for “Phenotypic plasticity and selection on leaf traits in response to snowmelt timing and summer precipitation”

Jocelyn Navarro, John M. Powers, Ayaka Paul, Diane R. Campbell

2022-02-16