Teasing apart the host-related, nutrient-related and temperature-related effects shaping the phenology and microbiome of the tropical seagrass Halophila stipulacea

Halophila stipulacea seagrass meadows are an ecologically important and threatened component of the ecosystem in the Gulf of Aqaba. Recent studies have demonstrated correlated geographic patterns for leaf epiphytic community composition and leaf morphology, also coinciding with different levels of water turbidity and nutrient concentrations. Based on these observations, workers have suggested an environmental microbial fingerprint, which may reflect various environmental stress factors seagrasses have experienced, and may add a holobiont level of plasticity to seagrasses, assisting their acclimation to changing environments and through range expansion. However, it is difficult to tease apart environmental effects from host-diversity dependent effects, which have covaried in field studies, although this is required in order to establish that differences in microbial community compositions among sites are driven by environmental conditions rather than by features governed by the host. In this study we carried out a mesocosm experiment, in which we studied the effects of warming and nutrient stress on the composition of epiphytic bacterial communities and on some phenological traits. We studied H. stipulacea collected from two different meadows in the Gulf of Aqaba, representing differences in the host and the environment alike. We found that the source site from which seagrasses were collected was the major factor governing seagrass phenology, although heat increased shoot mortality and nutrient loading delayed new shoot emergence. Bacterial diversity, however, mostly depended on the environmental conditions. The most prominent pattern was the increase in Rhodobacteraceae under nutrient stress without heat stress, along with an increase in Microtrichaceae. Together, the two taxa have the potential to maintain nitrate reduction followed by an anammox process, which can together buffer the increase in nutrient concentrations across the leaf surface. Our results thus corroborate the existence of environmental microbial fingerprints, which are independent from the host diversity, and support the notion of a holobiont level plasticity, both important to understand and monitor H. stipulacea ecology under the changing climate.


