---
title: "Data from: Meroplankton diversity, seasonality and life-history traits across the Barents Sea Polar Front revealed by high-throughput DNA barcoding"
author: "Raphaelle Descoteaux"
date: "09/05/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Loading packages
```{r}

options(tidyverse.quiet = TRUE)
library(tidyverse)
library(vegan)
library(ggpubr)
library(ggdendro)
library(ggrepel)
library(car)

```

## Loading data
```{r}
# Individual larvae data
Data_individuals <- read.csv("MeroInd_IndividualData.csv")

# Taxonomic assignments for each DNA sequence
Data_sequences <- read.csv("MeroInd_SequenceAssignments.csv")

# Decisions whether a DNA-based ID matches the visual ID
Data_match <- read.csv("MeroInd_IDMatch.csv")%>%
  select(-Number)

# Individual larvae data with DNA-based ID and match decision
Data_ind_match <- Data_individuals%>%
  left_join(Data_sequences, by = "Sequence")%>%
  left_join(Data_match, by = c("Visual.ID", "Best.ID", "Visual.ID.comments"))%>%
  mutate(ID = ifelse(Match == "Control", "Control", ifelse((Total.reads < 1000 & Percent.highest.seq < 75)| Total.reads <10, "Insufficient reads", ifelse(Match == "No", "Unknown", Best.ID)))) # Identify samples with insufficient reads (all samples with < 10 total reads as well as samples with less than 1000 reads AND less than 75% reads to most abundant sequence) and samples for which the visual ID and DNA-based ID do not match.
Data_ind_match$Total.reads <- as.numeric(Data_ind_match$Total.reads)
Data_ind_match$Percent.highest.seq <- as.numeric(Data_ind_match$Percent.highest.seq)
Data_ind_match$Best.ID.percent <- round(Data_ind_match$Best.ID.percent, digits = 1)

# Visualize distribution of samples according to number and percent of reads
# Used to make decision on cutoffs for considering a larva as having 'insufficient reads'
ggplot(Data_ind_match, aes(color = Match, shape = Meroplankton, y=Total.reads, x=Percent.highest.seq))+
  geom_point()+
  scale_y_continuous(trans='log10')

# Remove controls and include % match in the name of the larvae when <98% match to database
Data_nocontrols <- Data_ind_match%>%
  filter(ID!="Control")%>%
  mutate(IDpercent = ifelse(ID != "Unknown" & ID!= "Insufficient reads" & Best.ID.percent <98, Best.ID.percent, ""))%>%
  mutate(percent = ifelse(ID != "Unknown" & ID!= "Insufficient reads" & Best.ID.percent<98, "%", ""))%>%
  unite(Final.ID, ID, IDpercent, percent, remove = FALSE, sep = " ")

# Order months chronologically
Data_nocontrols$Month <- factor(Data_nocontrols$Month, levels = c("November", "January", "April", "June", "August"))

# Quantitative data based on visual ID 
Data_quant <- read.csv("MeroInd_Quantification.csv", na.strings= "n/a")
Data_quant$Total.subs <- as.numeric(Data_quant$Total.subs)
Data_quant<- Data_quant%>%
  mutate(Depth.range = Depth.max-Depth.min)%>%
  mutate(Volume.filtered = Depth.range * Net.opening)%>%
  mutate(Subs.ratio = (1-Previous.subs)*(Subs.volume/Sieved.volume))%>%
  mutate(Total.sample = Total.subs/Subs.ratio)%>%
  mutate(Density = Total.sample/Volume.filtered)%>%
  select(Month,Location, Depth.max, Depth.min, Visual.ID, Visual.ID.standardized, Density)

# Order months chronologically
Data_quant$Month <- factor(Data_quant$Month, levels = c("November", "January", "April", "June", "August"))

# Environmental variables
## MeroInd_Env.csv is not provided here due to copyright issues. The datasets may be available on request from the data owners. 
## These data are only necessary for producing the temperature-salinity (TS) plot as well as for constraining variables in the canonical correspondence analysis (CCA). The rest of the code should function well without it.
## An already-produced MeroInd_Env_Summary datasheet, without the mean temperature, salinity and fluorescence a values is provided here for reference. It can also be used to run the CCA analysis with only Location, Month, Layer and Water Mass as constraining variables.
## If you wish to recreate the MeroInd_Env.csv datasheet, make sure to include Location (North or South), Month (November, January, April, June, August), Layer (Surface, Intermediate1, Intermediate2, Deep), Depth.max (See provided MeroInd_Env_Summary for list of Depth.max per station) and well as Depth, Salinity, Temperature and Fluorescence obtained from the CTD files.

#Data_env <- read.csv("MeroInd_Env.csv")%>%
#  filter(Layer!="Unused")%>%
#  unite(Month, Layer, col = Month.layer, sep = "_", remove = FALSE)

#Data_env_summary <-  Data_env%>%
#  group_by(Location, Month, Layer, Depth.max)%>%
#  summarize(mean_sal=mean(Salinity), mean_temp=mean(Temperature), mean_fluo = mean(Fluorescence))%>%
#  arrange(Location, Month, Depth.max)%>%
#  mutate(WaterMass = if(mean_temp>3 & mean_sal>34.8){("AW")
#   } else if(mean_temp<0 & mean_sal<34.7){("ArW")
#   } else if(mean_temp<=2 & mean_sal>34.8){("BSW")
#   } else if(mean_sal<34.3){("SW")
#   } else if(mean_sal>34.8){("AW-BSW")
#   } else if(mean_sal<34.8){("ArW-BSW")
#   } else{"Mix"})
#Data_env_summary$Month <- dplyr::recode(Data_env_summary$Month, November = "Nov", January = "Jan", April = "Apr", June = "Jun", August = "Aug")

Data_env_summary <- read.csv("MeroInd_Env_Summary.csv") # If not using the code above to produce 
```

