## Set working directory
  setwd("//ZEUS/cos/users/toortwijn/Mauritania/Detox experiment/Final/Oecologia/Polished raw data")

  ## Loading libraries
  library(lme4)
  library(lmerTest)
  library(Rmisc)
  library(readxl)
  library(AICcmodavg)
  
  ## Data input
  Loripes_weight=read_xlsx('Loripes_weight.xlsx')
  Loripes_sulfur=read_xlsx('Loripes_sulfur.xlsx')
  Loripes_consumers=read_xlsx('Loripes_consumers.xlsx')
  Intake_data=read_xlsx('Intake_data.xlsx')
  Gizzard_data=read_xlsx('Gizzard_measurements.xlsx')
  
  ## Averaging sulfur percentage per treatment and location
	Sulfur_table=summarySE(Loripes_sulfur,measurevar='Percentage_Sulfur',groupvars=c('Treatment','Location'))
	
  ## Averaging intake rates on Loripes per treatment and location
	Intake_rate=summarySE(Intake_data,measurevar='Consumed_Loripes',groupvars=c('Treatment','Location'))

  ## Select only correctly weighed Loripes samples
  Loripes_drymass=subset(Loripes_weight,DM_wrong=='NA')
  Loripes_shellmass=subset(Loripes_weight,SM_wrong=='NA')
  
  ## Define Experiment_Day in drymass and shellmass dataframes
  Session_Day=aggregate(Experiment_Day~Session,Loripes_sulfur,max)
  Loripes_drymass=merge(Session_Day,Loripes_drymass)
  Loripes_shellmass=merge(Session_Day,Loripes_shellmass)
	
## Figure 1 ----
	par(mfrow=c(2,1))
	
## Figure 1a	
	par(mar=c(1,6,1.5,3))
	
	barplot(Sulfur_table$Percentage_Sulfur,space=c(0,0,0.5,0),ylab=expression(paste('Sulfur percentage',italic(' S'))),col=c('green3','brown','green3','brown'),cex.axis=1.2,cex.lab=1.5,ylim=c(1.5,3.5),xpd=F)
	legend(x=3,y=3.6,legend=c('A','B'),pch=15,cex=1.4,col=c('green3','brown'),bty='n',title="Location     ")  
	
	arrows(x0=0.5,x1=0.5,y0=Sulfur_table$Percentage_Sulfur[1]-Sulfur_table$se[1],y1=Sulfur_table$Percentage_Sulfur[1]+Sulfur_table$se[1],code=3,angle=90,length=0.1)
	arrows(x0=1.5,x1=1.5,y0=Sulfur_table$Percentage_Sulfur[2]-Sulfur_table$se[2],y1=Sulfur_table$Percentage_Sulfur[2]+Sulfur_table$se[2],code=3,angle=90,length=0.1)
	arrows(x0=3,x1=3,y0=Sulfur_table$Percentage_Sulfur[3]-Sulfur_table$se[3],y1=Sulfur_table$Percentage_Sulfur[3]+Sulfur_table$se[3],code=3,angle=90,length=0.1)
	arrows(x0=4,x1=4,y0=Sulfur_table$Percentage_Sulfur[4]-Sulfur_table$se[4],y1=Sulfur_table$Percentage_Sulfur[4]+Sulfur_table$se[4],code=3,angle=90,length=0.1)
	mtext("A)",3,at=-0.8,cex=2)
	text(0.5,1.6,'47',cex=0.9)
  text(1.5,1.6,'49',cex=0.9)
	text(3,1.6,'50',cex=0.9)
  text(4,1.6,'49',cex=0.9)
  
