#Adapt the working directory here
#(careful, there is one "setwd" line prior to each figure-associated section of code)
setwd("")

#calling the abundance data-------------
load("abunselectaxa.RData")
library(FD)

#retrieving the trait data---------
traits=list()
traits$Zoo=read.csv("Zoo_Traits_NSRS.csv",row.names = 1,na.strings="")
traits$Phyto=read.csv("Phyto_Traits_NSRS.csv",row.names = 1,na.strings="")
#log of size-related traits
traits$Zoo[,1:2]=log(traits$Zoo[,1:2])
traits$Phyto[,1:2]=log(traits$Phyto[,1:2])

#parameters for the computation of Gower distances, most importantly, weights attributed to each traits (also allows to parametrize fuzzy traits)
#CARFUL! THE NAMES MUST BE THE SAME AS IN TRAITS!!!
gowerparam=list()
gowerparam$Zoo$weights=c(0.5,0.5,1,0.5,0.5,0.5,0.5,1)
gowerparam$Phyto$weights=c(0.5,0.5,1,1/3,1/3,1/3,1,1,0.5,0.5)

#CHECK IF ALL TRAITS ARE AVAILABLE FOR ALL SPECIES
Zookept=read.csv("Zookept.csv",row.names=1)
traits$Zoo[Zookept$x,]
traits$Zoo=traits$Zoo[Zookept$x,]

Phytokept=read.csv("Phytokept.csv",row.names=1)
traits$Phyto[Phytokept$x,]
traits$Phyto=traits$Phyto[Phytokept$x,]

#Compute trait space, centroids and PCAs on centroids--------------
FTBA=list()

for (i in c("Phyto","Zoo")){
  FTBA[[i]]$Traits=traits[[i]]
  FTBA[[i]]$dist=gowdis(FTBA[[i]]$Traits,w=gowerparam[[i]]$weights)
  FTBA[[i]]$pcoa=pcoa(FTBA[[i]]$dist,correction="cailliez")
  FTBA[[i]]$TraitSpace=FTBA[[i]]$pcoa$vectors.cor
  
  for (j in names(abunselectaxa)){
    truc=dbFD(FTBA[[i]]$TraitSpace,abunselectaxa[[j]][[i]]$Months[,rownames(FTBA[[i]]$TraitSpace)],stand.x=F,corr="none",calc.FRic=F,calc.FDiv=F)
    FTBA[[i]][[j]]$centroids=truc$CWM
    
    FTBA[[i]][[j]]$PCA=dudi.pca(FTBA[[i]][[j]]$centroids,scale=F,center=T,scannf=F,nf=dim(FTBA[[i]][[j]]$centroids)[2])
  }
}

save(FTBA,file="FTBA NSRS.RData")#finished product of this section


###################################################################
#FIGURE 3:ecological meaning of Fcomp1 (i.e. PC1 from the PCA on centroids)------
###################################################################
#NOTE: the projection of taxa for ecological interpretation is performed here----

setwd("")
load("FTBA NSRS.RData")

#calling packages and defining some graphical characteristics---------
library(FD)
library(vioplot)
library(colourvalues)
library(pals)
library(Ternary)
library(png)

#some colors for the figure
RGBperso=function(a,b,c){
  cols=rgb((1-a)*0.9,(1-b)*0.9,(1-c)*0.9)
  return(cols)
}

bicol=c("dodgerblue3","darkorange2")
quadcol=c("dodgerblue3","darkorange2","forestgreen","darkorchid3")

#calling pictures to illustrate the axes
coppic=readPNG("PictureCop.png")
dinopic=readPNG("PictureDino.png")
diatopic=readPNG("PictureChaeto.png")
meropic=readPNG("PicturePluteus.png")

#Start-------------
pdf("Fig. 3 FComp1 and traits.pdf",width=6,height=8)
layout(rbind(matrix(c(8,2,3,4,5,6,7,8,8,1,1,1,1,1,1,8),8,2),matrix(c(7,2,3,4,5,6,7,7,1,1,1,1,1,7),7,2)+8),widths=c(1,2),heights=c(0.5,1,1,1,1,1,1,1.2,0.5,1,1,1,1,1,1.2))

#Phyto---------
par(mar=c(0,0,0,1))
plot(0,0,col="white",xlim=c(-0.1,1.05),ylim=c(0,6),xaxt="n",yaxt="n",ylab="",xlab="",yaxs="i",xaxs="i")
pos=c(0,0.25,0.5,0.75)
names(pos)=c("C02","C01","D02","D01")
lwd=2
cex=1.8

for (i in c("C02","C01","D02","D01")){
  #Here is the projection:
  proj=suprow(FTBA$Phyto[[i]]$PCA,FTBA$Phyto$TraitSpace)
  truc=proj$lisup[,1]
  names(truc)=rownames(proj$lisup)
  #"truc" now contains the coordinates of the taxa on PC1
  
  #This condition flips (or not) the axes of the different PCA (performed for the different standard areas)
  #to improve comparability, taking the axis 1 of the trait space as a reference
  #This does not impact the analysis as the direction of PCA axes (i.e. where is the negative side and shere is the positive side)
  #is defined randomly 
  if(cor(truc,FTBA$Phyto$TraitSpace[,1])<0){
    truc=-truc
  }
  truc=(truc-min(truc))/(max(truc)-min(truc))#this line scales the axis between 0 and 1 to improve visualisation (as described in the paper)
  truc2=FTBA$Phyto$Traits
  
  #THIS IS A SERIE OF COMMANDS TO PLOT EACH TRAIT:
  #size
  segments(x0=truc,y0=6-pos[i],y1=6-pos[i]-0.25,col=colour_values(apply(truc2[,1:2],1,mean),"viridis",alpha=0.5),lwd=lwd,lend=1)
  #cover
  segments(x0=truc,y0=5-pos[i],y1=5-pos[i]-0.25,col=rgb(t(col2rgb(bicol[as.factor(truc2[,3])])/255),alpha=0.5),lwd=lwd,lend=1)
  #shape
  segments(x0=truc,y0=4-pos[i],y1=4-pos[i]-0.25,col=rgb((t(255-col2rgb(rgb(truc2[4:6])))/255)*0.9,alpha=0.5),lwd=lwd,lend=1)
  #coloniality
  segments(x0=truc,y0=3-pos[i],y1=3-pos[i]-0.25,col=rgb(t(col2rgb(bicol[as.factor(truc2[,7])])/255),alpha=0.5),lwd=lwd,lend=1)
  #motility
  segments(x0=truc,y0=2-pos[i],y1=2-pos[i]-0.25,col=rgb(t(col2rgb(bicol[as.factor(truc2[,8])])/255),alpha=0.5),lwd=lwd,lend=1)
  #nutrition
  segments(x0=truc,y0=1-pos[i],y1=1-pos[i]-0.25,col=rgb(t(col2rgb(colour_values(truc2[,9],"magenta2green"))/255)*0.8,alpha=0.5),lwd=lwd,lend=1)
}
abline(h=seq(0,6,0.25),col="grey",lwd=1,xpd=F)
abline(h=0:6,lwd=1,xpd=F)
rect(ybottom=0,ytop=6,xleft=-0.1,xright=1.05,xpd=NA,lwd=1,col=NA)