## TS plot
```{r}
# Requires the Data_env object which is not provided here due to copyright issues on CTD data files (see above).
# Note that the published legend was produced separately

# ggplot(Data_env) + 
#   geom_point(aes(x=Salinity, y=Temperature, color = Month.layer, shape = Location), size = 3) +
#   scale_color_manual(values = c("chartreuse4", "chartreuse3", "chartreuse2", "chartreuse1", "orangered4", "orangered3", "orangered2", "orangered1", "deepskyblue4", "deepskyblue3", "deepskyblue1", "gold3", "gold2", "gold1", "mediumorchid4", "mediumorchid3", "mediumorchid1")) + 
#   theme_bw()+
#   theme(panel.grid= element_blank())+
#   annotate("rect", xmin = 34.8, xmax = 35.2, ymin = 3, ymax = 7, alpha=0.05, fill="red", color = "red")+
#   annotate("rect", xmin = 34, xmax = 34.7, ymin = -2, ymax = 0, alpha=0.05, fill="deepskyblue", color = "deepskyblue")+
#   annotate("rect", xmin = 34.8, xmax = 35.2, ymin = -2, ymax = 2, alpha=0.05, fill="black", color = "black")+
#   annotate("text", x = 35, y = 6.5, label = "Atlantic Water")+
#   annotate("text", x = 35, y = -1.4, label = "Barents Sea \nWater")+
#   annotate("text", x = 34.2, y = -1.6, label = "Arctic Water")+
#   theme(axis.text.x = element_text(size = 10))+
#   theme(axis.text.y = element_text(size = 10))+
#   theme(axis.title = element_text(size = 12))+
#   theme(legend.position = "right")+
#   scale_y_continuous(name = "Temperature (°C)")

```


