---
title: "Eelgrass pathogen synthesizes essential fatty acids"
author:
- Yoshioka RM^1^\*, Schram JB^1^, Galloway AWE^1^
- ^1^ University of Oregon Institute of Marine Biology, Charleston, Oregon, 97420, USA
- '*corresponding author: rmy@uoregon.edu'
date: '`r format(Sys.time(), "%d %B %Y")`'
output:
  word_document:
      reference_docx: RMarkdown_Word_Style_Draft.docx
      fig_caption: yes
  html_document: default
  pdf_document: default
  fig_caption: yes
csl: inter-research-science-center.csl
bibliography: ZoteroLibrary.bib
---

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

#### Keywords:
Eelgrass, pathogen, fatty acids, DHA, *Labyrinthula*

\  

### Abstract
Negative consequences of parasites and disease on hosts are usually better understood than their multifaceted ecosystem effects. The pathogen *Labyrinthula zosterae* (Lz) causes eelgrass wasting disease but has relatives that produce large quantities of nutritionally valuable long-chain polyunsaturated fatty acids (LCPUFA) such as docosahexaenoic acid (DHA). Here we quantify the fatty acids (FA) of Lz cultured on artificial media, eelgrass-based media, and eelgrass segments to investigate whether Lz may similarly produce LCPUFA. We also assess whether field-collected lesions show similar FA patterns to laboratory-inoculated eelgrass. We find that Lz produces DHA as its dominant FA along with other essential FA, even on an eelgrass-based media. DHA content was greater in both laboratory-inoculated and field-collected diseased eelgrass relative to their respective controls. If Lz's production scales *in situ*, it may present an unrecognized source of LCPUFA in eelgrass ecosystems.

\  

### 1. INTRODUCTION

By definition, parasites and disease negatively affect their hosts. However, their effects, either direct or via modification of host traits or abundance, are multifaceted and relatively understudied at community and ecosystem scales [@preston_disease_2016]. For example, parasites are recognized to modify biogeochemical processes [@breitbart_marine_2012] and account for substantial biomass [@kuris_ecosystem_2008]. Thus, it should be a priority to understand potential effects of parasites beyond their hosts.

Eelgrass wasting disease (EWD), afflicting the eelgrass *Zostera marina* L., caused die-offs of Atlantic eelgrass beds in the 1930s [@renn_wasting_1934] and 1980s [@short_eelgrass_1987]. It is caused by the parasitic protist *Labyrinthula zosterae* Muehlstein et al. [-@muehlstein_span_1991, hereafter "Lz"]. EWD and related *Labyrinthula*-caused diseases are concerning with seagrasses declining worldwide [@waycott_accelerating_2009], especially with EWD's wide geographic breadth [@sullivan_potential_2013] and diverse seagrass hosts infected by labyrinthulids [@vergeer_omnipresence_1994]. EWD can also be common within sites; In Washington, USA, Groner et al. [-@groner_plant_2016] found prevalences >40% at over half of their sites, with a 79% EWD prevalence at one site.

Lz belongs to the Labyrinthulomycetes, heterotrophic protists characterized by their ectoplasmic net [@raghukumar_ecology_2002]. The Labyrinthulomycetes include diverse pathogens [@raghukumar_ecology_2002] but are also known for abundantly producing long-chain polyunsaturated fatty acids (LCPUFA), namely eicosapentaenoic acid (EPA, 20:5ω-3) and docosahexaenoic acid (DHA, 22:6ω-3) [@kumon_new_2006;@armenta_single-cell_2013]. Such LCPUFA and some other fatty acids (FA), including linoleic acid (LIN, 18:2ω-6) and α-linolenic acid (ALA, 18:3ω-3), are important nutritional components for many consumers that rely on them for diverse physiological needs ranging from cell membrane fluidity to nervous system function [@parrish_essential_2009]. LIN and ALA, in particular, are considered "essential" FA, because most consumers poorly or cannot synthesize them in biologically relevant amounts [@parrish_essential_2009]. DHA and EPA are also often termed essential for the same reason, although ALA can serve as a precursor for them if the organism in consideration has elongation and desaturation abilities [@parrish_essential_2009]. Eelgrass has ALA as its dominant fatty acid but is poor in LCPUFA [@galloway_fatty_2012;@jaschinski_carbon_2008]. While increasing evidence shows animals are more capable of synthesizing LCPUFA than recognized [@kabeya_genes_2018], increased growth and reproduction associated with greater dietary long-chain essential FA [@winder_landsea_2017] indicate that exogenous sources of LCPUFA are still important. In eelgrass beds, red algae and epiphytic diatoms may supply EPA and some other LCPUFA [@galloway_fatty_2012;@jaschinski_carbon_2008]. Conversely, though, DHA is uncommon in macrophytes [@galloway_fatty_2012], suggesting DHA may instead enter eelgrass beds through the water column or microbial production.

Microbial bioconversion of biomolecules has been suggested as a mechanism for "trophic upgrading" [@klein_breteler_trophic_1999] in aquatic food webs. Specifically, microorganisms consume relatively nutritionally poor organisms or substrates and then synthesize more nutritionally valuable compounds for their predators. The ability for some planktonic heterotrophic protists to produce LCPUFA and other compounts is suggested to be one means of upgrading in pelagic food webs [@chu_quantitative_2008;@chu_species_2009]. Similarly, the capacity of Labyrinthulomycetes, which are typically saprotrophs, to synthesize LCPUFA is hypothesized to be a possible upgrading mechanism for detritus and an important ecological role for this taxon [@raghukumar_ecology_2002]. 

Here we investigate an intersection of Labyrinthulomycetes' pathogenicity and FA production. We quantify the FA of Lz on three substrates in the laboratory: serum seawater agar (SSA, a typical artificial medium), an eelgrass-based medium (biologically relevant substrate), and inoculated eelgrass segments (approximation of *in-situ* production). We use these substrates to determine how *in-vitro* FA production of Lz might vary with different resources and whether in-host Lz FA production might be detectable with potential host response. To assess whether substantial FA production may occur *in situ*, we also quantify FA of field-collected EWD-affected eelgrass. We focus particularly on ALA and DHA, hypothesizing that Lz presence will increase sample DHA at the expense of ALA, a key eelgrass FA and potential percursor for DHA.

```{r, include=FALSE}
# Clear workspace
rm(list=ls())
# Clear console
cat("\014")
# Set working directory appropriately
setwd("C:/Users/reynm/Dropbox/GRAD/01_EelgrassMediaProject")

# packages ----
library(reshape2)
library(vegan)
library(MASS)
library(ggplot2)
library(grid)
library(scales)
library(ggrepel)
library(dunn.test)
library(gridExtra)
library(plyr)
library(dplyr)
library(knitr)
library(flextable)
library(viridis)
library(ggpubr)
library(cowplot)

# 19:0 Standard Concentration calculation
C19_conc_ugperuL=C19_conc_mgpermL=(139.76-128.53)/10

```