text(c("C02","C01","D02","D01"),x=-0.05,y=seq(5.875,0.125,-0.25),font=4)
axis(side=1,font=2,las=1)
mtext("FComp1",1,xpd=NA,line=3,font=2,cex=0.8,at=0.5)
text(x=c(0,1),y=-0.85,pos=c(4,2),c("Non-motile\nautotrophs","Motile\nmixotrophs"),xpd=NA,font=4,col="grey30",cex=1)
polygon(x=c(0,0,-0.07),y=c(-1,-0.6,-0.8),col="black",xpd=NA)
polygon(x=c(1,1,1.07),y=c(-1,-0.6,-0.8),col="black",xpd=NA)
rasterImage(diatopic,xleft=0.2,xright=0.35,ybottom=-1,ytop=-0.6,xpd=NA)
rasterImage(dinopic,xleft=0.7,xright=0.8,ybottom=-1,ytop=-0.6,xpd=NA)

#legends phyto------------
text(c("Cell-size:","Cover:","Shape:","Coloniality:","Motility:","Nutritive strategy:"),x=-0.7,y=6:1-0.1,pos=4,xpd=NA,font=4,cex=1)

par(mar=c(1.5,6.5,1.5,5.5))
barplot(rep(1,100),col=colour_values(1:100,"viridis"),space=0,axes=F,xpd=T,border=NA,xaxs="i")
rect(xleft=0,xright=100,ybottom=0,ytop=1)
text(c(1,100),c(0.5,0.5),pos=c(2,4),c("Small","Large"),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(0.5,4,0.5,3))
plot(NA,NA,axes=F,xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",ylab="")
rect(ybottom=c(1.25,6.25),ytop=c(3.75,8.75),xleft=1,xright=3,col=bicol,xpd=NA)
text(c(3,3),c(2.5,7.5),c("Organic","Siliceous"),font=4,xpd=NA,col="grey30",pos=4,cex=0.9)

par(xpd=T,mar=c(1,3,1,2))
TernaryPlot(ylim=c(0,sqrt(3)/2),axis.tick=F,axis.labels=F,grid.lines=0,padding=0)
ColourTernary(TernaryPointValues(RGBperso),spectrum=NULL)
polygon(c(0,-0.5,0.5),c(sqrt(3)/2,0,0),border="black",col=NA)
text(c(0,-0.5,0.5),c(sqrt(3)/2,0,0),c("Round ","Irregular","Long"),pos=c(2,2,4),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(0.5,4,0.5,3))
plot(NA,NA,axes=F,xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",ylab="")
rect(ybottom=c(1.25,6.25),ytop=c(3.75,8.75),xleft=1,xright=3,col=bicol,xpd=NA)
text(c(3,3),c(2.5,7.5),c("Colonial","Solitary"),font=4,xpd=NA,col="grey30",pos=4,cex=0.9)

par(mar=c(0.5,4,0.5,3))
plot(NA,NA,axes=F,xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",ylab="")
rect(ybottom=c(1.25,6.25),ytop=c(3.75,8.75),xleft=1,xright=3,col=bicol,xpd=NA)
text(c(3,3),c(2.5,7.5),c("Motile","non-Motile"),font=4,xpd=NA,col="grey30",pos=4,cex=0.9)

par(mar=c(1.5,6.5,1.5,5.5))
barplot(rep(1,100),col=rgb(t(col2rgb(colour_values(1:100,"magenta2green"))/255)*0.8),space=0,axes=F,xpd=T,border=NA,xaxs="i")
rect(xleft=0,xright=100,ybottom=0,ytop=1)
text(c(1,100),c(0.5,0.5),pos=c(2,4),c("Heterotrophic","Phototrophic"),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(0,0,0,0))
plot(NA,NA,axes=F,ylim=c(0,8),xlim=c(0,10),ylab="",xlab="",xaxs="i",yaxs="i")
text(x=0,y=7.75,pos=4,"a) Phytoplankton",font=2,xpd=NA,cex=1.3)
#dev.off()

#Zoo--------------
par(mar=c(0,0,0,1))
plot(0,0,col="white",xlim=c(-0.1,1.05),ylim=c(0,5),xaxt="n",yaxt="n",ylab="",xlab="",yaxs="i",xaxs="i")
pos=c(0,0.25,0.5,0.75)
names(pos)=c("C02","C01","D02","D01")
lwd=2
cex=1.8

for (i in c("C02","C01","D02","D01")){
  #Here is the projection:
  proj=suprow(FTBA$Zoo[[i]]$PCA,FTBA$Zoo$TraitSpace)
  truc=proj$lisup[,1]
  names(truc)=rownames(proj$lisup)
  #"truc" now contains the coordinates of the taxa on PC1
  
  #This condition flips (or not) the axes of the different PCA (performed for the different standard areas)
  #to improve comparability, taking the axis 1 of the trait space as a reference
  #This does not impact the analysis as the direction of PCA axes (i.e. where is the negative side and shere is the positive side)
  #is defined randomly 
  if(cor(truc,FTBA$Zoo$TraitSpace[,1])>0){
    truc=-truc
  }
  truc=(truc-min(truc))/(max(truc)-min(truc))#this line scales the axis between 0 and 1 to improve visualisation (as described in the paper)
  truc2=FTBA$Zoo$Traits
  
  #THIS IS A SERIE OF COMMANDS TO PLOT EACH TRAIT:
  #size
  segments(x0=truc,y0=5-pos[i],y1=5-pos[i]-0.25,col=colour_values(apply(truc2[,1:2],1,mean),"viridis",alpha=0.5),lwd=lwd,lend=1)
  #life-cycle
  segments(x0=truc,y0=4-pos[i],y1=4-pos[i]-0.25,col=rgb(t(col2rgb(bicol[as.factor(truc2[,3])])/255),alpha=0.5),lwd=lwd,lend=1)
  #feeding mode
  segments(x0=truc,y0=3-pos[i],y1=3-pos[i]-0.25,col=rgb(t(col2rgb(colour_values(truc2[,4],"magenta2green"))/255)*0.8,alpha=0.5),lwd=lwd,lend=1)
  #trophic regime
  segments(x0=truc,y0=2-pos[i],y1=2-pos[i]-0.25,col=rgb(t(col2rgb(colour_values(c(truc2[,6],1,0),"magenta2green"))/255)*0.8,alpha=0.5)[1:dim(truc2)[1]],lwd=lwd,lend=1)
  #cover
  segments(x0=truc,y0=1-pos[i],y1=1-pos[i]-0.25,col=rgb(t(col2rgb(quadcol[as.factor(truc2[,8])])/255),alpha=0.5),lwd=lwd,lend=1)
}

abline(h=seq(0,5,0.25),col="grey",lwd=1,xpd=F)
abline(h=0:5,lwd=1,xpd=F)
rect(ybottom=0,ytop=5,xleft=-0.1,xright=1.05,xpd=NA,lwd=1,col=NA)

text(c("C02","C01","D02","D01"),x=-0.05,y=seq(4.875,0.125,-0.25),font=4)
axis(side=1,font=2,las=1)
mtext("FComp1",1,xpd=NA,line=3,font=2,at=0.5,cex=0.8)
text(x=c(0,1),y=-0.85,pos=c(4,2),c("Crustacean\nholoplankton","Meroplankton"),xpd=NA,font=4,col="grey30",cex=1)
polygon(x=c(0,0,-0.07),y=c(-1,-0.6,-0.8),col="black",xpd=NA)
polygon(x=c(1,1,1.07),y=c(-1,-0.6,-0.8),col="black",xpd=NA)
rasterImage(coppic,xleft=0.22,xright=0.32,ybottom=-1,ytop=-0.6,xpd=NA)
rasterImage(meropic,xleft=0.68,xright=0.75,ybottom=-1,ytop=-0.6,xpd=NA)

#legends zoo-------------------------
text(c("Body-size:","Life-cycle:","Feeding-mode:","Trophic regime:","Cover / Support:"),x=-0.7,y=5:1-0.1,pos=4,xpd=NA,font=4,cex=1)


par(mar=c(1.5,6.5,1.5,5.5))
barplot(rep(1,100),col=colour_values(1:100,"viridis"),space=0,axes=F,xpd=T,border=NA,xaxs="i")
rect(xleft=0,xright=100,ybottom=0,ytop=1)
text(c(1,100),c(0.5,0.5),pos=c(2,4),c("Small","Large"),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(0.5,4,0.5,3))
plot(NA,NA,axes=F,xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",ylab="")
rect(ybottom=c(1.25,6.25),ytop=c(3.75,8.75),xleft=1,xright=3,col=bicol,xpd=NA)
text(c(3,3),c(2.5,7.5),c("Holoplankton","Meroplankton"),font=4,xpd=NA,col="grey30",pos=4,cex=0.9)

par(mar=c(1.5,6.5,1.5,5.5))
barplot(rep(1,100),col=rgb(t(col2rgb(colour_values(1:100,"magenta2green"))/255)*0.8),space=0,axes=F,xpd=T,border=NA,xaxs="i")
rect(xleft=0,xright=100,ybottom=0,ytop=1)
text(c(1,100),c(0.5,0.5),pos=c(2,4),c("Active","Passive"),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(1.5,6.5,1.5,5.5))
barplot(rep(1,100),col=rgb(t(col2rgb(colour_values(1:100,"magenta2green"))/255)*0.8),space=0,axes=F,xpd=T,border=NA,xaxs="i")
rect(xleft=0,xright=100,ybottom=0,ytop=1)
text(c(1,100),c(0.5,0.5),pos=c(2,4),c("Carnivore","Herbivore"),font=4,xpd=NA,col="grey30",cex=0.9)

par(mar=c(0.5,4,0.5,3))
plot(NA,NA,axes=F,xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",ylab="")
rect(ybottom=c(1.25,6.25,1.25,6.25),ytop=c(3.75,8.75,3.75,8.75),xleft=c(-3.5,-3.5,7.5,7.5),xright=c(-1.5,-1.5,9.5,9.5),col=quadcol,xpd=NA)
text(c(-1.5,-1.5,9.5,9.5),c(2.5,7.5,2.5,7.5),c("Exoskeleton","Internal skeleton","Naked","Shell"),font=4,xpd=NA,col="grey30",pos=4,cex=0.9)

par(mar=c(0,0,0,0))
plot(NA,NA,axes=F,ylim=c(0,7),xlim=c(0,10),ylab="",xlab="",xaxs="i",yaxs="i")
text(x=0,y=6.75,pos=4,"b) Zooplankton",font=2,xpd=NA,cex=1.3)
#dev.off()

#end---------
dev.off()

##############################################################
#FIGURE 4: Seasonality (+ Supplementary Figure: seasonality in PC2)----------
##############################################################
setwd("")
load("FTBA NSRS.RData")

#calling packages, and defining some graphical characteristics-------
library(FD)
library(png)

Yearmin=1958#When do we start to plot?
Yearmax=2018#When do we stop plotting?

#This will be the plotting axis (time)
x=Yearmin:Yearmax
x=rep(x,each=12)+(((1:12)/12)-(1/24))

#pictures to illustrate the axes
coppic=readPNG("PictureCop.png")
dinopic=readPNG("PictureDino.png")
diatopic=readPNG("PictureChaeto.png")
meropic=readPNG("PicturePluteus.png")

#Periods (see paper for their definition)
begin=1958:1968
DEP1=1969:1977
RS1=1978:1988
DEP2=1989:1995
RS2=1996:2003
DEP3=2004:2013
end=2014:2018
periods=c(rep("margin",length(begin)),rep("DEP",length(DEP1)),rep("RS",length(RS1))
          ,rep("DEP",length(DEP2)),rep("RS",length(RS2)),rep("DEP",length(DEP3)),
          rep("margin",length(end)))
names(periods)=Yearmin:Yearmax

cols=c("grey30","dodgerblue3","darkorange2")
names(cols)=c("margin","DEP","RS")

#Figure 4:------------------

pdf(paste("Fig. 4 Seasonal results.pdf"),width=6,height=8)
layout(matrix(c(11,11,11,2,1,5,4,3,5,11,11,11,7,6,10,9,8,10),6,3,byrow=T),widths=c(1,1,0.25),heights=c(0.3,1,1,0.3,1,1))
par(cex=0.6)
for (j in c("Phyto","Zoo")){
  
  #this first loop gets the limits of variation in the different standard areas to provide 
  #a common color scale in the legend
  #the limits of that variation will be stored in "range"
  range=integer(0)
  
  for (i in c("C01","C02","D01","D02")){
    #the projection here is used to re-flip the axes (as in figure 2) and to bound the axes between 0 and 1
    #"Functional Composition axes" described in the paper
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,1]
    names(truc)=rownames(proj$lisup)
    
    #y will contain the PC1
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),1]
    #Conditions to flip axes (as in figure 2)
    if (j=="Phyto"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])<0){
        truc=-truc
        y=-y
      }
    }else if (j=="Zoo"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])>0){
        truc=-truc
        y=-y
      }
    }
    
    #bounding y between 0 and 1 by scaling it to the trait-space
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    range=c(range,y)
  }
  
  #this will be the limit of the color scale in the legend
  absmax=ceiling(max(range,na.rm=T)*100)/100
  absmin=floor(min(range,na.rm=T)*100)/100
  
  
  #Second loop, actually produces the graphs
  for (i in c("C01","C02","D01","D02")){
    #NOTE THAT WE START AS ABOVE, WE STILL NEED TO FLIP AND SCALE PC1 TO MAKE IT THE FCOMP1 DISPLAYED IN THE PAPER
    #the projection here is used to re-flip the axes (as in figure 2) and to bound the axes between 0 and 1
    #"Functional Composition axes" described in the paper
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,1]
    names(truc)=rownames(proj$lisup)
    
    #y will contain the PC1
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),1]
    #Conditions to flip axes (as in figure 2)
    if (j=="Phyto"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])<0){
        truc=-truc
        y=-y
      }
    }else if (j=="Zoo"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])>0){
        truc=-truc
        y=-y
      }
    }
    
    #bounding y between 0 and 1 by scaling it to the trait-space
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    
    
    #Start plotting:
    par(mar=c(3,2,2,1))
    plot(1,1,ylim=c(12.5,0.5),xlim=c(Yearmin-0.5,Yearmax+0.5),xlab="",ylab="",las=1,xaxs="i",yaxs="i",main="",yaxt="n",xaxt="n")
    mtext(i,3,0.5,at=Yearmin,font=4,cex=0.9)
    mtext(side=3,line=0.5,cex=0.7,at=1970,font=3,col="grey30",paste("(",sprintf("%.2f",round(100*FTBA[[j]][[i]]$PCA$eig[1]/sum(FTBA[[j]][[i]]$PCA$eig),2))," %)",sep=""))
    axis(side=2,at=1:12,labels=c("J","F","M","A","M","J","J","A","S","O","N","D"),las=1)
    
    rect(xleft=floor(x)-0.5,xright=floor(x)+0.5,ybottom=rep(1:12,length(Yearmin:Yearmax))-0.5,ytop=rep(1:12,length(Yearmin:Yearmax))+0.5,
         col=colourvalues::colour_values(c(y,absmin,absmax),"viridis",na_colour="#FFFFFF00")[1:length(y)],border=NA)
     
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=13.5,ytop=12.5,col=cols[periods],border=NA,xpd=NA)
    text(c("DEP1","RS1","DEP2","RS2","DEP3"),x=c(mean(DEP1),mean(RS1),mean(DEP2),mean(RS2),mean(DEP3)),y=13,col="white",xpd=NA)
    axis(side=1,pos=13.5,xpd=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=0.5,ytop=12.5,xpd=NA,col=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=12.5,ytop=13.5,xpd=NA,col=NA)
    
  }
  
  #additional code for the legend and some annotation
  par(mar=c(5,0,5,4))
  barplot(horiz=T,rep(1,100),col=c(rep("grey",length(1:(absmin*100))-1),colourvalues::colour_values((absmin*100):((absmax*100)-1)),rep("grey",length((absmax*100):100)-1)),space=0,border=NA,axes=F)
  axis(side=4,at=c(0,absmin*100,50,absmax*100,100),labels=c(0,absmin,0.5,absmax,1),las=1)
  rect(xleft=0,xright=1,ybottom=0,ytop=100,col=NA)
  mtext(side=4,line=2.5,"FComp1",cex=0.8,font=2)
  if(j=="Phyto"){
    text(x=c(-0.5,-0.5),y=c(-10,110),c("Non-motile\nautotrophs","Motile\nmixotrophs"),col="grey30",font=4,pos=4,xpd=NA,cex=1)
    rasterImage(diatopic,xleft=-0.2,xright=2.7,ybottom=-32,ytop=-15,xpd=NA)
    rasterImage(dinopic,xleft=0,xright=2,ybottom=120,ytop=135,xpd=NA)
    
  }else{
    text(x=c(-0.5,-0.5),y=c(-10,105),c("Crustacean\nholoplankton","Meroplankton"),col="grey30",font=4,pos=4,xpd=NA,cex=1)
    rasterImage(coppic,xleft=0.2,xright=2.2,ybottom=-30,ytop=-15,xpd=NA)
    rasterImage(meropic,xleft=0.2,xright=1.7,ybottom=110,ytop=122,xpd=NA)
    
  }
}