## Calculate abundance (Ind/m^3 and Ind/m^2) and richness
```{r Abundance and richness}

# Calculate how many larvae were DNA barcoded from each morphotype in each sample
Abund_visual_ID<- Data_nocontrols%>%
  group_by(Month, Location, Depth.max, Visual.ID)%>%
  summarize(No.visual= n())

# Calculate how many larvae belonged to each DNA barcode from each sample and morphotype
Abund_final_ID<- Data_nocontrols%>%
  group_by(Month, Location, Depth.max, Visual.ID, Final.ID, Meroplankton)%>%
  summarize(No.ind= n())

# Calculate the percentage of each morphotype made up of each DNA barcode in each sample
Abund_percent <- full_join(Abund_visual_ID,Abund_final_ID,by=c("Month","Location", "Depth.max", "Visual.ID"))%>%
  mutate(Percent = No.ind/No.visual)%>%
  arrange(Month, Location, Depth.max)

# Calculate the density (Ind/m^3) of each DNA barcode in each sample by combining quantification data with DNA barcoding data
Abund_sample <- full_join(Abund_percent,Data_quant,by=c("Month","Location", "Depth.max", "Visual.ID"))%>%
  arrange(Location, Month, Depth.max, Visual.ID)%>%
  mutate(Density.taxa = Percent*Density)%>%
  group_by(Month, Location, Depth.max, Depth.min, Final.ID, Meroplankton)%>%
  summarize(Density.combined = sum(Density.taxa))%>%
  filter(Meroplankton=="Yes")%>%
  filter(Final.ID != "Insufficient reads  " & Final.ID != "Unknown  ")%>%
  ungroup()%>%
  select(-Meroplankton)
Abund_sample$Month <- dplyr::recode(Abund_sample$Month, November = "Nov", January = "Jan", April = "Apr", June = "Jun", August = "Aug")
Abund_sample$Final.ID <- dplyr::recode(Abund_sample$Final.ID, "Pectinaria hyperborea  " = "Cistenides hyperborea", "Dendronotus niveus  " = "Dendronotus elegans") # Accepted name

# Calculate total abundance per sample
Abund_sample_total <- Abund_sample%>%
  group_by(Month, Location, Depth.max)%>%
  summarize(Density.total=sum(Density.combined))%>%
  arrange(Location, Month)%>%
  ungroup()%>%
  add_row(Month = "Jun", Location = "North", Depth.max = 15, Density.total = 0)

# Multiply density (Ind/m^3) by depth range (m) to get abundance integrated over whole water column (Ind/m^2)
Abund_station <- Abund_sample%>%
  group_by(Month, Location, Final.ID)%>%
  summarize(Abund=sum(Density.combined*(Depth.max-Depth.min)))

# Calculate total abundance per station
Abund_station_total <- Abund_station%>%
  group_by(Month, Location)%>%
  summarize(Abund.total=sum(Abund))%>%
  arrange(Location, Month)

# Add taxonomy
Taxo_mero <- Data_nocontrols%>%
  select(Final.ID, Phylum, Class, Order, Family, Genus, Species)%>%
  distinct()
Taxo_mero$Final.ID <- dplyr::recode(Taxo_mero$Final.ID, "Pectinaria hyperborea  " = "Cistenides hyperborea", "Dendronotus niveus  " = "Dendronotus elegans") # Accepted name

Abund_taxo <- Abund_station%>%
  left_join(Taxo_mero, by = "Final.ID")

# Table 
Abund_taxo_table <- Abund_taxo%>%
  unite(Sample, Location, Month, remove = TRUE)%>%
  spread(Sample, Abund, fill = 0)%>%
  relocate(Final.ID, Phylum, Class, Order, Family, Genus, Species, North_Nov, North_Jan, North_Apr, North_Jun, North_Aug, South_Jan, South_Apr, South_Jun, South_Aug)%>%
  arrange(Phylum, Class, Order, Family, Genus, Species)

# Calculate abundance per phylum
Abund_phylum <- Abund_taxo%>%
  group_by(Month, Location, Phylum)%>%
  summarize(Total=sum(Abund))

# Calculate abundance per taxon
Abund_taxa <- Abund_taxo%>%
  group_by(Month, Location, Phylum, Class, Final.ID)%>%
  summarize(Total=sum(Abund))%>%
  mutate(HigherTaxon = ifelse(Phylum == "Mollusca", Class, Phylum))

# Calculate number of taxa
Richness <- Abund_taxo%>%
  group_by(Month, Location, Phylum)%>%
  summarize(NoTaxa=n())

# Richness per sample
Richness_sample <- Abund_sample%>%
  group_by(Month, Location, Depth.max)%>%
  summarize(NoTaxa = n())

# Richness per site
Richness_site <- Abund_taxo%>%
  group_by(Location)%>%
  summarize(NoTaxa = n_distinct(Final.ID))

# Shared taxa
Richness_site_shared <- Abund_taxa%>%
  group_by(Final.ID, Location)%>%
  summarize(Abund = sum(Total))%>%
  spread(Location, Abund, fill = 0)%>%
  mutate(Location = ifelse(North > 0 & South > 0, "Shared", ifelse(North >0, "OnlyNorth", "OnlySouth")))%>%
  group_by(Location)%>%
  summarize(NoTaxa = n())

# Richness per site per month(integrated over whole water column)
Richness_site_month <- Richness%>%
  group_by(Month, Location)%>%
  summarize(No = sum(NoTaxa))%>%
  arrange(Location, Month)

# Appendix table
TableS1 <- Abund_taxo_table%>%
  select(-Order, -Family, -Genus, -Species)

# Richness per phylum
Richness_phylum <- TableS1 %>%
  group_by(Phylum, Class)%>%
  summarize(NoTaxa = n())

# Total number of taxa:
n_distinct(Abund_taxo$Final.ID)

```