```{r univariate, include=FALSE, message=FALSE}

# Read in data
# Each data set below is a sheet/table in Supplement 2, saved separately as a .csv
# note that concentrations in the data are ng per uL
data=read.csv("Yoshioka_etal_EelgrassPathogenFattyAcid_FA_full.csv") #Table S4
disease=read.csv("Yoshioka_etal_EelgrassPathogenFattyAcid_disease.csv") #Table S5
temp=read.csv("Yoshioka_etal_EelgrassPathogenFattyAcid_Temp.csv") #Table S6

# Temperature ----
temp$DT=strptime(temp$DT,format="%m/%d/%Y %H:%M")
expstart=strptime("11/4/2017 22:30",format="%m/%d/%Y %H:%M")
expend=strptime("1/24/2018 12:00",format="%m/%d/%Y %H:%M")
temp.exp=subset(temp,temp$DT>expstart&temp$DT<expend)
temp.mean=mean(temp.exp$TC)
temp.max=max(temp.exp$TC)
temp.min=min(temp.exp$TC)
temp.IQR=IQR(temp.exp$TC)
temp.med=median(temp.exp$TC)
temp.sd=sd(temp.exp$TC)

# Disease ----
# calculate lesion percent coverage based on lesion area and blade area
disease$sev_perc=disease$area_les/disease$area_seg*100
# Kruskal-Wallis test
kruskal.test(disease$sev_perc~disease$isolate)

# Dunn test posthoc with Bonferroni correction
dunn.test(disease$sev_perc, disease$isolate,method="bonferroni")

disease.med=aggregate(disease$sev_perc,by=list(disease$isolate),median)
disease.mean=aggregate(disease$sev_perc,by=list(disease$isolate),mean)

# Univariate fatty acid analyses ----

# "Metadata" length (# columns with descriptive information)
# This isn't actually metadata, but for matrices we only use numeric data
metadata=9

# Cleaning data
# Ensure numeric data are numeric (note that NAs will be introduced to replace "blank")
data$foilonly=as.numeric(as.character(data$foilonly))
data$foilsample=as.numeric(as.character(data$foilsample))

# Remove FA "measured" by GCMSSolutions that do not match MS identifications
data.sub=subset(data,select=-c(C8.0,C14.1w5))

# Set up working data
# note that current measures are in ng per uL, hence 1000 factor below
data.work=data.sub

# Calculate percent recovery
FAME_hexane = 1500 #final volume of FAME in hexane in uL
FAME_chloroform = 1000 #chloroform removed for transesterification in uL
lipid_chloroform = 1500 #total chloroform volume of lipid extracted in uL
ng_in_ug = 1000 #nanograms in a microgram

data.work$percentrecovery=data.work$C19.0/(C19_conc_ugperuL*data.work$C19volume/lipid_chloroform*FAME_chloroform/FAME_hexane*ng_in_ug)*100

# Scale to expected 19:0 standard
data.work$C19scale=(data.work$percentrecovery/100)^-1
# reorganize to bring recovery columns to front
data.work=cbind(data.work[c(1:metadata,length(data.work),(length(data.work)-1),(metadata+1):(length(data.work)-2))]) 
metadata=metadata+2 # added two metadata columns with "C19scale" and "percentrecovery"

# function:
C19scaling=function(x,w){
  x*w
}

# use mapply to apply function
data.work[,(metadata+1):length(data.work)]=mapply(C19scaling,
                                       x=data.work[,(metadata+1):length(data.work)],
                                       w=data.frame(data.work$C19scale)[,rep(1,length(data.work)-metadata)]
                                       )

# Subtract procedural blanks
# Need a cleaner method for this:
# Function to subset batches and subtract respective blanks
blank=function(x,y,z) {
  # x = dataframe, y = batch_extraction, z = metadata columns
cup1=subset(x, x$batch_extraction==y) # subset
cup2=cup1[,(z+1):length(cup1)] # take only data columns
cup3=cup1[cup1$sample_type=="blank",(z+1):length(cup1)] # dataframe of blanks
cup4=cup2-cup3[rep(1,length(cup1$C19.0)),] # subtraction
cup4[cup4<0]=0 # since concentrations cannot be negative, replace negatives with 0
cup1[,(z+1):length(cup1)]=cup4
cup1
}
# Use the function
data.2018.01.30=blank(data.work,"2018.01.30",metadata)
data.2018.02.01=blank(data.work,"2018.02.01",metadata)
data.2018.02.02=blank(data.work,"2018.02.02",metadata)
data.2018.02.04=blank(data.work,"2018.02.04",metadata)
data.2018.02.05=blank(data.work,"2018.02.05",metadata)
data.2019.01.30=blank(data.work,"2019.01.30",metadata)
data.2019.02.24=blank(data.work,"2019.02.24",metadata)

# rbind back into data.work
data.work=rbind(data.2018.01.30,data.2018.02.01,data.2018.02.02,data.2018.02.04,data.2018.02.05,data.2019.01.30,data.2019.02.24)

# Calculate sample weights
data.work$weight=data.work$foilsample-data.work$foilonly

# reorganize to bring "weight" to beginning
data.work=cbind(data.work[c(1:metadata,length(data.work),(metadata+1):(length(data.work)-1))]) 
metadata=metadata+1 # added another metadata column with "weight"

# Calculate FA weight per sample weights (ug per mg)
# function:
conc=function(x,Fh,Fc,lc,w){
  x*Fh/Fc*lc/w/1000
}
# use mapply to apply function:
data.work[,(metadata+1):length(data.work)]=mapply(conc,
                                         x=data.work[,(metadata+1):length(data.work)],
                                        Fh=FAME_hexane,
                                        Fc=FAME_chloroform,
                                        lc=lipid_chloroform,
                                         w=data.frame(data.work$weight)[,rep(1,length(data.work)-metadata)]
                                        )

# Second round of data cleaning
# Remove blanks
data.work=subset(data.work,data.work$sample_type!="blank")
# Drop 19:0 (standard) from further work
data.work=subset(data.work,select=-c(C19.0))

# Calculate proportions
# add FA concentrations for each sample
data.work$sum=apply(data.work[,(metadata+1):length(data.work)],1,sum)
# reorganize to bring sum column to front
data.work=cbind(data.work[c(1:metadata,length(data.work),(metadata+1):(length(data.work)-1))]) 
metadata=metadata+1
data.work.prop=data.work

# function
prop=function(x,y){
  x/y*100
}

# use mapply to apply function
data.work.prop[,(metadata+1):length(data.work.prop)]=mapply(prop,
                                                     x=data.work.prop[,(metadata+1):length(data.work.prop)],
                                                     y=data.frame(data.work.prop$sum)[,rep(1,
                                                       (length(data.work.prop)-metadata))]
                                                     ) 

# Further data cleaning and organization
# Melt data
data.m=melt(data.work,id.vars=c("filename","batch_GCMS","batch_extraction","isolate","substrate","sample_type","foilsample","foilonly","C19volume","C19scale","percentrecovery","weight","sum"),variable.name="FA",value.name="ug_per_mg")
data.prop.m=melt(data.work.prop,id.vars=c("filename","batch_GCMS","batch_extraction","isolate","substrate","sample_type","foilsample","foilonly","C19volume","C19scale","percentrecovery","weight","sum"),variable.name="FA",value.name="percent")

# Remove rare FA
# Remove if not on average at least 1% of at least 1 sample type
# aggregate to find means
mean.FA.prop=aggregate(data.prop.m$percent,by=list(data.prop.m$substrate,data.prop.m$sample_type,data.prop.m$FA),mean)
# rename columns
colnames(mean.FA.prop)=c("substrate","sample_type","FA","mean_percent")
# subset each sample type and determine which FA represent at least 1 %:
# SSA:
SSA_Lz_cells=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="Lz_cells" & mean.FA.prop$mean_percent>=1)
SSA_Lz_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=1)
SSA_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=1)
# EGA:
EGA_Lz_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="EGA" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=1)
EGA_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="EGA" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=1)
# Eelgrass Segments:
eelgrass_Lz_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=1)
eelgrass_substrate=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=1)
eelgrass_healthy=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="healthy" & mean.FA.prop$mean_percent>=1)
eelgrass_diseased=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="diseased" & mean.FA.prop$mean_percent>=1)
# subset each sample type and determine which FA represent at least 5 %:
# SSA:
SSA_Lz_cells5=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="Lz_cells" & mean.FA.prop$mean_percent>=5)
SSA_Lz_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=5)
SSA_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="SSA" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=5)
# EGA:
EGA_Lz_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="EGA" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=5)
EGA_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="EGA" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=5)
# Eelgrass Segments:
eelgrass_Lz_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="Lz_substrate" & mean.FA.prop$mean_percent>=5)
eelgrass_substrate5=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="substrate" & mean.FA.prop$mean_percent>=5)
eelgrass_healthy5=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="healthy" & mean.FA.prop$mean_percent>=5)
eelgrass_diseased5=subset(mean.FA.prop,mean.FA.prop$substrate=="eelgrass" & mean.FA.prop$sample_type=="diseased" & mean.FA.prop$mean_percent>=5)
# Create list of FA with over 1% representation (for further analyses)
FA_list=as.list(c(as.character(SSA_Lz_cells$FA),
                  as.character(SSA_Lz_substrate$FA),
                  as.character(SSA_substrate$FA),
                  as.character(EGA_Lz_substrate$FA),
                  as.character(EGA_substrate$FA),
                  as.character(eelgrass_Lz_substrate$FA),
                  as.character(eelgrass_substrate$FA)))
FA_list=unique(FA_list)

# Create list of FA with over 5% representation (determine arrows for NMDS, see below)
FA_list5=as.list(c(as.character(SSA_Lz_cells5$FA),
                  as.character(SSA_Lz_substrate5$FA),
                  as.character(SSA_substrate5$FA),
                  as.character(EGA_Lz_substrate5$FA),
                  as.character(EGA_substrate5$FA),
                  as.character(eelgrass_Lz_substrate5$FA),
                  as.character(eelgrass_substrate5$FA),
                  as.character(eelgrass_healthy5$FA),
                  as.character(eelgrass_diseased5$FA)))
FA_list5=unique(FA_list5)

# Subset concentration and proportion data to only include prominent FA
data.m.fin=data.m[data.m$FA %in% FA_list,]
data.prop.m.fin=data.prop.m[data.prop.m$FA %in% FA_list,]

# Univariate visualization/analyses (concentrations) ----

# All FA ----
# Just producing table of mean and sd FA , no actual analyses
FAall=aggregate(data.m.fin$ug_per_mg, by=list(data.m.fin$substrate,data.m.fin$sample_type,data.m.fin$FA),mean)
FAall.sd=aggregate(data.m.fin$ug_per_mg, by=list(data.m.fin$substrate,data.m.fin$sample_type,data.m.fin$FA),sd)
FAall=cbind(FAall,FAall.sd$x)
colnames(FAall)=c("substrate","sample_type","FA","mean","sd")

FAall.prop=aggregate(data.prop.m.fin$percent, by=list(data.prop.m.fin$substrate,data.prop.m.fin$sample_type,data.prop.m.fin$FA),mean)
FAall.prop.sd=aggregate(data.prop.m.fin$percent, by=list(data.prop.m.fin$substrate,data.prop.m.fin$sample_type,data.prop.m.fin$FA),sd)
FAall.prop=cbind(FAall.prop,FAall.prop.sd$x)
colnames(FAall.prop)=c("substrate","sample_type","FA","mean","sd")

# Plotting Options ----
panel_color="gray95"
dot_color="gray75"
color_palette=c("#c80000ff","#c88100ff","#69caffff")
fontfamily="sans"
viridisbegin=0.0
viridisend=0.95
pointstroke=1
viridisoption="viridis"

# ALA & DHA Plots ----

# subset data to exclude uninformative (2018) eelgrass segments
data.uni=data.m.fin[data.m.fin$batch_GCMS!="laby3_split5"|data.m.fin$substrate!="eelgrass",]
data.prop.uni=data.prop.m.fin[data.prop.m.fin$batch_GCMS!="laby3_split5"|data.prop.m.fin$substrate!="eelgrass",]

# Subset DHA and ALA
DHA=subset(data.uni,data.uni$FA=="C22.6w3")
ALA=subset(data.uni,data.uni$FA=="C18.3w3")
DHA.prop=subset(data.prop.uni,data.prop.uni$FA=="C22.6w3")
ALA.prop=subset(data.prop.uni,data.prop.uni$FA=="C18.3w3")

# Aggregate
DHA.median=aggregate(DHA$ug_per_mg, by=list(DHA$substrate,DHA$sample_type),median)
DHA.IQR=aggregate(DHA$ug_per_mg, by=list(DHA$substrate,DHA$sample_type),IQR)
DHA.median=cbind(DHA.median,DHA.IQR$x)
colnames(DHA.median)=c("substrate","sample_type","median","IQR")
DHA.median$substrate=factor(DHA.median$substrate,levels=c("blank","SSA","EGA","eelgrass"))

ALA.median=aggregate(ALA$ug_per_mg, by=list(ALA$substrate,ALA$sample_type),median)
ALA.IQR=aggregate(ALA$ug_per_mg, by=list(ALA$substrate,ALA$sample_type),IQR)
ALA.median=cbind(ALA.median,ALA.IQR$x)
colnames(ALA.median)=c("substrate","sample_type","median","IQR")
ALA.median$substrate=factor(ALA.median$substrate,levels=c("blank","SSA","EGA","eelgrass"))

DHA.prop.median=aggregate(DHA.prop$percent, by=list(DHA.prop$substrate,DHA.prop$sample_type),median)
DHA.prop.IQR=aggregate(DHA.prop$percent, by=list(DHA.prop$substrate,DHA.prop$sample_type),IQR)
DHA.prop.median=cbind(DHA.prop.median,DHA.prop.IQR$x)
colnames(DHA.prop.median)=c("substrate","sample_type","median","IQR")
DHA.prop.median$substrate=factor(DHA.prop.median$substrate,levels=c("blank","SSA","EGA","eelgrass"))

ALA.prop.median=aggregate(ALA.prop$percent, by=list(ALA.prop$substrate,ALA.prop$sample_type),median)
ALA.prop.IQR=aggregate(ALA.prop$percent, by=list(ALA.prop$substrate,ALA.prop$sample_type),IQR)
ALA.prop.median=cbind(ALA.prop.median,ALA.prop.IQR$x)
colnames(ALA.prop.median)=c("substrate","sample_type","median","IQR")
ALA.prop.median$substrate=factor(ALA.prop.median$substrate,levels=c("blank","SSA","EGA","eelgrass"))

# Plot ----
# ALA Concentration plot ----
ALA.plot=ggplot(data=ALA.median,aes(x=substrate,y=median,colour=sample_type,fill=sample_type))+
  geom_point(stat="identity",position=position_dodge(0.9),size=0, alpha=0)+
  geom_point(data=ALA,aes(x=substrate, y=ug_per_mg, group=sample_type), position=position_dodge(0.9),shape=16,colour=dot_color)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+ 
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  labs(x=NULL,y=expression("FA"~mu*"g"~"mg"^-1~"DW"))+
  geom_errorbar(aes(ymin=ifelse(median-IQR>0,median-IQR,10^-9),ymax=median+IQR),width=0,position=position_dodge(0.9))+ #set positive(non-0) limit for 0 values to allow plotting
  scale_y_log10(breaks = c(0.001,0.01,0.1,1,10,100),labels=c("0.001","0.01","0.1","1","10","100"))+
  geom_vline(xintercept=c(1.5,2.5), lwd=0.75, colour="white")+
  geom_point(stat="identity",position=position_dodge(0.9),size=3,shape=23,stroke=pointstroke)+
  theme_classic()+
  coord_cartesian(ylim=c(0.0005,200))+
  theme(text=element_text(family=fontfamily),
        axis.text.x = element_text(angle = 0,hjust=0.5,vjust=.5,size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=10),
        legend.title = element_text(size=10),
        legend.text = element_text(size=8),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = "none",
        panel.background = element_rect(fill = panel_color),
        strip.text = element_text(size=12),
        strip.background = element_blank())+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size=3, shape=23, width=0.1, linetype=0)))

# ALA Proportion plot ----
ALA.prop.plot=ggplot(data=ALA.prop.median,aes(x=substrate,y=median,colour=sample_type,fill=sample_type))+
  geom_point(stat="identity",position=position_dodge(0.9),size=0, alpha=0)+
  geom_point(data=ALA.prop,aes(x=substrate, y=percent, group=sample_type), position=position_dodge(0.9),shape=16,colour=dot_color)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+ 
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  labs(x=NULL,y="%TFA")+
  geom_errorbar(aes(ymin=ifelse(median-IQR>0,median-IQR,10^-9),ymax=median+IQR),width=0,position=position_dodge(0.9))+ #set positive(non-0) limit for 0 values to allow plotting
  geom_vline(xintercept=c(1.5,2.5), lwd=0.75, colour="white")+
  geom_point(stat="identity",position=position_dodge(0.9),size=3,shape=23,stroke=pointstroke)+
  theme_classic()+
  coord_cartesian(ylim=c(0,60))+
  theme(text=element_text(family=fontfamily),
        axis.text.x = element_text(angle = 0,hjust=0.5,vjust=.5,size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=10),
        legend.title = element_text(size=10),
        legend.text = element_text(size=8),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = "none",
        panel.background = element_rect(fill = panel_color),
        strip.text = element_text(size=12),
        strip.background = element_blank())+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size=3, shape=23, width=0.1, linetype=0)))

# DHA Concentration plot ----
DHA.plot=ggplot(data=DHA.median,aes(x=substrate,y=median,colour=sample_type,fill=sample_type))+
  geom_point(stat="identity",position=position_dodge(0.9),size=0, alpha=0)+
  geom_point(data=DHA,aes(x=substrate, y=ug_per_mg, group=sample_type), position=position_dodge(0.9),shape=16,colour=dot_color)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+ 
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  labs(x=NULL,y=NULL)+
  geom_errorbar(aes(ymin=ifelse(median-IQR>0,median-IQR,10^-9),ymax=median+IQR),width=0,position=position_dodge(0.9))+ #set positive(non-0) limit for 0 values to allow plotting
  scale_y_log10(breaks = c(0.001,0.01,0.1,1,10,100),labels=c("0.001","0.01","0.1","1","10","100"))+
  geom_vline(xintercept=c(1.5,2.5), lwd=0.75, colour="white")+
  geom_point(stat="identity",position=position_dodge(0.9),size=3,shape=23,stroke=pointstroke)+
  theme_classic()+
  coord_cartesian(ylim=c(0.0005,200))+
  theme(text=element_text(family=fontfamily),
        axis.text.x = element_text(angle = 0,hjust=0.5,vjust=.5,size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=10),
        legend.title = element_text(size=10),
        legend.text = element_text(size=8),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = "none",
        panel.background = element_rect(fill = panel_color),
        strip.text = element_text(size=12),
        strip.background = element_blank())+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size=3, shape=23, width=0.1, linetype=0)))

# DHA Proportion plot ----
DHA.prop.plot=ggplot(data=DHA.prop.median,aes(x=substrate,y=median,colour=sample_type,fill=sample_type))+
  geom_point(stat="identity",position=position_dodge(0.9),size=0, alpha=0)+
  geom_point(data=DHA.prop,aes(x=substrate, y=percent, group=sample_type), position=position_dodge(0.9),shape=16,colour=dot_color)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+ 
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  labs(x=NULL,y=NULL)+
  geom_errorbar(aes(ymin=ifelse(median-IQR>0,median-IQR,10^-9),ymax=median+IQR),width=0,position=position_dodge(0.9))+ #set positive(non-0) limit for 0 values to allow plotting
  geom_vline(xintercept=c(1.5,2.5), lwd=0.75, colour="white")+
  geom_point(stat="identity",position=position_dodge(0.9),size=3,shape=23,stroke=pointstroke)+
  theme_classic()+
  coord_cartesian(ylim=c(0,60))+
  theme(text=element_text(family=fontfamily),
        axis.text.x = element_text(angle = 0,hjust=0.5,vjust=0.5,size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=10),
        legend.title = element_text(size=10),
        legend.text = element_text(size=8),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = "none",
        panel.background = element_rect(fill = panel_color),
        strip.text = element_text(size=12),
        strip.background = element_blank())+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size=3, shape=23, width=0.1, linetype=0)))

# Extract legend

ALA.plot.legend=ggplot(data=ALA.median,aes(x=substrate,y=median,colour=sample_type,fill=sample_type))+
  geom_point(stat="identity",position=position_dodge(0.9),size=0, alpha=0)+
  geom_point(data=ALA,aes(x=substrate, y=ug_per_mg, group=sample_type), position=position_dodge(0.9),shape=16,colour=dot_color)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+ 
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,name = "sample type",labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,name = "sample type",alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  labs(x=NULL,y=NULL)+
  geom_errorbar(aes(ymin=ifelse(median-IQR>0,median-IQR,10^-9),ymax=median+IQR),width=0,position=position_dodge(0.9))+ #set positive(non-0) limit for 0 values to allow plotting
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x))+
  geom_vline(xintercept=c(1.5,2.5), lwd=0.75, colour="white")+
  geom_point(stat="identity",position=position_dodge(0.9),size=3,shape=23,stroke=pointstroke)+
  theme_classic()+
  coord_cartesian(ylim=c(0.00005,200))+
  theme(text=element_text(family=fontfamily),
        axis.text.x = element_text(angle = 90,hjust=1,vjust=.5,size=10),
        axis.text.y = element_text(size=8),
        axis.title = element_text(size=10),
        legend.title = element_text(size=10),
        legend.text = element_text(size=8),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = "bottom",
        panel.background = element_rect(fill = panel_color),
        strip.text = element_text(size=12),
        strip.background = element_blank())+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size=3, shape=23, width=0.1, linetype=0)))
figure1_legend=as_ggplot(get_legend(ALA.plot.legend))

# Put it all together ----
figure1ALA=plot_grid(ALA.prop.plot,ALA.plot,nrow=2,ncol=1,align = "v",labels=c("A","C"),label_size=12,label_x=0.183, label_y =0.98,label_fontface = "plain")
figure1ALA=annotate_figure(figure1ALA,top=text_grob(parse(text="~~~~~~~~~~~~~~~~ALA")))
figure1DHA=plot_grid(DHA.prop.plot,DHA.plot,nrow=2,ncol=1,align = "v", labels=c("B","D"),label_size=12,label_x=0.12,label_y =0.98,label_fontface = "plain")
figure1DHA=annotate_figure(figure1DHA,top=text_grob(parse(text="~~~~~~~~~~DHA")))
figure1=ggarrange(figure1ALA,figure1DHA,nrow=1,ncol=2,widths = c(1.075,1))
figure1=annotate_figure(figure1,bottom=text_grob("substrate",size=10))
figure1=ggarrange(figure1,figure1_legend,nrow=2,ncol=1,heights = c(9,1))
figure1

# Export figure 1 ----
pdf("Yoshioka_etal_EelgrassPathogenFattyAcid_figure_1.pdf",width=6.65,height=5)
figure1
dev.off()

png("Yoshioka_etal_EelgrassPathogenFattyAcid_figure_1.png",width=6.65*300,height=5*300,res=300)
figure1
dev.off()

# Analyze concentrations ----

# Kruskal-Wallis for SSA
KW.18.3w3.SSA=kruskal.test(ug_per_mg~sample_type,data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C18.3w3",])
KW.18.3w3.SSA.tab=data.frame(test="Kruskal-Wallis",comparison="total ALA by sample type",statistic=paste("chi2=",signif(KW.18.3w3.SSA$statistic,4),sep=""),df=KW.18.3w3.SSA$parameter,p=as.numeric(signif(KW.18.3w3.SSA$p.value,4)),substrate="SSA",FA="ALA")

KW.22.6w3.SSA=kruskal.test(ug_per_mg~sample_type,data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C22.6w3",])
KW.22.6w3.SSA.tab=data.frame(test="Kruskal-Wallis",comparison="total DHA by sample type",statistic=paste("chi2=",signif(KW.22.6w3.SSA$statistic,4),sep=""),df=KW.22.6w3.SSA$parameter,p=as.numeric(signif(KW.22.6w3.SSA$p.value,4)),substrate="SSA",FA="DHA")

# Posthoc Dunn test for SSA
dunn.18.3w3.SSA=dunn.test(data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C18.3w3",]$ug_per_mg, data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C18.3w3",]$sample_type,method="bonferroni")
dunn.18.3w3.SSA.tab=cbind(dunn.18.3w3.SSA$comparison,paste("Z=",signif(dunn.18.3w3.SSA$Z,4),sep=""),as.numeric(signif(dunn.18.3w3.SSA$P.adjusted,4)))
colnames(dunn.18.3w3.SSA.tab)=c("comparison","statistic","p")
dunn.18.3w3.SSA.tab=data.frame(dunn.18.3w3.SSA.tab)
dunn.18.3w3.SSA.tab$test="Dunn's"
dunn.18.3w3.SSA.tab$substrate="SSA"
dunn.18.3w3.SSA.tab$FA="ALA"
dunn.18.3w3.SSA.tab$p=as.numeric(as.character(dunn.18.3w3.SSA.tab$p))

dunn.22.6w3.SSA=dunn.test(data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C22.6w3",]$ug_per_mg, data.m.fin[data.m.fin$substrate=="SSA"&data.m.fin$FA=="C22.6w3",]$sample_type,method="bonferroni")
dunn.22.6w3.SSA.tab=cbind(dunn.22.6w3.SSA$comparison,paste("Z=",signif(dunn.22.6w3.SSA$Z,4),sep=""),as.numeric(signif(dunn.22.6w3.SSA$P.adjusted,4)))
colnames(dunn.22.6w3.SSA.tab)=c("comparison","statistic","p")
dunn.22.6w3.SSA.tab=data.frame(dunn.22.6w3.SSA.tab)
dunn.22.6w3.SSA.tab$test="Dunn's"
dunn.22.6w3.SSA.tab$substrate="SSA"
dunn.22.6w3.SSA.tab$FA="DHA"
dunn.22.6w3.SSA.tab$p=as.numeric(as.character(dunn.22.6w3.SSA.tab$p))

# Wilcoxon Signed Rank for EGA
data.m.fin.18.3w3.substrate=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="EGA"&data.m.fin$sample_type=="substrate")
data.m.fin.18.3w3.Lz_substrate=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="EGA"&data.m.fin$sample_type=="Lz_substrate")
MWU.18.3w3.EGA=wilcox.test(data.m.fin.18.3w3.substrate$ug_per_mg,data.m.fin.18.3w3.Lz_substrate$ug_per_mg,paired=T)
MWU.18.3w3.EGA.tab=data.frame(test="Wilcoxon Signed Rank",comparison="Lz-colonized substrate vs substrate alone",statistic=paste("V=",signif(MWU.18.3w3.EGA$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.EGA$p.value,4)),substrate="EGA",FA="ALA")

data.m.fin.22.6w3.substrate=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="EGA"&data.m.fin$sample_type=="substrate")
data.m.fin.22.6w3.Lz_substrate=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="EGA"&data.m.fin$sample_type=="Lz_substrate")
MWU.22.6w3.EGA=wilcox.test(data.m.fin.22.6w3.substrate$ug_per_mg,data.m.fin.22.6w3.Lz_substrate$ug_per_mg,paired=T)
MWU.22.6w3.EGA.tab=data.frame(test="Wilcoxon Signed Rank",comparison="Lz-colonized substrate vs substrate alone",statistic=paste("V=",signif(MWU.22.6w3.EGA$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.EGA$p.value,4)),substrate="EGA",FA="DHA")

# Wilcoxon Signed Ranks for Eelgrass Segments
# Field samples
data.m.fin.18.3w3.healthy=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="healthy")
data.m.fin.18.3w3.diseased=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="diseased")
MWU.18.3w3.eelgrass.field=wilcox.test(data.m.fin.18.3w3.healthy$ug_per_mg,data.m.fin.18.3w3.diseased$ug_per_mg,paired=T,alternative="greater")
MWU.18.3w3.eelgrass.field.tab=data.frame(test="Wilcoxon Signed Rank",comparison="healthy > diseased",statistic=paste("V=",signif(MWU.18.3w3.eelgrass.field$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.eelgrass.field$p.value,4)),substrate="eelgrass (field)",FA="ALA")

data.m.fin.22.6w3.healthy=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="healthy")
data.m.fin.22.6w3.diseased=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="diseased")
MWU.22.6w3.eelgrass.field=wilcox.test(data.m.fin.22.6w3.healthy$ug_per_mg,data.m.fin.22.6w3.diseased$ug_per_mg,paired=T,alternative="less")
MWU.22.6w3.eelgrass.field.tab=data.frame(test="Wilcoxon Signed Rank",comparison="healthy < diseased",statistic=paste("V=",signif(MWU.22.6w3.eelgrass.field$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.eelgrass.field$p.value,4)),substrate="eelgrass (field)",FA="DHA")

# Inoculated (2019 only)
data.m.fin.18.3w3.substrate=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="substrate"&data.m.fin$batch_GCMS=="laby4")
data.m.fin.18.3w3.Lz_substrate=subset(data.m.fin,data.m.fin$FA=="C18.3w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="Lz_substrate"&data.m.fin$batch_GCMS=="laby4")
MWU.18.3w3.eelgrass.inoc=wilcox.test(data.m.fin.18.3w3.substrate$ug_per_mg,data.m.fin.18.3w3.Lz_substrate$ug_per_mg,paired=T,alternative="greater")
MWU.18.3w3.eelgrass.inoc.tab=data.frame(test="Wilcoxon Signed Rank",comparison="control > inoculated",statistic=paste("V=",signif(MWU.18.3w3.eelgrass.inoc$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.eelgrass.inoc$p.value,4)),substrate="eelgrass (lab)",FA="ALA")

data.m.fin.22.6w3.substrate=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="substrate"&data.m.fin$batch_GCMS=="laby4")
data.m.fin.22.6w3.Lz_substrate=subset(data.m.fin,data.m.fin$FA=="C22.6w3"&data.m.fin$substrate=="eelgrass"&data.m.fin$sample_type=="Lz_substrate"&data.m.fin$batch_GCMS=="laby4")
MWU.22.6w3.eelgrass.inoc=wilcox.test(data.m.fin.22.6w3.substrate$ug_per_mg,data.m.fin.22.6w3.Lz_substrate$ug_per_mg,paired=T,alternative="less")
MWU.22.6w3.eelgrass.inoc.tab=data.frame(test="Wilcoxon Signed Rank",comparison="control < inoculated",statistic=paste("V=",signif(MWU.22.6w3.eelgrass.inoc$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.eelgrass.inoc$p.value,4)),substrate="eelgrass (lab)",FA="DHA")

# Combine test results ----
FAconc.tab=rbind.fill(KW.18.3w3.SSA.tab,dunn.18.3w3.SSA.tab,MWU.18.3w3.EGA.tab,MWU.18.3w3.eelgrass.inoc.tab,MWU.18.3w3.eelgrass.field.tab,KW.22.6w3.SSA.tab,dunn.22.6w3.SSA.tab,MWU.22.6w3.EGA.tab,MWU.22.6w3.eelgrass.inoc.tab,MWU.22.6w3.eelgrass.field.tab)

FAconc.tab$comparison=gsub(" - "," vs ",FAconc.tab$comparison)
FAconc.tab$comparison=gsub("vs substrate","vs substrate alone ",FAconc.tab$comparison)
FAconc.tab$comparison=gsub("Lz_substrate","Lz-colonized substrate",FAconc.tab$comparison)
FAconc.tab$comparison=gsub("Lz_cells","Lz cells",FAconc.tab$comparison)

# Analyze proportions ----

# Kruskal-Wallis for SSA
KW.18.3w3.prop.SSA=kruskal.test(percent~sample_type,data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C18.3w3",])
KW.18.3w3.prop.SSA.tab=data.frame(test="Kruskal-Wallis",comparison="percent ALA by sample type",statistic=paste("chi2=",signif(KW.18.3w3.prop.SSA$statistic,4),sep=""),df=KW.18.3w3.prop.SSA$parameter,p=as.numeric(signif(KW.18.3w3.prop.SSA$p.value,4)),substrate="SSA",FA="ALA")

KW.22.6w3.prop.SSA=kruskal.test(percent~sample_type,data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",])
KW.22.6w3.prop.SSA.tab=data.frame(test="Kruskal-Wallis",comparison="percent DHA by sample type",statistic=paste("chi2=",signif(KW.22.6w3.prop.SSA$statistic,4),sep=""),df=KW.22.6w3.prop.SSA$parameter,p=as.numeric(signif(KW.22.6w3.prop.SSA$p.value,4)),substrate="SSA",FA="DHA")

# Posthoc Dunn test for SSA
dunn.18.3w3.prop.SSA=dunn.test(data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C18.3w3",]$percent, data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C18.3w3",]$sample_type,method="bonferroni")
dunn.18.3w3.prop.SSA.tab=cbind(dunn.18.3w3.prop.SSA$comparison,paste("Z=",signif(dunn.18.3w3.prop.SSA$Z,3),sep=""),as.numeric(signif(dunn.18.3w3.prop.SSA$P.adjusted,3)))
colnames(dunn.18.3w3.prop.SSA.tab)=c("comparison","statistic","p")
dunn.18.3w3.prop.SSA.tab=data.frame(dunn.18.3w3.prop.SSA.tab)
dunn.18.3w3.prop.SSA.tab$test="Dunn's"
dunn.18.3w3.prop.SSA.tab$substrate="SSA"
dunn.18.3w3.prop.SSA.tab$FA="ALA"
dunn.18.3w3.prop.SSA.tab$p=as.numeric(as.character(dunn.18.3w3.prop.SSA.tab$p))

dunn.22.6w3.prop.SSA=dunn.test(data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",]$percent, data.prop.m.fin[data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",]$sample_type,method="bonferroni")
dunn.22.6w3.prop.SSA.tab=cbind(dunn.22.6w3.prop.SSA$comparison,paste("Z=",signif(dunn.22.6w3.prop.SSA$Z,3),sep=""),as.numeric(signif(dunn.22.6w3.prop.SSA$P.adjusted,3)))
colnames(dunn.22.6w3.prop.SSA.tab)=c("comparison","statistic","p")
dunn.22.6w3.prop.SSA.tab=data.frame(dunn.22.6w3.prop.SSA.tab)
dunn.22.6w3.prop.SSA.tab$test="Dunn's"
dunn.22.6w3.prop.SSA.tab$substrate="SSA"
dunn.22.6w3.prop.SSA.tab$FA="DHA"
dunn.22.6w3.prop.SSA.tab$p=as.numeric(as.character(dunn.22.6w3.prop.SSA.tab$p))

# Wilcoxon Signed Rank for EGA
data.prop.m.fin.18.3w3.prop.substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="EGA"&data.prop.m.fin$sample_type=="substrate")
data.prop.m.fin.18.3w3.prop.Lz_substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="EGA"&data.prop.m.fin$sample_type=="Lz_substrate")
MWU.18.3w3.prop.EGA=wilcox.test(data.prop.m.fin.18.3w3.prop.substrate$percent,data.prop.m.fin.18.3w3.prop.Lz_substrate$percent,paired=T)
MWU.18.3w3.prop.EGA.tab=data.frame(test="Wilcoxon Signed Rank",comparison="Lz-colonized substrate vs substrate alone",statistic=paste("V=",signif(MWU.18.3w3.prop.EGA$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.prop.EGA$p.value,4)),substrate="EGA",FA="ALA")

data.prop.m.fin.22.6w3.prop.substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="EGA"&data.prop.m.fin$sample_type=="substrate")
data.prop.m.fin.22.6w3.prop.Lz_substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="EGA"&data.prop.m.fin$sample_type=="Lz_substrate")
MWU.22.6w3.prop.EGA=wilcox.test(data.prop.m.fin.22.6w3.prop.substrate$percent,data.prop.m.fin.22.6w3.prop.Lz_substrate$percent,paired=T)
MWU.22.6w3.prop.EGA.tab=data.frame(test="Wilcoxon Signed Rank",comparison="Lz-colonized substrate vs substrate alone",statistic=paste("V=",signif(MWU.22.6w3.prop.EGA$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.prop.EGA$p.value,4)),substrate="EGA",FA="DHA")

# Wilcoxon Signed Ranks for Eelgrass Segments
# Field samples
data.prop.m.fin.18.3w3.prop.healthy=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="healthy")
data.prop.m.fin.18.3w3.prop.diseased=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="diseased")
MWU.18.3w3.prop.eelgrass.field=wilcox.test(data.prop.m.fin.18.3w3.prop.healthy$percent,data.prop.m.fin.18.3w3.prop.diseased$percent,paired=T,alternative="greater")
MWU.18.3w3.prop.eelgrass.field.tab=data.frame(test="Wilcoxon Signed Rank",comparison="healthy >  diseased",statistic=paste("V=",signif(MWU.18.3w3.prop.eelgrass.field$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.prop.eelgrass.field$p.value,4)),substrate="eelgrass (field)",FA="ALA")

data.prop.m.fin.22.6w3.prop.healthy=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="healthy")
data.prop.m.fin.22.6w3.prop.diseased=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="diseased")
MWU.22.6w3.prop.eelgrass.field=wilcox.test(data.prop.m.fin.22.6w3.prop.healthy$percent,data.prop.m.fin.22.6w3.prop.diseased$percent,paired=T,alternative="less")
MWU.22.6w3.prop.eelgrass.field.tab=data.frame(test="Wilcoxon Signed Rank",comparison="healthy < diseased",statistic=paste("V=",signif(MWU.22.6w3.prop.eelgrass.field$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.prop.eelgrass.field$p.value,4)),substrate="eelgrass (field)",FA="DHA")

# Inoculated (2019 only)
data.prop.m.fin.18.3w3.prop.substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="substrate"&data.prop.m.fin$batch_GCMS=="laby4")
data.prop.m.fin.18.3w3.prop.Lz_substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C18.3w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="Lz_substrate"&data.prop.m.fin$batch_GCMS=="laby4")
MWU.18.3w3.prop.eelgrass.inoc=wilcox.test(data.prop.m.fin.18.3w3.prop.substrate$percent,data.prop.m.fin.18.3w3.prop.Lz_substrate$percent,paired=T,alternative="greater")
MWU.18.3w3.prop.eelgrass.inoc.tab=data.frame(test="Wilcoxon Signed Rank",comparison="control > inoculated",statistic=paste("V=",signif(MWU.18.3w3.prop.eelgrass.inoc$statistic,4),sep=""),p=as.numeric(signif(MWU.18.3w3.prop.eelgrass.inoc$p.value,4)),substrate="eelgrass (lab)",FA="ALA")

data.prop.m.fin.22.6w3.prop.substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="substrate"&data.prop.m.fin$batch_GCMS=="laby4")
data.prop.m.fin.22.6w3.prop.Lz_substrate=subset(data.prop.m.fin,data.prop.m.fin$FA=="C22.6w3"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$sample_type=="Lz_substrate"&data.prop.m.fin$batch_GCMS=="laby4")
MWU.22.6w3.prop.eelgrass.inoc=wilcox.test(data.prop.m.fin.22.6w3.prop.substrate$percent,data.prop.m.fin.22.6w3.prop.Lz_substrate$percent,paired=T,alternative="less")
MWU.22.6w3.prop.eelgrass.inoc.tab=data.frame(test="Wilcoxon Signed Rank",comparison="control < inoculated",statistic=paste("V=",signif(MWU.22.6w3.prop.eelgrass.inoc$statistic,4),sep=""),p=as.numeric(signif(MWU.22.6w3.prop.eelgrass.inoc$p.value,4)),substrate="eelgrass (lab)",FA="DHA")

# Combine test results ----
FAprop.tab=rbind.fill(KW.18.3w3.prop.SSA.tab,dunn.18.3w3.prop.SSA.tab,MWU.18.3w3.prop.EGA.tab,MWU.18.3w3.prop.eelgrass.inoc.tab,MWU.18.3w3.prop.eelgrass.field.tab,KW.22.6w3.prop.SSA.tab,dunn.22.6w3.prop.SSA.tab,MWU.22.6w3.prop.EGA.tab,MWU.22.6w3.prop.eelgrass.inoc.tab,MWU.22.6w3.prop.eelgrass.field.tab)

FAprop.tab$comparison=gsub(" - "," vs ",FAprop.tab$comparison)
FAprop.tab$comparison=gsub("vs substrate","vs substrate alone ",FAprop.tab$comparison)
FAprop.tab$comparison=gsub("Lz_substrate","Lz-colonized substrate",FAprop.tab$comparison)
FAprop.tab$comparison=gsub("Lz_cells","Lz cells",FAprop.tab$comparison)

```

