Trait Plasticity based on climate transfer distance
Anoob Prakash
03/03/2021
Methods
family BLUP model
- Mixed effect model run to get the family level blups for use in the MCMC run.
- Garden is set as the fixed effect.
- Family and mBed (5 beds x 3 sites) set as random effect.
- Height growth for a season is used as the fitness proxy to compare the adaptive nature of trait plasticity.
- This helps in explaining whether the change in plasticity for a trait (positive/negative) is adaptive (top right of the plot) or maladaptive (bottom left of the plot).
# Creating mBed
fitness <- fitness %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
# lmer model for family blup
fitness_mod <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), fitness)
# extract the blups for the family from the model
fitness_blup <- ranef(fitness_mod)
fitness_blup <- fitness_blup$Family
fitness_blup <- cbind(Family=rownames(fitness_blup),fitness_blup)
rownames(fitness_blup) <- NULL
names(fitness_blup)[2] <- "fitness"Trait plasticity (as described in the paper)
To estimate phenotypic plasticity for each trait, we ran LMMs (Model IV, Table S1) that included the fixed effect of garden climate (gPC1, based on each garden’s eigenvector score along PC1 of the climate PCA for the year of measurement), and random effects of Bed (1|Bed) and Family with random intercepts and slopes (gPC1|Family). We used the ‘coef()’ function to extract the random slopes for each Family with gPC1, which provides a measure of plasticity to garden climate (Arnold et al. 2019). To test if the magnitude of plasticity reflected differences due to source climate variability, we modeled the absolute value of plasticity as a function of source climate (sPC1), with a random effect of Population (Model V, Table S1).
# setwd("Z:/Anoob/MCMCglmm")
suppressPackageStartupMessages({
# packages
require(MCMCglmm)
require(lmerTest)
require(dplyr)
require(tidyr)
require(data.table)
require(ggplot2)
require(cowplot)
require(plyr)
})
# data
budset_2019 <- read.csv("./trait_data/BudSet_2019.csv")
budset_2020 <- read.csv("./trait_data/BudSet_2020.csv")
budbreak_2020 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv")
heightGrowth_2019 <- read.csv("./trait_data/Growth_2019.csv")
heightGrowth_2020 <- read.csv("./trait_data/Growth_2020.csv")
# Plasticity <- read.csv("./trait_data/Plasticity.csv")
Meta <- read.table("./Exome/RS_Exome_metadata.txt", sep="\t",header=T)
PC_scores <- read.csv("./trait_data/PCA/PCA_score.csv", header=T)
PC_scores <- PC_scores[!is.na(PC_scores$Family),]
# meta data and PC
colnames(Meta)[colnames(Meta)=="Pop"] <- "Population"
Meta$Region <- as.factor(Meta$Region)
Meta$Region <- mapvalues(Meta$Region, from = c("C", "M", "E"), to = c("Core", "Margin", "Edge"));Meta$Region <- factor(Meta$Region,levels = c("Core", "Margin", "Edge"))
Meta <- merge(x=Meta,
y=PC_scores[,-c(1,2,3,5,17)],
by.x="Family",
by.y="Family")
Meta <- Meta[,c("Population","Family","Region","Latitude","Longitude","Elevation","PC1","PC2","PC3","PC4","PC5","PC6",
"PC7","PC8","PC9","PC10","PC11")]Garden climate
Garden climate
Garden_clim <- read.table("./ClimateNA/climNA_gardens_2019-20.txt", header=T)
Garden_clim## garden year MSP DD_0 DD5 DD18 PAS eFFP RH TD CMD MAR EXT PET
## 1 MD 2019 469 341 2649 368 94 297 66 25.4 212 15.2 36.4 866.3545
## 2 MD 2020 492 224 2548 362 44 299 67 23.7 131 14.8 36.4 872.7068
## 3 NC 2019 478 55 3873 803 5 324 64 20.4 277 17.1 37.7 1147.0480
## 4 NC 2020 723 52 3747 698 4 321 66 20.9 166 16.2 37.8 1094.4220
## 5 VT 2019 471 678 2451 407 172 293 67 31.9 147 13.9 39.6 746.0178
## 6 VT 2020 403 385 2633 483 74 294 67 28.8 226 14.2 38.9 793.8243selVars1 <- c('DD_0','DD18','MAR','PAS','MSP','RH','EXT','CMD','TD','eFFP','PET')
colnames(Garden_clim)[1] <- "Garden"
colnames(Garden_clim)[2] <- "Year"# create ID for garden means
# garden_identity <- Garden_climateNA1[,1]
# create ID for garden
garden_identity <- Garden_clim[,1:2]
# selected clim of gardens
data_p1 <- Garden_clim[,selVars1]garden_id <- garden_identity
colnames(garden_id)[1] <- "Garden"
colnames(garden_id)[2] <- "Year"
data_pca <- prcomp(Garden_clim[,-c(1:5)],scale. = T)
# Extract PC axes for plotting
PCAvalues <- data.frame(Garden = Garden_clim$Garden, Year = Garden_clim$Year,
data_pca$x)
PCAvalues$Year <- as.factor(PCAvalues$Year)
# Extract loadings of the variables
PCAloadings <- data.frame(Variables = rownames(data_pca$rotation), data_pca$rotation)
# write.csv(PCAvalues, row.names = F, "./trait_data/PCA/PCA_garden_score.csv")
# write.csv(PCAloadings, row.names = F, "./trait_data/PCA/PCA_garden_loadings.csv")
PCAvalues_mod <- PCAvalues
PCAvalues_mod$Garden[PCAvalues_mod$Garden=="VT"] <- "Vermont"
PCAvalues_mod$Garden[PCAvalues_mod$Garden=="MD"] <- "Maryland"
PCAvalues_mod$Garden[PCAvalues_mod$Garden=="NC"] <- "North_Carolina"Trait plasticity
Bud break 2020
- cGDD
budbreak_2020 <- merge(x=budbreak_2020,y=PCAvalues_mod[PCAvalues_mod$Year=="2020",], all.x=TRUE)
budbreak_2020$mBed <- paste0(budbreak_2020$Garden,"_",budbreak_2020$Bed)
#plasticity model
budbreak_PC1_mod <- lmer(data = budbreak_2020, na.action = na.omit,
cGDD ~ PC1 + (1|mBed) + (PC1|Family))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00453019 (tol = 0.002, component 1)summary(budbreak_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: cGDD ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: budbreak_2020
##
## REML criterion at convergence: 47653.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.6484 -0.5404 0.0786 0.6205 4.0111
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 305.24 17.471
## PC1 34.53 5.876 -0.97
## mBed (Intercept) 1081.09 32.880
## Residual 1144.54 33.831
## Number of obs: 4757, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 643.013 8.570 13.308 75.026 < 2e-16 ***
## PC1 -27.920 5.402 13.073 -5.169 0.000177 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.061
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00453019 (tol = 0.002, component 1)budbreak_PC1 <- ranef(budbreak_PC1_mod)
budbreak_PC1 <- budbreak_PC1$Family
budbreak_PC1 <- cbind(Family=rownames(budbreak_PC1),budbreak_PC1)
rownames(budbreak_PC1) <- NULL
budbreak_PC1$Family <- as.factor(budbreak_PC1$Family)
colnames(budbreak_PC1)[colnames(budbreak_PC1)=="PC1"] <- "BudBreak_2020"
budbreak_PC1 <- budbreak_PC1[,-2]
head(budbreak_PC1) ## Family BudBreak_2020
## 1 AB_05 -5.8545602
## 2 AB_08 2.5435443
## 3 AB_12 -3.8914145
## 4 AB_16 -4.1721421
## 5 AB_18 -0.8298989
## 6 ALB_01 -3.2717690# converting the negative values to absolute values to make plasticity values comparable across the traits
bb2020x_PC1 <- coef(budbreak_PC1_mod)
# hist(testx_PC1$PC1)
bb2020x_PC1 <- bb2020x_PC1$Family
bb2020x_PC1 <- cbind(Family=rownames(bb2020x_PC1),bb2020x_PC1)
rownames(bb2020x_PC1) <- NULL
bb2020x_PC1$abs <- abs(bb2020x_PC1$PC1)
hist(bb2020x_PC1$PC1)
hist(bb2020x_PC1$abs)
bb2020x_PC1 <- bb2020x_PC1[,-2]
colnames(bb2020x_PC1)[3] <- "BudBreak_cGDD"
bb2020x_PC1 <- bb2020x_PC1[,-2] - DOY
budbreak_2020_DOY <- merge(x=budbreak_2020,y=PCAvalues_mod[PCAvalues_mod$Year=="2020",], all.x=TRUE)
budbreak_2020_DOY$mBed <- paste0(budbreak_2020_DOY$Garden,"_",budbreak_2020_DOY$Bed)
#plasticity model
budbreak_DOY_PC1_mod <- lmer(data = budbreak_2020_DOY, na.action = na.omit,
BudBreak ~ PC1 + (1|mBed) + (PC1|Family))
summary(budbreak_DOY_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudBreak ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: budbreak_2020_DOY
##
## REML criterion at convergence: 30004.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.9704 -0.5442 0.0633 0.6057 4.0191
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 7.5674 2.7509
## PC1 0.5274 0.7262 -0.95
## mBed (Intercept) 0.9968 0.9984
## Residual 28.3259 5.3222
## Number of obs: 4757, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 130.7792 0.3089 22.3288 423.31 <2e-16 ***
## PC1 8.2513 0.1754 14.5242 47.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.155budbreak_DOY_PC1 <- ranef(budbreak_DOY_PC1_mod)
budbreak_DOY_PC1 <- budbreak_DOY_PC1$Family
budbreak_DOY_PC1 <- cbind(Family=rownames(budbreak_DOY_PC1),budbreak_DOY_PC1)
rownames(budbreak_DOY_PC1) <- NULL
budbreak_DOY_PC1$Family <- as.factor(budbreak_DOY_PC1$Family)
colnames(budbreak_DOY_PC1)[colnames(budbreak_DOY_PC1)=="PC1"] <- "BudBreak_2020_DOY"
budbreak_DOY_PC1 <- budbreak_DOY_PC1[,-2]
head(budbreak_DOY_PC1) ## Family BudBreak_2020_DOY
## 1 AB_05 -0.6581685
## 2 AB_08 0.2415947
## 3 AB_12 -0.4483200
## 4 AB_16 -0.4263305
## 5 AB_18 -0.1692816
## 6 ALB_01 -0.3407461# converting the negative values to absolute values to make plasticity values comparable across the traits
bb2020_DOYx_PC1 <- coef(budbreak_DOY_PC1_mod)
# hist(testx_PC1$PC1)
bb2020_DOYx_PC1 <- bb2020_DOYx_PC1$Family
bb2020_DOYx_PC1 <- cbind(Family=rownames(bb2020_DOYx_PC1),bb2020_DOYx_PC1)
rownames(bb2020_DOYx_PC1) <- NULL
bb2020_DOYx_PC1$abs <- abs(bb2020_DOYx_PC1$PC1)
hist(bb2020_DOYx_PC1$PC1)
hist(bb2020_DOYx_PC1$abs)
bb2020_DOYx_PC1 <- bb2020_DOYx_PC1[,-2]
colnames(bb2020_DOYx_PC1)[3] <- "BudBreak_DOY"
bb2020_DOYx_PC1 <- bb2020_DOYx_PC1[,-2] Bud set 2019
budset_2019 <- merge(x=budset_2019,y=PCAvalues_mod[PCAvalues_mod$Year=="2019",], all.x=TRUE)
budset_2019$mBed <- paste0(budset_2019$Garden,"_",budset_2019$Bed)
# plasticity model
budset_2019_PC1_mod <- lmer(data = budset_2019, na.action = na.omit,
BudSet ~ PC1 + (1|mBed) + (PC1|Family))
summary(budset_2019_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudSet ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: budset_2019
##
## REML criterion at convergence: 44149.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0560 -0.7264 0.2055 0.7190 2.3287
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 107.581 10.372
## PC1 2.326 1.525 -0.06
## mBed (Intercept) 20.846 4.566
## Residual 645.778 25.412
## Number of obs: 4685, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 256.4645 1.3623 19.1099 188.265 < 2e-16 ***
## PC1 3.1744 0.4383 14.2106 7.243 3.92e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 0.016budset_2019_PC1 <- ranef(budset_2019_PC1_mod)
budset_2019_PC1 <- budset_2019_PC1$Family
budset_2019_PC1 <- cbind(Family=rownames(budset_2019_PC1),budset_2019_PC1)
rownames(budset_2019_PC1) <- NULL
budset_2019_PC1$Family <- as.factor(budset_2019_PC1$Family)
colnames(budset_2019_PC1)[colnames(budset_2019_PC1)=="PC1"] <- "BudSet_2019"
budset_2019_PC1 <- budset_2019_PC1[,-2]
bs2019x_PC1 <- coef(budset_2019_PC1_mod)
# hist(testx_PC1$PC1)
bs2019x_PC1 <- bs2019x_PC1$Family
bs2019x_PC1 <- cbind(Family=rownames(bs2019x_PC1),bs2019x_PC1)
rownames(bs2019x_PC1) <- NULL
bs2019x_PC1$abs <- abs(bs2019x_PC1$PC1)
hist(bs2019x_PC1$PC1)
hist(bs2019x_PC1$abs)
bs2019x_PC1 <- bs2019x_PC1[,-c(2:3)]
colnames(bs2019x_PC1)[2] <- "BudSet_2019" Bud set 2020
budset_2020 <- merge(x=budset_2020,y=PCAvalues_mod[PCAvalues_mod$Year=="2020",], all.x=TRUE)
budset_2020$mBed <- paste0(budset_2020$Garden,"_",budset_2020$Bed)
# plasticity model
budset_2020_PC1_mod <- lmer(data = budset_2020, na.action = na.omit,
BudSet ~ PC1 + (1|mBed) + (PC1|Family))
summary(budset_2020_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: BudSet ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: budset_2020
##
## REML criterion at convergence: 37189.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1937 -0.6040 -0.2718 0.5537 3.2311
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 27.376 5.232
## PC1 1.165 1.080 0.14
## mBed (Intercept) 102.524 10.125
## Residual 319.734 17.881
## Number of obs: 4282, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 218.488 2.649 13.328 82.489 < 2e-16 ***
## PC1 5.832 1.668 13.061 3.495 0.00392 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.055budset_2020_PC1 <- ranef(budset_2020_PC1_mod)
budset_2020_PC1 <- budset_2020_PC1$Family
budset_2020_PC1 <- cbind(Family=rownames(budset_2020_PC1),budset_2020_PC1)
rownames(budset_2020_PC1) <- NULL
budset_2020_PC1$Family <- as.factor(budset_2020_PC1$Family)
colnames(budset_2020_PC1)[colnames(budset_2020_PC1)=="PC1"] <- "BudSet_2020"
budset_2020_PC1 <- budset_2020_PC1[,-2]
bs2020x_PC1 <- coef(budset_2020_PC1_mod)
# hist(testx_PC1$PC1)
bs2020x_PC1 <- bs2020x_PC1$Family
bs2020x_PC1 <- cbind(Family=rownames(bs2020x_PC1),bs2020x_PC1)
rownames(bs2020x_PC1) <- NULL
bs2020x_PC1$abs <- abs(bs2020x_PC1$PC1)
hist(bs2020x_PC1$PC1)
hist(bs2020x_PC1$abs)
bs2020x_PC1 <- bs2020x_PC1[,-c(2:3)]
colnames(bs2020x_PC1)[2] <- "BudSet_2020"Growth 2019
growth_2019 <- merge(x=heightGrowth_2019,y=PCAvalues_mod[PCAvalues_mod$Year=="2019",], all.x=TRUE)
growth_2019$mBed <- paste0(growth_2019$Garden,"_",growth_2019$Bed)
# plasticity model
growth_2019_PC1_mod <- lmer(data = growth_2019, na.action = na.omit,
Growth ~ PC1 + (1|mBed) + (PC1|Family))## boundary (singular) fit: see ?isSingularsummary(growth_2019_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: growth_2019
##
## REML criterion at convergence: 30471
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3461 -0.6843 -0.0734 0.5849 7.7062
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 7.36158 2.7132
## PC1 0.01984 0.1408 1.00
## mBed (Intercept) 2.96525 1.7220
## Residual 35.26015 5.9380
## Number of obs: 4681, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.0434 0.4772 15.9954 23.142 1e-13 ***
## PC1 0.1561 0.1575 13.1595 0.991 0.34
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 0.041
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular# ranef
growth_2019_PC1 <- ranef(growth_2019_PC1_mod)
growth_2019_PC1 <- growth_2019_PC1$Family
growth_2019_PC1 <- cbind(Family=rownames(growth_2019_PC1),growth_2019_PC1)
rownames(growth_2019_PC1) <- NULL
growth_2019_PC1$Family <- as.factor(growth_2019_PC1$Family)
colnames(growth_2019_PC1)[colnames(growth_2019_PC1)=="PC1"] <- "Growth_2019"
growth_2019_PC1 <- growth_2019_PC1[,-2]
# coef
gw2019x_PC1 <- coef(growth_2019_PC1_mod)
gw2019x_PC1 <- gw2019x_PC1$Family
gw2019x_PC1 <- cbind(Family=rownames(gw2019x_PC1),gw2019x_PC1)
rownames(gw2019x_PC1) <- NULL
gw2019x_PC1$abs <- abs(gw2019x_PC1$PC1)
hist(gw2019x_PC1$PC1)
hist(gw2019x_PC1$abs)
gw2019x_PC1 <- gw2019x_PC1[,-c(2:3)]
colnames(gw2019x_PC1)[2] <- "Growth_2019" Growth 2020
growth_2020 <- merge(x=heightGrowth_2020,y=PCAvalues_mod[PCAvalues_mod$Year=="2020",], all.x=TRUE)
growth_2020$mBed <- paste0(growth_2020$Garden,"_",growth_2020$Bed)
# plasticity model
growth_2020_PC1_mod <- lmer(data = growth_2020, na.action = na.omit,
Growth ~ PC1 + (1|mBed) + (PC1|Family))
summary(growth_2020_PC1_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ PC1 + (1 | mBed) + (PC1 | Family)
## Data: growth_2020
##
## REML criterion at convergence: 30330.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7661 -0.5754 -0.0754 0.4788 5.2978
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Family (Intercept) 11.0872 3.3297
## PC1 0.1973 0.4442 -0.86
## mBed (Intercept) 1.7268 1.3141
## Residual 45.6260 6.7547
## Number of obs: 4475, groups: Family, 334; mBed, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 16.2121 0.3989 20.5972 40.645 < 2e-16 ***
## PC1 0.8208 0.2257 13.2012 3.637 0.00294 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PC1 -0.089growth_2020_PC1 <- ranef(growth_2020_PC1_mod)
growth_2020_PC1 <- growth_2020_PC1$Family
growth_2020_PC1 <- cbind(Family=rownames(growth_2020_PC1),growth_2020_PC1)
rownames(growth_2020_PC1) <- NULL
growth_2020_PC1$Family <- as.factor(growth_2020_PC1$Family)
colnames(growth_2020_PC1)[colnames(growth_2020_PC1)=="PC1"] <- "Growth_2020"
growth_2020_PC1 <- growth_2020_PC1[,-2]
gw2020x_PC1 <- coef(growth_2020_PC1_mod)
# hist(testx_PC1$PC1)
gw2020x_PC1 <- gw2020x_PC1$Family
gw2020x_PC1 <- cbind(Family=rownames(gw2020x_PC1),gw2020x_PC1)
rownames(gw2020x_PC1) <- NULL
gw2020x_PC1$abs <- abs(gw2020x_PC1$PC1)
hist(gw2020x_PC1$PC1)
hist(gw2020x_PC1$abs)
gw2020x_PC1 <- gw2020x_PC1[,-c(2:3)]
colnames(gw2020x_PC1)[2] <- "Growth_2020"Plasticity
Plasticity <- Reduce(function(x, y) merge(x, y, all=TRUE),
list(Meta,bb2020x_PC1,bb2020_DOYx_PC1,bs2019x_PC1,bs2020x_PC1,gw2019x_PC1,gw2020x_PC1))
Plasticity <- Plasticity[!is.na(Plasticity$BudBreak_cGDD),]
# write.csv(Plasticity, "./trait_data/Plasticity.csv")