## Checking success rate of DNA barcoding
```{r}

# Number of individuals from each standardized visual ID that were successfully barcoded, mismatched (DNA ID not fitting the visual ID), or had insufficient DNA reads
DNA_success_category <-  Data_nocontrols%>%
  mutate(DNA.success = ifelse(Final.ID == "Insufficient reads  ", "Insufficient reads", ifelse(Best.ID == "No match", "No match", ifelse(Final.ID == "Unknown  ", "Mismatch", "Success"))))%>%
  group_by(Visual.ID.standardized, DNA.success)%>%
  summarize(No= n())

# Total number of individuals within each standardized visual ID category
DNA_success_total <-  Data_nocontrols%>%
  group_by(Visual.ID.standardized)%>%
  summarize(Total= n())

# Percentage of each DNA success category for each standardized visual ID category
DNA_success_percent <- merge(DNA_success_category,DNA_success_total,by=c("Visual.ID.standardized"))%>%
  mutate(Percent = No/Total)

# Dataframe above in table format
DNA_success_percent_wide <-DNA_success_percent%>%
  select(Visual.ID.standardized,DNA.success, Percent)%>%
  spread(DNA.success, Percent, fill= 0, drop=FALSE)%>%
  merge(DNA_success_total, by = "Visual.ID.standardized")%>%
  filter(Visual.ID.standardized %in% c("Bivalve", "Bryozoa", "Echinodermata", "Gastropoda", "Nemertea", "Polychaete", "Arthropoda"))%>%
  relocate("Visual.ID.standardized", "Total", "Success", "Insufficient reads", "No match", "Mismatch")

# DNA_Success_overall:
summarize(DNA_success_percent_wide, Overall = (sum(Total*Success)/sum(Total))*100)

# Number of individuals that were barcoded:
DNA_success_percent%>%
  summarize(Overall = sum(No))

# Number of barcoded individuals which were successful
DNA_success_percent%>%
  filter(DNA.success == "Success")%>%
  summarize(Overall = sum(No))

# Percentage of sequences with > 98%
DNA_match <- Data_nocontrols%>%
  mutate(DNA.success = ifelse(Final.ID == "Unknown  ", "Mismatch", ifelse(Final.ID == "Insufficient reads  ", "Insufficient reads", "Success")))%>%
  select(Final.ID, Best.ID.percent, Best.ID.source, Meroplankton, DNA.success)%>%
  filter(DNA.success == "Success" & Meroplankton == "Yes")%>%
  unique()
(nrow(DNA_match[DNA_match$Best.ID.percent >98,])/nrow(DNA_match))*100

# Calculate percent identification at each taxonomic resolution for MOTUs with >98% sequence match
Data_nocontrols%>%
  filter(Meroplankton == "Yes" & Match == "Yes" & Final.ID != "Insufficient reads  ")%>%
  select(Sequence, Final.ID, Best.ID.percent, Resolution)%>%
  unique()%>%
  left_join(Abund_taxo_table, by = "Final.ID")%>%
  filter(Best.ID.percent >= 98)%>%
  group_by(Resolution)%>%
  summarize(No = n())
```

## Make abundance and richness plots at phylum level
```{r}

# Abundance
Abundance_plot <- ggplot(Abund_phylum, aes(fill=Phylum, y=Total, x=Month)) + 
  geom_bar(position= position_dodge2(preserve = "single"), stat="identity")+
  facet_wrap(~Location, scales = "free")+
  scale_fill_manual(values = c("#009E73", "#56B4E9", "#E69F00", "#D55E00", "#CC79A7", "#0072B2", "#F0E442", "#999999"))+
  scale_y_continuous(name=bquote('Abundance'~(ind/m^2)), labels = scales::comma, trans = "sqrt", breaks = c(0, 100, 1000, 3000, 6000, 10000, 50000, 100000, 200000, 300000), limits = c(0, NA))+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Abundance")+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 10))+
  theme(legend.title = element_text(size = 12))+
  theme(strip.text.x = element_text(size = 12))

# Richness
Richness_plot <- ggplot(Richness, aes(fill=Phylum, y=NoTaxa, x=Month)) + 
  geom_bar(position= position_dodge2(preserve = "single"), stat="identity")+
  facet_wrap(~Location, scales = "free")+
  scale_fill_manual(values = c("#009E73", "#56B4E9", "#E69F00", "#D55E00", "#CC79A7", "#0072B2", "#F0E442", "#999999"))+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Richness")+
  scale_y_continuous(name = "Number of taxa", breaks = c(0, 5, 10), limits = c(0, 13))+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 10))+
  theme(legend.title = element_text(size = 12))+
  theme(strip.text.x = element_text(size = 12))

# Figure 4 - Phyla abundance and richness 
ggarrange(Abundance_plot, Richness_plot, labels = c("(A)", "(B)"), ncol = 1, nrow = 2, common.legend = TRUE, legend = "bottom", align = "v")

```