```{r flextable Build Table S3, echo=FALSE}
# decimal counter function
# Will count decimals except trailing 0s or sigfigs of integers
deccount=function(x){
  nchar(gsub("(.*\\.)|([0]*$)", "", as.character(x))) 
}

usedFA.median=aggregate(data.uni$ug_per_mg, by=list(data.uni$substrate,data.uni$sample_type,data.uni$FA),median)
usedFA.IQR=aggregate(data.uni$ug_per_mg, by=list(data.uni$substrate,data.uni$sample_type,data.uni$FA),IQR)
usedFA.median=cbind(usedFA.median,usedFA.IQR$x)
colnames(usedFA.median)=c("substrate","sample_type","FA","median","IQR")

usedFA.prop.median=aggregate(data.prop.uni$percent, by=list(data.prop.uni$substrate,data.prop.uni$sample_type,data.prop.uni$FA),median)
usedFA.prop.IQR=aggregate(data.prop.uni$percent, by=list(data.prop.uni$substrate,data.prop.uni$sample_type,data.prop.uni$FA),IQR)
usedFA.prop.median=cbind(usedFA.prop.median,usedFA.prop.IQR$x)
colnames(usedFA.prop.median)=c("substrate","sample_type","FA","median","IQR")

usedFA.median$medIQR=paste(ifelse((signif(usedFA.median$IQR,2)-floor(signif(usedFA.median$IQR,2)))!=0,round(usedFA.median$median,deccount(signif(usedFA.median$IQR,2))),round(usedFA.median$median,0)),"±",signif(usedFA.median$IQR,2),sep="")

usedFA.prop.median$medIQR=paste(ifelse((signif(usedFA.prop.median$IQR,2)-floor(signif(usedFA.prop.median$IQR,2)))!=0,round(usedFA.prop.median$median,deccount(signif(usedFA.prop.median$IQR,2))),round(usedFA.prop.median$median,0)),"±",signif(usedFA.prop.median$IQR,2),sep="")

usedFA.median$medIQR=paste(usedFA.median$medIQ," (",usedFA.prop.median$medIQ,"%)",sep="")

# Combine it all
usedFA.median.cast=dcast(usedFA.median,FA~substrate+sample_type,value.var="medIQR")

# Reorder columns: SSA then EGA then eelgrass
usedFA.median.cast=usedFA.median.cast[,c("FA",
                                         "SSA_Lz_cells",
                                         "SSA_Lz_substrate",
                                         "SSA_substrate",
                                         "EGA_Lz_substrate",
                                         "EGA_substrate",
                                         "eelgrass_diseased",
                                         "eelgrass_healthy",
                                         "eelgrass_Lz_substrate",
                                         "eelgrass_substrate")]

FA.medIQR.tab=usedFA.median.cast

FA.medIQR.tab$FA=sub("C","",FA.medIQR.tab$FA)
FA.medIQR.tab$FA=sub("c","",FA.medIQR.tab$FA)
FA.medIQR.tab$FA=sub("\\.",":",FA.medIQR.tab$FA)
FA.medIQR.tab$FA=sub("w0","n-x",FA.medIQR.tab$FA)
FA.medIQR.tab$FA=sub("w","n-",FA.medIQR.tab$FA)

FA.medIQR.flex.header=data.frame(col_keys=c("FA",
                                         "SSA_Lz_cells",
                                         "SSA_Lz_substrate",
                                         "SSA_substrate",
                                         "EGA_Lz_substrate",
                                         "EGA_substrate",
                                         "eelgrass_diseased",
                                         "eelgrass_healthy",
                                         "eelgrass_Lz_substrate",
                                         "eelgrass_substrate"), medium=c("FA","SSA","SSA","SSA","EGA","EGA","eelgrass","eelgrass","eelgrass","eelgrass"),
                                 treatment=c("FA","Lz cells (n=8)","Lz-colonized substrate (n=8)","susbtrate alone (n=4)","Lz-colonized substrate (n=8)","substrate alone (n=8)","field diseased (n=5)","field healthy (n=5)","Lz-colonized substrate (n=5)","substrate alone (n=5)"),
                                 stringsAsFactors = FALSE)

FA.medIQR.flex=regulartable(FA.medIQR.tab,col_keys=c("FA","break1","SSA_Lz_cells","SSA_Lz_substrate","SSA_substrate","break2","EGA_Lz_substrate","EGA_substrate","break3","eelgrass_diseased","eelgrass_healthy","eelgrass_Lz_substrate","eelgrass_substrate"))
FA.medIQR.flex=set_header_df(FA.medIQR.flex,mapping=FA.medIQR.flex.header,key="col_keys")
FA.medIQR.flex=merge_h(FA.medIQR.flex, part = "header")
FA.medIQR.flex=merge_v(FA.medIQR.flex, part = "header")
FA.medIQR.flex=fontsize(FA.medIQR.flex,size=5,part="body")
FA.medIQR.flex=fontsize(FA.medIQR.flex,size=5,part="header")
FA.medIQR.flex=flextable::font(FA.medIQR.flex,fontname="Times New Roman",part="all")
FA.medIQR.flex=bold(FA.medIQR.flex, part="header")
FA.medIQR.flex=theme_vanilla(FA.medIQR.flex)
FA.medIQR.flex=align(FA.medIQR.flex,align="center",part="all")
FA.medIQR.flex=autofit(FA.medIQR.flex)
FA.medIQR.flex=empty_blanks(FA.medIQR.flex)
```

