Characterizing host-pathogen interactions between Zostera marina and Labyrinthula zosterae

infected host and parasite. Methods: Eelgrass ( Zostera marina ) shoots were placed in one of two temperature treatments, 11° C or 18° C, acclimated for 10 days, and exposed to a waterborne inoculation containing infectious Labyrinthula zosterae ( Lz ) or sterile seawater. At two-and ﬁ ve-days post-exposure, pathogen load, visible disease signs, whole leaf phenolic content, and both host-and pathogen-transcriptomes were characterized. Results: Two days after exposure, more than 90% of plants had visible lesions and Lz DNA was detectable in 100% percent of sampled plants in the Lz exposed treatment. Concentrations of total phenolic compounds were lower after 5 days of combined exposure to warmer temperatures and Lz , but were unaffected in other treatments. Concentrations of condensed tannins were not affected by Lz or temperature, and did not change over time. Analysis of the eelgrass transcriptome revealed 540 differentially expressed genes in response to Lz exposure, but not temperature. Lz -exposed plants had gene expression patterns

Introduction: Seagrass meadows serve as an integral component of coastal ecosystems but are declining rapidly due to numerous anthropogenic stressors including climate change.Eelgrass wasting disease, caused by opportunistic Labyrinthula spp., is an increasing concern with rising seawater temperature.To better understand the host-pathogen interaction, we paired whole organism physiological assays with dual transcriptomic analysis of the infected host and parasite.
Methods: Eelgrass (Zostera marina) shoots were placed in one of two temperature treatments, 11°C or 18°C, acclimated for 10 days, and exposed to a waterborne inoculation containing infectious Labyrinthula zosterae (Lz) or sterile seawater.At two-and five-days post-exposure, pathogen load, visible disease signs, whole leaf phenolic content, and both host-and pathogentranscriptomes were characterized.
Results: Two days after exposure, more than 90% of plants had visible lesions and Lz DNA was detectable in 100% percent of sampled plants in the Lz exposed treatment.Concentrations of total phenolic compounds were lower after 5 days of combined exposure to warmer temperatures and Lz, but were unaffected in other treatments.Concentrations of condensed tannins were not affected by Lz or temperature, and did not change over time.Analysis of the eelgrass transcriptome revealed 540 differentially expressed genes in response to Lz exposure, but not temperature.Lz-exposed plants had gene expression patterns

Introduction
Despite occupying less than 0.1% of the ocean floor (primarily in estuaries and embayments), seagrasses provide numerous ecosystem services, including shoreline stabilization and protection from erosion, nutrient cycling, sediment retention, reduction of pathogens in the water column, and provision of critical nursery habitat for numerous invertebrate and fish species, many of which are of economic importance to the seafood industry (Nordlund et al., 2016;Lamb et al., 2017).Seagrass meadows are threatened by a wide range of global and local stressors associated with anthropogenic activity such as coastal development, increasing sea surface temperature, rising sea levels, sediment and nutrient runoff and disease (Orth et al., 2006;Sullivan et al., 2013).Consequently, over the past century, global seagrass beds have experienced approximately 19% loss in cover (Dunic et al., 2021).
Zostera marina, the most widely distributed seagrass species in the northern hemisphere (Short et al., 2007), is susceptible to infection with the protist Labyrinthula zosterae (Lz), which causes eelgrass wasting disease (EWD).While chronic Lz infections are nearly ubiquitous in distribution, acute epidemics of wasting disease can cause declines in eelgrass populations, most notably in the 1930s when Atlantic eelgrass populations were nearly eliminated along the coasts of North America and Europe (Short et al., 1987).Lz is a pathogenic Labyrinthulomycete and can be transmitted through direct contact between healthy and infected blades within an eelgrass bed (Muehlstein, 1992) or through the water column (Groner et al., 2014;Groner et al., 2018).Upon infection, Lz targets chloroplasts of eelgrass leaves, compromising photosynthesis, and ultimately causing necrosis of eelgrass leaves, and in some cases, entire plants (Muehlstein, 1992).Non-lethal impacts of EWD include reductions in growth and below-ground storage of carbohydrates, which may impact seasonal resilience of eelgrass meadows (Graham et al., 2021).Strains of Lz vary in virulence, raising questions about mechanisms of virulence employed by this pathogen and potential costs associated with those mechanisms (Martin et al., 2016).
Eelgrass has a robust defense system which includes proteins and surface-associated metabolites like fatty acids and phenolics, pcoumaric acid, rosmarinic acid, and zosteric acid (Papazian et al., 2019).These compounds play a wide-range of roles including but not limited to signaling molecules, antioxidant activity, free radical scavenging activity, and regulators of auxin transport (Ma et al., 2016).Phenolic compounds which are produced in response to stress and pathogen presence can be also considered as defenserelated secondary metabolites.Indeed, some phenolic compounds have been shown to inhibit growth of Lz in vitro (Vergeer and Develi, 1997) and concentrations of total phenols (% dry mass) are increased in diseased plants in the field (Groner et al., 2016).In terrestrial plants, phytohormones including abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA) play significant roles in regulating defense responses against a great number of biotic and abiotic factors (Thaler et al., 2004;Tamaoki et al., 2013), however these hormones have not been well explored in marine plants such as eelgrass.
As with many marine diseases, EWD is sensitive to temperature (Tracy et al., 2019).Both EWD prevalence and severity are correlated with sea surface temperatures in the Salish Sea (WA, USA), and in the Isles of Scilly (UK), (Bull et al., 2012;Groner et al., 2021).Increased temperatures enhance the growth-reducing effects of EWD on eelgrass (Bull et al., 2012).These changes are likely driven by both the host and pathogen.Lz is also sensitive to temperature, with faster in vitro growth documented at 18°C compared to 11°C for a strain isolated from the Salish Sea (WA, USA) (Dawkins et al., 2018); similarly, increased disease severity was noted in eelgrass exposed to Lz and held at 18˚C compared to 11˚C (Agnew et al., 2022).Recent studies have correlated heatwaves to substantial seagrass die-offs (Strydom et al., 2020;Groner et al., 2021), reduced restoration success (Aoki et al., 2020), and changes in fatty acid composition (Beca-Carretero et al., 2018), showing that seagrasses are sensitive to warming conditions.
To better understand the factors that contribute to Lz infection and provide insight on mechanisms of virulence and host defense, we characterized the production of defensive compounds in Z. marina and the transcriptional responses of both Z. marina and Lz following a controlled Lz challenge conducted at ambient and warm seawater temperatures.We hypothesize that both increased temperature and Lz infection will alter expression of genes related to immune function and stress response in Z. marina.More specifically, we hypothesize that phenolic and hormonal gene pathways will be enriched in response to presence of Lz, but we expected these responses may be reduced at warmer temperatures that could be stressful to Z. marina.Finally, we hypothesize that Lz will exhibit a change in gene expression in response to temperature, reflecting the association between EWD and warmer seawater temperatures found in nature (Brakel et al., 2019;Groner et al., 2021).