## Figure 1b
	par(mar=c(4,6,1.2,3))
	barplot(Intake_rate$Consumed_Loripes,space=c(0,0,0.5,0),ylab=expression(paste('Intake rate ',italic('I'),' (',italic('Loripes '),'per trial)')),cex.axis=1.2,cex.lab=1.5,col=c('green3','brown','green3','brown'),ylim=c(30,35),xpd=F)

	arrows(x0=0.5,x1=0.5,y0=Intake_rate$Consumed_Loripes[1]-Intake_rate$se[1],y1=Intake_rate$Consumed_Loripes[1]+Intake_rate$se[1],code=3,angle=90,length=0.1)
	arrows(x0=1.5,x1=1.5,y0=Intake_rate$Consumed_Loripes[2]-Intake_rate$se[2],y1=Intake_rate$Consumed_Loripes[2]+Intake_rate$se[2],code=3,angle=90,length=0.1)
	arrows(x0=3,x1=3,y0=Intake_rate$Consumed_Loripes[3]-Intake_rate$se[3],y1=Intake_rate$Consumed_Loripes[3]+Intake_rate$se[3],code=3,angle=90,length=0.1)
	arrows(x0=4,x1=4,y0=Intake_rate$Consumed_Loripes[4]-Intake_rate$se[4],y1=Intake_rate$Consumed_Loripes[4]+Intake_rate$se[4],code=3,angle=90,length=0.1)
	mtext("Control",1,at=1,cex=1.5,line=2)
	mtext("Treatment",1,at=3.5,cex=1.5,line=2)
	mtext("B)",3,at=-0.8,cex=2,line=1)
	text(0.5,30.3,'120',cex=0.9)
	text(1.5,30.3,'120',cex=0.9)
	text(3,30.3,'120',cex=0.9)
	text(4,30.3,'120',cex=0.9)
	

## Statistics Figure 1a, variation in sulfur percentage ------

  sulfmodelnames=c(1:15)
  sulfmodels=list()
  sulfmodels[[1]] <- lmer(Percentage_Sulfur~as.factor(Treatment)+Experiment_Day+Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[1]]='~T+D+L+(1|Ses)'
  sulfmodels[[2]] <- lmer(Percentage_Sulfur~as.factor(Treatment)+Experiment_Day+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[2]]='~T+D+(1|Ses)'
  sulfmodels[[3]] <- lmer(Percentage_Sulfur~as.factor(Treatment)+Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[3]]='~T+L+(1|Ses)'
  sulfmodels[[4]] <- lmer(Percentage_Sulfur~Experiment_Day+Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[4]]='~D+L+(1|Ses)'
  sulfmodels[[5]] <- lmer(Percentage_Sulfur~as.factor(Treatment)*Location+Experiment_Day+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[5]]='~T*L+D+(1|Ses)'
  sulfmodels[[6]] <- lmer(Percentage_Sulfur~as.factor(Treatment)*Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[6]]='~T*L+(1|Ses)'
  sulfmodels[[7]] <- lmer(Percentage_Sulfur~as.factor(Treatment)*Experiment_Day+Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[7]]='~T*D+L+(1|Ses)'
  sulfmodels[[8]] <- lmer(Percentage_Sulfur~as.factor(Treatment)*Experiment_Day+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[8]]='~T*D+(1|Ses)'
  sulfmodels[[9]] <- lmer(Percentage_Sulfur~as.factor(Treatment)+Experiment_Day*Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[9]]='~T+D*L+(1|Ses)'
  sulfmodels[[10]] <- lmer(Percentage_Sulfur~Experiment_Day*Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[10]]='~D*L+(1|Ses)'
  sulfmodels[[11]] <- lmer(Percentage_Sulfur~as.factor(Treatment)*Experiment_Day*Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[11]]='~T*D*L+(1|Ses)'
  sulfmodels[[12]] <- lmer(Percentage_Sulfur~as.factor(Treatment)+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[12]]='~T+(1|Ses)'
  sulfmodels[[13]] <- lmer(Percentage_Sulfur~Experiment_Day+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[13]]='~D+(1|Ses)'
  sulfmodels[[14]] <- lmer(Percentage_Sulfur~Location+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[14]]='~L+(1|Ses)'
  sulfmodels[[15]] <- lmer(Percentage_Sulfur~1+(1|Session),Loripes_sulfur,REML=FALSE)
  sulfmodelnames[[15]]='~1+(1|Ses)'
  
  ## Model selection
  Sulfmodel_table <- aictab(cand.set = sulfmodels, modnames = sulfmodelnames, sort = TRUE, second.ord=T)
  Sulfmodel_table
  
  ## Best approximating model
  summary(sulfmodels[[3]])
  
  ## Check if Sulfur is Length dependent by adding Length to the best approximating model and in a simple linear model only including length
  
  summary(lmer(Percentage_Sulfur~Location+as.factor(Treatment)+Length+(1|Session),Loripes_sulfur,REML=FALSE))
  summary(lm(Percentage_Sulfur~Length,Loripes_sulfur))
  