par(mar=c(0,0,0,0))
plot(axes=F,NA,NA,ylim=c(0,2.6),xlim=c(0,10),xaxs="i",yaxs="i")
text(x=c(0,0),y=c(0.15,2.45),c("b) Zooplankton","a) Phytoplankton"),xpd=NA,pos=4,font=2,cex=1.3)

dev.off()
#Supplementary Figure FComp2 Seasonal------------------

pdf(paste("Suppl Fig FComp2 Seasonal.pdf"),width=6,height=8)
layout(matrix(c(11,11,11,2,1,5,4,3,5,11,11,11,7,6,10,9,8,10),6,3,byrow=T),widths=c(1,1,0.25),heights=c(0.3,1,1,0.3,1,1))
par(cex=0.6)

#for comments see figure 3 above.
#Note that we do not flip the axes because, unlike axes 1,
#we haven't found a lot similarity in between them.
for (j in c("Phyto","Zoo")){
  range=integer(0)
  for (i in c("C01","C02","D01","D02")){
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,2]
    names(truc)=rownames(proj$lisup)
    
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),2]
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    range=c(range,y)
  }
  
  absmax=ceiling(max(range,na.rm=T)*100)/100
  absmin=floor(min(range,na.rm=T)*100)/100
  
  for (i in c("C01","C02","D01","D02")){
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,2]
    names(truc)=rownames(proj$lisup)
    
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),2]
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    
    par(mar=c(3,2,2,1))
    plot(1,1,ylim=c(12.5,0.5),xlim=c(Yearmin-0.5,Yearmax+0.5),xlab="",ylab="",las=1,xaxs="i",yaxs="i",main="",yaxt="n",xaxt="n")
    mtext(i,3,0.5,at=Yearmin,font=4,cex=0.9)
    mtext(side=3,line=0.5,cex=0.7,at=1970,font=3,col="grey30",paste("(",sprintf("%.2f",round(100*FTBA[[j]][[i]]$PCA$eig[2]/sum(FTBA[[j]][[i]]$PCA$eig),2))," %)",sep=""))
    axis(side=2,at=1:12,labels=c("J","F","M","A","M","J","J","A","S","O","N","D"),las=1)
    
    rect(xleft=floor(x)-0.5,xright=floor(x)+0.5,ybottom=rep(1:12,length(Yearmin:Yearmax))-0.5,ytop=rep(1:12,length(Yearmin:Yearmax))+0.5,
         col=colourvalues::colour_values(c(y,absmin,absmax),"viridis",na_colour="#FFFFFF00")[1:length(y)],border=NA)
    
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=13.5,ytop=12.5,col=cols[periods],border=NA,xpd=NA)
    text(c("DEP1","RS1","DEP2","RS2","DEP3"),x=c(mean(DEP1),mean(RS1),mean(DEP2),mean(RS2),mean(DEP3)),y=13,col="white",xpd=NA)
    axis(side=1,pos=13.5,xpd=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=0.5,ytop=12.5,xpd=NA,col=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=12.5,ytop=13.5,xpd=NA,col=NA)
  }
  
  par(mar=c(5,0,5,4))
  barplot(horiz=T,rep(1,100),col=c(rep("grey",length(0:(absmin*100))-1),colourvalues::colour_values((absmin*100):((absmax*100)-1)),rep("grey",length((absmax*100):100)-1)),space=0,border=NA,axes=F)
  axis(side=4,at=c(0,absmin*100,50,absmax*100,100),labels=c(0,absmin,0.5,absmax,1),las=1)
  rect(xleft=0,xright=1,ybottom=0,ytop=100,col=NA)
  mtext(side=4,line=2.5,"FComp2",cex=0.8,font=2)
}