## Make abundance plots at species level
```{r}

# Calculate total number of each higher taxon (Phylum or Class) at each sampling event
Abund_taxa_cumul <- Abund_taxa%>%
  group_by(Month, Location, HigherTaxon)%>%
  summarize(Cumul = sum(Total))

# Calculate percentage of each higher taxon made up of each species
Abund_taxa_percent <- Abund_taxa%>%
  unite(sample, Location, Month, remove = TRUE)%>%
  spread(sample, Total, fill=0)%>%
  gather(sample, Total, 5:13)%>%
  separate(sample, c("Location", "Month"), sep = "_")%>%
  left_join(Abund_taxa_cumul, by = c("Month", "Location", "HigherTaxon"))%>%
  replace_na(list(Cumul = 0))%>%
  mutate(Percent = ifelse(Cumul == 0, 0, (Total/Cumul)*100))%>%
  mutate(ID = ifelse((HigherTaxon == "Annelida" | HigherTaxon == "Gastropoda") & Percent < 10, "Other", Final.ID))
Abund_taxa_percent$Month <- factor(Abund_taxa_percent$Month, levels = c("Nov", "Jan", "Apr", "Jun", "Aug"))

# Echinodermata
Echinoderm <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Echinodermata")
Echinoderm_totals <- Echinoderm%>%
  group_by(Month, Location)%>%
  summarize(Abundance=sum(Total))
Echinoderm_Prop <- Echinoderm_totals%>%
  mutate(PropAbund=Abundance/(max(Echinoderm_totals$Abundance)))
Echinoderm$Month <- factor(Echinoderm$Month, levels = c("Nov", "Jan", "Apr", "Jun", "Aug"))
Echinoderm_max = max(Echinoderm_totals$Abundance)

Echinoderm_plot <- ggplot()+
  geom_bar(data=Echinoderm, aes(y=Total, x=Month, fill = ID), position="fill", stat="identity")+
  scale_fill_manual(values = c("#FE6100", "#31D969", "#E51031", "#92A4CC", "#648FFF", "#785EF0","#FFB000", "#A631BF", "#E1DFAF", "#F9AEBF", "#8AB17B"))+
  geom_line(data = Echinoderm_Prop, aes(x=Month, y = PropAbund, group = 1), size =1)+
  geom_point(data = Echinoderm_Prop, aes(x=Month, y = PropAbund), size = 3)+
  facet_wrap(~Location)+
  scale_y_continuous("Proportion of taxa", sec.axis = sec_axis(~ . * Echinoderm_max, name=bquote('Abundance'~(ind/m^2))), labels = scales::comma)+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Echinodermata")+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(axis.title.y.right = element_text(angle = 90))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 8))+
  theme(legend.title = element_text(size = 10))+
  theme(strip.text.x = element_text(size = 12))+
  theme(legend.key.size = unit(0.3, "cm"))+
  labs(fill = "Taxa")
 
# Annelida
Annelid <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Annelida")
Annelid_totals <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Annelida")%>%
  group_by(Month, Location)%>%
  summarize(Abundance=sum(Total))
Annelid_Prop <- Annelid_totals%>%
  mutate(PropAbund=Abundance/(max(Annelid_totals$Abundance)))
Annelid_max = max(Annelid_totals$Abundance)

Annelid_plot <- ggplot()+
  geom_bar(data=Annelid, aes(y=Total, x=Month, fill = ID), position="fill", stat="identity")+
  scale_fill_manual(values = c("#DC267F", "#648FFF", "#785EF0", "#92A4CC", "#1E5FE6", "#F9AEBF", "#FE6100", "#FFB000","#A9B69B", "#E51031", "#BC530D", "#8AB17B", "#31D969"))+
  geom_line(data = Annelid_Prop, aes(x=Month, y = PropAbund, group = 1), size =1)+
  geom_point(data = Annelid_Prop, aes(x=Month, y = PropAbund), size = 3)+
  facet_wrap(~Location)+
  scale_y_continuous("Proportion of taxa", sec.axis = sec_axis(~ . * Annelid_max, name=bquote('Abundance'~(ind/m^2))), labels = scales::comma)+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Annelida")+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(axis.title.y.right = element_text(angle = 90))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 8))+
  theme(legend.title = element_text(size = 10))+
  theme(strip.text.x = element_text(size = 12))+
  theme(legend.key.size = unit(0.3, "cm"))+
  labs(fill = "Taxa")

# Bivalvia
Bivalvia <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Bivalvia")
Bivalvia_totals <- Bivalvia%>%
  group_by(Month, Location)%>%
  summarize(Abundance=sum(Total))
Bivalvia_Prop <- Bivalvia_totals%>%
  mutate(PropAbund=Abundance/(max(Bivalvia_totals$Abundance)))
Bivalvia_max = max(Bivalvia_totals$Abundance)

Bivalvia_plot <- ggplot()+
  geom_bar(data=Bivalvia, aes(y=Total, x=Month, fill = ID), position="fill", stat="identity")+
  scale_fill_manual(values = c("#DC267F", "#FE6100", "#785EF0", "#648FFF", "#FFB000", "#A631BF", "#57E07D"))+
  geom_line(data = Bivalvia_Prop, aes(x=Month, y = PropAbund, group = 1), size =1)+
  geom_point(data = Bivalvia_Prop, aes(x=Month, y = PropAbund), size = 3)+
  facet_wrap(~Location)+
  scale_y_continuous("Proportion of taxa", sec.axis = sec_axis(~ . * Bivalvia_max, name=bquote('Abundance'~(ind/m^2))), labels = scales::comma)+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Bivalvia")+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(axis.title.y.right = element_text(angle = 90))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 8))+
  theme(legend.title = element_text(size = 10))+
  theme(strip.text.x = element_text(size = 12))+
  theme(legend.key.size = unit(0.3, "cm"))+
  labs(fill = "Taxa")

# Gastropoda
Gastropoda <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Gastropoda")
Gastropoda_totals <- Abund_taxa_percent%>%
  filter(HigherTaxon == "Gastropoda")%>%
  group_by(Month, Location)%>%
  summarize(Abundance=sum(Total))
Gastropoda_Prop <- Gastropoda_totals%>%
  mutate(PropAbund=Abundance/(max(Gastropoda_totals$Abundance)))
Gastropoda_max = max(Gastropoda_totals$Abundance)

Gastropoda_plot <- ggplot()+
  geom_bar(data=Gastropoda, aes(y=Total, x=Month, fill = ID), position="fill", stat="identity")+
  scale_fill_manual(values = c("#DC267F", "#648FFF", "#785EF0", "#92A4CC", "#1E5FE6", "#F9AEBF", "#FE6100", "#FFB000","#A9B69B", "#E51031", "#BC530D", "#8AB17B", "#31D969", "grey"))+
  geom_line(data = Gastropoda_Prop, aes(x=Month, y = PropAbund, group = 1), size =1)+
  geom_point(data = Gastropoda_Prop, aes(x=Month, y = PropAbund), size = 3)+
  facet_wrap(~Location)+
  scale_y_continuous("Proportion of taxa", sec.axis = sec_axis(~ . * Gastropoda_max, name=bquote('Abundance'~(ind/m^2))), labels = scales::comma)+
  theme_bw()+
  theme(panel.grid= element_blank())+
  ggtitle("Gastropoda")+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(axis.title.y.right = element_text(angle = 90))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 8))+
  theme(legend.title = element_text(size = 10))+
  theme(strip.text.x = element_text(size = 12))+
  theme(legend.key.size = unit(0.3, "cm"))+
  labs(fill = "Taxa")

# Combine into one Fig 
# Figure 5 - Taxon-level abundances
ggarrange(Bivalvia_plot, Gastropoda_plot, Echinoderm_plot, Annelid_plot, labels = c("(A)", "(B)", "(C)", "(D)"), font.label = list(size = 12), ncol = 1, nrow = 4, align = "v")

```