## Variation in Loripes Drymass & Shellmass ----
  
  par(mfrow=c(1,1))
  boxplot(DryMass~Location+Treatment,Loripes_drymass,at=c(1,2,3.4,4.4),names=NA,col=c('green3','brown','green3','brown'),xlab=NA,ylab=expression(paste('flesh weight per  ',italic('Loripes '),'(g DM)')))
  legend('topright',legend=c('A','B'),pch=15,cex=1.2,col=c('green3','brown'),bty='n',title="Location")  
  mtext("Control",1,at=1.5,line=2,cex=1.5)
  mtext("Treatment",1,at=3.9,line=2,cex=1.5)
  
  ## Drymass variation
  dmmodelnames <- c(1:15)
  dmmodels <- list()
  dmmodels[[1]] <- lmer(DryMass~as.factor(Treatment)+Experiment_Day+Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[1]]='~T+D+L+(1|Ses)'
  dmmodels[[2]] <- lmer(DryMass~as.factor(Treatment)+Experiment_Day+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[2]]='~T+D+(1|Ses)'
  dmmodels[[3]] <- lmer(DryMass~as.factor(Treatment)+Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[3]]='~T+L+(1|Ses)'
  dmmodels[[4]] <- lmer(DryMass~Experiment_Day+Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[4]]='~D+L+(1|Ses)'
  dmmodels[[5]] <- lmer(DryMass~as.factor(Treatment)*Location+Experiment_Day+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[5]]='~T*L+D+(1|Ses)'
  dmmodels[[6]] <- lmer(DryMass~as.factor(Treatment)*Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[6]]='~T*L+(1|Ses)'
  dmmodels[[7]] <- lmer(DryMass~as.factor(Treatment)*Experiment_Day+Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[7]]='~T*D+L+(1|Ses)'
  dmmodels[[8]] <- lmer(DryMass~as.factor(Treatment)*Experiment_Day+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[8]]='~T*D+(1|Ses)'
  dmmodels[[9]] <- lmer(DryMass~as.factor(Treatment)+Experiment_Day*Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[9]]='~T+D*L+(1|Ses)'
  dmmodels[[10]] <- lmer(DryMass~Experiment_Day*Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[10]]='~D*L+(1|Ses)'
  dmmodels[[11]] <- lmer(DryMass~as.factor(Treatment)*Experiment_Day*Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[11]]='~T*D*L+(1|Ses)'
  dmmodels[[12]] <- lmer(DryMass~as.factor(Treatment)+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[12]]='~T+(1|Ses)'
  dmmodels[[13]] <- lmer(DryMass~Experiment_Day+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[13]]='~D+(1|Ses)'
  dmmodels[[14]] <- lmer(DryMass~Location+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[14]]='~L+(1|Ses)'
  dmmodels[[15]] <- lmer(DryMass~1+(1|Session),Loripes_drymass,REML=FALSE)
  dmmodelnames[[15]]='~1+(1|Ses)'
  
  dmmodel_table <- aictab(cand.set = dmmodels, modnames = dmmodelnames, sort = TRUE, second.ord=T)
  dmmodel_table
  
  ## Best approximating model
  summary(dmmodels[[14]])
  
  ## Shellmass variation
  smmodelnames <- c(1:15)
  smmodels <- list()
  smmodels[[1]] <- lmer(ShellMass~Treatment+Experiment_Day+Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[1]]='~T+D+L+(1|Ses)'
  smmodels[[2]] <- lmer(ShellMass~Treatment+Experiment_Day+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[2]]='~T+D+(1|Ses)'
  smmodels[[3]] <- lmer(ShellMass~Treatment+Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[3]]='~T+L+(1|Ses)'
  smmodels[[4]] <- lmer(ShellMass~Experiment_Day+Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[4]]='~D+L+(1|Ses)'
  smmodels[[5]] <- lmer(ShellMass~Treatment*Location+Experiment_Day+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[5]]='~T*L+D+(1|Ses)'
  smmodels[[6]] <- lmer(ShellMass~Treatment*Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[6]]='~T*L+(1|Ses)'
  smmodels[[7]] <- lmer(ShellMass~Treatment*Experiment_Day+Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[7]]='~T*D+L+(1|Ses)'
  smmodels[[8]] <- lmer(ShellMass~Treatment*Experiment_Day+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[8]]='~T*D+(1|Ses)'
  smmodels[[9]] <- lmer(ShellMass~Treatment+Experiment_Day*Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[9]]='~T+D*L+(1|Ses)'
  smmodels[[10]] <- lmer(ShellMass~Experiment_Day*Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[10]]='~D*L+(1|Ses)'
  smmodels[[11]] <- lmer(ShellMass~Treatment*Experiment_Day*Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[11]]='~T*D*L+(1|Ses)'
  smmodels[[12]] <- lmer(ShellMass~Treatment+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[12]]='~T+(1|Ses)'
  smmodels[[13]] <- lmer(ShellMass~Experiment_Day+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[13]]='~D+(1|Ses)'
  smmodels[[14]] <- lmer(ShellMass~Location+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[14]]='~L+(1|Ses)'
  smmodels[[15]] <- lmer(ShellMass~1+(1|Session),Loripes_shellmass,REML=FALSE)
  smmodelnames[[15]]='~1+(1|Ses)'
  
  smmodel_table <- aictab(cand.set = smmodels, modnames = smmodelnames, sort = TRUE, second.ord=T)
  smmodel_table
  
  ## Best approximating model
  summary(smmodels[[15]])