par(mar=c(0,0,0,0))
plot(axes=F,NA,NA,ylim=c(0,2.6),xlim=c(0,10),xaxs="i",yaxs="i")
text(x=c(0,0),y=c(0.15,2.45),c("b) Zooplankton","a) Phytoplankton"),xpd=NA,pos=4,font=2,cex=1.3)

dev.off()


##############################################################
#FIGURE 5:Annual + sampling effort----------------------
##############################################################
#Note that the kruskal-Wallis test are realised here-----
setwd("")
load("FTBA NSRS.RData")

#calling packages, and defining some graphical characteristics----------
library(FD)
library(png)

Yearmin=1958#When do we start to plot?
Yearmax=2018#When do we stop ploting?
x=Yearmin:Yearmax
x=rep(x,each=12)+(((1:12)/12)-(1/24))

#pictures to illustrate axes
coppic=readPNG("PictureCop.png")
dinopic=readPNG("PictureDino.png")
diatopic=readPNG("PictureChaeto.png")
meropic=readPNG("PicturePluteus.png")

#Periods
begin=1958:1968
DEP1=1969:1977
RS1=1978:1988
DEP2=1989:1995
RS2=1996:2003
DEP3=2004:2013
end=2014:2018
periods=c(rep("margin",length(begin)),rep("DEP",length(DEP1)),rep("RS",length(RS1))
          ,rep("DEP",length(DEP2)),rep("RS",length(RS2)),rep("DEP",length(DEP3)),
          rep("margin",length(end)))
