Set up workspace

Load packages:

library(dplyr)
library(tidyr)
library(lme4)
library(lmerTest)
library(car)
library(emmeans)
library(pbkrtest)
library(parallel)
library(ggplot2)
library(geodist)

Read in data:

field<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-Field-2013-2016-archive.csv', header=T)

surv<-read.csv('~/Desktop/BCMA archive.R1/CFRs-2017surv.csv', header=T)

rp<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-RefPops-AllMerged.csv', header=T)

dry<-read.csv('~/Desktop/BCMA archive.R1/BCMA-final-DryDownCFR.csv', header=T)

nuc.div<-read.table("~/Desktop/BCMA archive.R1/Bs54.fold4.div.gene.chr.txt", header=T, na="-9")

Ecological drivers of selection

Variation in insect resistance across environments

Temporary arrays:

library(pbkrtest)

# Supplementary Table 1A
## subset data
use.a<-subset(field, Expt=="arrays")

## fit model
d4<-lmer(logLAR2 ~ Site + (1|Blk) + Geno + Category + (1|Fam) + Geno*Site + Geno*Category, data=use.a, 
         control=lmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5)), na.action = na.omit)
summary(d4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ Site + (1 | Blk) + Geno + Category + (1 | Fam) + Geno *  
##     Site + Geno * Category
##    Data: use.a
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 7778.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5676 -0.6813 -0.1647  0.5183  7.1082 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Blk      (Intercept) 0.0526425 0.22944 
##  Fam      (Intercept) 0.0004151 0.02037 
##  Residual             0.2504217 0.50042 
## Number of obs: 5193, groups:  Blk, 60; Fam, 10
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         0.82697    0.09276   60.09774   8.915 1.35e-12 ***
## SiteKC             -0.19488    0.10978   56.98266  -1.775  0.08120 .  
## SiteLA             -0.15946    0.10939   56.17866  -1.458  0.15047    
## SiteMM             -0.20490    0.11009   57.62952  -1.861  0.06783 .  
## SiteRM             -0.18783    0.10932   56.04588  -1.718  0.09130 .  
## SiteUA             -0.27554    0.10932   56.03514  -2.521  0.01459 *  
## GenoSS              0.15195    0.04653  284.11214   3.266  0.00122 ** 
## CategoryL           0.10875    0.08107   55.02170   1.341  0.18530    
## CategoryM           0.02283    0.07657   57.15397   0.298  0.76662    
## SiteKC:GenoSS      -0.12412    0.05285 5156.81455  -2.349  0.01888 *  
## SiteLA:GenoSS      -0.02542    0.05189 5156.57491  -0.490  0.62420    
## SiteMM:GenoSS      -0.05487    0.05406 5155.93256  -1.015  0.31017    
## SiteRM:GenoSS      -0.13429    0.05166 5156.76223  -2.600  0.00936 ** 
## SiteUA:GenoSS      -0.10267    0.05156 5157.47298  -1.991  0.04649 *  
## GenoSS:CategoryL   -0.04208    0.03755 5116.29865  -1.121  0.26249    
## GenoSS:CategoryM   -0.02329    0.03470 5119.50198  -0.671  0.50214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## random effects tests
ranova(d4)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## logLAR2 ~ Site + Geno + Category + (1 | Blk) + (1 | Fam) + Site:Geno + 
##     Geno:Category
##           npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>      19 -3889.1 7816.2                         
## (1 | Blk)   18 -4154.2 8344.4 530.21  1     <2e-16 ***
## (1 | Fam)   18 -3890.1 7816.2   1.98  1     0.1597    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## ANOVA with KR ddf
d4.a2<-anova(d4, ddf="Kenward-Roger") # takes a few minutes to run
d4.a2
## Type III Analysis of Variance Table with Kenward-Roger's method
##                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Site          2.84488 0.56898     5   52.0  2.2721 0.06077 .
## Geno          2.10329 2.10329     1    9.0  8.3990 0.01767 *
## Category      0.37910 0.18955     2   52.3  0.7569 0.47418  
## Site:Geno     2.95790 0.59158     5 5150.0  2.3623 0.03761 *
## Geno:Category 0.31615 0.15807     2 5121.5  0.6312 0.53198  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Figure 2A (left pane)
## extract least squares means
emmeans(d4, pairwise ~ Geno*Site)$emmeans
##  Geno Site emmean     SE  df asymp.LCL asymp.UCL
##  LL   401   0.871 0.0790 Inf     0.716     1.026
##  SS   401   1.001 0.0778 Inf     0.848     1.154
##  LL   KC    0.676 0.0775 Inf     0.524     0.828
##  SS   KC    0.682 0.0771 Inf     0.531     0.833
##  LL   LA    0.711 0.0770 Inf     0.561     0.862
##  SS   LA    0.816 0.0770 Inf     0.665     0.967
##  LL   MM    0.666 0.0779 Inf     0.513     0.819
##  SS   MM    0.741 0.0775 Inf     0.589     0.893
##  LL   RM    0.683 0.0769 Inf     0.532     0.834
##  SS   RM    0.679 0.0769 Inf     0.528     0.830
##  LL   UA    0.595 0.0769 Inf     0.445     0.746
##  SS   UA    0.623 0.0769 Inf     0.472     0.773
## 
## Results are averaged over the levels of: Category 
## Degrees-of-freedom method: asymptotic 
## Confidence level used: 0.95

## pairwise comparisons within environments
emmeans(d4, pairwise ~ Geno | Site)$contrasts
## Site = 401:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS  -0.13016 0.0414 Inf -3.146  0.0017 
## 
## Site = KC:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS  -0.00604 0.0377 Inf -0.160  0.8726 
## 
## Site = LA:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS  -0.10473 0.0363 Inf -2.885  0.0039 
## 
## Site = MM:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS  -0.07528 0.0393 Inf -1.914  0.0556 
## 
## Site = RM:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS   0.00414 0.0360 Inf  0.115  0.9084 
## 
## Site = UA:
##  contrast estimate     SE  df z.ratio p.value
##  LL - SS  -0.02748 0.0358 Inf -0.767  0.4430 
## 
## Results are averaged over the levels of: Category 
## Degrees-of-freedom method: asymptotic

Transplant experiments:

library(pbkrtest)
# Supplementary Table 1B

## Subset data
use.t<-subset(field, Expt=="transplants")

## Fit model
d3<-lmer(logLAR2 ~ Env + (1|Blk) + Geno + (1|Fam) + Geno*Env, data=use.t,
          control=lmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5)))
summary(d3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ Env + (1 | Blk) + Geno + (1 | Fam) + Geno * Env
##    Data: use.t
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 10267.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9669 -0.6982 -0.1579  0.6200  3.8882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Blk      (Intercept) 0.15584  0.3948  
##  Fam      (Intercept) 0.01019  0.1009  
##  Residual             0.88776  0.9422  
## Number of obs: 3674, groups:  Blk, 149; Fam, 19
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           1.9044     0.1290  191.5422  14.764  < 2e-16 ***
## EnvALD-2014          -0.6209     0.2019  167.6795  -3.075 0.002460 ** 
## EnvALD-2015           0.3895     0.1740  227.7899   2.238 0.026206 *  
## EnvGTH-2016          -0.8372     0.1806  358.4077  -4.637 4.97e-06 ***
## EnvMAH-2014          -0.7247     0.1951  201.2145  -3.715 0.000263 ***
## EnvMAH-2015          -0.9154     0.1877  200.7535  -4.876 2.20e-06 ***
## EnvSCH-2016          -0.2998     0.1516  201.2438  -1.978 0.049277 *  
## EnvSIL-2014          -0.9076     0.1936  188.1999  -4.689 5.26e-06 ***
## EnvSIL-2015          -1.2731     0.1759  187.2644  -7.238 1.14e-11 ***
## GenoSS                0.4168     0.1136  132.8703   3.670 0.000350 ***
## EnvALD-2014:GenoSS   -0.3611     0.1464 3538.4064  -2.466 0.013701 *  
## EnvALD-2015:GenoSS   -0.4212     0.1471 2728.9738  -2.863 0.004225 ** 
## EnvGTH-2016:GenoSS   -0.6919     0.1610 3287.7743  -4.297 1.78e-05 ***
## EnvMAH-2014:GenoSS   -0.2475     0.1548 3553.5025  -1.600 0.109797    
## EnvMAH-2015:GenoSS   -0.6139     0.1457 2600.6488  -4.212 2.61e-05 ***
## EnvSCH-2016:GenoSS   -0.4330     0.1190 2314.0075  -3.637 0.000282 ***
## EnvSIL-2014:GenoSS   -0.4764     0.1523 3551.2107  -3.129 0.001768 ** 
## EnvSIL-2015:GenoSS   -0.6775     0.1408 2489.5631  -4.813 1.58e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## random effects tests
ranova(d3)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## logLAR2 ~ Env + Geno + (1 | Blk) + (1 | Fam) + Env:Geno
##           npar  logLik   AIC    LRT Df Pr(>Chisq)    
## <none>      21 -5133.7 10309                         
## (1 | Blk)   20 -5312.2 10664 356.91  1  < 2.2e-16 ***
## (1 | Fam)   20 -5142.1 10324  16.80  1  4.148e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## ANOVA with KR ddf
d3.a2<-anova(d3, ddf="Kenward-Roger") # takes ~2.5 min to run
d3.a2 
## Type III Analysis of Variance Table with Kenward-Roger's method
##           Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## Env      199.208 24.9010     8  135.15 28.0462 < 2.2e-16 ***
## Geno       0.083  0.0830     1   15.96  0.0935    0.7638    
## Env:Geno  29.754  3.7193     8 2840.39  4.1884 5.333e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Figure 2A (center pane)
## extract least squares means
emmeans(d3, pairwise ~ Geno* Env)$emmeans
##  Geno Env      emmean     SE  df asymp.LCL asymp.UCL
##  LL   ALD-2013  1.904 0.1290 Inf     1.652     2.157
##  SS   ALD-2013  2.321 0.1260 Inf     2.074     2.568
##  LL   ALD-2014  1.284 0.1665 Inf     0.957     1.610
##  SS   ALD-2014  1.339 0.1630 Inf     1.020     1.659
##  LL   ALD-2015  2.294 0.1264 Inf     2.046     2.542
##  SS   ALD-2015  2.289 0.1221 Inf     2.050     2.529
##  LL   GTH-2016  1.067 0.1371 Inf     0.799     1.336
##  SS   GTH-2016  0.792 0.1117 Inf     0.573     1.011
##  LL   MAH-2014  1.180 0.1580 Inf     0.870     1.489
##  SS   MAH-2014  1.349 0.1543 Inf     1.047     1.652
##  LL   MAH-2015  0.989 0.1443 Inf     0.706     1.272
##  SS   MAH-2015  0.792 0.1395 Inf     0.519     1.065
##  LL   SCH-2016  1.605 0.0948 Inf     1.419     1.790
##  SS   SCH-2016  1.588 0.0905 Inf     1.411     1.766
##  LL   SIL-2014  0.997 0.1562 Inf     0.691     1.303
##  SS   SIL-2014  0.937 0.1542 Inf     0.635     1.239
##  LL   SIL-2015  0.631 0.1287 Inf     0.379     0.884
##  SS   SIL-2015  0.371 0.1312 Inf     0.113     0.628
## 
## Degrees-of-freedom method: asymptotic 
## Confidence level used: 0.95

## pairwise comparisons within environments
emmeans(d3, pairwise ~ Geno | Env)$contrasts
## Env = ALD-2013:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS  -0.41682 0.114 Inf -3.670  0.0002 
## 
## Env = ALD-2014:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS  -0.05573 0.123 Inf -0.452  0.6512 
## 
## Env = ALD-2015:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.00442 0.115 Inf  0.039  0.9693 
## 
## Env = GTH-2016:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.27507 0.136 Inf  2.016  0.0438 
## 
## Env = MAH-2014:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS  -0.16928 0.133 Inf -1.273  0.2030 
## 
## Env = MAH-2015:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.19711 0.113 Inf  1.749  0.0803 
## 
## Env = SCH-2016:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.01618 0.082 Inf  0.197  0.8436 
## 
## Env = SIL-2014:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.05963 0.130 Inf  0.458  0.6468 
## 
## Env = SIL-2015:
##  contrast estimate    SE  df z.ratio p.value
##  LL - SS   0.26066 0.106 Inf  2.455  0.0141 
## 
## Degrees-of-freedom method: asymptotic

Laboratory assay:

library(pbkrtest)
# Supplementary Table 1C
## subset data
use.l<-subset(field, Expt=="Lab_Tni")

## fit model
d5<-lmer(logLAR2 ~ (1|Blk) + Geno + (1|Fam), data=use.l, 
         control=lmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5)))