```{r multivariate, include=FALSE, message=FALSE} 
# Multivariate visualization/analyses (proportions)

# Cast data into wide format
data.fin=dcast(data.m.fin,substrate+sample_type+filename+isolate~FA,value.var="ug_per_mg")
data.prop.fin=dcast(data.prop.m.fin,substrate+sample_type+filename+isolate+batch_GCMS~FA,value.var="percent")

# subset data to exclude uninformative (2018) eelgrass segments
data.multi=data.prop.fin[data.prop.fin$batch_GCMS!="laby3_split5"|data.prop.fin$substrate!="eelgrass",]

# Set up matrices for analyses
data.mat=data.multi[,6:(length(data.multi))] # proportions

# arcsin square root transformation
data.mat=data.mat/100 #turn into proportions
data.mat=asin(sqrt(data.mat))

# set dataframe as matrix
data.mat=as.matrix(data.mat)

# Fit an NMDS! ----
FA_NMDS=metaMDS(data.mat,k=2,distance="euclidean",wascores=F,noshare=F,autotransform=F)
NMDS_stress=FA_NMDS$stress

# View Stressplot
stressplot(FA_NMDS)

# Extract scores from fit
FA_NMDS_plot=scores(FA_NMDS,display=c("sites")) # pull scores from MDS fit 
FA_NMDS_plot=data.frame(data.multi[1:5],FA_NMDS_plot)

# get arrows of FAs correlating to first two dimensions
FA_NMDS_envfit=envfit(FA_NMDS, data.mat, perm = 9999) # fit vector arrows of FAs correlating with dimensions

# pull scores vector arrows of FAs correlating with dimensions
FA_arrows=scores(FA_NMDS_envfit,display=c("vectors")) 

# calculate scaling multiplier based on min range of points in NMDS
scalefactor=0.3*min(max(FA_NMDS_plot$NMDS1)-min(FA_NMDS_plot$NMDS1), # makes scaling multiplier based on min range of points
                    max(FA_NMDS_plot$NMDS2)-min(FA_NMDS_plot$NMDS2))

#FA_arrows=FA_arrows*ordiArrowMul(FA_NMDS_plot) # scaling of arrows to fit plot area (with vegan's ordiArrowMul function)

FA_arrows=FA_arrows*scalefactor # scaling of arrows to fit plot area based on scalefactor above

FA_arrows=as.data.frame(FA_arrows) # save those as a datafullframe

# Use onlyFA with over 5% representation for arrows:
FA_arrows=FA_arrows[rownames(FA_arrows) %in% FA_list5,]

# FA arrow text formatting
arrow_angle=90-((atan2(FA_arrows$NMDS1, FA_arrows$NMDS2)* 180)/pi)
arrow_angle=ifelse(arrow_angle>90 & arrow_angle<270,arrow_angle+180,arrow_angle)
arrow_hjust=ifelse(arrow_angle>90,1,0)
FA_arrows$labs=as.character(rownames(FA_arrows))
FA_arrows$labs=gsub("C","",FA_arrows$labs)
FA_arrows$labs=gsub("B","",FA_arrows$labs)
FA_arrows$labs=sub("c","",FA_arrows$labs)
FA_arrows$labs=sub("\\.",":",FA_arrows$labs)
FA_arrows$labs=sub("w","*omega*-",FA_arrows$labs)

# Replace shorthand with "ALA" and "DHA" for labels
FA_arrows[FA_arrows$labs=="18:3*omega*-3",]$labs="ALA"
FA_arrows[FA_arrows$labs=="22:6*omega*-3",]$labs="DHA"

# Stress annotation
stressanot=grobTree(textGrob(paste("stress = ",round(NMDS_stress,3),sep=''),x=0.03,y=0.95,
                    hjust=0, gp=gpar(fontsize=10,fontfamily=fontfamily)))

# Re-order levels to make ordering consistent
FA_NMDS_plot$substrate=factor(FA_NMDS_plot$substrate,levels=c("blank","SSA","EGA","eelgrass"))

# Plot NMDS ----
Lz_NMDS=ggplot()+
   geom_point(data=FA_NMDS_plot,aes(NMDS1,NMDS2,colour=sample_type,shape=substrate,fill=sample_type),size=3,stroke=pointstroke)+
  geom_segment(data=FA_arrows, aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
               arrow = arrow(length = unit(0.2, "cm")),colour="gray50",alpha=1)+
  #scale_colour_manual(labels = c("Lz cells","Lz-colonized substrate","substrate alone"), values = color_palette)+
  scale_color_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,name = "sample type",labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  scale_fill_viridis(discrete=TRUE,begin=viridisbegin,end=viridisend,option=viridisoption,name = "sample type",alpha=0.5,labels = c("field diseased","field healthy","Lz cells","Lz-colonized substrate","substrate alone")) +  
  scale_shape_manual(labels = c("SSA","EGA","eelgrass"), values = c(21,23,24)) +
  geom_text(data=as.data.frame(FA_arrows[,c("NMDS1","NMDS2")]*1.1),
            aes(NMDS1, NMDS2, label =  FA_arrows$labs),alpha=0.5,size=3,parse=T,
            angle=arrow_angle, hjust=arrow_hjust,family=fontfamily) +
  theme_classic()+
  coord_fixed()+
  xlab("NMDS dimension 1")+
  ylab("NMDS dimension 2")+
  annotation_custom(stressanot)+
  theme(text = element_text(family=fontfamily),
        axis.ticks = element_blank(),
        axis.text = element_blank(),
        axis.title = element_text(size=10),
        legend.box = 'vertical',
        legend.box.just = 'left',
        legend.justification = c("left","top"),
        legend.margin = margin(0, 0, 0, 0),
        legend.position = c(0.03,0.9),
        #legend.position = "right",
        legend.title = element_text(size=10),
        legend.text = element_text(size=8,hjust=0))+
  guides(shape = guide_legend(override.aes = list(size = 3)),
         colour = guide_legend(override.aes = list(size = 3, shape=22)))

# Export figure 2 ----
pdf("Yoshioka_etal_EelgrassPathogenFattyAcid_figure_2.pdf",width=6.65,height=5)
Lz_NMDS
dev.off()

png("Yoshioka_etal_EelgrassPathogenFattyAcid_figure_2.png",width=6.65*300,height=5*300,res=300)
Lz_NMDS
dev.off()

# Set up for multivariate analysis
data.mat.head=data.multi[,1:5] # get just metadata
data.mat.head$combo=paste(data.mat.head$substrate,data.mat.head$sample_type,sep="_")

# PERMDISP2 with betadisper
data.dist=vegdist(data.mat,method="euclidean")
anova(betadisper(data.dist,factor(data.mat.head$combo)))

# PERMANOVA with adonis ----
mod.permanova=adonis(data.mat ~ sample_type*substrate, data=data.mat.head, method="euclidean", permutations=9999)

permanova_sample_type_df=mod.permanova$aov.tab$Df[1]
permanova_substrate_df=mod.permanova$aov.tab$Df[2]

permanova_sample_type_F=mod.permanova$aov.tab$F.Model[1]
permanova_substrate_F=mod.permanova$aov.tab$F.Model[2]

permanova_sample_type_R2=mod.permanova$aov.tab$R2[1]
permanova_substrate_R2=mod.permanova$aov.tab$R2[2]

```

