We chose to perform our survey in the area around the island of Hiddensee on the German Baltic Sea coast as it offers variable yet locally representative environments colonized by eelgrass meadows. The area is characterized by semi-enclosed, shallow lagoons (German: “Bodden”) which differ from the open coast with respect to depth distribution, salinity and nutrient concentrations. The surveyed lagoons (the Vitter Bodden and Schaproder Bodden) feature somewhat lower average salinity (8.8 ± 1.0‰, range 6.5–13.4‰) compared to the open coast (Libben bay, 9.4 ± 1.6‰, range 6.8–15.9‰). Nutrient levels are elevated in lagoons (total N: 38.2 ± 12.3, total P: 1.2 ± 0.62 μmol l^-1) compared to open coast waters (total N 19.9 ± 4.3, total P: 0.91 ± 0.4 μmol l^-1) mainly due to agricultural runoff and other human activities (Schiewer, 2008). Salinity and nutrient values are yearly averages from monitoring data 2005 – 2014 (Landesamt für Umwelt, Naturschutz und Geologie Mecklenburg-Vorpommern, unpublished data).

Zostera marina shoots were sampled via scuba diving at seven sites around the island of Hiddensee (Figure 1). At every site, five replicate shoots (above-ground parts) were collected within an area of approximately 1 m². Shoots were kept cool and in clean plastic bags until sample processing which was carried out within a few hours. Environmental variables such as water depth and parameters relating to the surrounding macrophyte vegetation (eelgrass and co-occurring macroalgae and aquatic plants) were determined as part of a simultaneous macrophyte inventory (Bühler, 2016).

The sampled eelgrass shoots were rinsed with sterile filtered seawater and handled under sterile conditions. Leaf areas with heavy macroscopic epibiosis and visibly degrading tissue were removed before biofilm sampling. Epibiotic biofilm from the entire leaf surface of each shoot was scraped off using sterile cotton swabs. The scraped off material was suspended in sterile-filtered seawater (10 ml) and kept cool. The suspensions were aliquoted into separate tubes and centrifuged (10000 × g, 10 min, 4°C) to pellet biofilm material. Pellets were frozen at -20°C until DNA and chlorophyll a extraction. Chlorophyll a was extracted overnight in acetone at +4°C and absorbance at 665 and 750 nm (A₆₆₅ and A₇₅₀) was measured spectrophotometrically (Genesys 20, Thermo Fisher Scientific). Biofilm chlorophyll a content (C_chla) was calculated as C_chla = 11.4 ^∗ (A₆₆₅ – A₇₅₀) and was expressed in μg chlorophyll a cm^-2 of biofilm.

DNA was extracted from biofilm pellets using the MoBio PowerSoil kit for DNA (MOBIO Laboratories). Mechanical lysis was achieved by bead beating in a FastPrep 24 5G (MP Biomedicals). Extracted DNA was amplified with primer pairs targeting the V4 region of the 16S rRNA gene [515f: 5′-GTGYCAGCMGCCGCGGTAA-3′, 806r: 5′-GGACTACNVGGGTWTCTAAT-3′ (Walters et al., 2016)] and the V7 region of the 18S rRNA gene [F-1183mod: 5′-AATTTGACTCAACRCGGG-3′, R-1443mod: 5′-GRGCATCACAGACCTG-3′ (Ray et al., 2016)] coupled to custom adaptor-barcode constructs. PCR amplification and Illumina MiSeq library preparation and sequencing (V3 chemistry) was carried out by LGC Genomics in Berlin. Sequences have been submitted to the NCBI short read archive under bioproject number PRJNA389390, and accession numbers SAMN07197184 – SAMN07197218 (16S rRNA gene sequences) and SAMN07197255 – SAMN07197289 (18S rRNA gene sequences).

Sequences clipped from adaptor and primer sequence remains were processed using the DADA2 package in R (version 1.2.0) (Callahan et al., 2016; R Development Core Team, 2017). Briefly, forward and reverse Illumina reads were truncated to 200 bp, filtered (maxEE = 2, truncQ = 2), dereplicated and error rates were estimated using the maximum possible error estimate from the data as an initial guess. Sample sequences were inferred, and paired forward and reverse reads were merged. Chimeric sequences were removed using the removeBimeraDenovo function. The resulting unique sequence variants (analogous to operational taxonomic units) were used to construct a table containing relative abundances of sequence variants across all samples. Sequence variants were taxonomically classified using a lowest common ancestor approach (Lanzén et al., 2012) based on the Silva database (Pruesse et al., 2007).