## Statistics figure 1b, variation in intake rates on Loripes ------
  irmodelnames <- c(1:14)
  irmodels <- list()
  irmodels[[1]] <- lmer(Consumed_Loripes~as.factor(Treatment)+Experiment_Day+Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[1]]='~T+D+L+(1|Ses)+(1|ID)'
  irmodels[[2]] <- lmer(Consumed_Loripes~as.factor(Treatment)+Experiment_Day*Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[2]]='~T+D*L+(1|Ses)+(1|ID)'
  irmodels[[3]] <- lmer(Consumed_Loripes~as.factor(Treatment)*Experiment_Day*Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[3]]='~T*D*L+(1|Ses)+(1|ID)'
  irmodels[[4]] <- lmer(Consumed_Loripes~as.factor(Treatment)*Location+Experiment_Day*Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[4]]='~T*L+D*L+(1|Ses)+(1|ID)'
  irmodels[[5]] <- lmer(Consumed_Loripes~as.factor(Treatment)+Experiment_Day+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[5]]='~T+D+(1|Ses)+(1|ID)'
  irmodels[[6]] <- lmer(Consumed_Loripes~Experiment_Day+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[6]]='~D+(1|Ses)+(1|ID)'
  irmodels[[7]] <- lmer(Consumed_Loripes~as.factor(Treatment)+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[7]]='~T+(1|Ses)+(1|ID)'
  irmodels[[8]] <- lmer(Consumed_Loripes~Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[8]]='~L+(1|Ses)+(1|ID)'
  irmodels[[9]] <- lmer(Consumed_Loripes~as.factor(Treatment)*Experiment_Day+Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[9]]='~T*D+L+(1|Ses)+(1|ID)'
  irmodels[[10]] <- lmer(Consumed_Loripes~as.factor(Treatment)*Location+Experiment_Day+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[10]]='~T*L+D+(1|Ses)+(1|ID)'
  irmodels[[11]] <- lmer(Consumed_Loripes~as.factor(Treatment)*Location+as.factor(Treatment)*Experiment_Day+Experiment_Day*Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[11]]='~T*L+T*D+D*L+(1|Ses)+(1|ID)'
  irmodels[[12]] <- lmer(Consumed_Loripes~as.factor(Treatment)+Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[12]]='~T+L+(1|Ses)+(1|ID)'
  irmodels[[13]] <- lmer(Consumed_Loripes~Experiment_Day+Location+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[13]]='~D+L+(1|Ses)+(1|ID)'
  irmodels[[14]] <- lmer(Consumed_Loripes~1+(1|Session)+(1|Bird_ID),Intake_data,REML=FALSE)
  irmodelnames[[14]]='~1+(1|Ses)+(1|ID)'
  
  irmodel_table <- aictab(cand.set = irmodels, modnames = irmodelnames, sort = TRUE, second.ord=T)
  irmodel_table
  
  ## Best approximating model
  summary(irmodels[[1]])
 