\  

### 2. MATERIALS AND METHODS

We collected EWD lesions for Lz isolation and eelgrass for experimental substrates in the South Slough in Charleston, Oregon, USA in 2017 and 2019. We collected only eelgrass leaves (i.e. not whole turions) under State of Oregon Parks and Recreation Department permit #008-16 and with permission of the South Slough National Estuarine Research Reserve. We performed laboratory work and FA analyses at the University of Oregon Institute of Marine Biology, Charleston, Oregon, USA.

\  

#### 2.1. Substrate and sample production
We prepared three types of substrates in 100-mm petri dishes: SSA, eelgrass agar (EGA), and eelgrass leaf segments (hereafter, "segments"), summarized in table 1 and detailed in Supplement 1. We used nine isolates (V17-V25) of Lz, cultured from EWD lesions. We aimed to capture a greater variety of Lz FA production using multiple isolates rather than a single isolate but did not sample to specifically test isolate differences. We inoculated four replicate plates of each substrate with 10 uL of ~10^6^ cells mL^-1^ of each Lz isolate in sterile seawater. Four replicate control plates for each substrate were treated the same but using sterile seawater as a sham inoculum. EGA did not initially yield sufficient Lz growth, so we plated pieces of Lz-colonized EGA onto new EGA plates. Isolate V22 still did not grow sufficiently and was omitted. We incubated all plates at `r signif(mean(temp.exp$TC),3)` ± `r signif(sd(temp.exp$TC),2)`°C (mean ± sd) with 12/12-hr light/dark fluorescent lighting. We allowed dark lesions typical of EWD to develop on inoculated segments before collection. Incubation durations varied with extremely different timings of disease/degradation on segments and substantial Lz growth on agar substrates (10 days for SSA, 42 days EGA, and 8-31 days segments). 