Experimental design
Similarly sized eelgrass ramets were collected from False Bay, San Juan Island, WA (48.550°N, 123.008°W) on 29 July 2015.Epiphytes, grazers and damaged leaves were removed from ramets and they were held in flow-through 14°C, 0.2 mm filter-sterilized seawater at Friday Harbor Labs, WA, USA.
Eelgrass ramets were placed in one of two temperature treatments (target temperatures 11°or 18°C, actual temperatures ranged ± 1°C from target) on 3 Aug 2015 and allowed to acclimate to these treatments for 10 days.These temperatures represent typical summer sea surface temperatures found locally along coastal eelgrass beds at high and low tides, respectively.A constant exposure to 18°C would be unusually warm in this region.Each treatment was replicated eight times for a total of 16 experimental units.Each replicate unit contained four ramets.Experimental units were 4-L containers with 3500 mL of 0.2 mm filter-sterilized seawater.Experimental units were split into groups of eight that were clustered in a common cooler.Each cooler received flow-thru 0.2 µm filtered and UV-irradiated seawater.Seawater temperatures were continuously monitored and adjusted using Honeywell UDA2182 pH controller and Honeywell Durafet III electrodes.After flowing into the cooler, water was then piped into the 8 experimental units.Water temperatures were maintained by circulating it through chilled water set 1°C below target temperatures and, if necessary, heating it with aquatic heaters placed in the coolers (i.e.O'Donnell et al., 2013).Two temperature loggers (iButtons, Whitewater, WI) set to record every 30 min were placed in each cooler for independent confirmation of temperature.Full-spectrum lights were kept to a 12: 12 hr light: dark schedule with a mean light output of 161 ± 3 (mean ± SE) µmol/m 2 /sec PAR below the water surface.To reduce algal blooms, which could clog plumbing and block light, tanks were treated daily with 0.67 ppm (final concentration) of Germanium dioxide (GeO 2 ; Markham and Hagmeier, 1982).
A 36-hr incubation with the pathogen inoculant was conducted after 10 days of temperature acclimation (Figure 1).Immediately prior to inoculation, all ramets were photographed, weighed, labelled with uniquely colored zip-ties, and pin-pricked at the sheath for later measurements of growth.All ramets within a replicate were placed in Whirlpaks with 98 mL of sterile seawater.Control treatments were inoculated with 2 mL of sterile seawater.Lz treatments were inoculated with 2 mL of inoculant for a final concentration of 2.5 x10 5 cells per mL (additional details in Supplementary Material).After adding the inoculate, the Whirlpaks were closed and floated in treatment containers to maintain target temperatures.36 hours later, ramets were removed from Whirlpaks and placed back into treatment containers.The Lz strain used in this experiment, 8.16.D, has been used in several previous experiments, (e.g., Groner et al., 2014;Groner et al., 2018) and was isolated in 2011 from nonflowering adult Zostera marina shoots that were collected at Picnic Cove, Shaw Island, 48°33.942'N, 122°55.448'W).Timeline of experimental treatments and sampling.After collection and cleaning of plants to remove grazers, epiphytes and diseased tissue, eelgrass ramets were held at 14°C for 1d and then transferred in groups of 4 ramets to 11°C or 18°C for a 10d acclimation period.Each group of 4 ramets was replicated 8 times during the acclimation period.After that, the groups were equally divided into control and Lzexposed treatments at each of the two temperatures.Lz or sham exposures occurred for 36 hr.Subsampling of plants occurred at days 2 and 5 post-Lz exposure.
Forty-eight hours after inoculation (experiment day 13), a random selection of half of the individuals (48 total) from each experimental unit were removed for sampling (Figure 1).Each of those individuals was photographed and weighed.Then the roots and shoots were frozen in liquid nitrogen for transcriptomic sampling.The oldest leaf was split lengthwise into two pieces.One half was frozen in liquid nitrogen for transcriptomic analysis, and the other was stored in 70% molecular grade ethanol at room temperature for subsequent measurement of Lz load using quantitative PCR (qPCR) as described in Groner et al. (2018).The second oldest leaf was frozen in liquid nitrogen for measurement of total phenolic compounds and condensed tannins.Five days after inoculation, the experiment was ended because the plants had severe disease signs.At this point, the remaining ramets were photographed and sampled (as described above) for measurement of visual disease signs (e.g., lesions with a distinct black border, sensu Groner et al., 2014) as well as phenol and tannin concentrations (Figure 1).