Introduction
Halophila stipulacea (Forsk) Ascherson is a small tropical seagrass, dominant in the Gulf of Aqaba [1][2][3], the northernmost edge of its natural range [4,5]. Here H. stipulacea forms large discontinuous meadows along a wide range of depths [1-3, 6, 7]. Ecosystem services and functions associated with H. stipulacea meadows in this region are numerous. H. stipulacea are attributed with high primary productivity, enriching the water with oxygen. They sequester blue carbon, partially mitigating ocean acidification for neighbouring reefs [8]. They also reduce pathogen loads in seawater [9] and uptake nutrients, improving water clarity for neighbouring ecosystems [reviewed in 10]. Importantly, they provide major nursery grounds for fish, crustaceans, gastropods and bivalves [1]. However, these meadows are severely threatened by anthropogenic pressures. Seawater in the Gulf of Aqaba has been warming 50% faster than the global mean coastal sea surface temperature trend of 0.17 ± 0.11 °C per decade [11,12]. In addition, they are faced with intensive coastal development (e.g., the Saraya laguna in Aqaba, and the Splash Park in Eilat). These [3threats to seagrasses and the ecological functions they provide are made worse by the long water residence time estimates in the Gulf [3-8 years; 13] caused by its semi-enclosed basin shape, and the effects of eutrophication are intensified several-fold [reviewed by 2].
With the ongoing decline of seagrasses worldwide [14] alongside the relatively slow responding community based indicators in most seagrass monitoring efforts [15][16][17], there is a growing need for fast and responsive indicators to changes in the ecological status of seagrasses [18]. In a study of three Gulf-of-Aqaba H. stipulacea meadows growing at a shared depth, Mejia et al. [19] was able to distinguish bacterial "environmental fingerprints" for the different sites, from a "core microbiome", which was shared among all sites. Bacterial differences among sites were attributed mostly to Rhodobacteraceae, which dominated particularly at the southernmost meadow, the site with highest nitrogen and phosphate concentration measurements, but with lowest turbidity. Additionally, plants from this meadow had the smallest leaf area, among the studied meadows. The existence of "environmental fingerprints" alongside a core microbiome was also corroborated by Rotini et al. [7], who compared the microbiomes in H. stipulacea meadows along a depth gradient of 4-28 m. The epiphytic leaf microbial communities of these meadows were found to be dominated by Gammaproteobacteria and Bacteroidetes in the light limiting conditions of the deeper sites, while Cyanobacteria and Rhodobacteraceae thrived in conditions of high light availability and hydrodynamics, in shallower sites. Based on the "environmental fingerprint" hypothesis, bacterial community shifts, and particularly Rhodobacteraceae relative abundance increase, may provide early indications to nutrient exposure [7,19,20], and possibly act as a buffer, protecting the seagrass host from excess nutrients through bacterial metabolism. However, field studies cannot tease apart the effects of environmental factors from those of possible differences in host diversity among seagrass meadows.
Outside of its natural range, H. stipulacea is a successful lessepsian migrant in the eastern Mediterranean Sea [21][22][23][24][25] and an invasive species in the Caribbean Sea [26,27]. The ability to migrate into new habitats could be related not only to its extensive morphological, biochemical and physiological plasticity [reviewed by 2,6,7,[28][29][30], but perhaps also to its bacterial "environmental fingerprints". Like many other organisms, seagrasses were suggested to establish symbiotic relationships with their associated microbial communities to form a functional unit (the "holobiont") that reacts as a whole to environmental changes [31][32][33]. It has been suggested that the epiphytic bacterial community of the seagrass leaves, the focus of this study, benefit from the organic carbon enriched microhabitat on the leaves [34] and comprises aerobic organotrophic bacterial species able to utilize the secreted polymers [7,34]. Diazotrophic bacteria such as cyanobacteria, which enhance nitrogen availability [35][36][37][38], may also be a part of this community, depending on nutrient concentrations in the water [39]. Alternatively, leaf cyanobacterial biofilms may reduce light availability [40].
In this study we aimed to understand whether the bacterial "environmental fingerprint" in H. stipulacea existed independently from host related diversity, which may exist among meadows in the field, in order to establish the utility of bacterial shifts as early warning signs for environmental stress. This would also provide evidence for holobiont level plasticity, supporting the rapid range extension of H. stipulacea. We also tested the interaction between phenotypic responses of H. stipulacea, including shoot production and mortality, with microbial community shifts, under temperature and nutrient stress. We carried out a mesocosm experiment in which seagrass phenology and leaf epiphyte community compositions Keywords: Seagrass, Halophila stipulacea, Phenology, Microbiome, Mesocosm experiment, Heatwave, Nutrient loading were the dependent variables, while the source seagrass meadow, water temperature and nutrient concentrations were the independent variables. We simulated shallow meadow conditions of 8-10 m depth, in which light was not a limiting factor. We hypothesised that both the environmental conditions and the source seagrass meadow would partially explain the variance in the dependent variables, and that the interaction of the independent factors would affect the seagrass performance and the microbiome differently than each factor separately.

Results
To study the effects of heat, nutrients and the source H. stipulacea meadow on the phenology and epiphytic microbiome of H. stipulacea in a controlled environment, we carried out a mesocosm experiment. We controlled for the source site by including seagrasses from two sites, Tur-Yam Beach (TY) and South Beach (SB), both located in the northern Gulf of Aqaba (Eilat, Israel). TY is offshore a small and active marina, close to a rarely used crude oil terminal, whereas SB is removed from any obvious anthropogenic pressures. While turbidity in the TY site is slightly higher than in the SB site, higher nutrient concentrations were measured in SB (further detailed in Mejia et al. 2016). To reduce the effect of microbial legacies, seagrasses were washed with freshwater and re-inoculated by exposure to natural seawater. The mesocosm aquaria system ( Fig. 1) included four temperature baths, each containing five aquaria filled with artificial seawater. Following acclimation, under 27 °C and no nutrient enrichment, the mesocosm conditions diverged to four different regimes, for the duration of the 40 days experimental phase: (i) control water temperature with no nutrient enrichment (CTCN), (ii) control water temperature with nutrient enrichment (CTN), (iii) increased water temperatures (31 °C) with no nutrient enrichment (TCN), and iv) the combination of the two stressors-increased water temperatures (31 °C) with nutrient enrichment (TN). The experimental phase was followed by a recovery step in which baseline temperatures were set in all the baths and nutrient addition was stopped.

Nutrient measurements, shoot formation and shoot mortality
The establishment of heatwave conditions (31 °C) in TCN and TN aquaria during the stress phase of the experiment (T1-T5) is demonstrated in Fig. 2A 4 3− concentrations ( Fig. 2E) were similar in all treatments at T0   (CTN, TN) treatments, respectively (0.0187 < p-value < 0.068), as expected. PO 4 3− concentrations at T8 were very high, albeit similar among treatments (3.21-4.40 µM on average), most likely due the the accumulation of PO 4 3− in all the aquaria. Shoot mortality appeared among TY plants only at T5, and only in heated aquaria (TCN, TN), in which a total of 15 dead shoots were counted (Fig. 2G). Shoot mortality in SB plants, however (Fig. 2F), was evident two weeks earlier at T3 with a total of 29 dead shoots across treatments and time points, including five from unheated aquaria (CTCN, CTN). Correspondingly, new shoot emergence ( Fig. 2H and I) was significantly higher in TY plants than SB plants (F = 16, p-value = 8 × 10 -5 ) and in non-enriched aquaria (CTCN, TCN) than in enriched aquaria (CTN,TN) with F = 4.6 and p-value = 0.003, in a type-3 factorial ANOVA. Time was a marginally significant factor (F = 2.3, p-value = 0.02). Therefore, according to phenology measurements, TY seagrass inherently performed better under thermal stress for both shoot death and emergence, temperature stress primarily promoted shoot death, more so in SB hosts more than in TY hosts, and nutrient loads suppressed the emergence of new shoots for both SB and TY hosts, down to a similar level.

Sequencing results
To study the temporal dynamics of the leaf epiphytic microbiota, we carried out a metabarcoding experiment. Following the exclusion of organelle sequences and chimeric sequences, as well as amplicon sequence variants (ASVs) with a frequency lower than 30 across all samples, the analysis included 142 mesocosm epiphyte samples, 20 mesocosm water samples and 10 field epiphyte samples, with a median read-pair count of 14,501 (i.e., 29,002 individual reads). Except for one of the retained mesocosm epiphyte samples, which had 5,851 reads, the filtered read count ranged from 8,312 to 33,256.

Bacterial diversity
Samples were rarified to 5000 read-pairs following the guidance of a rarefaction curve, which confirmed that most rare taxa were represented (Additional file 1: Fig.  S1), rarefaction curves were produced using filtered matrices. Water samples had significantly lower Faith PD values than epiphyte samples (ANOVA p-value < 10 -56 ; Additional file 2: Fig. S2). The first principal coordinate in a principal coordinate analysis (PCoA; Additional file 2: Fig. S2), explained 38.4% of the total variance and segregated the sample types, water vs. epiphytes, to distinct clusters. Water samples were more divergent from one another than epiphyte samples, and the main ASVs explaining the difference between the sample types belonged to the family Alteromonadaceae (ASV cc2916), the NS3a marine cluster (745c89) and the genus Winogradskyella (01f97d), which were more abundant in the water samples (Additional file 2: Fig. S2).
To study the bacterial diversity within the epiphytic communities, we excluded ASVs that were also found in the aquarium-water samples and regarded them as contamination. Ratification of the epiphytic samples was done following this step. Faith PD values were similar among time points and between the source sites (SB and TY). Treatments were a significant factor (ANOVA p-value = 0.01), due to significant differences between the TN and CTN treatments (q-value = 0.005; Fig. 3), both of which were enriched at T3 and T5, differing only in the water temperatures applied. The decrease in alpha diversity under CTN corresponded with the increase in NO 2 − under the same treatment (Fig. 2B).
The unweighted UniFrac PCoA analysis ( Fig. 4) revealed an increase in the divergence among samples from different treatments, following T0. As time progressed, the second axis (8-9% of the total variance) segregated the control samples (CTCN) from the rest of the treatments, and axis 1 (12-13%) separated CTN from the other treatments. The first and second axis in the Weighted UniFrac analysis (Additional file 3: Fig.  S3) explained a larger cumulative portion of the total variance (25-30%) with similar patterns but shorter distances among the treatments. According to the redundancy analysis ANOVA, a model accounting for time, , accounted for 21% of the total variance (p-value = 0.001), with small but significant effects of time (3.3%, q-value = 0.003), temperature (2.9%, q-value = 0.003), PO 4 3− concentration (2%, q-value = 0.007), (1.8%, q-value = 0.008) and NH 4 + (1.2%, q-value = 0.038). The source site had a borderline significant effect (q-value = 0.09) explaining 1% of the variance. The remaining 7.5% of the variance explained by the model can be attributed to interactions among these factors.
To further understand the importance of the source site in determining the response of biofilm communities to the treatment, we carried out a factorial ANOVA test using UniFrac pairwise distances. We tested whether within-SB, within-TY, and between-sites pairwise distances were significantly different (Fig. 5). We took care to only include within-treatment and within-time point pairwise distances and to test whether the selected distances differed among the treatments.
At T0, pairwise distances within TY were smaller than within SB or between source sites when considering weighted distances (p-value < 0.012) but not when using unweighted distances. In intermediate time points (Fig. 5C and G), particularly in T5, pairwise distance increased both within and between source sites, with a greater increase in the within-TY distances, (q-value = 0.047 and q-value = 0.024 for weighted and unweighted distances, respectively). Following recovery ( Fig. 5D and H), the within-TY distances increased further compared with the within SB distances (q-value = 0.005 for weighted distances) or the amongsite distances (q-value = 0.005 and q-value = 5 × 10 -4 ). It would therefore appear that the source site of the seagrass was one of the factors shaping the bacterial composition of the seagrass biofilm. Specifically, Time points T0 and T8 had baseline temperatures and no active nutrient enrichment in all baths. The percent total variance accounted for by each coordinate is indicated on the corresponding axis. The most important ASVs, following the importance definition by Legendre and Legendre [41], are represented by BiPlot analyses (gray arrows) and their taxonomic identifications are noted in the legend. A similar analysis, without exclusion of ASVs that were shared with water samples, yielded similar results (Additional file 12: File S1) TY epiphytic communities have diverged from each other at higher rates than SB epiphytic communities. This result coincides with the increased mortality and reduced generation of shoots that was observed in SB plants compared with TY plants (Fig. 2F-I; mentioned above).

Order dynamics
Relative abundances of the high abundance orders were similar among treatments and time points (Fig. 6), except for a few but important exceptions. The most abundant order (27-40%), Rhodobacterales (Alphaproteobacteria; Fig. 6M-P), had similar relative abundances in all treatments throughout the experiment, except for T8. By that time point, in which similar recovery conditions were already applied for three weeks in all the treatments, their relative abundance significantly increased (q-value < 0.009) under the legacy of control temperature conditions (27 °C) and particularly under the legacy of control temperatures with high nutrients concentration (TCN), in correspondence with the increased NO 2 − concentrations (Fig. 2). Flavobacteriales (Fg. 6I-L, Bacteroidetes; Flavobacteriia) diverged among the treatments, with increased relative abundances under nutrient enriching (CTN & TN) at time point T3 (q-value = 0.025) and under control temperatures (CTCN & CTN) at time point 5 (q-value = 0.048), but the relative abundance in the different treatments converged under recovery conditions at T8. Mean relative abundances of Flavobacteriales ranged from 6 to 14% across time points and treatments. Similarly, Rhizobiales ( Fig. 6Q-T, Alphaproteobacteria) diverged among treatments at time points T3 and T5 (q-value < 0.02), to converge at T8, after three weeks of recovery conditions. However, Rhizobiales consistently flourished only under heatwave conditions without nutrient loading (TCN). Their mean relative abundances ranged from 1.5% to 9% across time points and treatments. SBR1031-clade bacteria ( Fig. 6E-H, Chloroflexi, Anaerolineae), with mean relative abundances of 2-13%, had a higher relative abundance under control conditions (CTCN) at T3 (q-value < 0.03), but this increase did not persist, and this order was not among the eight most abundant orders following recovery. Cyanobacterales ( Fig. 6A-D) emerged as a high relative-abundance order (2-8%) at T5. They had higher relative-abundance under baseline temperatures with nutrient enriching (CTN) in T5 (q-value < 0.02).
ANCOM tests revealed additional differentially abundant orders, with low relative abundances (Additional file 4: Fig. S4). Alphaproteobacteria incertae sedis (an artificial group of several alphaproteobacteria genera with unclear placement) were more abundant in the control samples (CTCN) than in any of the treatment samples. Alteromonadales (Alteromonadaceae in particular, Additional file 4: Fig. S4), decreased with time under all the experimental and control regimes. BD7-11 (phylum Planctomycetota; Additional file 4: Fig. S4) and Myxococcales (Deltaproteobacteria) followed a similar pattern to one another, where, unlike Microtrichales and Rhodobacterales, they were suppressed under CTN. Caenarcaniphilales increased under stress compared to the T0 abundance in each treatment (CTN, TCN, TN), and flourished after recovery from heatwave and nutrient loading conditions (TN), JTB23 (Gammaproteobacteria;  Fig. S4), which flourished under nutrient loading with baseline temperatures (CTN), similarly to Rhodobacterales, the OM182 clade (Gammaproteobacteria), which increased under nutrient loading conditions (CTN, TN) and Sphingobacteriales (particularly the NS11-12 marine group, Additional file 4: Fig. S4), which flourished during recovery from heatwave and baseline nutrients conditions (TCN). Additional orders differed among treatments but were consistent with their relative abundance at T0, prior to the stress phase of the experiment.
A family level ANCOM analysis (Additional file 4: Fig. S4) largely mirrored the order level results, with the addition of the following families that did not belong to the above mentioned orders. Microscillaceae (order Cytophagales), which developed with time only under control conditions, Phormidiaceae (order Oscillatoriales), which developed best under nutrient loading conditions (CTN, TN) and flourished following recovery from the combined stress TN conditions and Rubinisphaeraceae (Planctomycetales), that increased under control temperature combined with nutrient enrichment conditions (CTN) and persisted even during the recovery.

Dynamics of key ASVs
The increase in the relative abundance of ASV 71a746 (Rhodobacterales), which is highly explanatory of the total variance (Fig. 4) seems to be responsible for the observed increase in the relative abundance of Rhodobacterales, which corresponds with the decrease in Faith PD observed under CTN (Fig. 3). ASV 71a746 is most closely related to AM691091, an extremophile from the East German Creek system in Canada [42] and to several Roseobacter sequences (Additional file 5: Fig. S5), a genus containing temperate and polar species [43], and was first described in seaweed [44]. Some Roseobacter reduce NO 3 − to NO 2 − [44]. The Microtrichales ASV 904c34 consistently developed only under the CTN conditions as well. Microtrichaceae are found in systems in which NH 4 + is oxidized with NO 2 − via the anammox system [45]. Methyloligellaceae ee5ce9 (Rhizobiales) are methanotrophs [46], which flourished most under TCN. Dynamic cyanobacterial ASVs included Rivularia 8c13b7, particularly successful under CTN, which appears to be related to extremophilic cyanobacteria from alkaline, saline and thermal environments (Additional file 6: Fig. S6). It was the most abundant cyanobacterial ASV, but it was displaced by T8. Pseudanabaena 2f6e75, another cyanobacterium, is closely related to an isolate from the sponge Axinella damicornis [KY744814 ; 47]. It emerged at T8 in all treatments, including CTN. Lyngbya bc8953, which is closely related to cyanobacterial isolates from the intestinal tract of herbivorous marine fish (HM630185) and black band disease coral tissue [DQ446127; 48], also emerged at T8, under TN. Saprospiraceae a2fde72, which is closely related to isolates from the surface of macroalgae (Additional file 7: Fig. S7; DQ269042) prevailed under all treatments, except for TN. Under the TN treatment, this ASV was displaced by another Saprospiraceae ASV (9ecde7), but reemerged following recovery. Interestingly, ASV 9ecde7 is related to isolates from a Guerrero Negro hypersaline microbial mat [49]. In the SBR1031 group, which was most successful under control conditions, ASV A4b f35075 appeared to be an exception. It emerged late under TN. The gammaproteobacterium Granulosicoccus ASV 9d7587, flourished only under CTN, similarly to Microtrichales and Rhodobacterales. This ASV emerged as important in Additional file 3: Fig. S3. Granulosicoccus was first identified in Antarctica, has low optimal temperatures, and can reduce nitrate [50]. The temporal dynamics of key ASVs and taxa are summarized in Additional file 8: Fig. S8.

Discussion
Studies of H. stipulacea meadows in the Gulf of Aqaba have revealed associations between the epiphytic aboveground microbiota composition and water nutrient concentrations [7,19,20]. Based on the observed associations, these studies have suggested that microbiota function within the seagrass holobiont and have pointed out their potential use as an ecological indicator of exposure to environmental stress by the seagrass host [20]. While in-situ studies might provide a better representation of realistic conditions to which seagrasses and their microbiota are exposed, it is often difficult to tease apart the relative importance of covarying factors, such as heat waves, nutrient loading and the host diversity. In addition, in field experiments, there is the possibility that legacy environmental conditions, rather than those measured during the experiment, are responsible for the observed bacterial community compositions [51]. This may lead to misinterpretation of the key factors, and of the meaning of bacterial dynamics as bioindicators of environmental stress. In this mesocosm experiment we set out to control for the temperature, nutrient concentrations and source meadow, to study their respective importance and to test the existence of "environmental fingerprint" properties previous workers attributed to seagrass microbiota [7,19,20]. A fundamental requirement for such a marker is that bacteria that were indicated as key players in-situ are also represented in the experimental system at the end of the acclimation step (T0), so that the observed consequent dynamics will bear relevance to the natural meadows. This equivalence between the experimental system and previous in-situ results are best reflected by the core ASVs. The mesocosm study system has 107 core ASVs shared between the natural (in situ) samples and the experimental samples, including ASVs belonging to taxa which responded to the experiment or that were specifically shown to respond to the experiment. The largest cohorts of core ASVs belonged to alpha-proteobacteria (48 ASVs, mainly Rhodobacterales) and Bacteroidia (25 ASVs, mostly Flavobacteriales), both key groups of the aboveground in situ H. stipulacea microbiota [19]. As a case in point, the influence of the core ASV 71a746 on the alpha diversity in the mesocosm was paramount. Two groups of dependent variables were quantified in this study, phenological and microbial. The two groups had fundamentally different governing factors. Phenological properties were largely dependent on the seagrass collection site, with overall better performance of the H. stipulacea from the TY. The two phenological traits, shoot mortality and shoot production, reacted differently to the two stress types we have simulated. Shoot death was accelerated by high temperatures, while the emergence of new shoots was slowed down by nutrient loading. Hence, both stressors are detrimental to the seagrass development, but operate through different mechanisms. The environment in SB is considered to be less disturbed than the environment in TY. However, Mejia et al. [19] measured higher nutrient concentrations in SB water and pore water than in TY water. It is therefore challenging to determine whether environmental conditions have driven TY seagrass to higher resilience, SB seagrass to higher fragility, or if phenological differences between the sites were at all shaped by environmental factors. Indeed, our experiment does not allow us to distinguish among several possible mechanisms such as plasticity, selection, gene flow or drift, which may have been responsible for the phenological differences between the two seagrass source-sites, because there is no genetic population structure data.
Our results also reveal a complex effect of temperature, nutrient concentrations and time on the epiphytic microbial community associated with H. stipulacea. Only 20% of the total compositional variance was attributed to factors which were accounted for or to interactions among them. Possibly, additional important independent factors, which were not accounted for, might have explained the remaining variance. For example, variability in leaf exudates within a seagrass population would likely be reflected in the microbial composition, and sediment microbiome dynamics may have directly modulated the nutrient measurements. Another possibility is that biotic interactions within the epiphytic communities were significant sources of variance. Depending on the initial epiphytic community compositions, such interactions would have had different consequences under similar conditions, particularly when taking the effect of ecological drift into account [52]. Under strong ecological drift, the resolution of abiotic perturbations in a microbial community is partially stochastic, especially if different taxa share functions and have equivalent fitnesses. The small size of the mesocosm, in comparison to the size of natural meadows and their environment, could have contributed to an increased ecological drift effect.
With the exception of time, the temperature had the largest effect on the explained variance of the bacterial composition, followed by the effects of the nutrient concentrations. The source site had a small and borderline significant effect. However, the unifrac distance distributions within and among source sites revealed that the TY communities, initially more conserved than the SB communities, diverged to a larger extent from one another than the SB epiphytic microbial communities, by the end of the experiment. This result coincided with lower mortality rates and higher growth rates in the TY seagrass and should thus be evaluated as a possible source of beneficial holobiont level plasticity of the TY seagrass in future studies. Still, the source site of the seagrass had a minor effect on the microbial community composition, in contrast with its effect on seagrass phenology, and the microbiota were mostly shaped by environmental factors. These results strongly support the "environmental fingerprint" hypothesis formulated by Mejia et al. [19], Rotini et al. [7] and Conte et al. [20] and highlight the early warning information that can be gained regarding exposure to stress, from monitoring the microbiota in wild meadows.
In terms of alpha diversity, only the combination of control temperatures (27 °C) with enriched nutrients (the CTN treatment) caused a reduction in Faith PD, in comparison with the control treatment. The reduction in alpha diversity was explained mainly by the increase in Rhodobacterales and particularly by ASV 71a746. Based on a phylogenetic analysis, ASV 71a746 is most closely related to genbank accession AM691091, an extremophile from the East German Creek system in Canada [42], within Roseobacter, a genus normally containing temperate and polar species [43], and was first described in seaweed [44]. It is therefore possible that this ASV represents the southern edge of its congeneric distribution, can only flourish under the baseline temperature of 27 °C, but effectively utilizes the excess nutrients. Interestingly, some Roseobacter spp. reduce NO 3 − to NO 2 − [44], which is consistent with the enriched NO 2 − measurements under the CTN scenario. Granulosicoccus ASV 9d7587, which is also related to cold water nitrate reducers, followed a similar pattern to that of the Rhodobacterales ASV 71a746. Concomitantly, Microtrchales, and Microtrichaceae in particular, also prevailed particularly under the CTN scenario. Microtrichaceae were found to be highly abundant in a partial nitrification-anammox system, where partial nitrification, such as that carried out by Roseobacter, produces NO 2 − , which is then used by Microtrichaceae to oxidize NH 4 + [45]. Therefore, attempts to isolate these microorganisms and experimentally study their contribution to N-cycles may shed light on the buffering that microorganisms can provide against nitrogen enriched environments. Interestingly, the highest seagrass mortality was observed under high temperature treatments, and almost never under nutrient enrichment alone (CTN). The presence of these bacteria on seagrass leaves, and their increase under high nutrient concentrations with baseline temperatures, may represent a resilience mechanism of H. stipulacea to nitrogen-enriched environments, possibly allowing them to outcompete other seagrass species in anthropogenically disturbed, or less oligotrophic areas than the Gulf of Aqaba. For example, it would be interesting to test whether similar dynamics are absent from Cymodocea nodosa, which is rapidly displaced by H. stipulacea in the Mediterranean Sea, particularly in a disturbed harbour on the Tunisian coast [24]. Previous studies have reported that Red Sea H. stipulacea grew faster, forming denser meadows in polluted areas [29,53], which is consistent with the low effect of nutrient loading on mortality that we report here. Such a buffering mechanism, however, would be sensitive to the increase in seawater temperatures, according to our results.
The experiment also sheds light on the response of cyanobacterial biofilms to each of the stressors. The relative abundance of Cyanobacteria increased by T8 under all treatments, even weeks after active loading ceased, probably due to the accumulation of nutrients in the aquaria. This was particularly true under TN. Cyanobacteria bloomed earliest under the CTN treatment, where temperatures were kept at baseline (27 °C) but nutrients were actively enriched. Rivularia spp., to which the early blooming ASV belonged, was the most abundant cyanobacterial ASV. Rivularia spp. was already shown to be sensitive to nutrient loading [54], and accordingly this ASV was displaced by an ASV belonging to Pseudanabaena sp., which emerged at T8 under all treatments. A Lyngbia sp. ASV which also emerged at T8 under TN is known to form biofilms on seagrass leaves and reduce light availability [40]. Interestingly, the relative abundance of Cyanobacteria was a major contributor to the difference between TY and SB in situ [19]. With higher abundances in the SB meadow, where nutrient concentrations were higher. Lastly, cyanobacteria, although considered to either increase nitrogen assimilation [35] and limit light availability [40], their increased abundance under CTN does not seem to be detrimental to the seagrass in the mesocosm. However, light was not a limiting factor in our experiment, and cyanobacteria fix nitrogen only when it is not otherwise available [55] so in fact, they may participate in the buffering of nutrient loading.

Conclusion
In this study, we were able to tease apart the impacts of environment dependent and host dependent factors on phenological and microbial properties of H. stipulacea, illustrating that the measured phenological properties are mostly host dependent while the epiphytic microbial composition is mostly environment dependent. This result supports the "environmental fingerprint" hypothesis raised in recent studies, and highlights the utility of microbiome shifts as bioindicators of nutrient exposure changes in H. stipulacea meadows. We further propose that bacterial community dynamics contribute to the holobiont's plasticity by buffering nutrient effects when high concentrations are encountered, which may facilitate the range extension of H. stipulacea into new habitats, including northerner latitudes with less oligotrophic waters than in their native range. Our experiment demonstrated that when exposed to either stressor, plants from the TY meadow, which is a site with medium anthropogenic impacts, performed better than plants from the SB meadow, which is a less disturbed site. However, understanding why this differentiation occurred is beyond the scope of our experiment. Our findings may have important implications concerning the future of H. stipulacea as global climate changes progress and the importance of the holobiont perspective in understanding them.

Plant collection
Intact and healthy H. stipulacea plants bearing [5][6] shoots were collected in July 2019 from 6-8 m depth (Irradiances of 250 μmol photons m −2 S −1 ) by scuba diving. Two meadows in the northern Gulf of Aqaba (Eilat, Israel) were visited: South Beach (SB; 29.497664°N 34.912737°E) and Tur-Yam Beach (TY; 29.516527°N 34.927205°E). At each site, 140 plants were collected from every 5-10 m to avoid pseudoreplicates. The plants were put in ziplock bags filled with seawater and were transported in a cooler box to the seagrass mesocosm facility at the Dead Sea Arava Science Center, where they were immersed in freshwater for 2-3 min and wiped off to remove organisms as much as possible with minimal damage to the plants. Six additional H. stipulacea were collected from each meadow on Dec 15 2019 and stored at − 80 °C until further processing to compare the post-acclimation mesocosm epiphytic communities (see below) with epiphytic communities in the field. As the field and mesocosm samples were not sampled at the same time, we used the field samples only to examine whether mesocosm 16S rRNA variants also occur in the field. The field samples were not used for diversity or differential abundance analyses.

Experimental setup
A mesocosm aquaria system (Fig. 1) was set up to simulate heatwave and eutrophication conditions, alongside control conditions. The system included four temperature baths, each controlled with a ProfiLux aquarium controller (GHL aquarium computers, Germany). Each bath contained five 60 L aquaria, layered with 6 cm of sieved and autoclaved natural coastal sediment, and filled with artificial seawater at a 40 practical salinity units concentration of Red Sea salt (www. redse afish. com). In each aquarium we planted 16-18 plants from each site (TY and SB) side by side, divided by a barrier between the two source sites. 100 ml of seawater from each site (TY and SB) were added to each aquarium once a week during the 4 weeks acclimation phase. Plants were acclimated for four weeks to baseline mesocosm conditions (27 °C, 250 μmol photons m −2 S −1 at the water surface, 12 h light/day), allowing them to recover from the potential stress inflicted during plant collection, transportation and transplantation. Seagrasses were then exposed to four different treatments in each of the four baths: i) control water temperature with no nutrient enrichment (CTCN), ii) control water temperature with nutrient enrichment (CTN), iii) increased water temperatures (31 °C) with no nutrient enrichment (TCN), and iv) the combination of the two stressors-increased water temperatures (31 °C) with nutrient enrichment (TN). To initiate these experimental conditions at the end of the 4-weeks acclimation period (T 0 ), using aquaria heaters, the temperature in two baths was gradually increased (0.8 °C/day) from 27 to 31 °C, ~ 4 °C above the average summer temperatures, simulating the increased temperature of the Red Sea at the end of the century [12,56,57]. The rest of the aquaria were left at control water temperatures of 27 °C. Similarly, nutrients were gradually loaded in the aquaria of two baths, simulating eutrophication, by adding crushed slow-release fertilizer pellets (Osmocote, 17:11:10 N:P:K) twice a week, until reaching a final nitrate concentration of 20 μm L −1 . Once reaching the target stress conditions of 31 °C and 100 μm nitrate (T1), these conditions were sustained for additional 5 weeks (T1-T4; the stress phase), before gradually reducing the temperature back to 27 °C at a 0.8 °C/ day rate, and nutrient enrichment was stopped (T4). Once the baseline temperature of 27 °C was regained in all the baths (T5), plants were allowed to recover from the stress conditions for 3 additional weeks (T5-T8). H. stipulacea acclimate to environmental changes within a two weeks window [3] and thus both the four-weeks acclimation phase, and the duration of the experiment, were expected to expose morphological and physiological changes in the seagrass. Throughout the experiment, water exchanges were made weekly (~ 10% of seawater volume) and light intensity and salinity was kept constant. Water temperatures were logged automatically every hour with GHL PT-1000 electrodes (GHL, Aquarium Computer, Germany; 2/bath), and manually daily (WTW 340i, WTW, Germany). Water samples were taken (50 ml, filtered through a 0.22 micron syringe filter, kept frozen) for future nutrient analysis. Nutrients were measured to confirm the establishment of nutrient concentration differences between the high and low nutrient treatments. It should be noted that the measured nutrients reflect not only the administered quantities but also the accumulation of nutrients and interactions with biotic factors.