From SSA plates, we collected Lz-colonized agar from one inoculated plate per isolate (Lz-colonized substrate, n=8) and agar from all control plates (substrate alone, n=4). We also scraped cells from all inoculated SSA plates and pooled them by isolate (Lz cells, n=8) to quantify FA per Lz weight and separate Lz FA from that in associated agar.  We collected paired colonized (Lz-colonized substrate, n=8) and uncolonized (substrate alone, n=8) agar from EGA plates due to described alternative culturing. Before collecting inoculated (Lz-colonized substrate) and control (substrate alone) segments, we photographed and quantified each segment's EWD severity as the percent segment area lesioned using ImageJ [@rasband_imagej_1997]. We stored segments separately but pooled them for FA analyses to ensure necessary sample masses. Because pooling reduced control sample size to n=1, we later produced a set of inoculated and control segments (n=5 each) using larger paired segments (each pair from the same leaf) and only the high-virulence isolate V20. Methods were otherwise the same. We stored and analyzed each of the later segments separately.

To assess whether FA of inoculated segments reflect those of EWD lesions *in situ*, we additionally collected five eelgrass leaves with eelgrass lesions. We excised apparently diseased (lesioned) tissue from healthy tissue, storing each separately (field diseased and field healthy, resp., n=5 each). We verified visual EWD diagnosis by plating small subsamples of diseased tissue on SSA.

