## Set working directory
  setwd("x")
  
  ## Loading libraries
  library(readxl)  
  library(ggplot2)
  library(gridExtra)
  library(Rmisc)
  library(SIBER)
  library(dplyr)
  library(lme4)
  library(MuMIn)

  ## Data input
  C=read.csv('C.csv')  
  N=read.csv('N.csv')
  Groups=read.csv('Groups.csv')  
  
  ## Enter diet per line to C and N files
  for(i in 1:386){
     if(C$Bird[i] != 'Food'){
        C$Diet[i]=subset(Groups,Bird==C$Bird[i])$Diet
       }
     }  
  for(i in 1:376){
    if(N$Bird[i] != 'Food'){
        N$Diet[i]=subset(Groups,Bird==N$Bird[i])$Diet
    }
  }  

  ## Remove samples 70Back_2 and 90Back_2, two old feathers with Waddensea isotope signature
  C=subset(C,Sample!='70Back_2'&Sample!='70Back_2_2'&Sample!='90Back_2'&Sample!='90Back_2_2')
  N=subset(N,Sample!='70Back_2'&Sample!='70Back_2_2'&Sample!='90Back_2'&Sample!='90Back_2_2')

  ## Select only food values for C and N
  FoodC=subset(C,Tissue=='Food')
  FoodC$Time=as.numeric(FoodC$Time)
  ## Throw out first Peringia values (old batch) and one outlier which probably contained shell material  
  FoodC=subset(FoodC,Diet=='Trouvit'|(Diet=='Peringia'&Time>6&d13C< -20))
  FoodC=aggregate(d13C~Time+Diet,FoodC,mean)
  FoodN=aggregate(d15N~Time+Diet,subset(N,Tissue=='Food'),mean)
  FoodN$Time=as.numeric(FoodN$Time)
  FoodN=subset(FoodN,Diet=='Trouvit'|(Diet=='Peringia'&Time>6))
  
  ## Explore food values
  boxplot(d13C~Diet,FoodC)
  boxplot(d15N~Diet,FoodN)
  
  ## Check if food values are stable over time
  par(mfrow=c(2,1))
  plot(d13C~Time,FoodC)
  plot(d15N~Time,FoodN)

  summary(lm(d13C~Time,subset(FoodC,Diet=='Peringia')))
  summary(lm(d13C~Time,subset(FoodC,Diet=='Trouvit')))
  summary(lm(d15N~Time,subset(FoodN,Diet=='Peringia')))
  summary(lm(d15N~Time,subset(FoodN,Diet=='Trouvit')))
  
  ## Calculate averages (±SD/SE) of d13C and d15N per food type
  Cfood=summarySE(FoodC,measurevar='d13C',groupvars='Diet')
  Nfood=summarySE(FoodN,measurevar='d15N',groupvars='Diet')
  
  ## Select bird tissue samples
  sampleC=subset(C,Tissue!='Food') 
  sampleN=subset(N,Tissue!='Food') 
    
  ### Exclude February and March tissue samples on Peringia diet
  sampleC=sampleC[-which(sampleC$Diet=='Peringia'&sampleC$Time=='February'),]
  sampleC=sampleC[-which(sampleC$Diet=='Peringia'&sampleC$Time=='March'),]
  sampleN=sampleN[-which(sampleN$Diet=='Peringia'&sampleN$Time=='February'),]
  sampleN=sampleN[-which(sampleN$Diet=='Peringia'&sampleN$Time=='March'),]

  ## Prepare for plotting
  sampleC$Type=factor(sampleC$Type,levels=c('Plasma','Cells','Primary covert','Breast','Back'))
  sampleN$Type=factor(sampleN$Type,levels=c('Plasma','Cells','Primary covert','Breast','Back'))

  plot1=merge(sampleC[,c(1:7)],sampleN[,c(1:7)])
  plot2=merge(FoodC,FoodN)
  names(plot2)[1]='Type'
  plot2$Type='Food'
  plotsiber=rbind(plot1[,c(5:8)],plot2)

  ## Figure 1
  ggplot(data = plotsiber, 
                 aes(x = d13C, 
                     y = d15N)) + 
    geom_point(aes(fill = Type, shape = Diet), size = 3) +
    ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
    xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
    theme(text = element_text(size=16)) + 
    scale_shape_manual(values=c(21,24))+
    theme_classic()+ 
    theme(legend.position = "none")+
    theme(text = element_text(size=22))+
    stat_ellipse(aes(group = interaction(Type, Diet), 
                     fill=Type,color=Type ),
                 size=1.3,
                 alpha = 0.2, 
                 level = 0.95,
                 type = "norm",
                 geom = "polygon")+
    scale_color_manual(values=c('yellow','red','grey55','grey80','grey30','blue'))+
    scale_fill_manual(values=c('yellow','red','grey55','grey80','grey30','blue'))

  ## Calculate discrimination per sample (difference between stable isotope value of sample and its corresponding diet)
  sampleC$Discrimination=0
  for(i in 1:291){
    if(sampleC$Diet[i]=='Peringia'){
      sampleC$Discrimination[i]=sampleC$d13C[i]-Cfood$d13C[1]
    }
    else{sampleC$Discrimination[i]=sampleC$d13C[i]-Cfood$d13C[2]
    }}
  
  sampleN$Discrimination=0
  for(i in 1:291){
    if(sampleN$Diet[i]=='Peringia'){
      sampleN$Discrimination[i]=sampleN$d15N[i]-Nfood$d15N[1]
    }
    else{sampleN$Discrimination[i]=sampleN$d15N[i]-Nfood$d15N[2]
    }}
  

  ## Figure 2: Boxplots of discrimination factors
  p1=ggplot(data=sampleC,aes(x=Type,y=Discrimination,fill=Diet))+
    geom_jitter(pch=1,cex=0.9,position=position_jitterdodge(jitter.width=0.25))+
    geom_boxplot(alpha=0.85)+
    theme_classic()+
    theme(legend.position="none")+
    ggtitle(expression(paste(''^{13},"C")))+
    theme(plot.title=element_text(hjust=0.5,size=30),
          axis.title.y = element_text(size=14),
          axis.text.y=element_text(size=14),
          axis.title.x=element_blank(),
          plot.margin = margin(5.5,5.5,20,5.5))+
    ylab(expression(paste("Discrimination factor  ",italic(Delta))))+
    annotate("text", x = 1.5, y = -0.55, label = "Blood",size=5) +
    annotate("text", x = 4, y = -0.55, label = "Feathers",size=5) +
    coord_cartesian(ylim = c(0, 4), clip = "off")
  
  p2=ggplot(data=sampleN,aes(x=Type,y=Discrimination,fill=Diet))+
    geom_jitter(pch=1,cex=0.9,position=position_jitterdodge(jitter.width=0.25))+
    geom_boxplot(alpha=0.85)+
    theme_classic()+
    theme(legend.position=c(0.8,0.15),legend.text=element_text(size=15),legend.title=element_text(size=15),
          legend.box.background = element_blank(),
          axis.title.y = element_blank(),axis.text.y = element_blank(),
          axis.title.x=element_blank())+
    ggtitle(expression(paste(''^{15},"N")))+
    theme(plot.title=element_text(hjust=0.5,size=30),
          plot.margin = margin(5.5,25.5,20,5.5))+
    annotate("text", x = 1.5, y = -0.55, label = "Blood",size=5) +
    annotate("text", x = 4, y = -0.55, label = "Feathers",size=5) +
    coord_cartesian(ylim = c(0, 4), clip = "off")
  
  grid.arrange(p1,p2,nrow=1)
 
  ## Check if adding 'individual bird' as random effect improves the model explaining variation in discrimination factors
  cmodel1=lm(Discrimination~Type*Diet,sampleC)
  cmodel2=lmer(Discrimination~Type*Diet+(1|Bird),sampleC)
  model.sel(cmodel1,cmodel2)
  nmodel1=lm(Discrimination~Type*Diet,sampleN)
  nmodel2=lmer(Discrimination~Type*Diet+(1|Bird),sampleN)
  model.sel(nmodel1,nmodel2)

  ## Calculate mean values per individual bird, time point,tissue type and diet
  avgC=summarySE(sampleC,measurevar='d13C',groupvars=c('Type','Diet'))
  avgN=summarySE(sampleN,measurevar='d15N',groupvars=c('Type','Diet'))
  ## Calculate mean discrimination factors per tissue type and diet  
  DisC=summarySE(sampleC,measurevar='Discrimination',groupvars=c('Type','Diet'))
  DisN=summarySE(sampleN,measurevar='Discrimination',groupvars=c('Type','Diet'))
  
  ## Test differences in discrimination factors between tissue types and diets
  anovaC=aov(Discrimination~Diet*Type,sampleC)
  summary(anovaC)
  TukeyHSD(anovaC)  
  
  anovaN=aov(Discrimination~Diet*Type,sampleN)
  summary(anovaN)
  TukeyHSD(anovaN)