Multivariate statistical analyses were carried out using functions in the vegan package [version 2.4-1, (Oksanen et al., 2016)] in R [version 3.3.1, (R Development Core Team, 2017)]. To visualize similarities in sequence variant composition (community composition) between sampling sites, non-metric multidimensional scaling (vegan function metaMDS) was performed on Hellinger-transformed sequence variant counts using Bray–Curtis distance. To explain the variation in community composition in response to the environment, abiotic and biotic variables were selected based on their explanatory power in PERMANOVA tests (vegan function adonis, variables were transformed (log_e) to achieve normal-distribution if needed). A final model including the two most relevant environmental variables was selected and distance-based redundancy analysis ordination (vegan function capscale, Bray–Curtis distance) constrained to these environmental variables was used to visualize variation.

Sequence variant richness was calculated by rarefying read counts to the lowest number of reads in a sample of each dataset (16S: 5861, 18S: 25780). Evenness was calculated as J = H/logS, where H is the Shannon diversity index and S is rarefied richness (Pielou, 1977). To identify sequence variants that were differently abundant in lagoon sites and open coast sites, differential abundance analysis as implemented in the DEseq2 R package (Love et al., 2014).

A partial mantel test was used to assess correlation between prokaryotic and eukaryotic community composition while correcting for environmental variables [vegan function mantel.partial, Bray–Curtis distance for community data, euclidean distance for environmental variables, (Legendre and Legendre, 2012)]. To analyze pairwise correlations between 16S and 18S sequence variants, co-occurrence networks were calculated in R using Pearson correlation. P-values were adjusted for multiple testing using the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). Only correlations with adjusted p-values < 0.01 were treated as significant and included in the network. The network was plotted and edited for clarity (symbols, colors etc.) in Cytoscape version 3.3.0 (Shannon et al., 2003).

This study represents the first effort to characterize both bacterial and eukaryotic epibiotic communities on eelgrass leaves and to relate their communities to potentially important environmental drivers and to each other. We found that the local abiotic environment, in the form of water depth, explained a great deal of variation for both communities. Water depth was one of the main differences between lagoon and open coast sites which likely in part explains the clear separation observed between these sites. However, water depth encompasses separate physical factors such as light penetration and wave shear forces which are difficult to further separate in this study. In addition, several other factors such as mean salinity, nutrient levels and water temperature differ on seasonal time scales in lagoons compared to the open coast (Schiewer, 2008; Blindow et al., 2016) and could also contribute to this variation.

Of the biotic environmental variables measured, eelgrass shoot leaf area (a measure of the size of each shoot) and biofilm chlorophyll a content explained the most variation in community composition of both bacteria and eukaryotes. The explanatory power of eelgrass leaf area, which is related to eelgrass productivity (Echavarría-Heras et al., 2010), suggests an influence of host physiology on epibiotic biofilm composition. The amount of variation explained by chlorophyll a content, a proxy for the contribution of photosynthetic organisms (algae and cyanobacteria) in the biofilms, further indicates that interactions within the biofilms, for example between eukaryotic algae and bacteria play an important role in shaping community composition.

Richness of bacterial sequence variants was correlated with leaf area and with eelgrass dry weight at the sampling sites. Larger leaf area corresponds directly to a larger area of biofilm sampled, since biofilm from entire eelgrass shoots was removed. This observed richness-area relationship agrees with similar biodiversity-area relationships in other microbial habitats (reviewed in Martiny et al., 2006) and among macroscopic organisms. However, an alternative explanation for this pattern relates to the classical notion that energy supply limits biodiversity (Currie, 1991). Since leaf area is a proxy for eelgrass productivity (Echavarría-Heras et al., 2010), its correlation with bacterial diversity could indicate that higher resource supply in the form of for example organic exudates promotes leaf surface bacterial diversity. In contrast, eukaryotic biodiversity (measured as sequence variant richness) seemed to be controlled by factors specific to lagoon vs. open coast habitats, rather than area.

Despite the small geographic area surveyed in this study, variability in community composition of both bacteria and eukaryotes was substantial among sites, which shows that eelgrass leaf biofilms are dynamic habitats on a regional scale. This observation is consistent with, and complementary to, recent results that revealed that eelgrass leaf surface microbiomes were more variable than root-associated microbiomes on a global scale (Fahimipour et al., 2017), as well as within an individual eelgrass patch (Ettinger et al., 2017). Nevertheless, we detected several sequence variants in all samples analyzed, supporting the possible existence of a core microbiome of eelgrass leaves. Several core sequence variants were classified to the family Rhodobacteraceae (Alphaproteobacteria), which are often surface-attached, chemoorganotrophic bacteria common in the marine environment (Dang and Lovell, 2002). Interestingly, this family was also consistently detected in the rhizospheres of eelgrass and its close relative Z. noltii (Cúcio et al., 2016; Ettinger et al., 2017) as well as on leaves and below-ground parts of the seagrass Halophila stipulacea in the red sea (Mejia et al., 2016) indicating that their members indeed play important roles in association with various seagrasses. The observation of a Methylophilaceae (Betaproteobacteria) sequence variant as an abundant core community member agrees with detection of this taxon on eelgrass and other co-occurring macrophytes in the brackish-water Chesapeake Bay (Crump and Koch, 2008) and suggests importance of one-carbon metabolism on eelgrass leaf surfaces (Lapidus et al., 2011).