We freeze-dried and weighed all samples, obtaining dry sample weights, prior to FA extraction. We extracted and analyzed FA using gas chromatography and mass spectrometry, modifying methods of Taipale et al. [-@taipale_suitability_2016;-@taipale_fatty_2013] as detailed in Supplement 1.

\  

#### 2.2. Data analysis
We analyzed data in R (version:`r paste(R.version$major,R.version$minor,sep=".")`). We standardized FA concentrations by sample dry weight (μg mg^-1^ DW) and calculated FA proportions on the total identified-FA concentrations within a sample (%TFA). We used nonparametric methods for univariate comparisons, as data were non-normal with small sample sizes.  We compared DHA and ALA concentrations and proportions among sample types within substrates. We used Wilcoxon Signed Rank ('Wilcoxon' hereafter) tests for EGA and eelgrass, as samples were paired, and Kruskal-Wallis ('KW' hereafter) tests for SSA. To address the specific hypothesis that Lz would increase DHA at the expense of ALA, we used one-tailed Wilcoxon tests for eelgrass samples but also present results two-tailed tests for transparency. We did not adjust these tests for multiple comparisons because each deliberately investigated a particular FA's patterns within a substrate. We used Dunn's tests with Bonferroni corrections for posthoc pairwise comparisons following KW tests [dunn.test in package dunn.test, @dinno_dunn.test_2017]. 

We used Euclidean distances of arcine-square-root-transformed proportions for multivariate analyses, including FA >1% on average of FA of any substrate-sample type combination (`r length(FA_list)` FA used). We visualized FA multivariate compositions using non-metric multidimensional scaling [NMDS, metaMDS in vegan package, @oksanen_vegan_2018] and compared them among substrates and sample types using permutational multivariate analysis of variance [PERMANOVA with 9999 permutations, adonis in vegan package, @oksanen_vegan_2018].

We do not display or analyze the earlier segment data as the control had n=1 due to pooling. Further, temporal and methodological differences between earlier and later segment sets discourage combining them, although they are qualitatively similar (Figure S1 in Supplement 1).

\  

### 3. RESULTS

All field-collected EWD samples produced Lz growth on SSA plates.

DHA was the dominant FA in SSA- and EGA-cultured Lz samples, followed by palmitic acid (16:0). In SSA-cultured Lz cells, DHA composed `r signif(mean(data.prop.m.fin[data.prop.m.fin$sample_type=="Lz_cells"&data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",]$percent),3)` ± `r signif(sd(data.prop.m.fin[data.prop.m.fin$sample_type=="Lz_cells"&data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",]$percent),3)`% (mean ± sd, n = `r length(data.prop.m.fin[data.prop.m.fin$sample_type=="Lz_cells"&data.prop.m.fin$substrate=="SSA"&data.prop.m.fin$FA=="C22.6w3",]$percent)`) of identified FA. Palmitic acid was the dominant FA in SSA and EGA substrate-alone samples. ALA was the dominant FA of all eelgrass samples.

On SSA, Lz cells had proportionally greater ALA than the substrate alone (Dunn's test, p<0.05, Fig. 1A). In contrast, ALA concentrations were greater in Lz cells than both Lz-colonized SSA or SSA alone (Dunn's test, p<0.05, Fig. 1C). For EGA, Lz-colonized substrate had greater ALA proportions and concentrations than the respective controls (Wilcoxon, p<0.05, Fig. 1A,C). In both laboratory-inoculated and field-collected eelgrass, ALA proportions and concentrations were lower in Lz-inoculated and diseased ones relative to respective controls (one-tailed Wilcoxon, p<0.05; two-tailed p=0.0625, Fig. 1A,C).

DHA proportions were not different between Lz cells or Lz-colonized substrate on SSA (Dunn's test, p=1.0), but both were greater than the SSA alone (Dunn's test, p<0.05, Fig. 1B). DHA concentrations were greater in Lz cells than both Lz-colonized SSA and SSA alone (Dunn's test, p<0.05, Fig. 1D). As with ALA, DHA proportions and concentrations were greater Lz-colonized EGA than EGA alone (Dunn's test, p<0.05, Fig. 1B,D). Eelgrass DHA showed an opposite pattern to ALA; Lz-inoculated segments and diseased tissue both had greater DHA proportions and concentrations than respective controls (one-tailed Wilcoxon, p<0.05; two-tailed p=0.0625, Fig. 1B,D).

We summarize univariate comparisons of proportions and concentrations in Tables S1 and S2, respectively, in Supplement 2.

The FA compositions of samples differed by substrate and sample type (PERMANOVA, substrate: F(`r permanova_substrate_df`) = `r signif(permanova_substrate_F,3)`, partial R^2^ = `r signif(permanova_substrate_R2,3)`; sample type: F(`r permanova_sample_type_df`) = `r signif(permanova_sample_type_F,3)`, partial R^2^ = `r signif(permanova_sample_type_R2,3)`, Fig. 2).

\  

### 4. DISCUSSION

Using both artificial (SSA) and eelgras-derived (EGA) substrates, we show that pathogenic Lz produces DHA as its dominant FA. Notably, as the agar substrates provided here lack significant quantities of even precursor FA, our results indicate Lz is capable of synthesizing and elongating many FA to produce DHA. Furthermore, while substrate was a significant predictor of FA composition in PERMANOVA, sample type explained much more of the variation (partial R^2^ = `r signif(permanova_substrate_R2,3)` versus `r signif(permanova_sample_type_R2,3)`, respectively) and suggests FA production of Labyrinthulids on artificial substrates may qualitatively reflect the abilities of those organisms *in situ*.

Lz's DHA content far exceeds that of most primary producers. In an analysis of 40 northeast Pacific marine macrophyte species spanning 21 orders in 4 phyla, no species exceeded on average 5 %TFA of DHA and only two had detectable DHA [@galloway_fatty_2012]. Field diseased eelgrass samples found here had `r signif(mean(data.prop.m.fin[data.prop.m.fin$sample_type=="diseased"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$FA=="C22.6w3",]$percent),3)` ± `r signif(sd(data.prop.m.fin[data.prop.m.fin$sample_type=="diseased"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$FA=="C22.6w3",]$percent),3)`% (mean ± sd, n = `r length(data.prop.m.fin[data.prop.m.fin$sample_type=="diseased"&data.prop.m.fin$substrate=="eelgrass"&data.prop.m.fin$FA=="C22.6w3",]$percent)`), placing it among the top macrophyte DHA sources in the northeast Pacific. For marine phytoplankton, Galloway & Winder [-@galloway_partitioning_2015] synthesized data for FA contents in 208 species and found the greatest DHA content in dinoflagellates (21±1 μg mg^-1^ DW, 17±9 %TFA, mean±sd), exceeded by Lz. Lz's DHA content does, expectedly, fall with many other heterotrophic protists researched for industrial LCPUFA production, including Lz's relatives [5 to over 50% lipid, @kumon_new_2006;@armenta_single-cell_2013]. With greater ecological relevance, five planktonic heterotrophic protists studied by Chu et al. [-@chu_quantitative_2008;-@chu_species_2009] fed biologically relevant diets (bacteria or phytoplankton) contained DHA at ~5-17 %TFA.