Measurement of total phenolic compounds and condensed tannins
Total phenolic compounds were measured in 96 well microplates using the methods described in Groner et al. (2018).Briefly, sampled eelgrass leaves were frozen (-20°C) and sent on dry ice to the Shannon Point Marine Center for analysis.Frozen tissues were lyophilized and ground to a fine powder in an SPEX mixer/mill (SPEX, Metuchen, NJ).Ground tissue (9-10 mg) was extracted overnight in 1 mL of HPLC-grade methanol, then diluted in ANSI Type I water (19 parts water to 1 part extract).Forty µL of 40% Folin-Ciocalteu's Reagent (Sigma F9252) was added to each 100 µL aliquot (n=3 per sample) of each of the diluted extracts, mixed for 5 min prior to the addition of 100 µL of 2N sodium carbonate.Samples were shaken for another 30 min in a 50°C chamber and then the absorbance of each cell was read at 765 nm.Concentrations were standardized using native standards from Z. marina collected from Ship Harbor WA, using caffeic acid as a secondary standard (Groner et al., 2016).The remaining methanolic extracts were used for measurements of condensed tannins, which were conducted with a sulfuric acid method (Bate-Smith and Rasper, 1969;Nitao et al., 2001) that was modified for use in a microplate reader.Condensed tannins (proanthocyanidins) were cleaved by the addition of 43% sulfuric acid in methanol (250 µL/well containing 100 µL of a methanol extract) for 30 min at 50°C for.Tannin cleavage produces red-colored chromophores (anthocyanins), which were then quantified spectrophotometrically at 550 nm using a Biotek Synergy multiplate reader.Cyanidin (Cayman Chemical, purity >98%) was used as a standard.All assessments were run in triplicate.

Quantification of Lz abundance
Samples (mean dry weight of leaf material = 7 mg) were preserved in 70% ethanol until DNA extractions were performed using a Qiagen DNeasy Plant Mini Kit following the manufacturer's instructions.Prior to DNA extractions, samples were removed from ethanol, patted dry, weighed and then transferred into a clean 1.5ml micro-centrifuge tube with one 3mm tungsten carbide bead (Qiagen) per tube.Samples were then placed on dry ice for 10-15 minutes and lysed using a Qiagen TissueLyser at a speed of 25 Hz for 5 minutes.
Presence of Lz DNA was detected and quantified using qPCR detection using primers and probe targeting the Lz ITS region (Bockelmann et al., 2013); Laby_ITS_Taq_f (5'-TTG AAC GTA ACA TTC GAC TTT CGT-3') and the reverse primer Laby_ITS_Taq_r (5' -ACG CAT GAA GCG GTC TTC TT -3') were used, in addition to the TaqMan probe Laby_ITS_probe (5' FAM-TGG ACG AGT GTG TTT TG -MGB-NFQ 3').qPCR reaction conditions and the standard curve (serial dilution between 1.37-1.37*10^4cells) were as described by Groner et al., 2018 on run on the Applied Biosystems 7500 Fast Real-Time PCR System.For each qPCR plate, we aimed for R2 = 0.999 and an efficiency of 90-110%; the mean reaction efficiency was 105.5% +/-1.5% (cross three total plates).We used a cut-off of 1.37 cells as our limit of detection (amplification of 1.37 cells was achieved in all reactions in each cell curve).
Statistical analyses of treatment effects on phenolic compounds, tannins, disease signs and Lz DNA Independent and interactive effects of temperature and Lz exposure on concentrations of phenols, concentrations of tannins, presence of suspected disease signs and log 10 transformed Lz DNA concentrations in leaf tissues were quantified for the time points when these measurements were taken (R v. 4.1.1).Linear regression was used to quantify treatment effects on total phenolic compounds, condensed tannins, and Lz DNA concentrations, while logistic regression was used to quantify treatment effects on the presence of disease signs (package lme4, Bates et al., 2015).In all models, the 4L container replicate was included as a random effect.

RNA extractions and transcriptome sequencing
Total RNA from four samples per treatment was extracted following the Qiagen RNeasy plant kit protocol; samples were chosen based on both quality of RNA and qPCR results.DNA was removed from extracted RNA using the Turbo DNA-free treatment according to the manufacturer's instructions (Ambion Inc., The Life Technologies Corporation ™ , Grand Island, NY).Removal of DNA was confirmed by using RNA (1 ml) as template in a qPCR targeting 18s ribosomal DNA as previously described (Burge and Friedman, 2012).RNA concentrations were quantified using the NanoDrop® ND-1000 (NanoDrop Technologies, Wilmington, DE).Samples were shipped to McGill University and Genome Quebec for library preparation and sequencing.
RNA quality was assessed prior to library preparation using an Agilent Bioanalyzer.The Illumina TruSeq stranded cDNA kit was used for library preparation and barcoding; quality assessment of libraries indicated that one sample was not of high enough quality for sequencing.The remaining samples were distributed across two sequencing lanes on an Illumina HiSeq 4000 platform (2x150 bp paired end reads).