## Making heatmap
```{r, fig.height= 18, fig.width = 11}

# Make abundance data into a matrix. Only keep taxa that make up at least 1% of at least one sample.
Dendro_df <- Abund_taxa%>%
  left_join(Abund_station_total, by = c("Month", "Location"))%>%
  mutate(Percent = (Total/Abund.total)*100)%>%
  group_by(Final.ID)%>%
  summarize(Max = max(Percent))%>%
  right_join(Abund_taxa, by = "Final.ID")%>%
  filter(Max >= 1)%>%
  mutate(ID = Final.ID)%>%
  ungroup()%>%
  select(Location, Month, Phylum, ID, Total)%>%
  unite(Station, Month, Location, sep = ".")%>%
  spread(Station, Total, fill = 0)
  
Dendro_matrix <- as.matrix(Dendro_df[,-c(1:2)])
rownames(Dendro_matrix) <- Dendro_df$ID

# Make a matrix that quantifies similarity between each pair of taxa
Dendro_matrix_sim <- vegdist(Dendro_matrix^0.25, method ="bray")

# Single linkage clustering
Dendro_clust<- hclust(Dendro_matrix_sim, method = "average")
summary(Dendro_clust)
Dendrogram <- as.dendrogram(Dendro_clust)

Dendro_model <- dendro_data(Dendrogram, type = "rectangle")

Fig6a <- ggplot(segment(Dendro_model)) + 
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) + 
  coord_flip() + 
  scale_y_reverse(expand = c(0.2, 0))+
  #geom_text(data = Dendro_data$labels, aes(x = x, y = y, label = label)) + 
  theme_dendro()

Dendro_order <- as.vector(Dendro_model$labels$label)  
Dendro_df$ID <- factor(x = Dendro_df$ID,
                               levels = (Dendro_order), 
                               ordered = TRUE)

# Making heatmap
Heatmap_df <- Dendro_df%>%
  gather(sample, abundance, 3:11)%>%
  separate(sample, c("month", "location"))%>%
  filter(abundance != 0)
Heatmap_df$month <- factor(Heatmap_df$month, levels = c("Nov", "Jan", "Apr", "Jun", "Aug"))

Fig6b <- ggplot(Heatmap_df, aes(x = month, y = ID))+
  geom_tile(aes(fill = abundance), color = "white") +
  facet_wrap(~location, scales = "free_x")+
  theme_bw()+
  theme(panel.grid= element_blank())+
  scale_fill_gradient(low = "#e9f0f6",high = "#1b3246", trans = "log", breaks = c(10, 100, 1000, 10000, 100000), labels = c(10, 100, 1000, 10000, 100000))+
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(axis.title.y = element_blank())+
  labs(fill = bquote('Abundance'~(ind/m^2)))+
  xlab("Month")+
  theme(legend.position="right")+
  theme(legend.text = element_text(size = 10))+
  theme(legend.title = element_text(size = 12))+
  theme(strip.text.x = element_text(size = 12))

# Phylum color bar
Fig6c <- ggplot(Heatmap_df, aes(x = 1, y = ID), width = 0.5)+  
  geom_tile(aes(fill = Phylum))+
  scale_fill_manual(values = c("#009E73", "#56B4E9", "#E69F00", "#D55E00", "#CC79A7", "#0072B2", "#F0E442", "#999999"))+
  xlab("Phylum")+
  theme_void()+
  theme(panel.grid= element_blank())+
  theme(legend.position="right")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank())

Fig6_legend_color_bar <- get_legend(Fig6c)
Fig6_legend_heatmap <-  get_legend(Fig6b)
Fig6_legends <- ggarrange(Fig6_legend_heatmap, Fig6_legend_color_bar, nrow = 2)

# Figure 6 - Heatmap
ggarrange(Fig6a, Fig6c, Fig6b, Fig6_legends,  ncol = 4, nrow = 1, widths = c(2, 0.5, 10, 2.5), legend = "none")

```