names(periods)=Yearmin:Yearmax

cols=c("grey30","dodgerblue3","darkorange2")
names(cols)=c("margin","DEP","RS")

#this will contain the result of the kruskal-wallis tests
testresults=list()
#Figure 5----------------------

#start------
pdf(paste("Fig. 5 results annual.pdf"),width=6,height=8)
layout(matrix(c(9,9,2,1,4,3,9,9,6,5,8,7),6,2,byrow=T),widths=c(1,1),heights=c(0.3,1,1,0.3,1,1))
par(cex=0.6)

for (j in c("Phyto","Zoo")){
  for (i in c("C01","C02","D01","D02")){
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,1]
    names(truc)=rownames(proj$lisup)
    
    #y will contain the PC1
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),1]
    #Conditions to flip axes (as in figure 2)
    if (j=="Phyto"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])<0){
        truc=-truc
        y=-y
      }
    }else if (j=="Zoo"){
      if(cor(truc,FTBA[[j]]$TraitSpace[,1])>0){
        truc=-truc
        y=-y
      }
    }
    
    #bounding y between 0 and 1 by scaling it to the trait-space
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    
    y=tapply(y,floor(x),mean,na.rm=T)
    
    samp=table(floor(as.numeric(rownames(FTBA[[j]][[i]]$centroids))))#this gives the number of months sampled
    
    if(i %in% c("C02","D02")){
      par(mar=c(3,8,2,0))
    }else{
      par(mar=c(3,4,2,4))
    }
    
    plot(NA,NA,xlim=c(Yearmin-0.5,Yearmax+0.5),ylim=c(0,12),ylab="",xlab="",axes=F,yaxs="i",xaxs="i")
    if(i %in% c("C01","D01")){
      axis(side=4,las=1,col.axis="grey40",col="grey40")
    }else{
      axis(side=4,las=1,col.axis="grey40",col="grey40",labels=F)
    }
    mtext(i,3,0.5,at=Yearmin,font=4,cex=0.9)
    mtext(side=3,line=0.5,cex=0.7,at=1970,font=3,col="grey30",paste("(",sprintf("%.2f",round(100*FTBA[[j]][[i]]$PCA$eig[1]/sum(FTBA[[j]][[i]]$PCA$eig),2))," %)",sep=""))
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=0,ytop=samp,col=cols[periods],border=NA)
    rect(xleft=c(min(DEP1)-0.5,min(DEP2)-0.5,min(DEP3)-0.5),xright=c(max(DEP1)+0.5,max(DEP2)+0.5,max(DEP3)+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    rect(xleft=c(min(RS1)-0.5,min(RS2)-0.5),xright=c(max(RS1)+0.5,max(RS2)+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    rect(xleft=c(Yearmin-0.5,max(DEP3)+0.5),xright=c(min(DEP1)-0.5,Yearmax+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=-1,ytop=0,col=cols[periods],border=NA,xpd=NA)
    text(c("DEP1","RS1","DEP2","RS2","DEP3"),x=c(mean(DEP1),mean(RS1),mean(DEP2),mean(RS2),mean(DEP3)),y=-0.5,col="white",xpd=NA)
    axis(side=1,pos=-1,xpd=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=-1,ytop=0,xpd=NA,col=NA)
    
    par(new=T)
    plot(Yearmin:Yearmax,y,xlab="",ylab="",las=1,xaxs="i",xlim=c(Yearmin-0.5,Yearmax+0.5),col=rgb(1,1,1,0),xaxt="n")
    
    #Comparisons in between the different DEP
    means=c(mean(y[as.character(DEP1)]),mean(y[as.character(DEP2)]),mean(y[as.character(DEP3)]))
    sd=c(sd(y[as.character(DEP1)]),sd(y[as.character(DEP2)]),sd(y[as.character(DEP3)]))
    
    #segments for means and sd of DEP
    segments(y0=means,
             x0=c(min(DEP1),min(DEP2),min(DEP3))-0.5,x1=c(max(DEP1),max(DEP2),max(DEP3))+0.5,
             col=cols["DEP"],lwd=3,lend=1)
    
    
    #plot the data
    points(Yearmin:Yearmax,y,type="b",pch=21,bg=cols[periods],lwd=1,cex=1)
    box()
    
    
    #Showing the result of the Kruskal test
    test1=kruskal.test(list(y[as.character(DEP1)],y[as.character(DEP2)]))
    test2=kruskal.test(list(y[as.character(DEP2)],y[as.character(DEP3)]))
    
    if (test1$p.value>0.05){
      mtext(side=3,line=0.25,at=mean(RS1),"NS",font=2,cex=0.6)
    }else{
      mtext(side=3,line=0.25,at=mean(RS1),sprintf("%.3f",round(abs(means[1]-means[2]),3)),font=2,cex=0.6)
    }
    if (test2$p.value>0.05){
      mtext(side=3,line=0.25,at=mean(RS2),"NS",font=2,cex=0.6)
    }else{
      mtext(side=3,line=0.25,at=mean(RS2),sprintf("%.3f",round(abs(means[3]-means[2]),3)),font=2,cex=0.6)
    }
    
    
    #keeping the test results on the side
    testresults[[j]][[i]]$effects=c(abs(means[1]-means[2]),abs(means[3]-means[2]))
    names(testresults[[j]][[i]]$effects)=c("RS1","RS2")
    testresults[[j]][[i]]$p=c(test1$p.value,test2$p.value)
    names(testresults[[j]][[i]]$p)=c("RS1","RS2")
  }
}
#additional code for anotation--------
par(mar=c(0,0,0,0))
plot(axes=F,NA,NA,ylim=c(0,2.6),xlim=c(0,10),xaxs="i",yaxs="i")
text(x=c(0,0),y=c(0.15,2.45),c("b) Zooplankton","a) Phytoplankton"),xpd=NA,pos=4,font=2,cex=1.3)
mtext(side=4,line=-1.5,at=c(-1,1.3),"Number of months sampled",col="grey40",font=4,cex=0.7)
mtext(side=2,line=-3,at=c(-1,1.3),"FComp1",font=2,cex=0.9)
text(x=0.55,y=c(0.4,2.2,-1.9,-0.1),c("Non-motile\nautotrophs","Motile\nmixotrophs","Crustacean\nholoplankton","Meroplankton"),col="grey30",font=4,pos=c(3,1,3,1),xpd=NA,cex=1)

polygon(x=c(0.55,0.3,0.8),y=c(-0.05,-0.1,-0.1),xpd=NA,col="black",border=NA)
polygon(x=c(0.55,0.3,0.8),y=c(-1.95,-1.9,-1.9),xpd=NA,col="black",border=NA)
polygon(x=c(0.55,0.3,0.8),y=c(2.25,2.2,2.2),xpd=NA,col="black",border=NA)
polygon(x=c(0.55,0.3,0.8),y=c(0.35,0.4,0.4),xpd=NA,col="black",border=NA)
rasterImage(coppic,xleft=0.15,xright=0.95,ybottom=-1.75,ytop=-1.55,xpd=NA)
rasterImage(meropic,xleft=0.2,xright=0.8,ybottom=-0.4,ytop=-0.2,xpd=NA)
rasterImage(diatopic,xleft=0.1,xright=1.1,ybottom=0.55,ytop=0.75,xpd=NA)
rasterImage(dinopic,xleft=0.1,xright=0.9,ybottom=1.85,ytop=2.05,xpd=NA)

#end---------
dev.off()

#Supplementary figure: Annual on FComp 2---------
#start------
pdf(paste("Suppl Fig FComp2 Annual.pdf"),width=6,height=8)
layout(matrix(c(9,9,2,1,4,3,9,9,6,5,8,7),6,2,byrow=T),widths=c(1,1),heights=c(0.3,1,1,0.3,1,1))
par(cex=0.6)

for (j in c("Phyto","Zoo")){
  for (i in c("C01","C02","D01","D02")){
    proj=suprow(FTBA[[j]][[i]]$PCA,FTBA[[j]]$TraitSpace)
    truc=proj$lisup[,1]
    names(truc)=rownames(proj$lisup)
    
    #y will contain the PC1
    y=FTBA[[j]][[i]]$PCA$li[as.character(x),2]
    
    
    #bounding y between 0 and 1 by scaling it to the trait-space
    y=(y-min(truc,na.rm=T))/(max(truc,na.rm=T)-min(truc,na.rm=T))
    
    y=tapply(y,floor(x),mean,na.rm=T)
    
    samp=table(floor(as.numeric(rownames(FTBA[[j]][[i]]$centroids))))#this gives the number of months sampled
    
    if(i %in% c("C02","D02")){
      par(mar=c(3,8,2,0))
    }else{
      par(mar=c(3,4,2,4))
    }
    
    plot(NA,NA,xlim=c(Yearmin-0.5,Yearmax+0.5),ylim=c(0,12),ylab="",xlab="",axes=F,yaxs="i",xaxs="i")
    if(i %in% c("C01","D01")){
      axis(side=4,las=1,col.axis="grey40",col="grey40")
    }else{
      axis(side=4,las=1,col.axis="grey40",col="grey40",labels=F)
    }
    mtext(i,3,0.5,at=Yearmin,font=4,cex=0.9)
    mtext(side=3,line=0.5,cex=0.7,at=1970,font=3,col="grey30",paste("(",sprintf("%.2f",round(100*FTBA[[j]][[i]]$PCA$eig[2]/sum(FTBA[[j]][[i]]$PCA$eig),2))," %)",sep=""))
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=0,ytop=samp,col=cols[periods],border=NA)
    rect(xleft=c(min(DEP1)-0.5,min(DEP2)-0.5,min(DEP3)-0.5),xright=c(max(DEP1)+0.5,max(DEP2)+0.5,max(DEP3)+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    rect(xleft=c(min(RS1)-0.5,min(RS2)-0.5),xright=c(max(RS1)+0.5,max(RS2)+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    rect(xleft=c(Yearmin-0.5,max(DEP3)+0.5),xright=c(min(DEP1)-0.5,Yearmax+0.5),ybottom=0,ytop=12,col=rgb(1,1,1,0.7),border=NA)
    
    rect(xleft=Yearmin:Yearmax-0.5,xright=Yearmin:Yearmax+0.5,ybottom=-1,ytop=0,col=cols[periods],border=NA,xpd=NA)
    text(c("DEP1","RS1","DEP2","RS2","DEP3"),x=c(mean(DEP1),mean(RS1),mean(DEP2),mean(RS2),mean(DEP3)),y=-0.5,col="white",xpd=NA)
    axis(side=1,pos=-1,xpd=NA)
    rect(xleft=Yearmin-0.5,xright=Yearmax+0.5,ybottom=-1,ytop=0,xpd=NA,col=NA)
    
    par(new=T)
    plot(Yearmin:Yearmax,y,xlab="",ylab="",las=1,xaxs="i",xlim=c(Yearmin-0.5,Yearmax+0.5),col=rgb(1,1,1,0),xaxt="n")
    
    #Comparisons in between the different DEP
    means=c(mean(y[as.character(DEP1)]),mean(y[as.character(DEP2)]),mean(y[as.character(DEP3)]))
    sd=c(sd(y[as.character(DEP1)]),sd(y[as.character(DEP2)]),sd(y[as.character(DEP3)]))
    
    #segments for means and sd of DEP
    segments(y0=means,
             x0=c(min(DEP1),min(DEP2),min(DEP3))-0.5,x1=c(max(DEP1),max(DEP2),max(DEP3))+0.5,
             col=cols["DEP"],lwd=3,lend=1)
    
    
    #plot the data
    points(Yearmin:Yearmax,y,type="b",pch=21,bg=cols[periods],lwd=1,cex=1)
    box()
    
    
    #Showing the result of the Kruskal test
    test1=kruskal.test(list(y[as.character(DEP1)],y[as.character(DEP2)]))
    test2=kruskal.test(list(y[as.character(DEP2)],y[as.character(DEP3)]))
    
    if (test1$p.value>0.05){
      mtext(side=3,line=0.25,at=mean(RS1),"NS",font=2,cex=0.6)
    }else{
      mtext(side=3,line=0.25,at=mean(RS1),sprintf("%.3f",round(abs(means[1]-means[2]),3)),font=2,cex=0.6)
    }
    if (test2$p.value>0.05){
      mtext(side=3,line=0.25,at=mean(RS2),"NS",font=2,cex=0.6)
    }else{
      mtext(side=3,line=0.25,at=mean(RS2),sprintf("%.3f",round(abs(means[3]-means[2]),3)),font=2,cex=0.6)
    }
    
    
    #keeping the test results on the side
    testresults[[j]][[i]]$effects=c(abs(means[1]-means[2]),abs(means[3]-means[2]))
    names(testresults[[j]][[i]]$effects)=c("RS1","RS2")
    testresults[[j]][[i]]$p=c(test1$p.value,test2$p.value)
    names(testresults[[j]][[i]]$p)=c("RS1","RS2")
  }
}
#additional code for anotation--------
par(mar=c(0,0,0,0))
plot(axes=F,NA,NA,ylim=c(0,2.6),xlim=c(0,10),xaxs="i",yaxs="i")
text(x=c(0,0),y=c(0.15,2.45),c("b) Zooplankton","a) Phytoplankton"),xpd=NA,pos=4,font=2,cex=1.3)
mtext(side=4,line=-1.5,at=c(-1,1.3),"Number of months sampled",col="grey40",font=4,cex=0.7)
mtext(side=2,line=-3,at=c(-1,1.3),"FComp2",font=2,cex=0.9)


#end---------
dev.off()


##############################################################
#Supplementary figure on Trait Space (results of PCoAs and projection of PC1 on centroids)-------
##############################################################
setwd("")
load("FTBA NSRS.RData")
library(FD)


for (i in c("Phyto","Zoo")){
  pdf(paste("Suppl Fig Trait Space characteristic and PC1 projection",i,".pdf"),width=8,height=7)
  layout(matrix(c(11,1,1,2,2,11,3,3,4,5,6,6,7,7,8,9,10,10),3,6,byrow=T),widths=c(1,1,2,2,1,1),heights=c(3,2,2))
  
  #first plot: description of the trait space with correlations with traits
  par(mar=c(4,4,2,2))
  plot(FTBA[[i]]$TraitSpace[,1],FTBA[[i]]$TraitSpace[,2],las=1,pch=16,col="dodgerblue3",xlab=paste("PCoA axis 1 (",round(FTBA[[i]]$pcoa$values$Rel_corr_eig[1]*100,2)," %)",sep=""),ylab=paste("PCoA axis 2 (",round(FTBA[[i]]$pcoa$values$Rel_corr_eig[2]*100,2)," %)",sep=""))
  corel=envfit(FTBA[[i]]$TraitSpace,FTBA[[i]]$Traits,choices=c(1,2,3,4))
  if (i =="Zoo"){
    plot(corel,col="grey30",font=2,xpd=NA,labels=list(factors=c("Holoplanktonic","Meroplanktonic","exoskeleton","internal skeleton","naked","shell")))
    mtext("Zooplankton",side=3,line=0.5,font=4,cex=1,at=min(FTBA[[i]]$TraitSpace[,1]))
  }else{
    plot(corel,col="grey30",font=2,xpd=NA,labels=list(factors=c("Organic","Siliceous","Colonial","Solitary","motile","non-motile")))
    mtext("Phytoplankton",side=3,line=0.5,font=4,cex=1,at=min(FTBA[[i]]$TraitSpace[,1]))
  }
  #simple barplots for eigenvalues of PCoAs
  par(mar=c(4,4,2,2))
  barplot(FTBA[[i]]$pcoa$values$Rel_corr_eig*100,las=1, ylab="PCoA axes relative eigenvalues (%)",xpd=NA)
  
  #set of plots (two by standard area) to display: 
  for(j in c("C02","C01","D02","D01")){
    #1) the centroids as well as the projection of the first axis of the PCA performed on the centroids 
    slope=FTBA[[i]][[j]]$PCA$c1[2,1]/FTBA[[i]][[j]]$PCA$c1[1,1]
    intercept=mean(FTBA[[i]][[j]]$centroids[,2])-mean(FTBA[[i]][[j]]$centroids[,1])*slope
    par(mar=c(1,1,2,2))
    plot(FTBA[[i]]$TraitSpace[,1],FTBA[[i]]$TraitSpace[,2],las=1,pch=16,col="dodgerblue3",xlab="",ylab="",axes=F)
    box()
    points(FTBA[[i]][[j]]$centroids[,1],FTBA[[i]][[j]]$centroids[,2],col=rgb(t(col2rgb("orangered2")/256),alpha=0.3),pch=3)
    abline(a=intercept,b=slope,col="grey30",lwd=3)
    mtext(j,side=3,line=0.5,at=min(FTBA[[i]]$TraitSpace[,1]),font=4,cex=1)
    
    #2) Key information on the first PCA axis
    par(mar=c(1,2.5,2,2))
    if(FTBA[[i]][[j]]$PCA$c1$CS1[1]<0){#this condition flips all the graphs in the same direction to ease comparison
      barplot(FTBA[[i]][[j]]$PCA$c1$CS1*-1,las=1,ylab="Principal axes of PC1",xpd=NA)
      legend(x="topright",bty="n",c(
        paste("Var Exp:",round(FTBA[[i]][[j]]$PCA$eig[1]/sum(FTBA[[i]][[j]]$PCA$eig)*100,2),"%"),
        paste("Range:",round(min(FTBA[[i]][[j]]$PCA$li$Axis1*-1),2),"to",round(max(FTBA[[i]][[j]]$PCA$li$Axis1*-1),2))))
    }else{
      barplot(FTBA[[i]][[j]]$PCA$c1$CS1,las=1,ylab="Principal axes of PC1",xpd=NA)
      legend(x="topright",bty="n",c(
        paste("Var Exp:",round(FTBA[[i]][[j]]$PCA$eig[1]/sum(FTBA[[i]][[j]]$PCA$eig)*100,2),"%"),
        paste("Range:",round(min(FTBA[[i]][[j]]$PCA$li$Axis1),2),"to",round(max(FTBA[[i]][[j]]$PCA$li$Axis1),2))))
    }
  }
  dev.off()
}