Transcriptome assembly and annotation
Prior to de novo assembly, reads were quality screened (FastQC Version 0.11.4,Babraham Bioinformatics), and sequences were screened for remaining Illumina adapter sequences (Trimmomatic, Bolger et al., 2014).De novo assembly was carried out using reads from all fifteen samples using Trinity (Grabherr et al., 2011;Haas et al., 2013, Trinity RNASeq Github).Following transcriptome assembly, reads for each sample were matched back to gene isoforms.
A multi-step annotation process was used to annotate the transcriptome in order to separate isoform consensus sequences (or contigs) matching to Z. marina and those non-matching contigs ("non-host").First, a database of annotations was created from the Z. marina genome (GCA_001185155.1)using a BLASTx search of the UniProtKB/Swiss-Prot database.The BLASTx algorithm (Altschul et al., 1990) was used to search the Z. marina database with a threshold E-value of 0.00001; each contig matching as Z. marina was filtered out.Similarly, the blastx algorithm was used to search the Swiss-Prot database with a threshold E-value of 0.00001 to form the non-Zostera annotation for subsequent analysis.Finally, to determine which non-host contigs have significant similarity to Labyrinthula and other closely related 'slime-mold' species, the blast algorithm was used to match non-host contigs against the recently added Labyrinthula sp.Ha genome (NCBI Accession # GCA_015227615.1) and also against the genomes of two other, better annotated genomes, a thraustochytrid, Hondaea fermentalgiana (NCBI Accession # GCA_002897355.1) and the amoeba model organism, Dictyostelium discoideum (NCBI Accession # GCF_000004695.1).
Gene ontology (GO) information was used to annotate both Z. marina and non-Zostera transcriptome data separately.Gene ontology IDs and associated GO Slim terms for biological processes, cellular components, and molecular function categories were downloaded from the UniProtKB/Swiss-Prot database.Swiss-Prot identifiers from BLASTx output were joined to gene ontology terms.Since assembly and annotation were conducted at the isoform level, transcriptome data was streamlined to conduct differential expression analysis on genes.Annotated sequences were sorted by contig name and e-value; contig names included a contig, gene, and isoform identifier.The isoform entry with the lowest e-value was retained for each gene.The final annotated dataset for both Z. marina and non-Zostera transcriptomes included only one entry per gene.Annotated genes were merged with count data to prepare for differential expression analysis.The data is deposited in the NCBI database in the SRA as Bioproject (PRJNA990835).

Analysis of Zostera and non-Zostera transcriptomes
Statistical differences in count data were conducted for Z. marina and non-Zostera annotated genes using EdgeR (Robinson et al., 2010;McCarthy et al., 2012;Chen et al., 2016).Briefly, we utilized EdgeR to filter and normalize the count data based on the library size and examine samples for outliers using MDS plots.Next, a multi-factorial model including estimated dispersion (i.e.estimate of biological variation) was used with multiple contrasts to calculate fold change, log CPM (counts per million), and raw and Benjamini-Hochberg adjusted p-values to correct for False Discovery Rate.
We tested for statistical differences using several contrasts followed by the primary questions of interest, respectively examining the relationships of the pathogen-by-temperature interaction, temperature alone, and pathogen exposure alone: 1) What genes are differentially expressed between Lz-exposed shoots at elevated temperatures?At ambient temperatures?
2) What genes are differentially expressed between elevated and ambient temperatures, irrespective of Lz-exposure?
3) What genes are differentially expressed between Lz-exposed and control shoots, irrespective of temperature?
Very few genes were differentially expressed in questions 1 and 2 (25 genes and 0 genes, respectively).Thus, for the remainder of the analysis we focused on the effect of Lz exposure alone.Differentially expressed genes were pooled into host and pathogen genes.Genes were considered differentially expressed based on FDR adjusted p-values from edgeR (p < 0.05).
To identify groups of Z. marina genes with similar expression patterns (eigengenes) and associate these patterns with experimental conditions, a Weighted Gene Co-expression Network Analysis (WGCNA) was conducted using the WGCNA package in R (version 1.69;Langfelder and Horvath, 2008).The WGCNA was only conducted for Z. marina genes, as transcriptome sequencing produced a high magnitude of information for Z. marina that was difficult to parse without additional analysis.After identifying significant modules, overrepresented biological process and molecular function GOterms were identified using the GO-MWU enrichment method (Wright et al., 2015).More details on WGCNA and GO-MWU analyses are available in the Supplemental Material.