## Increasing intake rates over time ----   
  
  ## Averaging intake rates per experiment day, treatment and location
  timeplot=summarySE(Intake_data,measurevar='Consumed_Loripes',groupvars=c('Experiment_Day','Treatment','Location'))
  timeplot$Experiment_Day=timeplot$Experiment_Day+1  # To start at Day=1
  
  # Figure 2 (third axis is added after the gizzard calculations)
  par(mar=c(4,6,1,5))
  plot(Consumed_Loripes~Experiment_Day,timeplot,type='n',ylim=c(19,46),cex.axis=2,cex.lab=2,xlab=expression(paste('Experiment Day ',italic('D'))),ylab=expression(paste('Intake rate ',italic('I'),' (',italic('Loripes '),'per trial)')))
  arrows(x0=timeplot$Experiment_Day, x1=timeplot$Experiment_Day, y0=timeplot$Consumed_Loripes-timeplot$se, y1=timeplot$Consumed_Loripes+timeplot$se, code=3, angle=90, length=0.1)
  points(Consumed_Loripes~Experiment_Day,subset(timeplot,Treatment==0&Location=='A'),cex=3,col='green3',pch=16)
  points(Consumed_Loripes+0.15~Experiment_Day,subset(timeplot,Treatment==10&Location=='A'),cex=3,col='green3',pch=-as.hexmode("25CE"))
  points(Consumed_Loripes~Experiment_Day,subset(timeplot,Treatment==0&Location=='B'),cex=3,col='brown',pch=16)
  points(Consumed_Loripes+0.15~Experiment_Day,subset(timeplot,Treatment==10&Location=='B'),cex=3,col='brown',pch=-as.hexmode("25CE"))
  text(timeplot$Experiment_Day,45.8,labels=timeplot$Location,cex=1)
  text(1.3,46.6,expression(paste(italic('Location:'))),cex=1)
  legend('bottomright',cex=1.4,legend=c('Control, Location A','Treatment, Location A','Control, Location B','Treatment, Location B','Observed gizzard mass'),pch=c(16,-as.hexmode("25CE"),16,-as.hexmode("25CE"),95),col=c('green3','green3','brown','brown','blue'))  

  ## Gizzard analysis
  # Averaging double measures on height and width
  Gizzard_data$H=rowMeans(Gizzard_data[c('Height1','Height2')],na.rm=TRUE)
  Gizzard_data$W=rowMeans(Gizzard_data[c('Width1','Width2')],na.rm=TRUE)
  
  # Calculating gizzard mass based on height and width
  Gizzard_data$mass=-1.09+3.78*Gizzard_data$H*Gizzard_data$W
  # Averaging gizzard mass per (measured) day
  Average_Gizzardmass=aggregate(mass~Experiment_Day,Gizzard_data,mean)
  # Calculate maximum shell processing rate based on average gizzard mass (in g/s)
  Average_Gizzardmass$maxshellint=10^-1.293*Average_Gizzardmass$mass^2
  # Calculate maximum shell processing rate per trial (30 minutes)
  Average_Gizzardmass$maxshellintptrial=Average_Gizzardmass$maxshellint*1800
  # Calculate maximum Loripes intake rate based on gizzard size and average shellmass per Loripes (same as intercept of best approximating model, which is the intercept-only model)
  Average_Gizzardmass$maxint=Average_Gizzardmass$maxshellintptrial/(mean(Loripes_shellmass$ShellMass)*1000)
  Average_Gizzardmass$Experiment_Day=c(2,6,9)
  # Adding the estimated maximum intake rates based on gizzard size
  points(maxint~Experiment_Day,Average_Gizzardmass,cex=5,pch=95,col='blue')
  # Third axis with the required gizzard mass
  y=c(5.5,6.0,6.5,7.0,7.5,8)
  ad=numeric(6)
  axis2=data.frame(y,ad)
  axis2$ad=(10^-1.293*axis2$y^2)*1800/(mean(Loripes_shellmass$ShellMass)*1000)
  axis(4,at=axis2$ad,labels=axis2$y,cex.axis=2)
  mtext("Required gizzard mass (g)", side=4, line=3,cex=2)

  ## Averaging drymass per treatment and location
  Drymass_table=summarySE(Loripes_drymass,measurevar='DryMass',groupvars=c('Treatment','Location'))
  Drymass_table2=data.frame(Drymass_table[,3:5])            # select useful columns
  Sulfur_table2=data.frame(Sulfur_table[,c(1:5)])           # select useful columns
  Sulfur_Drymass=cbind(Sulfur_table2,Drymass_table2)        # merge both
  names(Sulfur_Drymass)[c(3,5,6,8)]=c('NS','sdS','NDM','sdDM')
  
  # Calculate Variance of sulfur and drymass
  Sulfur_Drymass$varS=Sulfur_Drymass$sdS^2
  Sulfur_Drymass$varDM=Sulfur_Drymass$sdDM^2
  # Multiply percentage sulfur with drymass
  Sulfur_Drymass$DMS=Sulfur_Drymass$Percentage_Sulfur*Sulfur_Drymass$DryMass
  # Calculate variance of the product
  Sulfur_Drymass$varDMS=Sulfur_Drymass$varDM*Sulfur_Drymass$varS + Sulfur_Drymass$varDM*Sulfur_Drymass$Percentage_Sulfur^2 + Sulfur_Drymass$varS*Sulfur_Drymass$DryMass^2
  # Standard deviation out of the variance
  Sulfur_Drymass$sdDMS=sqrt(Sulfur_Drymass$varDMS)
  # Product in mg sulfur per Loripes:
  Sulfur_Drymass$DMS=Sulfur_Drymass$DMS*10 
  Sulfur_Drymass$sdDMS=Sulfur_Drymass$sdDMS*10
  # Calculate standard error (using N of sulfur measurements because it is the lowest)
  Sulfur_Drymass$seS=Sulfur_Drymass$sdDMS/sqrt(Sulfur_Drymass$NS)
  
  ## Figure 3, Intake rate against mg sulfur per Loripes
  par(mar=c(5,5,1,2))
  plot(Consumed_Loripes~sd,Intake_rate,cex=2,type='n',pch=16,xlim=c(0.34,0.63),ylim=c(31,35),
       ylab=expression(paste('Intake rate ',italic('I'),' (',italic('Loripes '),'per trial)')),
       xlab=expression(paste('mg sulfur per ',italic('Loripes'))),cex.lab=1.5,cex.axis=1.5)
  curve(12/(x),add=T,col='grey80')
  curve(15/(x),add=T,col='grey80')
  curve(18/(x),add=T,col='grey80')
  arrows(x0=Sulfur_Drymass$DMS,x1=Sulfur_Drymass$DMS,y0=Intake_rate$Consumed_Loripes-Intake_rate$se,y1=Intake_rate$Consumed_Loripes+Intake_rate$se,code=3,angle=90,length=0.1)
  arrows(x0=Sulfur_Drymass$DMS-Sulfur_Drymass$seS,x1=Sulfur_Drymass$DMS+Sulfur_Drymass$seS,y0=Intake_rate$Consumed_Loripes,y1=Intake_rate$Consumed_Loripes,code=3,angle=90,length=0.1)
  points(Sulfur_Drymass$DMS[1],Intake_rate$Consumed_Loripes[1],cex=3,pch=16,col='green3')
  points(Sulfur_Drymass$DMS[3],Intake_rate$Consumed_Loripes[3]+0.028,cex=3,pch=-as.hexmode("25CE"),col='green3')  # +0.028 to get the center of the circle exactly at the intersection of the arrow bars
  points(Sulfur_Drymass$DMS[2],Intake_rate$Consumed_Loripes[2],cex=3,pch=16,col='brown')
  points(Sulfur_Drymass$DMS[4],Intake_rate$Consumed_Loripes[4]+0.028,cex=3,pch=-as.hexmode("25CE"),col='brown')
  legend('topright',legend=c('Control, Location A','Treatment, Location A','Control, Location B','Treatment, Location B'),pch=c(16,-as.hexmode("25CE"),16,-as.hexmode("25CE")),col=c('green3','green3','brown','brown'),cex=1.2,pt.cex=1.8)
  text(12/31+0.01,31,12,col='grey80',cex=1.5)
  text(15/31+0.01,31,15,col='grey80',cex=1.5)
  text(18/31+0.01,31,18,col='grey80',cex=1.5)
  arrows(x0=Sulfur_Drymass$DMS[1]-0.007,x1=Sulfur_Drymass$DMS[3]+0.009,y0=Intake_rate$Consumed_Loripes[1]+0.11,y1=Intake_rate$Consumed_Loripes[3]-0.12,lwd=2,length=0.3,col='slateblue2')
  arrows(x0=Sulfur_Drymass$DMS[2]-0.007,x1=Sulfur_Drymass$DMS[4]+0.01,y0=Intake_rate$Consumed_Loripes[2]+0.11,y1=Intake_rate$Consumed_Loripes[4]-0.1,lwd=2,length=0.3,col='slateblue2')
  
  ## Figure 4, Loripes consumption at start of captivity
  par(mfrow=c(1,3))
  par(oma=c(5,5,1,1))
  par(mar=c(0.5,0.5,0.5,0.5))
  
  boxplot(Bill_length~Consuming,subset(Loripes_consumers,Year==2015),ylim=c(31,39.5),col=c('slategray3','slategray3'),xlab='',ylab='',cex.axis=2)
  text(1.5,39.2,'2015',cex=2)
  text(1.5,38.75,paste('N =',length(subset(Loripes_consumers,Year==2015)$Bill_length)),cex=1.5)
  boxplot(Bill_length~Consuming,subset(Loripes_consumers,Year==2017),ylim=c(31,39.5),col=c('slategray3','slategray3'),yaxt='n',cex.axis=2)
  text(1.5,39.2,'2017',cex=2)
  text(1.5,38.75,paste('N =',length(subset(Loripes_consumers,Year==2017)$Bill_length)),cex=1.5)
  boxplot(Bill_length~Consuming,subset(Loripes_consumers,Year==2018),ylim=c(31,39.5),col=c('slategray3','slategray3'),yaxt='n',cex.axis=2)
  text(1.5,39.2,'2018',cex=2)
  text(1.5,38.75,paste('N =',length(subset(Loripes_consumers,Year==2018)$Bill_length)),cex=1.5)
  mtext(expression(paste(italic('Loripes'),' consumption')), side=1, outer=T, cex=2, line=3.3)
  mtext('Bill length (mm)',side=2, outer=T, cex=2, line=3)

  with(Loripes_consumers,shapiro.test(Bill_length[Consuming=='Yes']))
  with(Loripes_consumers,shapiro.test(Bill_length[Consuming=='No']))
  var.test(Bill_length~Consuming,Loripes_consumers)
  t.test(Bill_length~Consuming,Loripes_consumers,var.equal=TRUE)