Phenological seagrass population descriptors
Seagrass phenology was a dependent variable in this experiment and was accounted for by quantifying the production of new shoots and the death of shoots. At the end of each time point (T 0 to T 8 ) the number of newly produced shoots and number of dead shoots since the previous time point were counted across each bath. The effect of temperature and nutrient stress treatments, source seagrass site and time on seagrass phenology was tested with a type-3 factorial anova, as implemented in StatModels [58].

Nutrient concentration measurements
Mesocosm seawater samples for nutrient analysis were collected with plastic syringes (~ 250 mL) above the H. stipulacea shoots of each aquarium (n = 4). Seawater samples were instantly filtered using sterile syringe filters (cellulose acetate; 0. 45

Epiphytic community samples collection
Epiphytic community samples were collected from both the mesocosm aquaria and from the H. stipulacea meadows at TY and SB, to study the bacterial dynamics in the experiment and to evaluate the relevance of the experimental results to the microbiota of the natural H. stipulacea population. In the mesocosm, two leaves from the third shoot of one SB and one TY plant were collected from four tanks per treatment (n = 4) in each of the following time points: following acclimation at the beginning of the experiment (T0), during the induced stress period (T3), at the end of the induced stress period (T5) and following the three weeks recovery period (T8). At each of these time points, one water sample was collected from each treatment. Third shoot leaves from wild and mesocosm plants were placed in DNeasy PowerSoil C1 solution (Qiagen) and sonicated for 3 min at 1.2 kHz in a sonicator bath. The C1 solution containing the sheared epiphytes underwent immediate DNA extraction. Water samples were filtered onto mixed cellulose esters 0.22-μm-pore-size filters, which were stored at − 80 °C until further processing.

16S rRNA metabarcoding
DNA was extracted from the epiphyte C1 solution using the DNeasy PowerSoil kit (Qiagen), and from water microbiome filters using the PowerWater DNA extraction kit (Qiagen), following the manufacturer's instructions. Metabarcoding libraries were prepared with a two step PCR protocol. For the first PCR reaction (PCR1), the V4 16S rRNA region was amplified, following [60], with the forward primer 515f 5'-tcg tcg gca gcg tca gat gtg tat aag aga cag GGT GCC AGC MGC CGC GGT AA-3' and the reverse primer 806R 5′-gtc tcg tgg gct cgg aga tgt gta taa gag aca gGA CTA CHV GGG TWT CTA AT-3′, along with artificial overhang sequences (lowercase). In the second PCR reaction (PCR2), sample specific barcode sequences and Illumina flow cell adapters were attached, using the forward primer ′5-AAT GAT ACG GCG ACC ACC GAG ATC TAC ACt cgt cgg cag cgt cag atg tgt ata aga gac ag-′3 and the reverse primer ′5-CAA GCA GAA GAC GGC ATA CGA GAT XXX XXX XXg tct cgt ggg ctc gg-′3′, including Illumina adapters (uppercase), overhang complementary sequences (lowercase), and sample specific DNA barcodes ('X' sequence). The PCR reactions were carried out in triplicate, with the KAPA HiFi HotStart ReadyMix PCR Kit (KAPA Biosystems), in a volume of 25 µl, including 2 µl of DNA template and following the manufacturer's instructions. PCR1 started with a denaturation step of 3 min at 95 °C, followed by 30 cycles of 20 s denaturation at 98 °C, 15 s of annealing at 55 °C and 7 s polymerization at 72 °C. The reaction was finalized with another minute-long polymerization step. PCR2 was carried out in a volume of 25 µl as well, with 2 µl of the PCR1 product as DNA template.
It started with a 3-min denaturation step at 95 °C, followed by 8 cycles of 20 s denaturation at 98 °C, 15 s of annealing at 55 °C and 7 s polymerization at 72 °C. PCR2 was also finalized with another 60-s polymerization step. Products of PCR1 and PCR2 were purified using AMPure XP PCR product cleanup and size selection kit (Beckman Coulter), following the manufacturer's instructions, and normalised based on Quant-iT PicoGreen (Invitrogen) quantifications. The fragment size distribution in the pooled libraries was examined on a TapeStation 4200 (Agilent) and the libraries were sequenced on an iSeq-100 Illumina platform, producing 150 bp paired end reads. Sequence data was deposited in the National Center for Biotechnology Information (NCBI) data bank, under BioProject accession PRJNA750596.

Amplicon sequence variants (ASVs) and taxonomy assignment
All the analyses carried out for this study are available as a Jupyter notebook in a github repository (GitHub: https:// git. io/ JBBeV, Zenodo https:// doi. org/ 10. 5281/ zenodo. 52172 77), along with the sequence data, intermediate and output files. The bioinformatics analysis was carried out within the Qiime2 [61] framework. DADA2 [62] was used to trim PCR primers, quality-filter, error correct, dereplicate and merge the read pairs, and to remove chimeric sequences, to produce the ASVs. Throughout the text, specific ASVs are referred to by the first 6 to 8 characters of their MD5 digests, which correspond with the biom table headers and the sequence IDs in the ASV fasta file. For taxonomic assignment, a naive Bayes classifier was trained using taxonomically identified reference sequences from the Silva 138 SSU-rRNA database [63] for the V4 fragments. All ASVs that were identified as mitochondrial or chloroplast sequences were filtered out from the feature table along with ASVs that had less than 30 occurrences across the dataset. An ASV phylogenetic tree was built with MAFFT 7.3 [64] for sequence alignment, and FastTree 2.1 [65], with the default masking options of the q2-phylogeny Qiime2 plugin. The phylogenetic tree was built denovo, including the observed ASV sequences only.

Microbial diversity analyses
ASVs shared between the epiphytic and planktonic communities were excluded from the biodiversity and differential abundance analyses of epiphytes to avoid artifact between-site similarity. It is impossible to know whether ASVs that occur both in the water and on the plant's surface are residents or water contamination. Possibly, such ASVs that have much increased densities on the leaf than in the water are residents, but we did not encounter such bacteria. Furthermore, diversity indices were not affected by this exclusion in preliminary analyses. Microbial diversity was estimated based on Faith's phylogenetic diversity [Faith PD;66] for alpha diversity and weighted and unweighted UniFrac distance [67] matrices for beta diversity, all of which consider the phylogenetic relationships among ASVs. Ordination of the beta-diversity pairwise distances was carried out with a principal coordinates analysis [PCoA; 41,68]. The contribution to bacterial compositional variance by time, temperature, nutrient concentrations and the source seagrass site was computed with the redundancy analysis anova procedure [69], as implemented in Vegan 2.5 [70]. Their contribution to alpha and beta diversity was additionally tested with a factorial ANOVA using the q2-longitudinal plugin [71]. P-values were corrected for multiple testing using the Benjamini-Hochberg procedure [72]. Corrected p-values are referred to as q-values throughout the text. All microbial diversity analyses were carried out using rarified tables. ANOVA tests were similarly carried out to test for significant differences between "within source site" and "among source site" pairwise UniFrac distance distributions.

Differentially abundant and explanatory taxa and ASVs
Two groups of order level taxa were considered when testing for differentially abundant orders, including the eight most relatively-abundant orders at each time point as well as low relative-abundance orders. The Kruskal Wallis test [73] was used to test for differential abundances among the highly abundant orders, at each time point separately, and additional differently abundant orders were identified with the Analysis of composition of microbiomes (ANOCM) procedure [74]. ANCOM was used to test for differential abundance at the family level as well, and BiPlot [41] analyses were used to highlight "important" ASVs, following the importance definition of Legendre and Legendre [41], which best explained the weighted and unweighted UniFrac pairwise distance matrices.

Core ASVs
Core ASVs were defined as ASVs found in all the wild samples and in all the samples of at least one of the treatments at time point zero. Similar core communities were identified for SB and TY epiphyte samples separately. To ensure that the regarded core communities bear relevance to the field, ASVs that were not present in field samples from the SB and TY sites were excluded from the core microbiome.