Disease progression and pathogen load
The inoculations were successful in transmitting Lz to the exposed plants.Two days after exposure, Lz concentrations in eelgrass tissue were 1.1 orders of magnitude greater in the 18°C treatment as compared to the 11°C treatment (with concentrations of 4.2 ± 4.3 x10 4 and 3.7 ± 3.6 x10 3 (mean ± SE) Lz cells/mg of eelgrass tissue, respectively) (Figure 2).These pathogen concentrations were significantly higher than the control treatments, where no Lz was detected using qPCR (t value = 12.3, p < 0.0001) No other factors affected Lz abundance in our linear regression (all p > 0.05).
Visible disease signs were recorded in all treatments (Figure 3).Pathogen exposure was a significant predictor of disease signs on day 2 and a marginally significant predictor of disease signs on day 5. On day 2, the disease signs were recorded in 93.7 ± 8.8% (mean +/-SE) of the exposed plants and 28.8 ± 19.7% of the control plants (t = 2.1, p = 0.04).Neither temperature nor the temperature by exposure interaction were predictors of visual disease signs.These values increased on day 5 where disease signs were recorded in 100 ± 0% (mean +/-SE) of the exposed plants and 37.5 ± 17.7% of the control plants (t = 1.99, p = 0.07).
The concentrations of total phenolic compounds and condensed tannins were not impacted by temperature or Lz treatments on day 2 (Figure 4).By day 5, total phenolic compounds were significantly lower in the warmer exposed treatment, where they were 0.70 ± 0.17% of dry mass as opposed to ranging between 1.15 and 1.39% of dry mass in the other three treatments (t value = -3.298,p < 0.001 for the temperature by exposed interaction).There was no significant impact of temperature or Lz exposure on the concentration of condensed tannins in eelgrass leaves on day 5.

Transcriptome assembly and annotation
Trinity de novo assembly yielded 1,065,804 isoforms (767,296 potential genes) with an average read length of 654.71 bp (median of 349) and an N50 value of 1075.BLAST searches of the de novo assembly yielded 281,600 contigs matching the Z. marina genome.784,204 contigs did not match the Z. marina genome and thus were categorized as 'non-host' contigs.

Differential expression analysis -host contigs
Of the 3,096 contigs identified as host genes and annotated, 540 genes were differentially expressed, with 338 genes with significantly higher expression levels in exposed blades and 202 genes with lower expression levels (Supplementary Table 1).The WGCNA identified thirteen module eigengenes comprised of genes with similar expression patterns for 3,096 Z. marina contigs (Figure 5A; Supplementary Table 2).Of these eigengenes, four -named "brown", "yellow", "blue", and "magenta" -were significantly  associated with Lz exposure (Figure 5B).These modules, labeled as color categories, are described in detail below.
The module eigengene "yellow" was positively correlated with Lz exposure (r 2 0.52, p= 0.049) and contained five differentially expressed genes (Figures 5B, C: Supplementary Table 2).While there were no significantly enriched GO terms in this module, common biological processes included proteolysis, cellular response to oxidative stress, response to osmotic stress, regulation of jasmonic acid metabolic process, and response to abscisic acid, and common molecular functions were binding of ATP, RNA, and metal ions (Supplementary Table 2).

Discussion
Exposure of eelgrass ramets to Lz resulted in rapid progression of EWD, with necrotic lesions forming within two days postexposure and plant death occurring after five days.The temperature did not alter the pathogen load or progression of EWD in the Lz exposed plants.This may be due to the high concentration of pathogen in the exposure, which was higher than in previous studies with this Lz strain (Groner et al., 2014;Groner et al., 2018) and/or pre-acclimation of the plants to the temperature treatments prior to exposure.Analysis of the eelgrass transcriptome revealed changes in gene expression, with patterns consistent with increased defensive responses through altered regulation of genes associated with phytohormone biosynthesis, stress response, and immune function.Analysis of non-Zostera genes (likely Lz genes), revealed expression patterns suggestive of Lz disrupting host immune responses and undergoing phagocytosis.This study reveals the importance of specific pathways related to plant defense against pathogen presence, and provides molecular evidence to support previous phenotypic observations of hostpathogen interactions of eelgrass wasting disease.