In eelgrass, patterns of laboratory segments and field-collected tissue were similar (Fig. 1,2). The culturing of Lz from field-collected lesioned eelgrass confirms our identification of EWD in the field, supporting the aggreement between laboratory and field samples. The eelgrass results imply a tradeoff between ALA and DHA with infection, with ALA being reduced and DHA being produced. A straightfoward explanation could be that Lz converts ALA to DHA, but because we did not specifically trace ALA's fate, we cannot definitively conclude this mechanism. Also, we emphasize that this tradeoff is not one-for-one, as ALA is likely being converted to other FA or being lost in degradation accompanying disease. Disease-associated degradation also likely explains the difference in ALA patterns between agar substrates and eelgrass: while ALA could be increased in ALA-less agar substrates with Lz presence, the ALA could be lost in ALA-rich eelgrass though degradation.  

While Lz-associated DHA increases found in both laboratory and field eelgrass were detectable, they were modest (Fig. 1B,D), which may indicate minimal relevance *in situ*. However, considering eelgrass's abundance as the foundation species in its ecosystems, the ubiquity of EWD [@sullivan_potential_2013], and the general scarcity of DHA in eelgrass bed primary producers [@galloway_fatty_2012;@kharlamenko_identification_2001], even this small production may be significant for eelgrass beds. Overall, our results highlight a potential role for Lz in secondary production in its communities, but ecosystem-level effects must still be confirmed. If abundant enough, Lz could be an important upgrader [@klein_breteler_trophic_1999] of LCPUFA-poor eelgrass.

Beyond considering FA production at community or ecosystem scales, investigating the potential relationship between disease severity and FA production should be a priority for future study. Specifically, investigations into the utility of FA production for pathogenic Lz may provide clues to how strain-specific virulence is modulated. Thraustochytrids, relatives of Lz, use FA like DHA for storage and the ectoplasmic net [@jain_docosahexaenoic_2007], an important component of infection for Lz [@muehlstein_host_1992]. Thus, FA production may be linked to disease in mechanism, not simply pathogen load.

We here demonstrated that a marine parasite produces exeptionally large amounts of DHA, finding that Lz may be an underrecognized source of LCPUFA in eelgrass ecosystems. While unsurprising in light of its relatives [@kumon_new_2006;@armenta_single-cell_2013], Lz's FA production is remarkable in expanding parasitism's potential ecological effects. If Lz's FA production scales *in situ*, this parasite may trophically upgrade food quality to consumers of its host via disease. This encourages studies on FA production in relation to disease in living plants to elucidate FA consequences of EWD. Future research on the trophic pathways of diseased tissue are needed to show whether Lz FA production enters higher consumers via enhanced eelgrass nutritional quality. This work demonstrates that EWD may not just be a destructive force: through producing valuable biomolecules that can regulate organismal and ecosystem production [@winder_landsea_2017], Lz could be a creator as well.

\  

### Acknowledgements

We thank the staff of the South Slough National Estuarine Research Reserve, Zofia Knorek, Nicole Nakata, Isaiah-Peacott Ricardos, and Alexa Romersa for field and laboratory assistance. We also thank Morgan Eisenlord for fantastic discussion on this subject and manuscript. We thank the three anonymous reviewers for constructive comments and suggestions, substantially improving our manuscript. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1309047. AWEG and JBS were supported by startup funds granted to AWEG by the University of Oregon.

\  

##### .

### Tables and Figures

Table 1. Summary of Lz substrates and sample types.
```{r flextable Table 1, echo=FALSE}
tab1=data.frame(
  "Substrate"=c("Serum seawater agar (SSA)","Eelgrass agar (EGA)","Eelgrass","Eelgrass"),
  "Sample.types"=c("Lz cells, Lz-colonized substrate, substrate alone","Lz-colonized substrate, substrate alone","Laboratory-inoculated segments: Lz-colonized substrate, substrate alone","Field-collected tissue: field diseased, field healthy"),
  "Eelgrass"=c("None","2nd/3rd leaves of turion, cleaned of epibionts and blended","7-cm piece of eelgrass from 2nd/3rd leaves of turion, cleaned of epibionts","Leaves cleaned of epibionts, portions of diseased and healthy tissue separated"),
  "Other.components"=c("Horse serum, peptones, yeast extract, glucose, germanium dioxide, noble agar","Filtered seawater, noble agar","Sterile seawater, sterile seawater agar plates","None"),
  "Sterilization"=c("Heating (dissolve agar), autoclaving","Heating (dissolve agar), autoclaving","None","None"),
  "Relevance"=c("Traditional artificial medium used to culture Lz","Medium plate with only eelgrass-derived resources without eelgrass host response","Controlled but realistic production of EWD-lesioned tissue","in-situ FA production")
)

tab1.flex=regulartable(tab1)
tab1.flex=set_header_labels(tab1.flex,Eelgrass="Eelgrass ?",Sample.types="Sample types",Other.components="Other components")
tab1.flex=merge_v(tab1.flex,j="Substrate")
tab1.flex=fontsize(tab1.flex,size=8,part="body")
tab1.flex=fontsize(tab1.flex,size=8,part="header")
tab1.flex=flextable::font(tab1.flex,fontname="Times New Roman",part="all")
tab1.flex=set_formatter( tab1.flex, 
  df = function(x) ifelse(is.na(x), "", x )  )
tab1.flex=bold(tab1.flex, part="header")
tab1.flex=align(tab1.flex,align="left",part="all")
tab1.flex=width(tab1.flex,width=1)
tab1.flex=width(tab1.flex,j="Other.components",width=1.5)
tab1.flex

```

##### .

``` {r fig1, echo=FALSE, warning=FALSE, fig.height = 5, fig.width = 6.65, dpi=300}
figure1
```

Fig. 1. (A,C) ALA (α-linolenic acid, 18:3ω-3) and (B,D) DHA (docosahexaenoic acid, 22:6ω-3) by substrate and sample type. Median (A,B) proportions (% total identified fatty acids) and (C,D) concentrations (μg mg^-1^ sample dry weight) are shown as diamonds with error bars (± 1 interquartile range). Substrates: serum seawater agar (SSA), eelgrass agar (EGA), and (eelgrass) leaf segments or tissue. Sample types: field-collected diseased eelgrass, field-collected healthy eelgrass, isolated *Labyrinthula zosterae* (Lz) cells, Lz-colonized substrate, and substrate alone. Grey points are individual analyzed samples. Sample sizes for each substrate-sample type combination are n=5 for each eelgrass-sample type combination, n=4 for SSA-substrate alone, and n=8 for all other combinations. Data analyses are summarized in Tables S1 and S2 in Supplement 2.

##### .

``` {r fig2, echo=FALSE, warning=FALSE, fig.height = 5, fig.width = 6.65, dpi=300}
Lz_NMDS
```

Fig. 2. Non-metric multidimensional scaling (NMDS; on Euclidean distances of arcsine-transformed FA proportions) plot of fatty acid (FA) compositions by substrate and sample type as in Fig. 1. Each point is an individual analyzed sample. Vectors show relative correlations of FA (>5% total identified FA of at least substrate-sample type combination) with NMDS dimensions.

##### .

Table S1: Univariate statistical analyses for fatty acid (FA) proportions. Kruskal-Wallis tests compared FA concentrations among SSA treatments followed by Dunn's test for posthoc pairwise comparisons.
```{r flextable Table S1, echo=FALSE}
FAprop.tab=FAprop.tab[,c(7,6,2,1,3:5)]
FAprop.flex=regulartable(FAprop.tab)
FAprop.flex=merge_v(FAprop.flex,j=c("FA","substrate"))
FAprop.flex=fontsize(FAprop.flex,size=8,part="body")
FAprop.flex=fontsize(FAprop.flex,size=8,part="header")
FAprop.flex=flextable::font(FAprop.flex,fontname="Times New Roman",part="all")
FAprop.flex=set_formatter(FAprop.flex, 
  df = function(x) ifelse(is.na(x), "", x )  )
FAprop.flex=bold(FAprop.flex, part="header")
FAprop.flex=align(FAprop.flex,align="left",part="all")
FAprop.flex=autofit(FAprop.flex)
FAprop.flex

```

##### .

Table S2: Univariate statistical analyses for fatty acid (FA) concentrations. 
```{r flextable Table S2, echo=FALSE}
FAconc.tab=FAconc.tab[,c(7,6,2,1,3:5)]
FAconc.flex=regulartable(FAconc.tab)
FAconc.flex=merge_v(FAconc.flex,j=c("FA","substrate"))
FAconc.flex=fontsize(FAconc.flex,size=8,part="body")
FAconc.flex=fontsize(FAconc.flex,size=8,part="header")
FAconc.flex=flextable::font(FAconc.flex,fontname="Times New Roman",part="all")
FAconc.flex=set_formatter( FAconc.flex, 
  df = function(x) ifelse(is.na(x), "", x )  )
FAconc.flex=bold(FAconc.flex, part="header")
FAconc.flex=align(FAconc.flex,align="left",part="all")
FAconc.flex=autofit(FAconc.flex)
FAconc.flex

```

##### .

Table S3: Summary of fatty acids (FA) used in this study. FA contributing to >1% of the identified FA for a given substrate-sample type combination are listed. Concentrations (as μg mg^-1^ dry weight) and proportions (as percentages in parentheses) of each FA are show for each substrate-sample type combination as medianinterquartile range. For FA, “x” denotes an unknown ω double bond position and “B” or “C” distinguish separately retained FA of the same notation.

```{r flextable Table S3, echo=FALSE}
FA.medIQR.flex
```

##### .
### LITERATURE CITED