summary(d5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logLAR2 ~ (1 | Blk) + Geno + (1 | Fam)
##    Data: use.l
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 681.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6669 -0.4081  0.2268  0.6597  1.6656 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Fam      (Intercept) 0.0000   0.0000  
##  Blk      (Intercept) 0.1957   0.4424  
##  Residual             1.3450   1.1598  
## Number of obs: 214, groups:  Fam, 9; Blk, 5
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   3.5028     0.2296   5.3897  15.254 1.22e-05 ***
## GenoSS       -0.9513     0.1590 208.0067  -5.982 9.48e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## GenoSS -0.372
## convergence code: 0
## boundary (singular) fit: see ?isSingular

## random effects tests
ranova(d5)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## logLAR2 ~ Geno + (1 | Blk) + (1 | Fam)
##           npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>       5 -340.86 691.72                         
## (1 | Blk)    4 -348.70 705.40 15.683  1   7.49e-05 ***
## (1 | Fam)    4 -340.86 689.72  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## ANOVA with KR ddf
d5.a2<-anova(d5, ddf="Kenward-Roger") 
d5.a2
## Type III Analysis of Variance Table with Kenward-Roger's method
##      Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
## Geno 48.114  48.114     1 6.8801  35.772 0.00059 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Figure 2A (right pane)
## extract least squares means
emmeans(d5, pairwise ~ Geno)$emmeans
##  Geno emmean    SE   df lower.CL upper.CL
##  LL     3.50 0.230 5.10     2.92     4.09
##  SS     2.55 0.226 4.84     1.97     3.14
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

## pairwise comparisons within environments
emmeans(d5, pairwise ~ Geno)$contrasts
##  contrast estimate    SE   df t.ratio p.value
##  LL - SS     0.951 0.159 6.88 5.981   0.0006 
## 
## Degrees-of-freedom method: kenward-roger

Fecundity selection by herbivores

Effect of damage on the probability of reproduction:

library(car)
# Supplementary Table 2A
## filter field data to vernalized adult transplants
adult<-filter(field, Cohort != "s2013" & Expt=="transplants" & CensusYr < 2017)

## fit model
prrep2<-glmer(didrep ~ Env*LAR2 + (1|Blk), data=adult, family=binomial(link="logit"),
             control=glmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(prrep2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: didrep ~ Env * LAR2 + (1 | Blk)
##    Data: adult
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##   3144.7   3247.3  -1555.3   3110.7     3077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6727 -0.7020  0.2494  0.5977  5.4780 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Blk    (Intercept) 0.6283   0.7927  
## Number of obs: 3094, groups:  Blk, 132
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       0.739506   0.318324   2.323  0.02017 *  
## EnvALD-2015      -0.305176   0.411330  -0.742  0.45813    
## EnvGTH-2016      -0.527352   0.389808  -1.353  0.17610    
## EnvMAH-2014       0.733268   0.448010   1.637  0.10169    
## EnvMAH-2015      -0.378877   0.425719  -0.890  0.37348    
## EnvSCH-2016       2.460741   0.397835   6.185 6.20e-10 ***
## EnvSIL-2014      -1.174447   0.436319  -2.692  0.00711 ** 
## EnvSIL-2015      -1.151725   0.400879  -2.873  0.00407 ** 
## LAR2             -0.033620   0.013593  -2.473  0.01339 *  
## EnvALD-2015:LAR2  0.006939   0.016984   0.409  0.68284    
## EnvGTH-2016:LAR2  0.004065   0.035227   0.115  0.90813    
## EnvMAH-2014:LAR2  0.017932   0.018926   0.947  0.34341    
## EnvMAH-2015:LAR2 -0.010229   0.021379  -0.478  0.63234    
## EnvSCH-2016:LAR2 -0.026026   0.018573  -1.401  0.16112    
## EnvSIL-2014:LAR2 -0.073384   0.043195  -1.699  0.08934 .  
## EnvSIL-2015:LAR2  0.153389   0.034531   4.442 8.91e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(prrep2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: didrep
##                Chisq Df Pr(>Chisq)    
## (Intercept)   5.3969  1    0.02017 *  
## Env         157.1284  7  < 2.2e-16 ***
## LAR2          6.1172  1    0.01339 *  
## Env:LAR2     33.2727  7  2.355e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(prrep2)
## [1] 3094

## random effect test
prrep2b<-glm(didrep ~ Env*LAR2, data=adult, family=binomial(link="logit"))
anova(prrep2, prrep2b)
## Data: adult
## Models:
## prrep2b: didrep ~ Env * LAR2
## prrep2: didrep ~ Env * LAR2 + (1 | Blk)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## prrep2b   16 3281.1 3377.7 -1624.5   3249.1                         
## prrep2    17 3144.7 3247.3 -1555.3   3110.7 138.37  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect of damage on total fruit size:

library(car)

# Supplementary Table 2B
## filter field data to plants in 2016 Schofield Garden (best replication) that did reproduce
use.sch<-filter(adult, Env=="SCH-2016" & didrep==1 & !is.na(TotalRepro) & TotalRepro > 0)

## fit model
repamt2<-lmer(log(TotalRepro) ~ LAR2 + (1|Blk), data=use.sch)
summary(repamt2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TotalRepro) ~ LAR2 + (1 | Blk)
##    Data: use.sch
## 
## REML criterion at convergence: 1696.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9464 -0.5097  0.1398  0.6902  2.2962 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Blk      (Intercept) 0.03459  0.1860  
##  Residual             0.40868  0.6393  
## Number of obs: 848, groups:  Blk, 30
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.435306   0.044865  42.236673  54.281   <2e-16 ***
## LAR2         -0.007115   0.003025 845.345277  -2.352   0.0189 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## LAR2 -0.415
Anova(repamt2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(TotalRepro)
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 2946.4000  1    < 2e-16 ***
## LAR2           5.5318  1    0.01867 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(repamt2)
## [1] 848
ranova(repamt2)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log(TotalRepro) ~ LAR2 + (1 | Blk)
##           npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>       4 -848.04 1704.1                         
## (1 | Blk)    3 -864.20 1734.4 32.327  1  1.303e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Viability selection by herbivores

Effect of damage on survival:

library(car)

# Supplementary Table 4
s1<-glmer(X2017_Surv0 ~ X2016_LAR2 + (1|Blk), data=surv, family=binomial(link="logit"))
summary(s1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: X2017_Surv0 ~ X2016_LAR2 + (1 | Blk)
##    Data: surv
## 
##      AIC      BIC   logLik deviance df.resid 
##    686.5    701.8   -340.2    680.5     1216 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0877  0.1869  0.2483  0.2983  0.9724 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Blk    (Intercept) 0.8843   0.9404  
## Number of obs: 1219, groups:  Blk, 30
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.76623    0.23575  11.734   <2e-16 ***
## X2016_LAR2  -0.01366    0.01199  -1.139    0.255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## X2016_LAR2 -0.385
Anova(s1, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: X2017_Surv0
##                Chisq Df Pr(>Chisq)    
## (Intercept) 137.6811  1     <2e-16 ***
## X2016_LAR2    1.2975  1     0.2547    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s1)
## [1] 1219

## random effect test
s1b<-glm(X2017_Surv0 ~ X2016_LAR2, data=surv, family=binomial(link="logit"))
anova(s1, s1b)
## Data: surv
## Models:
## s1b: X2017_Surv0 ~ X2016_LAR2
## s1: X2017_Surv0 ~ X2016_LAR2 + (1 | Blk)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## s1b    2 737.21 747.43 -366.61   733.21                         
## s1     3 686.49 701.81 -340.25   680.49 52.721  1  3.845e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Variation in survival across environments

Temporary arrays:

library(car)

# Supplementary Table 3A
## define array dataset
a<-use.a

## fit full model
s2<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Blk) + (1|Fam), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 
summary(s2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category +  
##     (1 | Blk) + (1 | Fam)
##    Data: a
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##   3246.8   3367.0  -1605.4   3210.8     5862 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -11.5696   0.0990   0.1702   0.3416   2.3645 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Blk    (Intercept) 0.8219   0.9066  
##  Fam    (Intercept) 0.1814   0.4259  
## Number of obs: 5880, groups:  Blk, 60; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       1.11661    0.43191   2.585 0.009730 ** 
## SiteKC            1.10565    0.44723   2.472 0.013427 *  
## SiteLA            2.62562    0.48901   5.369 7.91e-08 ***
## SiteMM            0.32161    0.43655   0.737 0.461299    
## SiteRM            4.17569    0.64793   6.445 1.16e-10 ***
## SiteUA            3.72388    0.58413   6.375 1.83e-10 ***
## GenoSS            0.82496    0.35940   2.295 0.021709 *  
## CategoryL        -0.12367    0.37978  -0.326 0.744701    
## CategoryM        -0.43696    0.35515  -1.230 0.218566    
## SiteKC:GenoSS    -0.03145    0.28434  -0.111 0.911929    
## SiteLA:GenoSS    -1.12606    0.35326  -3.188 0.001434 ** 
## SiteMM:GenoSS    -0.42347    0.24258  -1.746 0.080858 .  
## SiteRM:GenoSS    -2.18724    0.57477  -3.805 0.000142 ***
## SiteUA:GenoSS    -1.19828    0.53969  -2.220 0.026398 *  
## GenoSS:CategoryL  0.16940    0.26739   0.634 0.526381    
## GenoSS:CategoryM  0.15592    0.23746   0.657 0.511431    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s2)
## [1] 5880

## Wald test
Anova(s2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Surv2
##                 Chisq Df Pr(>Chisq)    
## (Intercept)    6.6837  1    0.00973 ** 
## Site          89.2772  5  < 2.2e-16 ***
## Geno           5.2690  1    0.02171 *  
## Category       1.7201  2    0.42313    
## Site:Geno     26.6144  5   6.78e-05 ***
## Geno:Category  0.5318  2    0.76653    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## random effects tests
s2b<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Fam), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

s2c<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + Geno*Category + (1|Blk), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

anova(s2, s2b)
## Data: a
## Models:
## s2b: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category + 
## s2b:     (1 | Fam)
## s2: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category + 
## s2:     (1 | Blk) + (1 | Fam)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## s2b   17 3510.2 3623.7 -1738.1   3476.2                         
## s2    18 3246.8 3367.0 -1605.4   3210.8 265.42  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(s2, s2c)
## Data: a
## Models:
## s2c: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category + 
## s2c:     (1 | Blk)
## s2: Surv2 ~ Site + Geno + Category + Geno * Site + Geno * Category + 
## s2:     (1 | Blk) + (1 | Fam)
##     npar    AIC  BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## s2c   17 3308.5 3422 -1637.2   3274.5                         
## s2    18 3246.8 3367 -1605.4   3210.8 63.704  1  1.446e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## parametric bootstrap
    # parametric bootstrap#
    # method is very slow,#
    # i.e. hours, so they #
    # are commented out. ##
    # delete #s to run it #

### define reduced models eliminating effects

#### drop geno*site interaction
s2.gs<-glmer(Surv2 ~ Site + Geno + Category + Geno*Category + (1|Blk) + (1|Fam), data=a,
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop genotype (main effect + interactions)
s2.g<-glmer(Surv2 ~ Site + Category + (1|Blk) + (1|Fam), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop site (main effect + interactions)
s2.s<-glmer(Surv2 ~ Geno + Category + Geno*Category + (1|Blk) + (1|Fam), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop geno*category interaction
s2.gc<-glmer(Surv2 ~ Site + Geno + Category + Geno*Site + (1|Blk) + (1|Fam), data=a,
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop category (main effect + interactions)
s2.c<-glmer(Surv2 ~ Site + Geno + Geno*Site +(1|Blk) + (1|Fam), data=a, 
          family=binomial(link="logit"),  
          control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

### run bootstrapping in parallel
library(parallel)
nc<-detectCores()
cl<-makeCluster(rep("localhost", nc))

#### test geno*site
# s2.pb.gs<-PBmodcomp(s2, s2.gs, cl=cl)
# s2.pb.gs

#### test geno
# s2.pb.g<-PBmodcomp(s2, s2.g, cl=cl)
# s2.pb.g

#### test site
# s2.pb.s<-PBmodcomp(s2, s2.s, cl=cl)
# s2.pb.s

#### test category
# s2.pb.c<-PBmodcomp(s2, s2.c, cl=cl)
# s2.pb.c

#### test geno*category
# s2.pb.gc<-PBmodcomp(s2, s2.gc, cl=cl)
# s2.pb.gc

### turn off cluster
stopCluster(cl)

# Figure 2C (left pane)
## extract least squares means
emmeans(s2, pairwise ~ Geno* Site, type="response")$emmeans
##  Geno Site  prob      SE  df asymp.LCL asymp.UCL
##  LL   401  0.717 0.07394 Inf     0.554     0.838
##  SS   401  0.866 0.04337 Inf     0.756     0.930
##  LL   KC   0.884 0.03848 Inf     0.785     0.941
##  SS   KC   0.950 0.01891 Inf     0.897     0.976
##  LL   LA   0.972 0.01152 Inf     0.938     0.988
##  SS   LA   0.967 0.01343 Inf     0.928     0.985
##  LL   MM   0.778 0.06282 Inf     0.632     0.877
##  SS   MM   0.853 0.04616 Inf     0.739     0.923
##  LL   RM   0.994 0.00361 Inf     0.981     0.998
##  SS   RM   0.979 0.00905 Inf     0.952     0.991
##  LL   UA   0.991 0.00497 Inf     0.974     0.997
##  SS   UA   0.988 0.00590 Inf     0.969     0.995
## 
## Results are averaged over the levels of: Category 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale

## pairwise comparisons within environments
emmeans(s2, pairwise ~ Geno | Site, type="response")$contrasts
## Site = 401:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.393 0.127 Inf -2.896  0.0038 
## 
## Site = KC:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.406 0.143 Inf -2.567  0.0103 
## 
## Site = LA:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.212 0.497 Inf  0.470  0.6384 
## 
## Site = MM:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.601 0.192 Inf -1.592  0.1115 
## 
## Site = RM:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       3.504 2.147 Inf  2.046  0.0407 
## 
## Site = UA:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.303 0.755 Inf  0.457  0.6474 
## 
## Results are averaged over the levels of: Category 
## Tests are performed on the log odds ratio scale

Transplant experiments:

library(car)

# Supplementary Table 3B
## define transplant dataset
t<-use.t

## fit full model
s3<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 
summary(s3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
##    Data: use.t
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
##      AIC      BIC   logLik deviance df.resid 
##   7134.1   7270.7  -3547.0   7094.1     6840 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7314 -0.6663  0.2340  0.6305  4.9096 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Blk    (Intercept) 1.0591   1.0291  
##  Fam    (Intercept) 0.0546   0.2337  
## Number of obs: 6860, groups:  Blk, 150; Fam, 19
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.63717    0.30195   2.110 0.034847 *  
## EnvALD-2014         1.60151    0.53629   2.986 0.002824 ** 
## EnvALD-2015        -1.16282    0.39404  -2.951 0.003168 ** 
## EnvGTH-2016        -2.70750    0.37217  -7.275 3.47e-13 ***
## EnvMAH-2014        -0.94173    0.43968  -2.142 0.032206 *  
## EnvMAH-2015        -0.01805    0.44683  -0.040 0.967778    
## EnvSCH-2016         1.24843    0.37153   3.360 0.000779 ***
## EnvSIL-2014        -0.15597    0.45916  -0.340 0.734089    
## EnvSIL-2015         0.61895    0.42558   1.454 0.145847    
## GenoSS             -0.26820    0.20970  -1.279 0.200893    
## EnvALD-2014:GenoSS -0.07596    0.35969  -0.211 0.832754    
## EnvALD-2015:GenoSS  0.34721    0.23224   1.495 0.134895    
## EnvGTH-2016:GenoSS  1.05173    0.23880   4.404 1.06e-05 ***
## EnvMAH-2014:GenoSS  0.10735    0.25520   0.421 0.674017    
## EnvMAH-2015:GenoSS  0.50360    0.28238   1.783 0.074516 .  
## EnvSCH-2016:GenoSS  0.75182    0.24825   3.028 0.002458 ** 
## EnvSIL-2014:GenoSS -0.05341    0.26036  -0.205 0.837475    
## EnvSIL-2015:GenoSS -0.68525    0.25734  -2.663 0.007749 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(s3)
## [1] 6860

## random effects tests
s3b<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5)))
s3c<-glmer(Surv2 ~ Env + Geno + Geno*Env + (1|Blk), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5)))

anova(s3, s3b)
## Data: use.t
## Models:
## s3b: Surv2 ~ Env + Geno + Geno * Env + (1 | Fam)
## s3: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## s3b   19 7847.7 7977.5 -3904.8   7809.7                         
## s3    20 7134.1 7270.7 -3547.0   7094.1 715.59  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(s3, s3c)
## Data: use.t
## Models:
## s3c: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk)
## s3: Surv2 ~ Env + Geno + Geno * Env + (1 | Blk) + (1 | Fam)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## s3c   19 7159.1 7288.9 -3560.5   7121.1                         
## s3    20 7134.1 7270.7 -3547.0   7094.1 27.033  1      2e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Wald test
Anova(s3, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Surv2
##                Chisq Df Pr(>Chisq)    
## (Intercept)   4.4527  1    0.03485 *  
## Env         198.1527  8  < 2.2e-16 ***
## Geno          1.6359  1    0.20089    
## Env:Geno     61.4704  8  2.397e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## parametric bootstrap
    # parametric bootstrap#
    # method is very slow,#
    # i.e. hours, so they #
    # are commented out. ##
    # delete #s to run it #

### define reduced models eliminating effects

#### drop geno*site interaction
s3.gs<-glmer(Surv2 ~ Env + Geno + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop genotype (main effect + interaction)
s3.g<-glmer(Surv2 ~ Env + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

#### drop site (main effect + interaction)
s3.s<-glmer(Surv2 ~ Geno + (1|Blk) + (1|Fam), data=use.t, family=binomial(link="logit"), control=glmerControl(optimizer = "bobyqa", 
                               optCtrl = list(maxfun=2e5))) 

### run bootstrapping in parallel
nc<-detectCores()
cl<-makeCluster(rep("localhost", nc))

#### test geno*site
# s3.pb.gs<-PBmodcomp(s3, s3.gs, cl=cl)
# s3.pb.gs

#### test effect of geno
# s3.pb.g<-PBmodcomp(s3, s3.g, cl=cl)
# s3.pb.g

#### test effect of site
# s3.pb.s<-PBmodcomp(s3, s3.s, cl=cl)
# s3.pb.s

### turn off cluster
stopCluster(cl)

# Figure 2C (right pane)
## extract least squares means
emmeans(s3, pairwise ~ Geno* Env, type="response")$emmeans
##  Geno Env       prob     SE  df asymp.LCL asymp.UCL
##  LL   ALD-2013 0.654 0.0683 Inf    0.5113     0.774
##  SS   ALD-2013 0.591 0.0709 Inf    0.4488     0.720
##  LL   ALD-2014 0.904 0.0404 Inf    0.7907     0.959
##  SS   ALD-2014 0.869 0.0495 Inf    0.7391     0.940
##  LL   ALD-2015 0.372 0.0645 Inf    0.2559     0.504
##  SS   ALD-2015 0.390 0.0650 Inf    0.2726     0.522
##  LL   GTH-2016 0.112 0.0246 Inf    0.0721     0.170
##  SS   GTH-2016 0.216 0.0388 Inf    0.1499     0.302
##  LL   MAH-2014 0.424 0.0848 Inf    0.2719     0.593
##  SS   MAH-2014 0.386 0.0804 Inf    0.2441     0.550
##  LL   MAH-2015 0.650 0.0790 Inf    0.4846     0.786
##  SS   MAH-2015 0.702 0.0722 Inf    0.5446     0.822
##  LL   SCH-2016 0.868 0.0279 Inf    0.8033     0.914
##  SS   SCH-2016 0.914 0.0193 Inf    0.8682     0.945
##  LL   SIL-2014 0.618 0.0877 Inf    0.4385     0.770
##  SS   SIL-2014 0.540 0.0899 Inf    0.3660     0.704
##  LL   SIL-2015 0.778 0.0551 Inf    0.6525     0.868
##  SS   SIL-2015 0.575 0.0754 Inf    0.4250     0.712
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale

## pairwise contrasts within environments
emmeans(s3, pairwise ~ Geno | Env, type="response")$contrasts
## Env = ALD-2013:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.308 0.274 Inf  1.279  0.2009 
## 
## Env = ALD-2014:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.411 0.485 Inf  1.000  0.3172 
## 
## Env = ALD-2015:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.924 0.168 Inf -0.434  0.6643 
## 
## Env = GTH-2016:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.457 0.092 Inf -3.889  0.0001 
## 
## Env = MAH-2014:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.175 0.273 Inf  0.692  0.4889 
## 
## Env = MAH-2015:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.790 0.191 Inf -0.972  0.3311 
## 
## Env = SCH-2016:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       0.617 0.130 Inf -2.293  0.0219 
## 
## Env = SIL-2014:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       1.379 0.328 Inf  1.351  0.1767 
## 
## Env = SIL-2015:
##  contrast odds.ratio    SE  df z.ratio p.value
##  LL / SS       2.595 0.551 Inf  4.491  <.0001 
## 
## Tests are performed on the log odds ratio scale

Drought reponse

Population-level response in the field

These analyses (Supplementary Table 5) were conducted using standard least squares linear models in JMP.

Experimental dry-down with CFR-HIFs

Genetic variation in morphology in response to drought:

library(car)
dry$Rack<-as.character(dry$Rack)
dry$Fam<-as.character(dry$Fam)

# Supplementary Table 6

## separate statistical models for each morphological response

### leaf length 2 (final height)

len<-lmer(LL_Length2_mm ~ Geno + Trmt + Geno*Trmt + LL1_mm + (1|Rack) + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(len)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LL_Length2_mm ~ Geno + Trmt + Geno * Trmt + LL1_mm + (1 | Rack) +  
##     (1 | Fam)
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 7363.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5996 -0.4916  0.0870  0.6153  3.5471 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept)  23.040   4.80   
##  Fam      (Intercept)   6.551   2.56   
##  Residual             122.681  11.08   
## Number of obs: 954, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     28.96138    1.89243  36.53828  15.304  < 2e-16 ***
## GenoMM          -1.70151    1.91204  10.66410  -0.890  0.39316    
## TrmtWet         -3.15346    1.04385 922.77176  -3.021  0.00259 ** 
## LL1_mm           0.53371    0.03178 926.63889  16.794  < 2e-16 ***
## GenoMM:TrmtWet   7.73979    1.48079 917.95306   5.227 2.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.499                     
## TrmtWet     -0.181  0.261              
## LL1_mm      -0.537 -0.012 -0.151       
## GnMM:TrmtWt  0.083 -0.367 -0.711  0.189
Anova(len, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LL_Length2_mm
##                Chisq Df Pr(>Chisq)    
## (Intercept) 234.2056  1  < 2.2e-16 ***
## Geno          0.7919  1   0.373525    
## Trmt          9.1264  1   0.002519 ** 
## LL1_mm      282.0298  1  < 2.2e-16 ***
## Geno:Trmt    27.3194  1  1.725e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(len)
## [1] 954
ranova(len)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LL_Length2_mm ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + 
##     Geno:Trmt
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        8 -3681.6 7379.3                         
## (1 | Rack)    7 -3727.6 7469.2 91.868  1  < 2.2e-16 ***
## (1 | Fam)     7 -3693.5 7400.9 23.615  1  1.177e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### total leaves
tlf<-lmer(TotalLvs ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(tlf)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TotalLvs ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +  
##     LL1_mm
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 4667.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5549 -0.5734 -0.0095  0.4870 12.1697 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.38291  0.6188  
##  Fam      (Intercept) 0.04476  0.2116  
##  Residual             7.32034  2.7056  
## Number of obs: 958, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     10.384840   0.326204  98.442935  31.835   <2e-16 ***
## GenoMM          -0.272180   0.281167  19.964449  -0.968   0.3446    
## TrmtWet          0.505336   0.252187 941.935067   2.004   0.0454 *  
## LL1_mm           0.150553   0.007403 719.815659  20.338   <2e-16 ***
## GenoMM:TrmtWet   1.048174   0.359182 937.492029   2.918   0.0036 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.409                     
## TrmtWet     -0.269  0.433              
## LL1_mm      -0.719 -0.029 -0.142       
## GnMM:TrmtWt  0.130 -0.610 -0.710  0.182
Anova(tlf, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: TotalLvs
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 1013.4948  1    < 2e-16 ***
## Geno           0.9371  1    0.33303    
## Trmt           4.0153  1    0.04509 *  
## LL1_mm       413.6266  1    < 2e-16 ***
## Geno:Trmt      8.5160  1    0.00352 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(tlf)
## [1] 958
ranova(tlf)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## TotalLvs ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
##            npar  logLik    AIC     LRT Df Pr(>Chisq)    
## <none>        8 -2333.9 4683.9                          
## (1 | Rack)    7 -2342.7 4699.3 17.4105  1  3.012e-05 ***
## (1 | Fam)     7 -2334.4 4682.8  0.9143  1      0.339    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### LMA
lma<-lmer(LMAest ~ Geno + Trmt + Geno*Trmt + (1|Rack)+ (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(lma)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LMAest ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +  
##     LL1_mm
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: -3393.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1426 -0.3309 -0.1003  0.1852 13.6523 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Rack     (Intercept) 4.061e-06 0.002015
##  Fam      (Intercept) 1.103e-06 0.001050
##  Residual             1.300e-04 0.011401
## Number of obs: 568, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     1.979e-02  1.559e-03  6.274e+01  12.694  < 2e-16 ***
## GenoMM         -8.462e-04  1.244e-03  9.808e+00  -0.680  0.51226    
## TrmtWet         3.215e-03  1.852e-03  4.252e+02   1.736  0.08327 .  
## LL1_mm          1.053e-04  3.876e-05  4.256e+02   2.716  0.00688 ** 
## GenoMM:TrmtWet  3.319e-03  2.704e-03  5.472e+02   1.227  0.22025    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.366                     
## TrmtWet     -0.104  0.247              
## LL1_mm      -0.794 -0.041 -0.109       
## GnMM:TrmtWt  0.036 -0.334 -0.639  0.122
Anova(lma, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LMAest
##                Chisq Df Pr(>Chisq)    
## (Intercept) 161.1296  1  < 2.2e-16 ***
## Geno          0.4624  1   0.496526    
## Trmt          3.0141  1   0.082545 .  
## LL1_mm        7.3763  1   0.006609 ** 
## Geno:Trmt     1.5062  1   0.219726    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(lma)
## [1] 568
ranova(lma)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LMAest ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
##            npar logLik     AIC    LRT Df Pr(>Chisq)  
## <none>        8 1697.0 -3377.9                       
## (1 | Rack)    7 1695.3 -3376.7 3.1922  1    0.07399 .
## (1 | Fam)     7 1696.6 -3379.2 0.6450  1    0.42191  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### leaf width 
wid<-lmer(LL_Width2_mm ~ Geno + Trmt + Geno*Trmt + LL1_mm + (1|Rack) + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(wid)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LL_Width2_mm ~ Geno + Trmt + Geno * Trmt + LL1_mm + (1 | Rack) +  
##     (1 | Fam)
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 3762.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8023 -0.5186  0.0553  0.6007  3.0446 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.5363   0.7323  
##  Fam      (Intercept) 0.1952   0.4418  
##  Residual             2.7882   1.6698  
## Number of obs: 951, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      6.339402   0.302091  30.180346  20.985  < 2e-16 ***
## GenoMM          -0.250119   0.318963  10.134307  -0.784   0.4509    
## TrmtWet         -0.390530   0.157776 918.670530  -2.475   0.0135 *  
## LL1_mm           0.061854   0.004812 928.831192  12.855  < 2e-16 ***
## GenoMM:TrmtWet   1.128147   0.223651 913.548309   5.044 5.49e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.522                     
## TrmtWet     -0.171  0.237              
## LL1_mm      -0.508 -0.012 -0.154       
## GnMM:TrmtWt  0.079 -0.333 -0.712  0.189
Anova(wid, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LL_Width2_mm
##                Chisq Df Pr(>Chisq)    
## (Intercept) 440.3730  1  < 2.2e-16 ***
## Geno          0.6149  1    0.43294    
## Trmt          6.1267  1    0.01332 *  
## LL1_mm      165.2502  1  < 2.2e-16 ***
## Geno:Trmt    25.4443  1  4.553e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wid)
## [1] 951
ranova(wid)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LL_Width2_mm ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + 
##     Geno:Trmt
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        8 -1881.4 3778.9                         
## (1 | Rack)    7 -1923.5 3860.9 84.066  1  < 2.2e-16 ***
## (1 | Fam)     7 -1898.2 3810.4 33.562  1  6.904e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### growth rate
gr2<-lmer(GrowthRt ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(gr2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthRt ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +  
##     LL1_mm
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 1877.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5996 -0.4916  0.0870  0.6153  3.5471 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.07111  0.2667  
##  Fam      (Intercept) 0.02022  0.1422  
##  Residual             0.37865  0.6153  
## Number of obs: 954, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      1.608965   0.105135  36.538313  15.304  < 2e-16 ***
## GenoMM          -0.094528   0.106225  10.664104  -0.890  0.39316    
## TrmtWet         -0.175192   0.057991 922.771759  -3.021  0.00259 ** 
## LL1_mm          -0.025905   0.001766 926.638887 -14.672  < 2e-16 ***
## GenoMM:TrmtWet   0.429988   0.082266 917.953056   5.227 2.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.499                     
## TrmtWet     -0.181  0.261              
## LL1_mm      -0.537 -0.012 -0.151       
## GnMM:TrmtWt  0.083 -0.367 -0.711  0.189
Anova(gr2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: GrowthRt
##                Chisq Df Pr(>Chisq)    
## (Intercept) 234.2057  1  < 2.2e-16 ***
## Geno          0.7919  1   0.373525    
## Trmt          9.1264  1   0.002519 ** 
## LL1_mm      215.2767  1  < 2.2e-16 ***
## Geno:Trmt    27.3194  1  1.725e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gr2)
## [1] 954
ranova(gr2)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## GrowthRt ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        8 -938.68 1893.4                         
## (1 | Rack)    7 -984.62 1983.2 91.868  1  < 2.2e-16 ***
## (1 | Fam)     7 -950.49 1915.0 23.615  1  1.177e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### green leaves
glf<-lmer(NumLvs_Grn ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(glf)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NumLvs_Grn ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +  
##     LL1_mm
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 4437.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6275 -0.5629 -0.0521  0.5536  4.5457 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.20639  0.4543  
##  Fam      (Intercept) 0.08516  0.2918  
##  Residual             5.74036  2.3959  
## Number of obs: 959, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      9.771777   0.300931  61.771991  32.472  < 2e-16 ***
## GenoMM          -0.920397   0.285998  15.092490  -3.218  0.00571 ** 
## TrmtWet         -0.485014   0.222697 943.990125  -2.178  0.02966 *  
## LL1_mm           0.047072   0.006592 841.281137   7.141 2.01e-12 ***
## GenoMM:TrmtWet   1.516285   0.317653 939.977577   4.773 2.10e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.455                     
## TrmtWet     -0.253  0.376              
## LL1_mm      -0.692 -0.027 -0.146       
## GnMM:TrmtWt  0.120 -0.529 -0.710  0.186
Anova(glf, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: NumLvs_Grn
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 1054.4218  1  < 2.2e-16 ***
## Geno          10.3568  1    0.00129 ** 
## Trmt           4.7433  1    0.02941 *  
## LL1_mm        50.9936  1  9.267e-13 ***
## Geno:Trmt     22.7854  1  1.811e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(glf)
## [1] 959
ranova(glf)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## NumLvs_Grn ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + 
##     Geno:Trmt
##            npar  logLik    AIC     LRT Df Pr(>Chisq)   
## <none>        8 -2218.8 4453.7                         
## (1 | Rack)    7 -2224.0 4461.9 10.2554  1   0.001363 **
## (1 | Fam)     7 -2220.8 4455.5  3.8712  1   0.049120 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### LWC
lwc<-lmer(LWC_Ros ~ Geno + Trmt + Geno*Trmt + (1|Rack) + (1|Fam) + LL1_mm, data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(lwc)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LWC_Ros ~ Geno + Trmt + Geno * Trmt + (1 | Rack) + (1 | Fam) +  
##     LL1_mm
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: -1275
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.4728 -0.4046  0.0925  0.5552  3.0604 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  Rack     (Intercept) 3.073e-03 5.543e-02
##  Fam      (Intercept) 1.002e-18 1.001e-09
##  Residual             5.077e-03 7.125e-02
## Number of obs: 569, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     6.049e-01  1.347e-02  8.356e+01  44.911   <2e-16 ***
## GenoMM         -6.411e-03  6.663e-03  5.372e+02  -0.962    0.336    
## TrmtWet         1.220e-01  1.246e-02  5.527e+02   9.796   <2e-16 ***
## LL1_mm         -2.655e-03  2.478e-04  5.467e+02 -10.714   <2e-16 ***
## GenoMM:TrmtWet  3.015e-02  1.702e-02  5.344e+02   1.772    0.077 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TrmtWt LL1_mm
## GenoMM      -0.220                     
## TrmtWet     -0.084  0.278              
## LL1_mm      -0.591 -0.043 -0.075       
## GnMM:TrmtWt  0.027 -0.397 -0.608  0.120
## convergence code: 0
## boundary (singular) fit: see ?isSingular
Anova(lwc, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LWC_Ros
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 2017.0092  1    < 2e-16 ***
## Geno           0.9257  1    0.33598    
## Trmt          95.9684  1    < 2e-16 ***
## LL1_mm       114.7917  1    < 2e-16 ***
## Geno:Trmt      3.1395  1    0.07642 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(lwc)
## [1] 569
ranova(lwc)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LWC_Ros ~ Geno + Trmt + LL1_mm + (1 | Rack) + (1 | Fam) + Geno:Trmt
##            npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>        8 637.50 -1259.0                         
## (1 | Rack)    7 565.18 -1116.3 144.66  1     <2e-16 ***
## (1 | Fam)     7 637.50 -1261.0   0.00  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Morphological predictors of water use under drought:

library(car)

# Supplementary Table 7

## Supp Table 7A: size predictors
wu<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LL_Length2_mm + Geno*LL_Length2_mm + LL1_mm  + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5))) 
summary(wu) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LL_Length2_mm + Geno *  
##     LL_Length2_mm + LL1_mm + (1 | Fam)
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 467.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7835 -0.6379 -0.0872  0.5889  4.0839 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.051055 0.2260  
##  Fam      (Intercept) 0.009389 0.0969  
##  Residual             0.121463 0.3485  
## Number of obs: 486, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)            5.681253   0.105799 116.568927  53.699  < 2e-16 ***
## GenoMM                -0.039813   0.129640  78.073446  -0.307 0.759585    
## LL_Length2_mm          0.007380   0.002061 471.062504   3.582 0.000377 ***
## LL1_mm                 0.012136   0.001668 458.176040   7.274 1.52e-12 ***
## GenoMM:LL_Length2_mm  -0.001726   0.002429 455.324507  -0.711 0.477632    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM LL_L2_ LL1_mm
## GenoMM      -0.646                     
## LL_Lngth2_m -0.665  0.537              
## LL1_mm      -0.091 -0.018 -0.460       
## GnMM:LL_L2_  0.549 -0.845 -0.615  0.001
Anova(wu, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LSM_daily_norack
##                        Chisq Df Pr(>Chisq)    
## (Intercept)        2883.5450  1  < 2.2e-16 ***
## Geno                  0.0943  1  0.7587674    
## LL_Length2_mm        12.8290  1  0.0003413 ***
## LL1_mm               52.9137  1  3.485e-13 ***
## Geno:LL_Length2_mm    0.5051  1  0.4772678    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu)
## [1] 486
ranova(wu)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LSM_daily_norack ~ Geno + LL_Length2_mm + LL1_mm + (1 | Rack) + 
##     (1 | Fam) + Geno:LL_Length2_mm
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        8 -233.57 483.15                         
## (1 | Rack)    7 -274.35 562.71 81.558  1  < 2.2e-16 ***
## (1 | Fam)     7 -240.04 494.08 12.934  1  0.0003227 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Supp Table 7B: leaf number
wu2<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + TotalLvs + Geno*TotalLvs + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5))) 
summary(wu2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + TotalLvs + Geno * TotalLvs +  
##     (1 | Fam)
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 521.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0183 -0.6387 -0.0592  0.5886  3.2906 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.05144  0.2268  
##  Fam      (Intercept) 0.01954  0.1398  
##  Residual             0.14007  0.3743  
## Number of obs: 487, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       5.677153   0.135205 106.297637  41.989  < 2e-16 ***
## GenoMM           -0.133408   0.198389 115.703617  -0.672    0.503    
## TotalLvs          0.048018   0.007232 462.263441   6.640 8.81e-11 ***
## GenoMM:TotalLvs   0.001593   0.011543 457.947492   0.138    0.890    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM TtlLvs
## GenoMM      -0.594              
## TotalLvs    -0.816  0.522       
## GnMM:TtlLvs  0.480 -0.878 -0.588
Anova(wu2, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LSM_daily_norack
##                   Chisq Df Pr(>Chisq)    
## (Intercept)   1763.0923  1  < 2.2e-16 ***
## Geno             0.4522  1     0.5013    
## TotalLvs        44.0913  1  3.134e-11 ***
## Geno:TotalLvs    0.0190  1     0.8902    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu2)
## [1] 487
ranova(wu2)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LSM_daily_norack ~ Geno + TotalLvs + (1 | Rack) + (1 | Fam) + 
##     Geno:TotalLvs
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        7 -260.68 535.37                         
## (1 | Rack)    6 -296.46 604.93 71.564  1  < 2.2e-16 ***
## (1 | Fam)     6 -275.65 563.30 29.934  1   4.47e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Supp Table 7C: LMA
### exclude LMA extreme outliers
dry.use<-filter(dry, LMAest<.1)

### fit model
wu3<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LMAest + Geno*LMAest  + (1|Fam), data=dry.use, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(wu3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LMAest + Geno * LMAest +  
##     (1 | Fam)
##    Data: dry.use
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 523.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6013 -0.6387 -0.0511  0.5754  3.3613 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.06547  0.2559  
##  Fam      (Intercept) 0.02452  0.1566  
##  Residual             0.14977  0.3870  
## Number of obs: 477, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     5.94063    0.13275  70.79656  44.752  < 2e-16 ***
## GenoMM          0.05163    0.18322  65.84865   0.282    0.779    
## LMAest         21.59883    4.55808 445.24385   4.739  2.9e-06 ***
## GenoMM:LMAest  -8.08351    6.74018 438.38550  -1.199    0.231    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM LMAest
## GenoMM      -0.615              
## LMAest      -0.754  0.511       
## GenMM:LMAst  0.478 -0.817 -0.633
Anova(wu3, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LSM_daily_norack
##                 Chisq Df Pr(>Chisq)    
## (Intercept) 2002.7414  1  < 2.2e-16 ***
## Geno           0.0794  1     0.7781    
## LMAest        22.4542  1  2.152e-06 ***
## Geno:LMAest    1.4383  1     0.2304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu3)
## [1] 477
ranova(wu3)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LSM_daily_norack ~ Geno + LMAest + (1 | Rack) + (1 | Fam) + Geno:LMAest
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        7 -261.95 537.90                         
## (1 | Rack)    6 -300.86 613.71 77.804  1  < 2.2e-16 ***
## (1 | Fam)     6 -279.55 571.10 35.192  1  2.987e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Supp Table 7D: all predictors together
wu4<-lmer(LSM_daily_norack ~ Geno + (1|Rack) + LL_Length2_mm + TotalLvs + LMAest + LL1_mm  + (1|Fam), data=dry, control=lmerControl(optimizer = "bobyqa",
                                  optCtrl = list(maxfun=2e5)))
summary(wu4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LSM_daily_norack ~ Geno + (1 | Rack) + LL_Length2_mm + TotalLvs +  
##     LMAest + LL1_mm + (1 | Fam)
##    Data: dry
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
## 
## REML criterion at convergence: 450.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8906 -0.6486 -0.0963  0.5863  4.3035 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Rack     (Intercept) 0.050236 0.22413 
##  Fam      (Intercept) 0.009146 0.09563 
##  Residual             0.120494 0.34712 
## Number of obs: 478, groups:  Rack, 32; Fam, 10
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)     5.527581   0.121694 174.944551  45.422  < 2e-16 ***
## GenoMM         -0.116958   0.068845   7.616420  -1.699    0.130    
## LL_Length2_mm   0.007819   0.001787 458.644204   4.375 1.50e-05 ***
## TotalLvs        0.007588   0.007518 458.892479   1.009    0.313    
## LMAest          3.494983   1.770551 449.056322   1.974    0.049 *  
## LL1_mm          0.010221   0.002008 451.196513   5.091 5.24e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) GenoMM LL_L2_ TtlLvs LMAest
## GenoMM      -0.322                            
## LL_Lngth2_m -0.393  0.039                     
## TotalLvs    -0.551  0.033 -0.165              
## LMAest      -0.300  0.036  0.350 -0.150       
## LL1_mm       0.296 -0.056 -0.467 -0.445 -0.255
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
Anova(wu4, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: LSM_daily_norack
##                   Chisq Df Pr(>Chisq)    
## (Intercept)   2063.1495  1  < 2.2e-16 ***
## Geno             2.8861  1    0.08935 .  
## LL_Length2_mm   19.1444  1  1.212e-05 ***
## TotalLvs         1.0186  1    0.31285    
## LMAest           3.8965  1    0.04839 *  
## LL1_mm          25.9175  1  3.563e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(wu4)
## [1] 478
ranova(wu4)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LSM_daily_norack ~ Geno + LL_Length2_mm + TotalLvs + LMAest + 
##     LL1_mm + (1 | Rack) + (1 | Fam)
##            npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>        9 -225.15 468.30                         
## (1 | Rack)    8 -264.98 545.96 79.661  1  < 2.2e-16 ***
## (1 | Fam)     8 -231.50 478.99 12.696  1  0.0003664 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Experimental dry-down with broad panel of accessions

Genetic correlation between LWC and BC-ratio:

# Supplementary Table 8A
gwas.drt.reg<-lm(LWC_Dry_LSM ~ BCratio + Group, data=rp)
summary(gwas.drt.reg)
## 
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + Group, data = rp)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0100583 -0.0030995 -0.0005596  0.0025897  0.0163209 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.237e-01  1.330e-03 469.139  < 2e-16 ***
## BCratio     -2.625e-05  9.181e-06  -2.858  0.00465 ** 
## GroupCOL     9.843e-04  1.386e-03   0.710  0.47848    
## GroupNOR     1.650e-03  1.332e-03   1.238  0.21685    
## GroupUTA     1.490e-03  1.314e-03   1.134  0.25815    
## GroupWES     2.611e-03  1.366e-03   1.912  0.05716 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004347 on 231 degrees of freedom
##   (325 observations deleted due to missingness)
## Multiple R-squared:  0.04238,    Adjusted R-squared:  0.02165 
## F-statistic: 2.044 on 5 and 231 DF,  p-value: 0.07337
Anova(gwas.drt.reg, type="3")
## Anova Table (Type III tests)
## 
## Response: LWC_Dry_LSM
##             Sum Sq  Df    F value    Pr(>F)    
## (Intercept) 4.1580   1 2.2009e+05 < 2.2e-16 ***
## BCratio     0.0002   1 8.1709e+00  0.004646 ** 
## Group       0.0001   4 1.2381e+00  0.295475    
## Residuals   0.0044 231                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg)
## [1] 237

# Supplementary Table 8B
gwas.drt.reg.2<-lm(LWC_Dry_LSM ~ BCratio + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(gwas.drt.reg.2)
## 
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + gPCA1 + gPCA2 + gPCA3 + 
##     gPCA4 + gPCA5, data = rp)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0099812 -0.0029282 -0.0005765  0.0027871  0.0171346 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.250e-01  6.275e-04 996.054   <2e-16 ***
## BCratio     -2.180e-05  8.412e-06  -2.591   0.0101 *  
## gPCA1       -1.794e-03  8.976e-03  -0.200   0.8418    
## gPCA2       -4.061e-03  5.621e-03  -0.722   0.4707    
## gPCA3       -8.980e-03  5.768e-03  -1.557   0.1208    
## gPCA4       -2.559e-03  6.387e-03  -0.401   0.6890    
## gPCA5       -5.023e-03  5.542e-03  -0.906   0.3656    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00432 on 248 degrees of freedom
##   (307 observations deleted due to missingness)
## Multiple R-squared:  0.03447,    Adjusted R-squared:  0.01111 
## F-statistic: 1.475 on 6 and 248 DF,  p-value: 0.187
Anova(gwas.drt.reg.2, type="3")
## Anova Table (Type III tests)
## 
## Response: LWC_Dry_LSM
##              Sum Sq  Df    F value  Pr(>F)    
## (Intercept) 18.5115   1 9.9212e+05 < 2e-16 ***
## BCratio      0.0001   1 6.7133e+00 0.01014 *  
## gPCA1        0.0000   1 3.9900e-02 0.84176    
## gPCA2        0.0000   1 5.2180e-01 0.47074    
## gPCA3        0.0000   1 2.4239e+00 0.12077    
## gPCA4        0.0000   1 1.6060e-01 0.68898    
## gPCA5        0.0000   1 8.2160e-01 0.36559    
## Residuals    0.0046 248                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg.2)
## [1] 255

# Supplementary Table 8C
gwas.drt.reg.3<-lm(LWC_Dry_LSM ~ BCratio + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(gwas.drt.reg.3)
## 
## Call:
## lm(formula = LWC_Dry_LSM ~ BCratio + colutaPCA1 + colutaPCA2 + 
##     colutaPCA3 + colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.009775 -0.003266 -0.000705  0.002878  0.014638 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.249e-01  8.644e-04 722.962   <2e-16 ***
## BCratio     -2.456e-05  1.277e-05  -1.924   0.0572 .  
## colutaPCA1   4.516e-03  5.258e-03   0.859   0.3924    
## colutaPCA2   1.020e-03  5.113e-03   0.200   0.8422    
## colutaPCA3   9.411e-03  5.049e-03   1.864   0.0651 .  
## colutaPCA4   3.561e-03  5.147e-03   0.692   0.4906    
## colutaPCA5  -5.994e-03  5.016e-03  -1.195   0.2349    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.004719 on 103 degrees of freedom
##   (452 observations deleted due to missingness)
## Multiple R-squared:  0.1116, Adjusted R-squared:  0.05985 
## F-statistic: 2.156 on 6 and 103 DF,  p-value: 0.0532
Anova(gwas.drt.reg.3, type="3")
## Anova Table (Type III tests)
## 
## Response: LWC_Dry_LSM
##              Sum Sq  Df    F value  Pr(>F)    
## (Intercept) 11.6405   1 5.2267e+05 < 2e-16 ***
## BCratio      0.0001   1 3.7000e+00 0.05717 .  
## colutaPCA1   0.0000   1 7.3760e-01 0.39242    
## colutaPCA2   0.0000   1 3.9800e-02 0.84222    
## colutaPCA3   0.0001   1 3.4751e+00 0.06515 .  
## colutaPCA4   0.0000   1 4.7860e-01 0.49061    
## colutaPCA5   0.0000   1 1.4276e+00 0.23490    
## Residuals    0.0023 103                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(gwas.drt.reg.3)
## [1] 110

Molecular signatures of selection

Phenotype-environment correlations

Effects of climate PCs on BC-ratio:

library(car)

# Supplementary Table 9A
clim.use<-filter(rp, !is.na(Group))
clim.pca<-lm(BCratio ~ Group + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=clim.use)
summary(clim.pca)
## 
## Call:
## lm(formula = BCratio ~ Group + Latitude + Longitude + Elev_m + 
##     Prin1 + Prin2 + Prin3 + Prin4 + Prin5, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -69.774 -23.979   9.857  23.473  54.352 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 767.37533  317.83140   2.414 0.016426 *  
## GroupCOL    -22.88750   13.80455  -1.658 0.098485 .  
## GroupNOR     22.43580   11.47151   1.956 0.051522 .  
## GroupUTA      2.69566    9.97427   0.270 0.787165    
## GroupWES     34.34359   12.21289   2.812 0.005283 ** 
## Latitude     -3.20089    2.81812  -1.136 0.257037    
## Longitude     3.88218    2.24420   1.730 0.084796 .  
## Elev_m       -0.05521    0.02316  -2.384 0.017829 *  
## Prin1         6.04188    1.55452   3.887 0.000128 ***
## Prin2         0.43460    1.67307   0.260 0.795247    
## Prin3        -2.20075    2.02478  -1.087 0.278047    
## Prin4        -4.37423    2.46018  -1.778 0.076528 .  
## Prin5        -1.38109    2.58700  -0.534 0.593879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.94 on 270 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.249,  Adjusted R-squared:  0.2156 
## F-statistic: 7.459 on 12 and 270 DF,  p-value: 7.317e-12
Anova(clim.pca, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)    
## (Intercept)   5579   1  5.8294 0.016426 *  
## Group        11399   4  2.9778 0.019735 *  
## Latitude      1235   1  1.2901 0.257037    
## Longitude     2864   1  2.9925 0.084796 .  
## Elev_m        5438   1  5.6823 0.017829 *  
## Prin1        14456   1 15.1060 0.000128 ***
## Prin2           65   1  0.0675 0.795247    
## Prin3         1131   1  1.1814 0.278047    
## Prin4         3025   1  3.1613 0.076528 .  
## Prin5          273   1  0.2850 0.593879    
## Residuals   258389 270                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca)
## [1] 283

# Supplementary Table 9B
clim.pca.2<-lm(BCratio ~ gPCA1+gPCA2+gPCA3+gPCA4+gPCA5 + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=rp)
summary(clim.pca.2)
## 
## Call:
## lm(formula = BCratio ~ gPCA1 + gPCA2 + gPCA3 + gPCA4 + gPCA5 + 
##     Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 + Prin4 + 
##     Prin5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -79.67 -22.80   9.02  22.39  58.01 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 280.22771  216.93295   1.292 0.197476    
## gPCA1       -68.89008   72.28922  -0.953 0.341403    
## gPCA2        -9.50791   56.25149  -0.169 0.865896    
## gPCA3       -36.25943   55.71415  -0.651 0.515689    
## gPCA4       -75.94926   46.63565  -1.629 0.104502    
## gPCA5       -32.33279   37.33405  -0.866 0.387192    
## Latitude     -0.02286    2.35447  -0.010 0.992261    
## Longitude     0.75566    1.66746   0.453 0.650762    
## Elev_m       -0.05055    0.02170  -2.330 0.020523 *  
## Prin1         5.41360    1.47403   3.673 0.000286 ***
## Prin2         1.96863    1.70342   1.156 0.248768    
## Prin3        -3.89549    1.77502  -2.195 0.028992 *  
## Prin4        -3.25136    2.39034  -1.360 0.174831    
## Prin5        -1.87265    2.72882  -0.686 0.493110    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.32 on 287 degrees of freedom
##   (261 observations deleted due to missingness)
## Multiple R-squared:  0.2199, Adjusted R-squared:  0.1846 
## F-statistic: 6.223 on 13 and 287 DF,  p-value: 2.65e-10
Anova(clim.pca.2, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   1636   1  1.6687 0.1974757    
## gPCA1          891   1  0.9082 0.3414028    
## gPCA2           28   1  0.0286 0.8658960    
## gPCA3          415   1  0.4236 0.5156886    
## gPCA4         2601   1  2.6522 0.1045020    
## gPCA5          736   1  0.7500 0.3871915    
## Latitude         0   1  0.0001 0.9922615    
## Longitude      201   1  0.2054 0.6507621    
## Elev_m        5322   1  5.4269 0.0205230 *  
## Prin1        13227   1 13.4883 0.0002863 ***
## Prin2         1310   1  1.3356 0.2487680    
## Prin3         4723   1  4.8164 0.0289919 *  
## Prin4         1814   1  1.8502 0.1748309    
## Prin5          462   1  0.4709 0.4931099    
## Residuals   281444 287                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca.2)
## [1] 301

# Supplementary Table 9C
clim.pca.3<-lm(BCratio ~ colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5 + Latitude + Longitude + Elev_m + Prin1 + Prin2 + Prin3 +Prin4 + Prin5, data=rp)
summary(clim.pca.3)
## 
## Call:
## lm(formula = BCratio ~ colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5 + Latitude + Longitude + Elev_m + 
##     Prin1 + Prin2 + Prin3 + Prin4 + Prin5, data = rp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.414 -27.121   4.959  24.451  63.798 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1056.9635   467.1186   2.263  0.02548 * 
## colutaPCA1   111.6806    80.2442   1.392  0.16661   
## colutaPCA2    84.9290    96.2540   0.882  0.37939   
## colutaPCA3   -75.0512    45.2579  -1.658  0.09991 . 
## colutaPCA4   -42.8002    46.0758  -0.929  0.35483   
## colutaPCA5    66.0869    47.0720   1.404  0.16296   
## Latitude      -6.8074     4.9414  -1.378  0.17093   
## Longitude      4.0521     3.9138   1.035  0.30262   
## Elev_m        -0.1031     0.0431  -2.392  0.01835 * 
## Prin1         10.0592     3.7850   2.658  0.00896 **
## Prin2         -0.2182     3.2787  -0.067  0.94704   
## Prin3         -6.9895     4.1252  -1.694  0.09284 . 
## Prin4        -10.1140     5.1457  -1.966  0.05170 . 
## Prin5         -4.4414     5.6652  -0.784  0.43462   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.74 on 118 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.3266, Adjusted R-squared:  0.2525 
## F-statistic: 4.403 on 13 and 118 DF,  p-value: 4.879e-06
Anova(clim.pca.3, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)   
## (Intercept)   5830   1  5.1200 0.025481 * 
## colutaPCA1    2205   1  1.9370 0.166613   
## colutaPCA2     886   1  0.7785 0.379386   
## colutaPCA3    3131   1  2.7500 0.099913 . 
## colutaPCA4     982   1  0.8629 0.354833   
## colutaPCA5    2244   1  1.9711 0.162960   
## Latitude      2161   1  1.8978 0.170930   
## Longitude     1221   1  1.0720 0.302620   
## Elev_m        6513   1  5.7204 0.018349 * 
## Prin1         8042   1  7.0631 0.008959 **
## Prin2            5   1  0.0044 0.947040   
## Prin3         3269   1  2.8708 0.092840 . 
## Prin4         4399   1  3.8634 0.051700 . 
## Prin5          700   1  0.6146 0.434617   
## Residuals   134357 118                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(clim.pca.3)
## [1] 132

Axis loading in Supplementary Table 10 analyzed in JMP.

Effects of raw climate predicors on BC-ratio:

library(car)

# Supplementary Table 11

## max temp. in warmest month
### model A
bc5<-lm(BCratio ~ bio5_12 + Group, data=clim.use)
summary(bc5)
## 
## Call:
## lm(formula = BCratio ~ bio5_12 + Group, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.632 -25.666   6.684  24.537  49.753 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  86.48899   22.77335   3.798 0.000179 ***
## bio5_12      -0.11655    0.09834  -1.185 0.236964    
## GroupCOL    -14.60503    9.21810  -1.584 0.114247    
## GroupNOR      5.50130    8.86755   0.620 0.535514    
## GroupUTA     -0.02333    8.81616  -0.003 0.997891    
## GroupWES     30.75290    9.20246   3.342 0.000947 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.1 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1704, Adjusted R-squared:  0.1554 
## F-statistic: 11.38 on 5 and 277 DF,  p-value: 5.361e-10
Anova(bc5, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  14862   1 14.4234 0.0001794 ***
## bio5_12       1447   1  1.4046 0.2369645    
## Group        58590   4 14.2150 1.447e-10 ***
## Residuals   285429 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5)
## [1] 283

### model B
bc5b<-lm(BCratio ~ bio5_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5 , data=rp)
summary(bc5b)
## 
## Call:
## lm(formula = BCratio ~ bio5_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -82.12 -13.50  13.75  23.36  52.56 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   71.03179   22.80899   3.114 0.002025 ** 
## bio5_12       -0.03121    0.10072  -0.310 0.756861    
## gPCA1       -167.34469   61.77015  -2.709 0.007139 ** 
## gPCA2        -78.02524   39.83432  -1.959 0.051081 .  
## gPCA3       -141.42467   41.65034  -3.396 0.000779 ***
## gPCA4         -9.20119   46.11283  -0.200 0.841980    
## gPCA5        -50.53263   39.64434  -1.275 0.203433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.71 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.06937,    Adjusted R-squared:  0.0505 
## F-statistic: 3.677 on 6 and 296 DF,  p-value: 0.00154
Anova(bc5b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  11022   1  9.6983 0.0020250 ** 
## bio5_12        109   1  0.0960 0.7568614    
## gPCA1         8341   1  7.3395 0.0071385 ** 
## gPCA2         4360   1  3.8367 0.0510815 .  
## gPCA3        13103   1 11.5296 0.0007786 ***
## gPCA4           45   1  0.0398 0.8419802    
## gPCA5         1847   1  1.6247 0.2034331    
## Residuals   336407 296                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5b)
## [1] 303

### model C
bc5c<-lm(BCratio ~ bio5_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5 , data=rp)
summary(bc5c)
## 
## Call:
## lm(formula = BCratio ~ bio5_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -69.47 -36.99  11.97  30.74  59.31 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  101.3733    45.2398   2.241   0.0268 *
## bio5_12       -0.2083     0.2015  -1.034   0.3033  
## colutaPCA1    91.4868    37.1469   2.463   0.0151 *
## colutaPCA2    50.0404    36.7856   1.360   0.1762  
## colutaPCA3  -103.6384    43.0082  -2.410   0.0174 *
## colutaPCA4   -61.0502    37.4362  -1.631   0.1055  
## colutaPCA5    -6.9976    36.7991  -0.190   0.8495  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.79 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.1523, Adjusted R-squared:  0.1116 
## F-statistic: 3.743 on 6 and 125 DF,  p-value: 0.001853
Anova(bc5c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value  Pr(>F)  
## (Intercept)   6795   1  5.0212 0.02680 *
## bio5_12       1446   1  1.0683 0.30332  
## colutaPCA1    8208   1  6.0656 0.01515 *
## colutaPCA2    2504   1  1.8505 0.17618  
## colutaPCA3    7858   1  5.8068 0.01742 *
## colutaPCA4    3599   1  2.6594 0.10545  
## colutaPCA5      49   1  0.0362 0.84950  
## Residuals   169147 125                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc5c)
## [1] 132

## annual precipitation
### model A
bc12<-lm(BCratio ~ bio12_12 + Group, data=clim.use)
summary(bc12)
## 
## Call:
## lm(formula = BCratio ~ bio12_12 + Group, data = clim.use)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -68.14 -24.23   9.05  24.26  55.12 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  21.98641   13.54149   1.624 0.105592    
## bio12_12      0.07421    0.02084   3.561 0.000435 ***
## GroupCOL    -14.71369    9.01905  -1.631 0.103942    
## GroupNOR     10.83660    8.83265   1.227 0.220910    
## GroupUTA     -0.64720    8.60981  -0.075 0.940134    
## GroupWES     32.58201    8.94739   3.642 0.000323 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.47 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.2027, Adjusted R-squared:  0.1883 
## F-statistic: 14.08 on 5 and 277 DF,  p-value: 2.821e-12
Anova(bc12, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   2611   1  2.6362 0.1055916    
## bio12_12     12555   1 12.6778 0.0004354 ***
## Group        65541   4 16.5453 3.598e-12 ***
## Residuals   274322 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12)
## [1] 283

### model B
bc12b<-lm(BCratio ~ bio12_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc12b)
## 
## Call:
## lm(formula = BCratio ~ bio12_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -76.86 -17.50  12.58  23.50  51.82 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   30.91277   11.02639   2.804  0.00539 ** 
## bio12_12       0.06638    0.02176   3.051  0.00249 ** 
## gPCA1       -181.64707   60.89438  -2.983  0.00309 ** 
## gPCA2       -124.36306   42.03333  -2.959  0.00334 ** 
## gPCA3       -164.22040   41.42907  -3.964 9.25e-05 ***
## gPCA4        -30.01382   45.89501  -0.654  0.51364    
## gPCA5        -48.59144   38.96171  -1.247  0.21333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.2 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.09744,    Adjusted R-squared:  0.07915 
## F-statistic: 5.326 on 6 and 296 DF,  p-value: 3.09e-05
Anova(bc12b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)    
## (Intercept)   8663   1  7.8598 0.005389 ** 
## bio12_12     10257   1  9.3057 0.002491 ** 
## gPCA1         9808   1  8.8982 0.003092 ** 
## gPCA2         9649   1  8.7538 0.003339 ** 
## gPCA3        17319   1 15.7124 9.25e-05 ***
## gPCA4          471   1  0.4277 0.513641    
## gPCA5         1714   1  1.5554 0.213325    
## Residuals   326259 296                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12b)
## [1] 303

### model C
bc12c<-lm(BCratio ~ bio12_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc12c)
## 
## Call:
## lm(formula = BCratio ~ bio12_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.899 -33.183   8.784  31.353  58.354 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -18.16546   25.23954  -0.720  0.47304   
## bio12_12      0.13934    0.04788   2.910  0.00427 **
## colutaPCA1   90.27632   35.77871   2.523  0.01288 * 
## colutaPCA2   24.14510   36.86889   0.655  0.51374   
## colutaPCA3  -59.25154   42.59940  -1.391  0.16673   
## colutaPCA4  -19.98069   37.59934  -0.531  0.59608   
## colutaPCA5  -26.24634   36.42508  -0.721  0.47253   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.75 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.1993, Adjusted R-squared:  0.1609 
## F-statistic: 5.186 on 6 and 125 DF,  p-value: 8.547e-05
Anova(bc12c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)   
## (Intercept)    662   1  0.5180 0.473040   
## bio12_12     10826   1  8.4704 0.004275 **
## colutaPCA1    8137   1  6.3665 0.012884 * 
## colutaPCA2     548   1  0.4289 0.513741   
## colutaPCA3    2473   1  1.9346 0.166727   
## colutaPCA4     361   1  0.2824 0.596077   
## colutaPCA5     664   1  0.5192 0.472528   
## Residuals   159766 125                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc12c)
## [1] 132

##  precip. in wettest month
### model A
bc13<-lm(BCratio ~ bio13_12 + Group, data=clim.use)
summary(bc13)
## 
## Call:
## lm(formula = BCratio ~ bio13_12 + Group, data = clim.use)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -64.21 -29.50   9.78  23.94  49.34 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  31.4261    12.1107   2.595  0.00997 **
## bio13_12      0.5100     0.1571   3.246  0.00131 **
## GroupCOL    -14.3482     9.0559  -1.584  0.11424   
## GroupNOR      7.5600     8.7325   0.866  0.38739   
## GroupUTA     -0.2106     8.6444  -0.024  0.98058   
## GroupWES     26.1274     8.9596   2.916  0.00383 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.59 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1967, Adjusted R-squared:  0.1822 
## F-statistic: 13.57 on 5 and 277 DF,  p-value: 7.555e-12
Anova(bc13, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   6718   1  6.7336  0.009966 ** 
## bio13_12     10511   1 10.5355  0.001315 ** 
## Group        46721   4 11.7070 8.467e-09 ***
## Residuals   276365 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13)
## [1] 283

### model B
bc13b<-lm(BCratio ~ bio13_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc13b)
## 
## Call:
## lm(formula = BCratio ~ bio13_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -77.96 -21.38  10.20  22.89  55.10 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.5987     9.4592   2.072  0.03914 *  
## bio13_12       0.7694     0.1604   4.797 2.56e-06 ***
## gPCA1       -176.3879    59.3996  -2.970  0.00323 ** 
## gPCA2       -121.1120    39.2882  -3.083  0.00224 ** 
## gPCA3       -122.3305    39.9222  -3.064  0.00238 ** 
## gPCA4        -59.2207    45.6079  -1.298  0.19513    
## gPCA5        -51.7694    38.1163  -1.358  0.17544    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.48 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.1362, Adjusted R-squared:  0.1187 
## F-statistic: 7.779 on 6 and 296 DF,  p-value: 8.75e-08
Anova(bc13b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)    
## (Intercept)   4528   1  4.2929 0.039138 *  
## bio13_12     24271   1 23.0087 2.56e-06 ***
## gPCA1         9302   1  8.8180 0.003227 ** 
## gPCA2        10024   1  9.5028 0.002245 ** 
## gPCA3         9905   1  9.3895 0.002383 ** 
## gPCA4         1779   1  1.6860 0.195135    
## gPCA5         1946   1  1.8447 0.175436    
## Residuals   312245 296                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13b)
## [1] 303

### model C
bc13c<-lm(BCratio ~ bio13_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc13c)
## 
## Call:
## lm(formula = BCratio ~ bio13_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -68.963 -30.347   5.277  30.042  59.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -18.2133    19.7722  -0.921 0.358742    
## bio13_12      1.2838     0.3438   3.734 0.000285 ***
## colutaPCA1   71.5247    35.2574   2.029 0.044619 *  
## colutaPCA2   60.5407    35.1439   1.723 0.087425 .  
## colutaPCA3  -63.1587    38.9497  -1.622 0.107421    
## colutaPCA4  -14.7675    36.5713  -0.404 0.687049    
## colutaPCA5  -26.4076    35.4656  -0.745 0.457911    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.04 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.2308, Adjusted R-squared:  0.1939 
## F-statistic: 6.252 on 6 and 125 DF,  p-value: 9.11e-06
Anova(bc13c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   1042   1  0.8485 0.3587423    
## bio13_12     17118   1 13.9423 0.0002851 ***
## colutaPCA1    5053   1  4.1154 0.0446191 *  
## colutaPCA2    3644   1  2.9675 0.0874249 .  
## colutaPCA3    3228   1  2.6294 0.1074207    
## colutaPCA4     200   1  0.1631 0.6870488    
## colutaPCA5     681   1  0.5544 0.4579112    
## Residuals   153474 125                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc13c)
## [1] 132

## precip. in driest month
### model A
bc14<-lm(BCratio ~ bio14_12 + Group, data=clim.use)
summary(bc14)
## 
## Call:
## lm(formula = BCratio ~ bio14_12 + Group, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.044 -25.250   7.081  24.648  56.509 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  38.4368    13.6267   2.821 0.005138 ** 
## bio14_12      0.6967     0.3375   2.064 0.039938 *  
## GroupCOL    -12.9648     9.2224  -1.406 0.160906    
## GroupNOR      9.5391     9.1198   1.046 0.296485    
## GroupUTA      1.3543     8.8079   0.154 0.877915    
## GroupWES     35.7158     9.6086   3.717 0.000244 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.94 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1788, Adjusted R-squared:  0.164 
## F-statistic: 12.06 on 5 and 277 DF,  p-value: 1.398e-10
Anova(bc14, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   8115   1  7.9564  0.005138 ** 
## bio14_12      4346   1  4.2606  0.039938 *  
## Group        60517   4 14.8330 5.388e-11 ***
## Residuals   282531 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14)
## [1] 283

### model B
bc14b<-lm(BCratio ~ bio14_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc14b)
## 
## Call:
## lm(formula = BCratio ~ bio14_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -82.06 -13.58  13.93  23.23  52.08 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   62.53316    9.46275   6.608 1.81e-10 ***
## bio14_12       0.05315    0.33684   0.158   0.8747    
## gPCA1       -166.87764   61.88091  -2.697   0.0074 ** 
## gPCA2        -78.61588   41.28409  -1.904   0.0578 .  
## gPCA3       -142.25518   44.38452  -3.205   0.0015 ** 
## gPCA4         -7.87691   46.18717  -0.171   0.8647    
## gPCA5        -49.29633   39.66441  -1.243   0.2149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.72 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.06914,    Adjusted R-squared:  0.05028 
## F-statistic: 3.664 on 6 and 296 DF,  p-value: 0.001587
Anova(bc14b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  49644   1 43.6703 1.806e-10 ***
## bio14_12        28   1  0.0249  0.874729    
## gPCA1         8267   1  7.2725  0.007403 ** 
## gPCA2         4122   1  3.6262  0.057845 .  
## gPCA3        11678   1 10.2724  0.001498 ** 
## gPCA4           33   1  0.0291  0.864700    
## gPCA5         1756   1  1.5446  0.214912    
## Residuals   336488 296                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14b)
## [1] 303

### model C
bc14c<-lm(BCratio ~ bio14_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc14c)
## 
## Call:
## lm(formula = BCratio ~ bio14_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -68.62 -37.15  11.88  29.97  63.92 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   25.2221    25.1162   1.004  0.31722   
## bio14_12       1.0042     0.8477   1.185  0.23845   
## colutaPCA1   101.4081    38.9316   2.605  0.01031 * 
## colutaPCA2    10.5721    49.7907   0.212  0.83219   
## colutaPCA3  -107.3120    40.2074  -2.669  0.00862 **
## colutaPCA4   -57.9553    36.8986  -1.571  0.11879   
## colutaPCA5    -5.1293    36.7422  -0.140  0.88920   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.74 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.1545, Adjusted R-squared:  0.114 
## F-statistic: 3.808 on 6 and 125 DF,  p-value: 0.001611
Anova(bc14c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value   Pr(>F)   
## (Intercept)   1361   1  1.0084 0.317215   
## bio14_12      1894   1  1.4031 0.238445   
## colutaPCA1    9157   1  6.7849 0.010308 * 
## colutaPCA2      61   1  0.0451 0.832195   
## colutaPCA3    9614   1  7.1233 0.008618 **
## colutaPCA4    3329   1  2.4670 0.118789   
## colutaPCA5      26   1  0.0195 0.889200   
## Residuals   168699 125                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc14c)
## [1] 132

##  precip. in wettest quarter
### model A
bc16<-lm(BCratio ~ bio16_12 + Group, data=clim.use)
summary(bc16)
## 
## Call:
## lm(formula = BCratio ~ bio16_12 + Group, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.664 -29.912   9.722  23.467  48.310 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  32.00080   11.68921   2.738  0.00659 ** 
## bio16_12      0.18621    0.05498   3.387  0.00081 ***
## GroupCOL    -15.24902    9.03630  -1.688  0.09263 .  
## GroupNOR      7.58900    8.71502   0.871  0.38462    
## GroupUTA     -0.47207    8.62853  -0.055  0.95641    
## GroupWES     26.45400    8.93369   2.961  0.00333 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.54 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1993, Adjusted R-squared:  0.1849 
## F-statistic: 13.79 on 5 and 277 DF,  p-value: 4.912e-12
Anova(bc16, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   7453   1  7.4946 0.0065889 ** 
## bio16_12     11407   1 11.4703 0.0008096 ***
## Group        50361   4 12.6602 1.783e-09 ***
## Residuals   275470 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16)
## [1] 283

### model B
bc16b<-lm(BCratio ~ bio16_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc16b)
## 
## Call:
## lm(formula = BCratio ~ bio16_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -77.55 -21.46  10.57  23.10  54.71 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   21.7671     9.0143   2.415 0.016353 *  
## bio16_12       0.2702     0.0563   4.798 2.54e-06 ***
## gPCA1       -181.1240    59.4423  -3.047 0.002519 ** 
## gPCA2       -120.6882    39.2658  -3.074 0.002312 ** 
## gPCA3       -143.0147    39.7630  -3.597 0.000378 ***
## gPCA4        -62.7583    45.7826  -1.371 0.171479    
## gPCA5        -52.5443    38.1174  -1.378 0.169095    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.48 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.1363, Adjusted R-squared:  0.1187 
## F-statistic: 7.782 on 6 and 296 DF,  p-value: 8.689e-08
Anova(bc16b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   6151   1  5.8309 0.0163529 *  
## bio16_12     24288   1 23.0253 2.539e-06 ***
## gPCA1         9794   1  9.2846 0.0025193 ** 
## gPCA2         9965   1  9.4471 0.0023118 ** 
## gPCA3        13645   1 12.9361 0.0003776 ***
## gPCA4         1982   1  1.8791 0.1714793    
## gPCA5         2004   1  1.9002 0.1690949    
## Residuals   312229 296                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16b)
## [1] 303

### model C
bc16c<-lm(BCratio ~ bio16_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc16c)
## 
## Call:
## lm(formula = BCratio ~ bio16_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -69.695 -30.620   5.397  29.917  60.002 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -20.1810    19.3294  -1.044 0.298475    
## bio16_12      0.4825     0.1229   3.924 0.000143 ***
## colutaPCA1   89.4113    34.8671   2.564 0.011520 *  
## colutaPCA2   33.6537    35.1160   0.958 0.339733    
## colutaPCA3  -68.4244    37.8858  -1.806 0.073313 .  
## colutaPCA4   -3.0381    37.1866  -0.082 0.935016    
## colutaPCA5  -22.5663    35.1136  -0.643 0.521618    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 34.86 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.2388, Adjusted R-squared:  0.2023 
## F-statistic: 6.536 on 6 and 125 DF,  p-value: 5.057e-06
Anova(bc16c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   1324   1  1.0901 0.2984745    
## bio16_12     18711   1 15.3993 0.0001428 ***
## colutaPCA1    7990   1  6.5759 0.0115202 *  
## colutaPCA2    1116   1  0.9184 0.3397328    
## colutaPCA3    3963   1  3.2619 0.0733133 .  
## colutaPCA4       8   1  0.0067 0.9350160    
## colutaPCA5     502   1  0.4130 0.5216177    
## Residuals   151881 125                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc16c)
## [1] 132

## precip. in driest quarter
### model A
bc17<-lm(BCratio ~ bio17_12 + Group, data=clim.use)
summary(bc17)
## 
## Call:
## lm(formula = BCratio ~ bio17_12 + Group, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.214 -23.017   7.181  24.505  60.344 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.0806    14.2954   2.384   0.0178 *  
## bio17_12      0.2483     0.1085   2.288   0.0229 *  
## GroupCOL    -13.4026     9.1740  -1.461   0.1452    
## GroupNOR      9.6519     9.0543   1.066   0.2874    
## GroupUTA      0.5041     8.7451   0.058   0.9541    
## GroupWES     35.8184     9.5011   3.770   0.0002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.88 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1816, Adjusted R-squared:  0.1669 
## F-statistic:  12.3 on 5 and 277 DF,  p-value: 8.863e-11
Anova(bc17, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)   5777   1  5.6836    0.0178 *  
## bio17_12      5320   1  5.2339    0.0229 *  
## Group        61719   4 15.1799 3.103e-11 ***
## Residuals   281557 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17)
## [1] 283

### model B
bc17b<-lm(BCratio ~ bio17_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc17b)
## 
## Call:
## lm(formula = BCratio ~ bio17_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -82.26 -13.52  13.91  23.16  52.18 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.314e+01  1.031e+01   6.127 2.85e-09 ***
## bio17_12     8.892e-03  1.059e-01   0.084  0.93315    
## gPCA1       -1.664e+02  6.180e+01  -2.692  0.00750 ** 
## gPCA2       -7.790e+01  4.177e+01  -1.865  0.06316 .  
## gPCA3       -1.407e+02  4.308e+01  -3.267  0.00122 ** 
## gPCA4       -8.145e+00  4.618e+01  -0.176  0.86013    
## gPCA5       -4.960e+01  3.960e+01  -1.252  0.21144    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.72 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.06909,    Adjusted R-squared:  0.05022 
## F-statistic: 3.661 on 6 and 296 DF,  p-value: 0.001598
Anova(bc17b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  42678   1 37.5408 2.848e-09 ***
## bio17_12         8   1  0.0070  0.933146    
## gPCA1         8240   1  7.2477  0.007503 ** 
## gPCA2         3955   1  3.4786  0.063157 .  
## gPCA3        12131   1 10.6703  0.001217 ** 
## gPCA4           35   1  0.0311  0.860128    
## gPCA5         1783   1  1.5683  0.211441    
## Residuals   336508 296                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17b)
## [1] 303

### model C
bc17c<-lm(BCratio ~ bio17_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc17c)
## 
## Call:
## lm(formula = BCratio ~ bio17_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -69.30 -37.29  12.32  31.80  58.77 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   33.7665    25.9872   1.299   0.1962  
## bio17_12       0.2040     0.2509   0.813   0.4178  
## colutaPCA1    89.8805    37.1306   2.421   0.0169 *
## colutaPCA2    36.4796    40.6205   0.898   0.3709  
## colutaPCA3  -106.9157    44.1366  -2.422   0.0169 *
## colutaPCA4   -51.0097    37.0126  -1.378   0.1706  
## colutaPCA5    -8.2230    36.9518  -0.223   0.8243  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.85 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.1495, Adjusted R-squared:  0.1087 
## F-statistic: 3.663 on 6 and 125 DF,  p-value: 0.002195
Anova(bc17c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value  Pr(>F)  
## (Intercept)   2292   1  1.6883 0.19622  
## bio17_12       897   1  0.6609 0.41779  
## colutaPCA1    7955   1  5.8596 0.01693 *
## colutaPCA2    1095   1  0.8065 0.37088  
## colutaPCA3    7966   1  5.8679 0.01685 *
## colutaPCA4    2578   1  1.8993 0.17061  
## colutaPCA5      67   1  0.0495 0.82426  
## Residuals   169695 125                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc17c)
## [1] 132

## precip. in coldest quarter
### model A
bc19<-lm(BCratio ~ bio19_12 + Group, data=clim.use)
summary(bc19)
## 
## Call:
## lm(formula = BCratio ~ bio19_12 + Group, data = clim.use)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.060 -26.411   8.437  24.632  54.141 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  49.74063    9.95284   4.998 1.03e-06 ***
## bio19_12      0.08309    0.04290   1.937  0.05380 .  
## GroupCOL    -14.30358    9.17478  -1.559  0.12014    
## GroupNOR      6.31626    8.83765   0.715  0.47540    
## GroupUTA     -1.21043    8.74631  -0.138  0.89003    
## GroupWES     29.01108    9.02699   3.214  0.00146 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.97 on 277 degrees of freedom
##   (201 observations deleted due to missingness)
## Multiple R-squared:  0.1773, Adjusted R-squared:  0.1625 
## F-statistic: 11.94 on 5 and 277 DF,  p-value: 1.776e-10
Anova(bc19, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  25521   1 24.9763 1.030e-06 ***
## bio19_12      3832   1  3.7506    0.0538 .  
## Group        56035   4 13.7095 3.261e-10 ***
## Residuals   283044 277                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19)
## [1] 283

### model B
bc19b<-lm(BCratio ~ bio19_12 + gPCA1+gPCA2+gPCA3+gPCA4+gPCA5, data=rp)
summary(bc19b)
## 
## Call:
## lm(formula = BCratio ~ bio19_12 + gPCA1 + gPCA2 + gPCA3 + gPCA4 + 
##     gPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -77.15 -19.29  11.87  23.98  54.12 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   47.88087    6.21147   7.708 1.94e-13 ***
## bio19_12       0.12179    0.04444   2.740 0.006510 ** 
## gPCA1       -168.83018   60.86565  -2.774 0.005892 ** 
## gPCA2       -101.70364   40.19061  -2.531 0.011907 *  
## gPCA3       -160.96162   41.49379  -3.879 0.000129 ***
## gPCA4        -35.38125   46.53140  -0.760 0.447636    
## gPCA5        -49.58843   39.07503  -1.269 0.205418    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.3 on 296 degrees of freedom
##   (259 observations deleted due to missingness)
## Multiple R-squared:  0.0921, Adjusted R-squared:  0.0737 
## F-statistic: 5.004 on 6 and 296 DF,  p-value: 6.664e-05
Anova(bc19b, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  65882   1 59.4203 1.944e-13 ***
## bio19_12      8326   1  7.5097 0.0065096 ** 
## gPCA1         8531   1  7.6941 0.0058922 ** 
## gPCA2         7100   1  6.4036 0.0119072 *  
## gPCA3        16684   1 15.0480 0.0001292 ***
## gPCA4          641   1  0.5782 0.4476365    
## gPCA5         1786   1  1.6105 0.2054178    
## Residuals   328190 296                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19b)
## [1] 303

### model C
bc19c<-lm(BCratio ~ bio19_12 + colutaPCA1+colutaPCA2+colutaPCA3+colutaPCA4+colutaPCA5, data=rp)
summary(bc19c)
## 
## Call:
## lm(formula = BCratio ~ bio19_12 + colutaPCA1 + colutaPCA2 + colutaPCA3 + 
##     colutaPCA4 + colutaPCA5, data = rp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -68.11 -37.24  11.64  30.79  57.92 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   38.8541    18.6604   2.082   0.0394 *
## bio19_12       0.1172     0.1356   0.864   0.3894  
## colutaPCA1    82.2193    37.1110   2.215   0.0285 *
## colutaPCA2    35.7888    40.5232   0.883   0.3788  
## colutaPCA3  -108.1497    42.6199  -2.538   0.0124 *
## colutaPCA4   -39.1684    40.5749  -0.965   0.3362  
## colutaPCA5   -21.9217    41.2226  -0.532   0.5958  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 36.83 on 125 degrees of freedom
##   (430 observations deleted due to missingness)
## Multiple R-squared:  0.1501, Adjusted R-squared:  0.1093 
## F-statistic:  3.68 on 6 and 125 DF,  p-value: 0.002119
Anova(bc19c, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq  Df F value  Pr(>F)  
## (Intercept)   5882   1  4.3354 0.03937 *
## bio19_12      1012   1  0.7460 0.38941  
## colutaPCA1    6659   1  4.9084 0.02854 *
## colutaPCA2    1058   1  0.7800 0.37884  
## colutaPCA3    8736   1  6.4391 0.01239 *
## colutaPCA4    1264   1  0.9319 0.33624  
## colutaPCA5     384   1  0.2828 0.59582  
## Residuals   169580 125                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(bc19c)
## [1] 132

Validate p-values from linear models using permutation tests:

# Supplementary Table 11, sections A
set.seed(1)
nperm<-10000

var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5, bc12, bc13, bc14, bc16, bc17, bc19)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')

i<-1
j<-1

start<-Sys.time()

for (i in 1:length(var.list)){
  perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
  colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
  perm.out.group<-perm.out.climate
  for(j in 1:nperm){
    perm.BC<-sample(clim.use$BCratio)
    tmp<-cbind(clim.use, perm.BC)
    model<-lm(paste('perm.BC ~', var.list[i], ' + Group'), data=tmp)
    test<-data.frame(Anova(model, type="3"))
    perm.out.climate[j,]<-test[2,]
    perm.out.group[j,]<-test[3,]
  }
  
  # calculate empirical cdfs from permuted test statistics 
  pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
  pct.group<-ecdf(perm.out.group$F.value)
  
  # calculate observed test statistics from linear models
  obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=clim.use)), type=3)", sep="")))[2,3]
  obs.F.group<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=clim.use)), type=3)", sep="")))[3,3]
  
  # calculate p-values for climate and group from permutation tests
  p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
  p.group<-(1-pct.group(obs.F.group)) # 1 tailed test
  
  # calculate location of 95% tail in climate F distribution
  p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]

  # plot observed F onto null distribution
  p<-ggplot(data=perm.out.climate) + 
    geom_histogram(aes(x = F.value), 
                   color="lightgray", fill="lightgray", bins = 50) +
  labs(title=paste(name.list[i], " (P = ", 
                   sprintf("%0.4f", p.clim),")", sep=""), 
       x = "Permuted F-value", y = "Frequency") +
    geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50), 
                  lineend = 'butt', linejoin = 'mitre', 
                  arrow=arrow(length=unit(.1, "inches")), 
                  size=.5, color="red") +
    annotate("text", x=obs.F.clim, y=1000, 
             label="Observed F", hjust=0.5, color="red", size=4)+
    geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
    theme(plot.title = element_text(size=6))+
    theme_bw()
  
  # save plot
  svg(filename = paste("SupFig4-PermTest-", var.list[i], ".svg", sep=""),
         width = 4, height = 2, pointsize = 6)
  print(p)
  dev.off()
  
  # print output for table
  print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
  print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
  print(paste("Group: F = ", sprintf("%0.4f", obs.F.group), ", P = ", sprintf("%0.6f", p.group), sep=""))
  print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 1.4046, P = 0.230600"
## [1] "Group: F = 14.2150, P = 0.000000"
## [1] "Elapsed time: 1.59"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 12.6778, P = 0.000900"
## [1] "Group: F = 16.5453, P = 0.000000"
## [1] "Elapsed time: 3.2"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 10.5355, P = 0.001400"
## [1] "Group: F = 11.7070, P = 0.000000"
## [1] "Elapsed time: 4.81"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 4.2606, P = 0.038800"
## [1] "Group: F = 14.8330, P = 0.000000"
## [1] "Elapsed time: 6.29"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 11.4703, P = 0.000500"
## [1] "Group: F = 12.6602, P = 0.000000"
## [1] "Elapsed time: 7.77"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 5.2339, P = 0.022200"
## [1] "Group: F = 15.1799, P = 0.000000"
## [1] "Elapsed time: 9.25"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 3.7506, P = 0.054000"
## [1] "Group: F = 13.7095, P = 0.000000"
## [1] "Elapsed time: 10.72"
# Supplementary Table 11, sections B
set.seed(1)

nperm<-10000

var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5b, bc12b, bc13b, bc14b, bc16b, bc17b, bc19b)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')

i<-1
j<-1

start<-Sys.time()

for (i in 1:length(var.list)){
  perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
  colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
  perm.out.pca1<-perm.out.climate
  perm.out.pca2<-perm.out.climate
  perm.out.pca3<-perm.out.climate
  perm.out.pca4<-perm.out.climate
  perm.out.pca5<-perm.out.climate

  for(j in 1:nperm){
    perm.BC<-sample(rp$BCratio)
    tmp<-cbind(rp, perm.BC)
    model<-lm(paste('perm.BC ~', var.list[i], ' + gPCA1 + gPCA2 + gPCA3 + gPCA4 + gPCA5'), data=tmp)
    test<-data.frame(Anova(model, type="3"))
    perm.out.climate[j,]<-test[2,]
    perm.out.pca1[j,]<-test[3,]
    perm.out.pca2[j,]<-test[4,]
    perm.out.pca3[j,]<-test[5,]
    perm.out.pca4[j,]<-test[6,]
    perm.out.pca5[j,]<-test[7,]
  }
  
  # calculate empirical cdfs from permuted test statistics 
  pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
  pct.pca1<-ecdf(perm.out.pca1$F.value)
  pct.pca2<-ecdf(perm.out.pca2$F.value)
  pct.pca3<-ecdf(perm.out.pca3$F.value)
  pct.pca4<-ecdf(perm.out.pca4$F.value)
  pct.pca5<-ecdf(perm.out.pca5$F.value)
  
  # calculate observed test statistics from linear models
  obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[2,3]
  
  obs.F.pca1<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[3,3]
  obs.F.pca2<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[4,3]
  obs.F.pca3<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[5,3]
  obs.F.pca4<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[6,3]
  obs.F.pca5<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[7,3]
  
  # calculate p-values for climate and group from permutation tests
  p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
  p.pca1<-(1-pct.pca1(obs.F.pca1)) # 1 tailed test
  p.pca2<-(1-pct.pca2(obs.F.pca2)) # 1 tailed test
  p.pca3<-(1-pct.pca3(obs.F.pca3)) # 1 tailed test
  p.pca4<-(1-pct.pca4(obs.F.pca4)) # 1 tailed test
  p.pca5<-(1-pct.pca5(obs.F.pca5)) # 1 tailed test
  
  # calculate location of 95% tail in climate F distribution
  p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]

  # plot observed F onto null distribution
  p<-ggplot(data=perm.out.climate) + 
    geom_histogram(aes(x = F.value), 
                   color="lightgray", fill="lightgray", bins = 50) +
  labs(title=paste(name.list[i], " (P = ", 
                   sprintf("%0.4f", p.clim),")", sep=""), 
       x = "Permuted F-value", y = "Frequency") +
    geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50), 
                  lineend = 'butt', linejoin = 'mitre', 
                  arrow=arrow(length=unit(.1, "inches")), 
                  size=.5, color="red") +
    annotate("text", x=obs.F.clim, y=1000, 
             label="Observed F", hjust=0.5, color="red", size=4)+
    geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
    theme(plot.title = element_text(size=6))+
    theme_bw()
  
  # save plot
  svg(filename = paste("SupFig4-PermTest-", var.list[i], "-GeneticPCA-all.svg", sep=""),
         width = 4, height = 2, pointsize = 6)
  print(p)
  dev.off()
  
  # print output for table
  print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
  print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
  print(paste("Genetic PCA1: F = ", sprintf("%0.4f", obs.F.pca1), ", P = ", sprintf("%0.6f", p.pca1), sep=""))
  print(paste("Genetic PCA2: F = ", sprintf("%0.4f", obs.F.pca2), ", P = ", sprintf("%0.6f", p.pca2), sep=""))
  print(paste("Genetic PCA3: F = ", sprintf("%0.4f", obs.F.pca3), ", P = ", sprintf("%0.6f", p.pca3), sep=""))
  print(paste("Genetic PCA4: F = ", sprintf("%0.4f", obs.F.pca4), ", P = ", sprintf("%0.6f", p.pca4), sep=""))
  print(paste("Genetic PCA5: F = ", sprintf("%0.4f", obs.F.pca5), ", P = ", sprintf("%0.6f", p.pca5), sep=""))
  print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 0.0960, P = 0.759100"
## [1] "Genetic PCA1: F = 7.3395, P = 0.006800"
## [1] "Genetic PCA2: F = 3.8367, P = 0.049900"
## [1] "Genetic PCA3: F = 11.5296, P = 0.000700"
## [1] "Genetic PCA4: F = 0.0398, P = 0.840900"
## [1] "Genetic PCA5: F = 1.6247, P = 0.197400"
## [1] "Elapsed time: 2.58"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 9.3057, P = 0.002600"
## [1] "Genetic PCA1: F = 8.8982, P = 0.002700"
## [1] "Genetic PCA2: F = 8.7538, P = 0.003100"
## [1] "Genetic PCA3: F = 15.7124, P = 0.000200"
## [1] "Genetic PCA4: F = 0.4277, P = 0.510800"
## [1] "Genetic PCA5: F = 1.5554, P = 0.210900"
## [1] "Elapsed time: 5.21"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 23.0087, P = 0.000000"
## [1] "Genetic PCA1: F = 8.8180, P = 0.003200"
## [1] "Genetic PCA2: F = 9.5028, P = 0.001900"
## [1] "Genetic PCA3: F = 9.3895, P = 0.002500"
## [1] "Genetic PCA4: F = 1.6860, P = 0.197700"
## [1] "Genetic PCA5: F = 1.8447, P = 0.178900"
## [1] "Elapsed time: 7.88"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 0.0249, P = 0.872400"
## [1] "Genetic PCA1: F = 7.2725, P = 0.007800"
## [1] "Genetic PCA2: F = 3.6262, P = 0.063600"
## [1] "Genetic PCA3: F = 10.2724, P = 0.002000"
## [1] "Genetic PCA4: F = 0.0291, P = 0.866900"
## [1] "Genetic PCA5: F = 1.5446, P = 0.214500"
## [1] "Elapsed time: 10.48"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 23.0253, P = 0.000000"
## [1] "Genetic PCA1: F = 9.2846, P = 0.002700"
## [1] "Genetic PCA2: F = 9.4471, P = 0.002100"
## [1] "Genetic PCA3: F = 12.9361, P = 0.000500"
## [1] "Genetic PCA4: F = 1.8791, P = 0.168600"
## [1] "Genetic PCA5: F = 1.9002, P = 0.162300"
## [1] "Elapsed time: 13.16"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 0.0070, P = 0.932700"
## [1] "Genetic PCA1: F = 7.2477, P = 0.008100"
## [1] "Genetic PCA2: F = 3.4786, P = 0.065300"
## [1] "Genetic PCA3: F = 10.6703, P = 0.001400"
## [1] "Genetic PCA4: F = 0.0311, P = 0.866100"
## [1] "Genetic PCA5: F = 1.5683, P = 0.210900"
## [1] "Elapsed time: 15.79"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 7.5097, P = 0.007100"
## [1] "Genetic PCA1: F = 7.6941, P = 0.005300"
## [1] "Genetic PCA2: F = 6.4036, P = 0.011400"
## [1] "Genetic PCA3: F = 15.0480, P = 0.000200"
## [1] "Genetic PCA4: F = 0.5782, P = 0.455600"
## [1] "Genetic PCA5: F = 1.6105, P = 0.214400"
## [1] "Elapsed time: 18.38"
# Supplementary Table 11, sections C
set.seed(1)

nperm<-10000

var.list<-c('bio5_12', 'bio12_12', 'bio13_12', 'bio14_12', 'bio16_12', 'bio17_12', 'bio19_12')
mod.list<-list(bc5c, bc12c, bc13c, bc14c, bc16c, bc17c, bc19c)
name.list<-c('Max. T in warmest month', 'Annual precipitation', 'Precip. in wettest month', 'Precip. in driest month', 'Precip. in wettest quarter', 'Precip. in driest quarter', 'Precip. in coldest quarter')

i<-1
j<-1

start<-Sys.time()

for (i in 1:length(var.list)){
  perm.out.climate<-data.frame(matrix(NA, nrow=nperm, ncol=4))
  colnames(perm.out.climate)<-c("Sum.Sq", "Df","F.value", "P")
  perm.out.pca1<-perm.out.climate
  perm.out.pca2<-perm.out.climate
  perm.out.pca3<-perm.out.climate
  perm.out.pca4<-perm.out.climate
  perm.out.pca5<-perm.out.climate

  for(j in 1:nperm){
    perm.BC<-sample(rp$BCratio)
    tmp<-cbind(rp, perm.BC)
    model<-lm(paste('perm.BC ~', var.list[i], ' + colutaPCA1 + colutaPCA2 + colutaPCA3 + colutaPCA4 + colutaPCA5'), data=tmp)
    test<-data.frame(Anova(model, type="3"))
    perm.out.climate[j,]<-test[2,]
    perm.out.pca1[j,]<-test[3,]
    perm.out.pca2[j,]<-test[4,]
    perm.out.pca3[j,]<-test[5,]
    perm.out.pca4[j,]<-test[6,]
    perm.out.pca5[j,]<-test[7,]
  }
  
  # calculate empirical cdfs from permuted test statistics 
  pct.clim <- ecdf(perm.out.climate$F.value) # get CDF within null distribution
  pct.pca1<-ecdf(perm.out.pca1$F.value)
  pct.pca2<-ecdf(perm.out.pca2$F.value)
  pct.pca3<-ecdf(perm.out.pca3$F.value)
  pct.pca4<-ecdf(perm.out.pca4$F.value)
  pct.pca5<-ecdf(perm.out.pca5$F.value)
  
  # calculate observed test statistics from linear models
  obs.F.clim<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[2,3]
  
  obs.F.pca1<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[3,3]
  obs.F.pca2<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[4,3]
  obs.F.pca3<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[5,3]
  obs.F.pca4<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[6,3]
  obs.F.pca5<-eval(parse(text=paste("Anova((lm(", mod.list[[i]]$call[2], ", data=rp)), type=3)", sep="")))[7,3]
  
  # calculate p-values for climate and group from permutation tests
  p.clim <-(1-pct.clim(obs.F.clim)) # 1 tailed test
  p.pca1<-(1-pct.pca1(obs.F.pca1)) # 1 tailed test
  p.pca2<-(1-pct.pca2(obs.F.pca2)) # 1 tailed test
  p.pca3<-(1-pct.pca3(obs.F.pca3)) # 1 tailed test
  p.pca4<-(1-pct.pca4(obs.F.pca4)) # 1 tailed test
  p.pca5<-(1-pct.pca5(obs.F.pca5)) # 1 tailed test
  
  # calculate location of 95% tail in climate F distribution
  p.95.clim<-arrange(perm.out.climate, F.value)$F.value[0.95*nrow(perm.out.climate)]

  # plot observed F onto null distribution
  p<-ggplot(data=perm.out.climate) + 
    geom_histogram(aes(x = F.value), 
                   color="lightgray", fill="lightgray", bins = 50) +
  labs(title=paste(name.list[i], " (P = ", 
                   sprintf("%0.4f", p.clim),")", sep=""), 
       x = "Permuted F-value", y = "Frequency") +
    geom_segment(mapping=aes(x=obs.F.clim, y=900, xend=obs.F.clim, yend=50), 
                  lineend = 'butt', linejoin = 'mitre', 
                  arrow=arrow(length=unit(.1, "inches")), 
                  size=.5, color="red") +
    annotate("text", x=obs.F.clim, y=1000, 
             label="Observed F", hjust=0.5, color="red", size=4)+
    geom_vline(xintercept=p.95.clim, colour="black", linetype="dashed")+
    theme(plot.title = element_text(size=6))+
    theme_bw()
  
  # save plot
  svg(filename = paste("SupFig4-PermTest-", var.list[i], "-GeneticPCA-coluta.svg", sep=""),
         width = 4, height = 2, pointsize = 6)
  print(p)
  dev.off()
  
  # print output for table
  print(paste("Climate predictor: ", var.list[i], " (", name.list[i], ")", sep=""))
  print(paste("Climate: F = ", sprintf("%0.4f", obs.F.clim), ", P = ", sprintf("%0.6f", p.clim), sep=""))
  print(paste("Genetic PCA1: F = ", sprintf("%0.4f", obs.F.pca1), ", P = ", sprintf("%0.6f", p.pca1), sep=""))
  print(paste("Genetic PCA2: F = ", sprintf("%0.4f", obs.F.pca2), ", P = ", sprintf("%0.6f", p.pca2), sep=""))
  print(paste("Genetic PCA3: F = ", sprintf("%0.4f", obs.F.pca3), ", P = ", sprintf("%0.6f", p.pca3), sep=""))
  print(paste("Genetic PCA4: F = ", sprintf("%0.4f", obs.F.pca4), ", P = ", sprintf("%0.6f", p.pca4), sep=""))
  print(paste("Genetic PCA5: F = ", sprintf("%0.4f", obs.F.pca5), ", P = ", sprintf("%0.6f", p.pca5), sep=""))
  print(paste("Elapsed time:", round((Sys.time()-start), digits=2)))
}
## [1] "Climate predictor: bio5_12 (Max. T in warmest month)"
## [1] "Climate: F = 1.0683, P = 0.303800"
## [1] "Genetic PCA1: F = 6.0656, P = 0.017600"
## [1] "Genetic PCA2: F = 1.8505, P = 0.180500"
## [1] "Genetic PCA3: F = 5.8068, P = 0.015100"
## [1] "Genetic PCA4: F = 2.6594, P = 0.102000"
## [1] "Genetic PCA5: F = 0.0362, P = 0.846100"
## [1] "Elapsed time: 2.6"
## [1] "Climate predictor: bio12_12 (Annual precipitation)"
## [1] "Climate: F = 8.4704, P = 0.004800"
## [1] "Genetic PCA1: F = 6.3665, P = 0.012700"
## [1] "Genetic PCA2: F = 0.4289, P = 0.521200"
## [1] "Genetic PCA3: F = 1.9346, P = 0.166400"
## [1] "Genetic PCA4: F = 0.2824, P = 0.599600"
## [1] "Genetic PCA5: F = 0.5192, P = 0.469600"
## [1] "Elapsed time: 5.24"
## [1] "Climate predictor: bio13_12 (Precip. in wettest month)"
## [1] "Climate: F = 13.9423, P = 0.000200"
## [1] "Genetic PCA1: F = 4.1154, P = 0.047200"
## [1] "Genetic PCA2: F = 2.9675, P = 0.091100"
## [1] "Genetic PCA3: F = 2.6294, P = 0.102600"
## [1] "Genetic PCA4: F = 0.1631, P = 0.681700"
## [1] "Genetic PCA5: F = 0.5544, P = 0.454700"
## [1] "Elapsed time: 7.84"
## [1] "Climate predictor: bio14_12 (Precip. in driest month)"
## [1] "Climate: F = 1.4031, P = 0.242600"
## [1] "Genetic PCA1: F = 6.7849, P = 0.011000"
## [1] "Genetic PCA2: F = 0.0451, P = 0.829800"
## [1] "Genetic PCA3: F = 7.1233, P = 0.007000"
## [1] "Genetic PCA4: F = 2.4670, P = 0.120800"
## [1] "Genetic PCA5: F = 0.0195, P = 0.890000"
## [1] "Elapsed time: 10.43"
## [1] "Climate predictor: bio16_12 (Precip. in wettest quarter)"
## [1] "Climate: F = 15.3993, P = 0.000100"
## [1] "Genetic PCA1: F = 6.5759, P = 0.013100"
## [1] "Genetic PCA2: F = 0.9184, P = 0.347400"
## [1] "Genetic PCA3: F = 3.2619, P = 0.067700"
## [1] "Genetic PCA4: F = 0.0067, P = 0.935700"
## [1] "Genetic PCA5: F = 0.4130, P = 0.528100"
## [1] "Elapsed time: 13.11"
## [1] "Climate predictor: bio17_12 (Precip. in driest quarter)"
## [1] "Climate: F = 0.6609, P = 0.423000"
## [1] "Genetic PCA1: F = 5.8596, P = 0.016600"
## [1] "Genetic PCA2: F = 0.8065, P = 0.367500"
## [1] "Genetic PCA3: F = 5.8679, P = 0.017100"
## [1] "Genetic PCA4: F = 1.8993, P = 0.165300"
## [1] "Genetic PCA5: F = 0.0495, P = 0.823700"
## [1] "Elapsed time: 15.7"
## [1] "Climate predictor: bio19_12 (Precip. in coldest quarter)"
## [1] "Climate: F = 0.7460, P = 0.394400"
## [1] "Genetic PCA1: F = 4.9084, P = 0.029900"
## [1] "Genetic PCA2: F = 0.7800, P = 0.384100"
## [1] "Genetic PCA3: F = 6.4391, P = 0.009700"
## [1] "Genetic PCA4: F = 0.9319, P = 0.342800"
## [1] "Genetic PCA5: F = 0.2828, P = 0.585900"
## [1] "Elapsed time: 18.3"

Allelic variation in BCMA

Structural variants in BCMA3 and transcribed amino acid sequences were determined from .fasta sequence data using scripts in ‘/BCMA archive.R1/Sanger data/Sequence analysis/sanger_seq_analysis-v2.Rmd’.

Distribution of and covariance among alleles:

# Contingency tables of alleles
table(rp$Genotype, rp$Group)
##                                
##                                 Admixed COL NOR UTA WES
##   Complete exons + No deletion        5  18  18  14  22
##   Premature stop + 8bp deletion       1  16   0   3   0
##   Premature stop + No deletion        0   2   0   1   0
table(rp$aa.concat, rp$Group)
##     
##      Admixed COL NOR UTA WES
##   **       1  16   0   3   0
##   LV       2  10   3   5   0
##   VM       0   0   2   0  14
##   VV       3  10  13  10   8
table(rp$Genotype, rp$aa.concat, rp$Group)
## , ,  = Admixed
## 
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0  2  0  3
##   Premature stop + 8bp deletion  1  0  0  0
##   Premature stop + No deletion   0  0  0  0
## 
## , ,  = COL
## 
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0  8  0 10
##   Premature stop + 8bp deletion 16  0  0  0
##   Premature stop + No deletion   0  2  0  0
## 
## , ,  = NOR
## 
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0  3  2 13
##   Premature stop + 8bp deletion  0  0  0  0
##   Premature stop + No deletion   0  0  0  0
## 
## , ,  = UTA
## 
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0  4  0 10
##   Premature stop + 8bp deletion  3  0  0  0
##   Premature stop + No deletion   0  1  0  0
## 
## , ,  = WES
## 
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0  0 14  8
##   Premature stop + 8bp deletion  0  0  0  0
##   Premature stop + No deletion   0  0  0  0
table(rp$Genotype, rp$aa.concat)
##                                
##                                 ** LV VM VV
##   Complete exons + No deletion   0 18 18 51
##   Premature stop + 8bp deletion 20  0  0  0
##   Premature stop + No deletion   0  3  0  0
table(rp$Deletion, rp$Full_Length)
##               
##                Complete exons Premature stop
##   8bp deletion              0             20
##   No deletion              87              3

# Measure geographic distances of alleles
full.length<-rp %>% filter(Genotype=="Complete exons + No deletion") %>%dplyr::select(Longitude, Latitude)

stop.only<-rp %>% filter(Genotype =="Premature stop + No deletion") %>% dplyr::select(Latitude, Longitude)

stop.del<-rp%>% filter(Genotype=="Premature stop + 8bp deletion") %>% dplyr::select(Latitude, Longitude)

georange(full.length)/1000
##      minimum      maximum 
##    0.3648679 1216.6627968
georange(stop.only)/1000
##  minimum  maximum 
## 104.7254 328.2988
georange(stop.del)/1000
##     minimum     maximum 
##   0.4168024 729.0933943

Effects of BCMA3 variants on BC-ratio:

library(car)

# Data subset with only full-length alleles
final.fl<-filter(rp, Full_Length=="Complete exons")

# Supplementary Table 12A:
met.sep.1<-lm(BCratio~Group+Full_Length+Deletion, data=rp)
summary(met.sep.1)
## 
## Call:
## lm(formula = BCratio ~ Group + Full_Length + Deletion, data = rp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.522  -2.136   3.575  15.882  35.852 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                  70.27      30.34   2.316  0.02327 * 
## GroupCOL                     10.92      12.57   0.869  0.38777   
## GroupNOR                     11.44      12.64   0.905  0.36858   
## GroupUTA                     16.58      12.98   1.277  0.20558   
## GroupWES                     35.37      12.36   2.862  0.00546 **
## Full_LengthPremature stop   -71.98      27.57  -2.610  0.01091 * 
## DeletionNo deletion         -13.29      28.74  -0.463  0.64504   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.61 on 75 degrees of freedom
##   (480 observations deleted due to missingness)
## Multiple R-squared:  0.5283, Adjusted R-squared:  0.4905 
## F-statistic:    14 on 6 and 75 DF,  p-value: 1.288e-10
Anova(met.sep.1, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq Df F value  Pr(>F)  
## (Intercept)   3798  1  5.3655 0.02327 *
## Group         9763  4  3.4483 0.01215 *
## Full_Length   4823  1  6.8141 0.01091 *
## Deletion       151  1  0.2139 0.64504  
## Residuals    53089 75                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.1)
## [1] 82

# Supplementary Table 12B:
met.sep.2<-lm(BCratio~Group+aa.concat, data=final.fl)
summary(met.sep.2)
## 
## Call:
## lm(formula = BCratio ~ Group + aa.concat, data = final.fl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.740  -1.448   4.764  14.395  43.028 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   48.123     13.942   3.452  0.00101 **
## GroupCOL      12.243     15.225   0.804  0.42440   
## GroupNOR       8.441     14.638   0.577  0.56626   
## GroupUTA      17.877     15.209   1.175  0.24433   
## GroupWES      26.352     16.124   1.634  0.10725   
## aa.concatVM   20.502     14.056   1.459  0.14973   
## aa.concatVV   13.251      9.336   1.419  0.16081   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.55 on 62 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1913, Adjusted R-squared:  0.1131 
## F-statistic: 2.445 on 6 and 62 DF,  p-value: 0.0347
Anova(met.sep.2, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq Df F value   Pr(>F)   
## (Intercept)   9710  1 11.9137 0.001009 **
## Group         3368  4  1.0331 0.397388   
## aa.concat     2070  2  1.2698 0.288089   
## Residuals    50532 62                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.2)
## [1] 69

# Supplementary Table 12C:
met.sep.3<-lm(BCratio~Group+AA148+AA268, data=final.fl)
summary(met.sep.3)
## 
## Call:
## lm(formula = BCratio ~ Group + AA148 + AA268, data = final.fl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.740  -1.448   4.764  14.395  43.028 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   55.373     17.842   3.104  0.00288 **
## GroupCOL      12.243     15.225   0.804  0.42440   
## GroupNOR       8.441     14.638   0.577  0.56626   
## GroupUTA      17.877     15.209   1.175  0.24433   
## GroupWES      26.352     16.124   1.634  0.10725   
## AA148V        13.251      9.336   1.419  0.16081   
## AA268V        -7.250     10.902  -0.665  0.50851   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.55 on 62 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1913, Adjusted R-squared:  0.1131 
## F-statistic: 2.445 on 6 and 62 DF,  p-value: 0.0347
Anova(met.sep.3, type="3")
## Anova Table (Type III tests)
## 
## Response: BCratio
##             Sum Sq Df F value   Pr(>F)   
## (Intercept)   7851  1  9.6324 0.002879 **
## Group         3368  4  1.0331 0.397388   
## AA148         1642  1  2.0145 0.160808   
## AA268          360  1  0.4422 0.508505   
## Residuals    50532 62                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
nobs(met.sep.3)
## [1] 69

FST

FST was calculated using methods and scripts described in Wang et al. 2019.

Linkage disequilibrium

Pairwise LD among SNPs in the BCMA CFR-HIF region was calculated using the software PLINK (Purcell et al. 2007). PLINK output was modified for plotting heatmaps using a custom script (ld2matrix.py).

Nucleotide diversity

Pi for BCMA3 and other genes in the B. stricta genome were calculated using DnaSP, after correcting for possible sequence errors in regions supported by a single read (Python scripts: calculate_per_site_div.lable_singleton.py and calculate_per_window_diversity.masked_singleton.py).

Filter pi data:

# designate observed pi BCMA3 using DnaSP
pi.BCMA3 <- 0.00591 

# Only compare pi to genes with a similar length to BCMA3, since length is correlated with nucleotide diversity (Supplementary Figure 9). Also, filter out genes with pi = 0.
min.num.sites <- 360
nuc.div.OK <- filter(nuc.div, Pi >= 0 & site >= min.num.sites)

# Three main regions of the B. stricta genome have previously been identified as having elevated nucleotide diversity due to historical balancing selection (Wang et al. 2019). We also exclude these regions from comparisons testing for balancing selection on BCMA3, since we already know they have experienced balancing selection.

## by default, assign all genes to the region "background"
nuc.div.OK$reg <- "background"

## flag the regions of elevated polymorphism on LG1, LG5, and LG7 as identified in Wang et al. 2019
nuc.div.OK$reg[nuc.div.OK$start >= 14423807 & 
              nuc.div.OK$end <= 24375392 & nuc.div.OK$LG ==1] <- "LG1"
nuc.div.OK$reg[nuc.div.OK$start >= 471982 &
              nuc.div.OK$end <= 3400076 & nuc.div.OK$LG == 5] <- "LG5"
nuc.div.OK$reg[nuc.div.OK$start >= 451506 &
              nuc.div.OK$end <= 5083168 & nuc.div.OK$LG == 7] <- "LG7"

## filter to focus only on background regions
nuc.div.OK.2 <- filter(nuc.div.OK, reg == "background")
nrow(nuc.div.OK.2) # this is the number of comparable genes retained
## [1] 1689

Compare BCMA to the genome-wide distribution of pi:

# calculate the proportion of comparable genes for which pi > BCMA3
x <- nuc.div.OK.2$Pi > pi.BCMA3
p.val <- mean(x) 
p.val
## [1] 0.02901125

Average pi in comparable genes:

mean(nuc.div.OK.2$Pi)
## [1] 0.001690409

Effect size of BCMA on insect resistance

Get mean and SD of herbivore damage in each transplant environment:

# subset to transplant experiments
field.use<-filter(field, Expt=="transplants")

# mean herbivory for each geno in each env
h.means<-field.use %>% dplyr::select(Geno, Env, LAR2) %>% drop_na() %>% 
  group_by(Geno, Env) %>%
  summarise(Mean_dam=mean(LAR2))

# reshape means so there's one column for each geno
h.means<-spread(h.means, Geno, Mean_dam)

# SD of herbivory for each geno in each env
h.sd<-field.use%>%dplyr::select(Geno, Env, LAR2) %>% drop_na() %>%
  group_by(Env) %>%
  summarise(SD_dam=sd(LAR2))

Calculate effect size:

# vector of differences in genotypic means/SD per environment
## initialize vector
diff.vec<-rep(NA, nrow(h.sd))

## calculate effect size in each environment
i<-1
for(i in 1:length(diff.vec)){
  mean.LL<-h.means$LL[i]
  mean.SS<-h.means$SS[i]
  stdev<-h.sd$SD_dam[i]
  
  std.diff<-abs(mean.LL-mean.SS)/stdev
  diff.vec[i]<-std.diff
}

# calculate mean effect size across environments
mean(diff.vec)
## [1] 0.1715397

# calculate min and max effect sizes
min(diff.vec)
## [1] 0.01260256
max(diff.vec)
## [1] 0.4329623

# view differences by environment
data.frame(cbind(h.means$Env, as.numeric(diff.vec)))
##         X1                 X2
## 1 ALD-2013    0.3853977028103
## 2 ALD-2014 0.0594271726528837
## 3 ALD-2015 0.0854734188355042
## 4 GTH-2016  0.432962264183009
## 5 MAH-2014 0.0126025590252014
## 6 MAH-2015   0.15228075876984
## 7 SCH-2016  0.144273660559002
## 8 SIL-2014 0.0771454958018908
## 9 SIL-2015  0.194294400846806

Record session info

Package versions, etc.

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] geodist_0.0.4    ggplot2_3.3.2    pbkrtest_0.4-8.6 emmeans_1.5.0   
##  [5] car_3.0-9        carData_3.0-4    lmerTest_3.1-2   lme4_1.1-23     
##  [9] Matrix_1.2-18    tidyr_1.1.1      dplyr_1.0.2     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5          mvtnorm_1.1-1       lattice_0.20-41    
##  [4] zoo_1.8-8           digest_0.6.25       plyr_1.8.6         
##  [7] R6_2.4.1            cellranger_1.1.0    evaluate_0.14      
## [10] coda_0.19-3         pillar_1.4.6        rlang_0.4.7        
## [13] curl_4.3            multcomp_1.4-13     readxl_1.3.1       
## [16] rstudioapi_0.11     minqa_1.2.4         data.table_1.13.0  
## [19] nloptr_1.2.2.2      rmarkdown_2.3       labeling_0.3       
## [22] splines_4.0.2       statmod_1.4.34      stringr_1.4.0      
## [25] foreign_0.8-80      munsell_0.5.0       compiler_4.0.2     
## [28] numDeriv_2016.8-1.1 xfun_0.16           pkgconfig_2.0.3    
## [31] htmltools_0.5.0     tidyselect_1.1.0    tibble_3.0.3       
## [34] rio_0.5.16          codetools_0.2-16    crayon_1.3.4       
## [37] withr_2.2.0         MASS_7.3-52         grid_4.0.2         
## [40] nlme_3.1-148        xtable_1.8-4        gtable_0.3.0       
## [43] lifecycle_0.2.0     magrittr_1.5        scales_1.1.1       
## [46] zip_2.0.4           estimability_1.3    stringi_1.4.6      
## [49] farver_2.0.3        ellipsis_0.3.1      generics_0.0.2     
## [52] vctrs_0.3.2         boot_1.3-25         sandwich_2.5-1     
## [55] openxlsx_4.1.5      TH.data_1.0-10      tools_4.0.2        
## [58] forcats_0.5.0       glue_1.4.1          purrr_0.3.4        
## [61] hms_0.5.3           abind_1.4-5         survival_3.2-3     
## [64] yaml_2.2.1          colorspace_1.4-1    knitr_1.29         
## [67] haven_2.3.1