Host responses
The genome of Z. marina allowed us to identify 3,096 host contigs, and we detected 540 differentially expressed host genes responding to Lz infection as opposed to increased temperature (Supplementary Table 1).We found expression changes in genes and/or pathways involved with immunity such as pathogen detection, defense-related cell-cell signaling, defense-related metabolite production, and apoptosis.We also detected genes potentially involved in salt stress.The occurrence of many differentially expressed defense-response related GO terms, including genes associated with microbial detection and defense response to bacterial or fungal pathogens, indicates that Z. marina is responding to Lz infection.Plant innate immune system pathways are highly conserved and are initiated by recognition of pathogen-associated molecular patterns (PAMPs) and damage-associated molecular patterns (DAMPs) by membrane localized receptors (Choi and Klessig, 2016).Calcium signaling is a major pathway in pattern-triggered immunity, where binding of PAMPs and DAMPs to membrane receptors initiates the rapid release of calcium ions (Ca 2+ ) by receptor like kinases and their Ca 2+ dependent protein kinases and mitogen activated protein kinases (Hou et al., 2018).For example, calcium/calmodulin-dependent protein kinase (CcaMK) genes have been found to play a role in maintaining endosymbionts and in pathogen resistance in tomato, by promoting accumulation of hydrogen peroxide (Wang et al., 2015).Genes related to calciumand protein-kinases show up eight and 43 times, respectively, in the blue module.In particular, upregulation of the calcium/calmodulindependent protein kinase (CcaMK) type 1 (CaM kinase I) suggests the role of Lz-related PAMP/DAMP detection and may warrant further exploration.
Pathogen detection leads to the induction of signaling pathways to induce immune and defense responses.We found that the cAMP-dependent protein kinase catalytic subunit gamma is significantly upregulated with Lz-exposure.Changes in the abundance of cAMPs can affect the activity of other signaling pathways and cellular processes, such as protein kinases and transcription factors (S ́wiezȧwska et al., 2018), demonstrating that Z. marina may respond to Lz infection by altering upstream signal transduction.
Phytohormone production in plants is important for defenserelated signaling and is turned on after PAMP/DAMP detection.For example, compounds such as jasmonic acid (JA), salicylic acid (SA), and abscisic acid (ABA) are involved in general activation of plant defense systems against pathogens (Thaler et al., 2004;Tamaoki et al., 2013).In particular, SA is used in defense against biotic factors such as insects, mites, fungi, and bacteria and JA is used in defense pathways against saprophytic microbes like Lz (Tamaoki et al., 2013;Zhang et al., 2020).ABA is triggered in response to both abiotic and biotic stressors, and has complex interactions with both JA and SA responses to infection (Fan et al., 2009).We found GO terms for SA, JA, and ABA phytohormones [GO:0080140, GO:0009695, GO:0009753, and GO:0009751] were highly represented in the blades experimentally exposed to Lz, although these gene ontologies were not significantly enriched.Gene expression patterns in this study were indicative of phytohormone signaling via enrichment of serine/threonineprotein kinase PCRK2 gene which is involved in the resistance to bacterial pathogens and salicylate biosynthesis during pathogenic infection (Sreekanta et al., 2015).Methylesterase 1 gene, which is required during the conversion of methyl salicylate into SA was downregulated in exposed blades.These data suggest that Z. marina is altering production of these defense-related compounds.We found gene expression changes in two genes related to JA biosynthesis: chloroplastic rhomboid-like protein 11 (regulation of jasmonic acid metabolic process [GO:0080140]) and phospholipase A (defense response to fungus [GO:0050832]; jasmonic acid biosynthetic process [GO:0009695]).Phospholipase A is responsible for the release of linolenic acid from the chloroplast membrane during defensive signaling, which induces JA production (Yang et al., 2007;Mata-Peŕez et al., 2015).The detection of differentially regulated genes in the JA pathway is particularly interesting, as the role of JA in seagrasses is unclear.A congener of Z. marina, Z. muelleri, lost genes encoding jasmonate methyltransferase, which converts JA into a volatile compound, methyl jasmonate, but maintain other genes along the JA pathway, including for jasmonate synthesis, and signaling (Lee et al., 2016).Although the JA-associated genes were not differentially expressed in our study, changes to expression in these genes suggests a role for the jasmonic acid pathway in responding to Lz infection.
Gene expression changes also indicated down-stream effects of phytohormone production.In maize roots, expression of CHC1 gene increases in response to SA or ABA, suggesting involvement of CHC1 in the SA signaling pathway in maize defense responses (Zeng et al., 2013).Upregulation of the probable clathrin heavy chain 1 gene (CHC1) ([GO:0005198]) in Lz-infected blades suggests a similar response in eelgrass.We detected gene expression changes in factors involved in stress signaling which either require ABA or interact with ABA to trigger defense responses.For example, the Stype anion channel SLAH3 is an anion channel that can be triggered by ABA during a stress response (Roelfsema et al., 2012), and was downregulated in response to Lz infection.The B3 domaincontaining protein was also downregulated.This protein is a transcription factor for ABA and inhibits the production of ABA (Brady et al., 2003), so a decrease in expression of B3 could result in an increase in ABA production, indicating a greater stress response to Lz infection.One specific pathway in which ABA affects plant defense is in the callose pathway, which is involved in defense of fungal pathogens in plants.Specifically, ABA mutants displayed impaired resistance to necrotrophic fungi, and ABA addition to the infection site mimicked callose deposition and increased resistance to necrotrophic fungi (Mauch-Mani and Mauch, 2005).Therefore, it is not surprising to observe altered ABA expression in Z. marina infected with Lz, which has necrotic capabilities.
Additional immune responses in plants include the production of antimicrobial compounds and apoptosis of infected cells, the latter of which would be important for intracellular infections, such as with Lz.Upregulation of peroxisomal targeting signal 1 receptor indicates potential defenses at sites of pathogen entry.This gene is associated with concentrating antifungal glucosinolate derivatives in the peroxisome for eventual transport to sites of fungal infection entry (Bednarek et al., 2009).It would be valuable to assess whether these compounds can inhibit Lz as well.Differential expression of GDP-L-fucose synthase supports the role of apoptosis as an immune response to Lz. GDP-L-fucose plays a role in stomatal and apoptosis-related defense in Arabidopsis, indicating a similar role in Z. marina immunity (Zhang et al., 2019).Differential expression of EKC/KEOPS subunit genes, which are involved in apoptotic processes via telomere capping and elongation (Downey et al., 2006), also support the hypothesis that apoptosis is an important defense against Lz.
The production of phenolic compounds is another major immune response of plants.Many of phenolic metabolites detected in Z. marina have been implicated in Z. marina defense.These include flavonoids in sulfated and unsulfated form, as well as acids like caffeic, p-coumaric, ferulic, and zosteric, and rosmarinic acids, all of which are hydroxy, sulfoxy, or ester forms of cinnamic acid (Papazian et al., 2019).Resistance of Z. marina to Lz is hypothesized to depend upon phenolic acids, especially caffeic acid which has demonstrated activity against Lz in bioassays (Vergeer and Develi, 1997;Trevathan-Tackett et al., 2015).Lz has also been shown to induce production of total phenols in some cases (McKone and Tanner, 2009), and inhibit them in others, potentially as a mechanism to circumvent this defense.In addition, production of phenols in seagrasses are reduced at warmer temperatures (Vergeer et al., 1995).This is reflected in the reduced concentrations of phenols found five days post Lz-exposure in the warmer temperature treatment.
The shikimate pathway is the mechanism through which plants produce aromatic phenolic compounds in response to stress (Santos-Sańchez et al., 2019).Four genes -caffeic acid 3-O-methyltransferase, caffeoylshikimate esterase, 4-coumarate-CoA ligase-like 7, and cinnamoyl-CoA reductase -responsible for coding enzymes involved in the biosynthesis of different forms of cinnamic acids were differentially expressed in the exposed treatment, indicating a role for the shikimate pathway in Lz infection response.Caffeoylshikimate esterase (CSE), converts caffeoyl shikimate into caffeic acid, which has been shown to inhibit Lz growth in vitro (Vergeer and Develi, 1997;Vanholme et al., 2013).Another differentially expressed gene, caffeic acid 3-O-methyltransferase, is associated with converting caffeic acid into products associated with the lipid biosynthesis pathway.4coumarate-CoA is also associated with production of additional secondary compounds including isoflavonoids and furanocoumarins as phytoalexins (i.e., inhibitors of pathogen growth) (Hamberger and Hahlbrock, 2004).Cinnamoyl-CoA reductase is involved in the formation of phenolic compounds associated with the hypersensitive response (Lauvergeat et al., 2001).The hypersensitive response in plants is akin to the innate immune response in animals and is responsible for isolated apoptosis surrounding an infection in order to prevent further spread (Morel and Dangl, 1997).
Apart from individual roles in the shikimate pathway, all four differentially expressed genes in this pathway are also associated with lignin biosynthesis.The role of lignin is unclear in aquatic plants such as Z. marina.While it may be protective against microbial attack, it is typically concentrated in the slower growing rhizomes, and not in the blades (Klap et al., 2000), where gene expression was measured in this study.
In addition to genes involved in immune activation, hormone pathways, and phenolic production, the transcriptional patterns of Lz exposed plants and gene enrichment results suggest responses to environmental stress, including salt stress.Stress-related genes were significantly upregulated in the exposed treatment, indicating that Lz infection may cause damage to cells in a way that triggers this response.The mitochondrial phosphate carrier protein 3 (MPT3) is induced under high salinity conditions (Zhu et al., 2012).Increased expression of the MPT3 gene suggests that Z. marina may be experiencing salt stress as Lz breaks through host physical barriers.Similarly, the sorting nexin 1 protein regulates the process of accumulating nitrous oxide, which is an important signaling molecule in responses to abiotic stressors, including salt stress (Li et al., 2018).Higher sorting nexin 1 protein gene expression implies that Z. marina may accumulate more nitrous oxide in an effort to respond to abiotic stress.Finally, the upregulation of the calcium/ calmodulin-dependent protein kinase type 1 (CaM kinase I) gene suggests changing cytosol calcium concentrations, which is a known plant response to stressors such as anoxia, increases temperature, and increased salinity (White and Broadley, 2003).Changes in expression of genes associated with stress response in Z. marina may inform our understanding of how molecular functioning is altered upon exposure to environmental changes.This information can also shed light on how Lz infection cause stressful conditions for Z. marina.