## Community analysis
```{r}
# Note that the CCA analysis relies on the environmental variables which will differ slightly depending on whether you are using the provided Data_env_summary or whether you made your own using CTD data provided separately by the data users. Overall, however, results remain the same whether mean_temp, mean_sal and mean_fluo are included or not. The code below is built to exclude mean_temp, mean_sal and mean_fluo but these just need to be added to two lines of code highlighted below to adapt to the complete Data_env dataframe.


# Make dataframe with combined larval abundances and environmental factors
Commun_BioEnv <- Abund_sample%>%
  spread(Final.ID, Density.combined, fill = 0)%>%
  mutate(Depth = as.numeric(ifelse(Depth.max == "321", "285", Depth.max)))%>% # Correcting mismatch between Depth.max in Abund_sample and in Data_env_summary
  select(-Depth.max)%>%
  rename(Depth.max = Depth)%>%
  left_join(Data_env_summary, by = c('Location', 'Month', 'Depth.max'))%>%
  filter(Location!='North'|Month!='Jan'|Depth.max!='191')%>% # This station only has Bohuslania matsmichaeli and is throwing off the analysis
  mutate(Loc = substr(Location, 1,1), Mon = substr(Month, 1, 3))%>%
  select(-"Bohuslania matsmichaeli  ")%>%
  unite(Name, Mon, Loc, Depth.max, sep = ".")
rownames(Commun_BioEnv) <- Commun_BioEnv$Name  

# Only environmental factors
Commun_Env <- Commun_BioEnv%>%
  select(Location, Month, Layer, WaterMass) # Add mean_temp, mean_sal and mean_fluo if working with the CTD data (see notes in "Loading data" chunk)
Commun_Env$Month <- factor(Commun_Env$Month, levels = c("Nov", "Jan", "Apr", "Jun", "Aug"))
Commun_Env$Layer <- factor(Commun_Env$Layer, levels = c("Surface", "Intermediate1", "Intermediate2", "Deep"))
Commun_Env$Location <- factor(Commun_Env$Location, levels = c("North", "South"))

# Only larval abundances
Bio <- Commun_BioEnv%>%
  ungroup()%>%
  select(-"Month", -"Location", - "Name", -"Depth.min", -"Layer", -"WaterMass") # Include mean_temp, mean_sal and mean_fluo if working with the CTD data
rownames(Bio) <- Commun_BioEnv$Name 


# Figure 3 - Species accumulation curve
accum <- specaccum(Bio)
accum_data <- data.frame(Sites=accum$sites, Richness=accum$richness, SD=accum$sd)

ggplot() +
  geom_point(data=accum_data, aes(x=Sites, y=Richness)) + 
  geom_line(data=accum_data, aes(x=Sites, y=Richness)) +
  geom_ribbon(data=accum_data ,aes(x=Sites, ymin=(Richness-2*SD), ymax=(Richness+2*SD)), alpha=0.2)+
  theme(axis.text = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme_bw()+
  theme(panel.grid= element_blank())+
  xlab("Number of samples")+
  ylab("Taxon richness")


# CCA stepwise selection
cca_all <- cca(Bio^0.25~., data=Commun_Env)
anova (cca_all) # P < 0.01 - yes, significant!
cca_0 <- cca(Bio^0.25~ 1, data=Commun_Env) # empty model only with intercept
cca_step_select <- ordistep(cca_0, scope = formula (cca_all), direction = 'both', permutations = 49999, trace = T) #stepwise selection 
cca_new <- cca_step_select
summary(cca_new)
cca_new$anova$'Pr(>F)' <- p.adjust(cca_new$anova$'Pr(>F)', method = 'holm') #adjusted p-value
cca_new$anova #still significant
summary(cca_new)
anova(cca_new)
plot(cca_new)

cca_scores<-scores(cca_new,display=c("sp","wa","lc","bp","cn"))

# Select only the species that are significantly driving the ordination to show on the plot
species_fit<-envfit(cca_scores$sites, Bio^0.25)
species_pvals <- as.data.frame(species_fit$vectors$pvals)%>%
  rownames_to_column(var = "Species")%>%
  rename(pvals = "species_fit$vectors$pvals")
species_list<-as.data.frame(cca_scores$species)%>%
  rownames_to_column(var = "Species")%>%
  left_join(species_pvals, by = "Species")%>%
  filter(pvals < 0.05)

#Plotting
sites <- as.data.frame(cca_scores$sites)%>%
  select(CCA1, CCA2)%>%
  rownames_to_column()%>%
  separate(rowname, c("Month", "Location", "Depth"))%>%
  mutate(WaterMass = Commun_Env$WaterMass, Layer = Commun_Env$Layer)
sites$Month <- factor(sites$Month, levels = c("Nov", "Jan", "Apr", "Jun", "Aug"))
sites$Location <- dplyr::recode(sites$Location, N = "North", S = "South")
sites$Month <- dplyr::recode(sites$Month, Nov = "November", Jan = "January", Apr = "April", Jun = "June", Aug = "August")

# Centroids for months
centroids <- as.data.frame(cca_scores$centroids)
centroids$Month_short <- as.factor(c("Nov", "Jan", "Apr", "Jun", "Aug"))
summary(cca_new)

# Figure 8 - CCA 
ggplot(sites, aes(x = CCA1, y = CCA2)) +
  geom_point(aes(shape=Location, color = Month), size = 4) +
  geom_text_repel(data=species_list, aes(label=Species, fontface=3), position=position_jitter(width=0.2,height=0.25), size=2)+
  geom_text(data=centroids, aes(label=Month_short, fontface=2), size=3, colour="blue")+
  theme_bw()+
  theme(panel.grid= element_blank())+
  scale_color_manual(values = c("mediumorchid3", "deepskyblue3", "chartreuse3", "gold2", "orangered3")) +
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 10))+
  theme(legend.title = element_text(size = 12))+
  labs(x="CCA1 (12.3%)", y="CCA2 (9.6%)")

# Size component - did not meet assumptions of ANOVA so use non-parametric test instead
kruskal.test(Density.total ~ Month, data = Abund_sample_total)
kruskal.test(Density.total ~ Location, data = Abund_sample_total)

# Richness
richness.aov <- aov(NoTaxa ~ Month + Location, data = Richness_sample)
Anova(richness.aov, type = "III") # For unbalanced designs
plot(richness.aov, 2)
residuals_richness <- residuals(richness.aov)
shapiro.test(residuals_richness)
plot(richness.aov, 1) # Visualize homogeneity of variance
leveneTest(No ~ Month * Location, data = Richness_site_month) # Levene's test for homogeneity of variance

```