The strong correlation between the bacterial and eukaryotic communities observed in this study has at least two possible causes. First, the abiotic and biotic environment experienced by both communities could influence them in similar ways. This explanation is partially supported by the similar amounts of variation explained by water depth, leaf area and chlorophyll a content for both communities. Second, bacterial and eukaryotic taxa may interact directly, resulting in a bacterial community that is driven by the composition of the eukaryotic community, or vice versa. The topology of the co-occurrence network (Figure 5) suggests some importance of the latter scenario, where certain eukaryotic taxa (e.g., classified as Hydrozoa, Mytilus sp. and Cocconeidaceae) feature a disproportionate number of significant correlations with bacterial taxa. This indicates that bacterial communities in mature biofilms on eelgrass leaves are to some extent a reflection of the combined microbiomes of several epibiotic eukaryotes such as hydrozoans, mussels and diatoms in addition to any eelgrass-specific microbiome that may exist. Conversely, some bacterial taxa correlate with several eukaryotic taxa. For example, a cyanobacterial sequence variant classified as Rivularia sp. correlates with several ciliate taxa which may for example be attracted to this cyanobacterium as grazers. In a previous study, bacterial and eukaryotic communities of the tropical seagrass Enhalus acoroides were found to respond similarly to environmental CO₂ exposure, suggesting that seagrass bacterial and eukaryotic epibiotic communities are highly correlated also in other seagrass species and environments. However, potential direct interactions between bacterial and eukaryotic taxa were not addressed in this study (Hassenrück et al., 2015).

To date, many studies addressing the microbiomes of marine organisms have focused exclusively on prokaryote communities (see e.g., Hentschel et al., 2012; Wahl et al., 2012; Bourne et al., 2016 and references therein). In contrast, free-living microbial eukaryotes are receiving increasing attention and molecular community analysis has led to several recent discoveries that emphasize their underexplored diversity and crucial roles in the environment (Liu et al., 2009; Lima-Mendez et al., 2015). It is likely that insights of a similar caliber can be gained by applying molecular tools to eukaryotes associated to macroscopic hosts (Andersen et al., 2013), beyond relatively well-described eukaryotic symbionts such as zooxanthella in corals (Muller-Parker et al., 2015). We have shown that the bacterial microbiome is strongly correlated with epibiotic eukaryote communities, including microbial and macroscopic taxa, indicating that interactions between eukaryotes and bacteria drive community assembly on eelgrass surfaces to some extent. By ignoring eukaryotes in bacterial microbiome studies, a large part of the variation in bacterial communities could remain unexplained or misinterpreted. Indeed, the high diversity of both eukaryotic primary producers such as diatoms and heterotrophic taxa such as ciliates, rhizarians and many macroscopic and microscopic animals paint a picture of complex microbial ecosystems encompassing several trophic levels on eelgrass leaves. In these “microbial jungles,” bacteria are bound to play important roles which may be especially relevant in interaction with their eukaryote neighbors.

In addition to considering the understudied eukaryotes in microbiome studies, gaining a mechanistic understanding into microbiome assembly requires means to infer causative relationships between microbiome dynamics and abiotic and biotic factors (Widder et al., 2016). The present study included a limited number of samples from a small geographic area, providing a valuable first assessment of environmental drivers of community assembly and potential microbial interactions. Further insight could be gained by manipulative experiments, whereby for example, seagrass surfaces are exposed to colonization by microorganisms under controlled conditions. In addition, large observational microbiome datasets from a range of eelgrass meadows encompassing wider environmental gradients would allow more statistical power to tease apart different drivers of community assembly. Modeling approaches such as structural equation modeling (Pugesek et al., 2003) could allow assessment of causative relationships between the environment, the seagrass host and its microbiome. A major challenge lies in linking microbiome dynamics to processes that are relevant on the ecosystem level. Manipulative experiments that operate on an ecosystem scale, i.e., mesocosm experiments mimicking seagrass meadows, offer a promising approach to manipulate the environment while tracking microbiome dynamics and their links to critical ecosystem fluxes such as primary production and respiration. Ultimately, understanding assembly and function of seagrass microbiomes could contribute to tackling some of the challenges faced by seagrass meadow ecosystems under global climate change and other anthropogenic stressors.