Pathogen responses
By separating host from non-host contigs, we aimed to identify genes or pathways that may contribute to the virulence of L. zosterae during experimental infection of Z. marina.We found up-regulation of genes potentially involved in breakdown of host defense, chemotaxis and phagocytosis, and metabolism.Although an unannotated genome of another species of Labyrinthula (Labyrinthula sp.Ha) is newly available, very few contigs had significant matches to the current version of the genome, and only five of the differentially expressed genes from the non-host dataset had significant BLAST hits to this genome.However, the majority of the non-host contigs had significant matches to related slime-mold species including Dictyostelium species and Hondaea fermentalgiana (a thraustochytrid closely related to the genus Labyrinthula), suggesting that the identified non-host genes were likely from Lz inoculated in the infection treatment.
Zostera species (including Z. marina) produce high concentrations of a sulfate ester, zosteric acid, as a defense mechanism that protects blades from microbial biofouling (Grignon-Dubois et al., 2012;Vilas-Boas et al., 2017).We found a 24-fold up-regulation of a non-host contig that matched to the arylsulfatase gene (GBG30795.1) of H. fermentalgiana.Arylsulfatases are enzymes that break down sulfatides and organic sulfate esters.In other protists, arylsulfatases are produced in response to sulfate deprivation (Niedermeyer et al., 1987).However, for the fungal plant pathogen, Colletotrichum gloeosporiodes, arylsulfatase expression increases during penetration of host leaf tissue after experimental inoculation (Goodwin et al., 2000).In Lz, an increase in arylsulfatase expression may facilitate attachment by breaking down the anti-biofouling zosteric acid produced by the host.However, the potential role of arylsulfatases to facilitate protist pathogen attachment and colonization has yet to be investigated.Initial studies describing EWD observed Lz directly penetrating mesophyll cell walls, damage suggestive of feeding on plant cell organelles such as chloroplasts, and moving rapidly through the host tissues (Muehlstein, 1992;Raghukumar, 2002).In support of these microscopic observations and supporting the role of plant consumption as a mechanism of pathogenesis, we identified several up-regulated genes that are likely involved in movement, i.e. chemotaxis and cytoskeleton rearrangement and organization, as well as phagocytosis.For example, contigs homologous to shkB, sgcA and roco5 genes were all up-regulated after experimental inoculation of Lz.In D. discoideum, null mutants for shkB and sgcA have reduced chemotaxis and phagocytosis ability (Moniakis et al., 2001;Veltman and Van Haastert, 2006;Veltman et al., 2008), and null mutants for the roco5 genes have slow or no motility (Sawai et al., 2007).Six up-regulated non-host contigs were homologous to genes that are up-regulated during phagocytosis and/or have gene products which are a part of the 'macropinocytosis proteome' of D. discoideum (Sillo et al., 2008;Journet et al., 2012).Additionally, in D. discoideum, dhcA proteins cluster to phagosomes and promote phagolysosome fusion (Rai et al., 2016), and lysosomal aspartic protease CatD, a ubiquitous hydrolytic enzyme with a protein of homology in D. discoideum, is known to be a lysosomal pepsin in the macropinosome and is involved with phagosome maturation (Vines and King, 2019).
Genes involved in amino acid and lipid metabolism were upregulated among the non-host contigs.M42 peptidase, a cocatalytic metallopepsidase (cytosolic enzyme) described in bacteria and archaea (reviewed by Appolaire et al., 2016), was the most expressed non-host contig.M42 may play a role in pathogenesis in bacteria by using or modifying exogenous proteins including immunological antibacterial peptides.Other peptidases were up-regulated, such as Carboxypeptidase Y and an otubain-like ubiquitin thioesterase, and may be involved with the catabolism of both large and small peptides in Lz.Ethylmalony-CoA decarboxylase, polyketide synthase, beta-ketoadipyl-CoA thiolase, and patatin-like phosphatase domain-containing protein are all involved with fatty acid synthesis and may support membrane formation and increase membrane fluidity during cell division or for the extension of ectoplasmic nets of Lz.Upregulation of these metabolism related contigs along with the upregulation of genes associated with transcription and translation support increased metabolism and/or cellular proliferation.
Surprisingly, no genes potentially linked to the production of zoospores were up-regulated in this study.Other mechanisms of microbial virulence, such as toxin production, iron sequestration, and host immune invasion (Finlay and Falkow, 1997), were also not up-regulated in this study.Together the non-host transcriptome suggests the possible enzymatic degradation of plant cell defenses, supports the rapid dissemination of Lz cells within host tissue, and supports the voracious consumption of plant tissue by invading Lz cells as mechanisms of virulence in seagrass wasting disease.Time course experiments, with lower exposure concentrations (sensu Bower et al., 1989) would be useful in elucidating both the genes and virulence mechanisms involved in penetration vs. pathogenesis during this disease process.