## Size histograms
```{r, fig.height= 5}

# Filter for our two species of interest
Size <- filter(Data_nocontrols, Final.ID %in% c("Hiatella sp. K  ", "Ophiocten gracilis  ") & Length > 1)%>%
  select(Location, Month, Visual.ID, Final.ID, Length, Feature.measured)

# Calculate mean size at each sampling point
Size_mean <- Size%>%
  group_by(Location, Month, Final.ID)%>%
  summarize(Mean = mean(Length))

# Calculate how many larvae were measured at each location
Size_n <- Size%>%
  group_by(Location, Final.ID)%>%
  summarize(Observations=n())%>%
  mutate(N = "n")%>%
  unite(Obs, N, Observations, sep = " = ", remove = TRUE)%>%
  mutate(x = ifelse(Final.ID == "Hiatella sp. K  ", 330, 2500))%>%
  mutate(y= ifelse(Location == "North", 10, 40))%>%
  mutate(Month = "November")

# Figure 7 - Histograms
ggplot(Size, aes(x = Length, fill = Month)) + 
  geom_histogram(position= position_dodge2(preserve = "single"), bins = 15)+
  facet_grid(Location~Final.ID, scales = "free")+
  geom_vline(data=Size_mean, aes(xintercept=Mean, color=Month),linetype="dashed", size = 1)+
  scale_fill_manual(values = c("mediumorchid3", "deepskyblue3", "chartreuse3", "gold2", "orangered3"))+
  scale_color_manual(values = c("mediumorchid3", "deepskyblue3", "chartreuse3", "gold2", "orangered3"))+
  theme_bw()+
  theme(panel.grid= element_blank())+
  theme(legend.position="bottom")+
  scale_x_continuous(name="Length (µm)") +
  scale_y_continuous(name="Number of individuals")+  
  theme(axis.text.x = element_text(size = 10))+
  theme(axis.text.y = element_text(size = 10))+
  theme(axis.title = element_text(size = 12))+
  theme(plot.title = element_text(size = 12))+
  theme(legend.text = element_text(size = 10))+
  theme(legend.title = element_text(size = 12))+
  theme(strip.text.x = element_text(size = 12))+
  theme(strip.text.y = element_text(size = 12, angle = 0))+
  geom_text(data = Size_n, mapping = aes(x = x, y = y, label = Obs))

```