###################################################
#### Do species from colder regions have higher acclimation potential than species from warmer environments? 
###################################################

#Import dataset
# Import datatable and phylogenetic tree

sp_means= read.csv("sp_means_bio6_min_temp_se.csv")
tree2=readRDS("tree_freeze_2.rds")



######################
#### Phylogeny #######
######################


##### Matching phylogenetic tree with datatable

name.check(tree2,sp_means, data.names=sp_means$species)# check that name in the tree matches with the data

matches<-match(tree2$tip.label,sp_means$species, nomatch=0)
not<-subset(tree2$tip.label, matches ==0)
t2 <-drop.tip(tree2, not)


name.check(t2,sp_means, data.names=sp_means$species)# check that name in the tree matches with the data


## Check whether they are lined up correctly
sp_means$species == t2$tip.label

#Unsurprisingly, not all of the species are in their correct places; we will have to reorder the data frame. Here is one way to do it.

sp_means <- sp_means[match(t2$tip.label, sp_means$species), ]
sp_means$species == t2$tip.label #check again if they are linedup



################################
#### Phylogenetic regression ###
################################

#Run Phylogenetic regression
# min min temperature
d=gls(diff_mean~min_min_temp, data=sp_means,
      correlation=corBrownian(phy=t2), method="ML")

summary(d) #p-value= 0.001

# bio 6 
d=gls(diff_mean~mean_bio6c, data=sp_means,
      correlation=corBrownian(phy=t2), method="ML")

summary(d) #p-value= 0.0009

###### PLOT ######


#Plot
colors= c("Deciduous"="tan1","Evergreen"="darkorchid4")
d=ggplot(data=sp_means, aes(x=min_min_temp, y=diff_mean, color= LP)) + geom_point(size = 3.5, aes(shape=Section))+
  scale_color_manual(values=colors)+
  geom_smooth(method = "lm", se = TRUE, color='black')+
  theme_bw()+
  theme(plot.background = element_blank()
        ,panel.grid.major = element_blank()
        ,panel.grid.minor = element_blank())+
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=16,face="bold"))+
  xlab("Minumum temperature of coldest month (ºC)") + ylab("Acclimation potential")+
  annotate(geom="text", x=-5, y=43, label="p-value= 0.001",fontface="italic",
           color="black", size = 4.5)

print(d + scale_shape_manual(values = c(17,16,15))) # To set the symbols manually.