Conclusion
With increased outbreaks of marine disease, continued influence of anthropogenic stressors, and rapidly changing oceans contributing to significant declines in seagrass meadows across the globe, there is an increasing need to understand host-pathogen dynamics across environmental scales.Given the rapid decline of Z. marina globally coupled with active restoration, we have increased need for understanding of host-pathogen dynamics and how this might impact future restoration.This study demonstrates the value of dual host-pathogen transcriptomics for understanding physiological consequences of disease across environmental gradients.While the Lz pathogen exhibited a limited repertoire of virulence approaches, including degradation of cell walls, consumption of intracellular materials and movement within the host, the eelgrass host exhibited a wide range of responses to infection, from cascading phytohormone signals, to altered production of phenols, increased PAMP/DAMP signaling and apoptosis.As our knowledge of Lz increases, biologically realistic inoculation experiments can be used to further understand the role of changing nearshore conditions in facilitating or repressing disease in this important coastal habitat-forming species.

FIGURE 2 Lz
FIGURE 2 Lz cell equivalents (per mg eelgrass tissue) from qPCR analysis of eelgrass leaf tissue 2 days post-inoculation.Bars indicate standard errors.Letters and bars indicate significantly different treatments.

FIGURE 3
FIGURE 3 Prevalence of observed disease signs on eelgrass blades 2 and 5 days post-inoculation.Error bars indicate ± 1 standard error of the replicate mean.Letters and bars indicate significantly different treatments.

FIGURE 4
FIGURE 4 Effects of temperature and Lz treatment on condensed tannins and total phenolic compounds (as a percentage of eelgrass dry mass) 2 and 5 days post-inoculation.Data are means (± 1 SE).Letters and bars indicate significantly different treatments.
FIGURE 5Results of Z. marina WCGNA analysis.(A) Cluster dendrogram of 3,096 genes used in WGCNA analysis.Modules were identified (Dynamic Tree Cut), then merged based on a cut height of 0.34 (Merged dynamic).All modules contain at least 30 genes.(B) Significant correlations between module eigengenes (rows) and Lz infection status (column).The R 2 value is provided for each significant correlation, with P-value in parentheses.Positive correlations are represented by red boxes, while negative correlations are represented by blue boxes, with the intensity of the color increasing with correlation strength.(C) Expression of module eigengenes significantly correlated Lz exposure for exposed and control eelgrass.

TABLE 1
Biological process and molecular function GOterms enriched within significant module eigengenes for Z. marina.