Responses of Invasive and Native Populations of the Seagrass Halophila stipulacea to Simulated Climate Change

Introduction

Seagrass meadows thrive worldwide in shallow sedimentary shorelines, where they fulfill important ecological services (including high primary productivity, production, and burial of organic carbon, nutrient cycling, and sediment stabilization) estimated at US$ 2.8 × 10⁶ km^–2 yr^–1 (Costanza et al., 2014). Yet, seagrass ecosystems are facing a global crisis due to both direct (reduced water quality, coastal development, and poor land use) and indirect (climate change) anthropogenic perturbations (Orth et al., 2006; Waycott et al., 2009; Short et al., 2011). Loss of seagrass ecosystems entails devastation of the associated biodiversity and primary productivity, reduction of local fishing grounds, enhancement of coastal erosion, and loss of major carbon sinks culminating in ecological and socio-economic dislocations (Wyllie-Echeverria and Cox, 1999; Orth et al., 2006; Fourqurean et al., 2012; Costanza et al., 2014; Arias-Ortiz et al., 2018). Increasing seawater temperatures associated with global warming directly affect seagrasses in several ways, including changes in growth rates, physiological functions, patterns of sexual reproduction, phenology, and seed germination (Marbà et al., 1996; Orth et al., 2006; Marbà and Duarte, 2010; Jordà et al., 2012). In addition, ocean warming and the associated increase in the frequency, duration and severity of heat waves (Easterling et al., 2000; Schär et al., 2004; Oliver et al., 2018) act indirectly by facilitating the establishment and spreading of invasive species (Stachowicz et al., 2002; Bradley et al., 2010) alongside reducing the biotic resistance of native communities to invader establishment (Diez et al., 2012).

The seagrass Halophila stipulacea (Forsskål) Ascherson (order Alismatales, family Hydrocharitaceae) is a dioecious, small tropical species native to the Red Sea, Persian Gulf, and Indian Ocean (den Hartog, 1970; Green and Short, 2003). Within the northern Gulf of Aqaba (GoA), one of the focal regions of this study, H. stipulacea is the most dominant, and sometimes only seagrass species (El Shaffai, 2011; Mejia et al., 2016), forming extensive mono-specific meadows found in both shallow and deep environments (1–50 m depth; Sharon et al., 2011; Winters et al., 2017), in soft sediments ranging from fine sand to gravel (Mejia et al., 2016; Rotini et al., 2017).

Following the opening of the Suez Canal in 1869, H. stipulacea soon became a Lessepsian migrant (den Hartog, 1972; Lipkin, 1975a) and has since become established in many parts of the eastern Mediterranean Sea (Lipkin, 1975a, b; van der Velde and den Hartog, 1992; Gambi et al., 2009; Katsanevakis and Tsiamis, 2009; Sghaier et al., 2011). Surprisingly, in 2002 H. stipulacea was also reported from the Caribbean Sea (Hector and David, 2004), and merely a decade later it was found in most Eastern Caribbean island nations as well as along the coast of the continental South American (Willette and Ambrose, 2009, 2012; Maréchal et al., 2013; Vera et al., 2014; Willette et al., 2014; Steiner and Willette, 2015; van Tussenbroek et al., 2016; Scheibling et al., 2018). In the Caribbean, H. stipulacea has been shown to displace native Caribbean seagrass species (e.g., Syringodium filiforme, Halodule wrightii, and Halophila decipiens) by physically displacing local Caribbean seagrass species by monopolizing their spaces (Willette and Ambrose, 2012), leading to a major change in the Caribbean’s seagrass landscape (Steiner and Willette, 2015).

Although H. stipulacea was included in the “100 Worst Invasive Alien Species in the Mediterranean” (Streftaris and Zenetos, 2006), in this basin, evidence for its “invasive” characteristics are scarce. Observations point toward a relatively limited “invasion success” in this region, as can also be inferred from the limited number of studies reporting competitive displacement of native Mediterranean seagrasses by the invasive H. stipulacea (Georgiou et al., 2016). However, some evidence of the potential competitive displacement by the invasive or “alien” H. stipulacea does exist. Sghaier et al. (2014) showed that a small patch (0.2 ha) of H. stipulacea in Cap Monastir Marina (eastern Tunisian coast), grew to cover more than 2.2 ha in only 4 years, and, in the process, displaced 50% of the native Cymodocea nodosa. More evidence of the future potential of H. stipulacea in the Mediterranean, comes from its slow (12 km yr^–1; Chiquillo et al., 2018) but continuous spread and expansion in this basin, moving westwards and northwards (Lipkin, 1975a; Gambi et al., 2009, 2018; Sghaier et al., 2011).

As H. stipulacea seems to be highly adapted to a wide range of environmental conditions (Gambi et al., 2009; Sharon et al., 2011) it has been predicted that the ongoing tropicalization of the Mediterranean Sea (i.e., becoming warmer and saltier; Bianchi and Morri, 2003; Borghini et al., 2014; Ozer et al., 2017) accompanied by the recent doubling of the Suez Canal (Galil et al., 2015) can facilitate H. stipulacea to outcompete local seagrass species (Sghaier et al., 2014) or to occupy newly available habitat following extirpation of local Mediterranean temperate seagrass species (Marbà and Duarte, 2010; Jordà et al., 2012; Chefaoui et al., 2018; Marín-Guirao et al., 2018; Savva et al., 2018).

In addition to the changing Mediterranean, H. stipulacea is also facing climate change within its native range, e.g., the northern GoA (Fine et al., 2013). Surrounded by deserts, the GoA receives very little rain and practically no runoffs from neighboring arid lands, creating an oligotrophic and warm tropical sea (Winters et al., 2006; Reiss and Hottinger, 2012). While studies have highlighted the potential effects of ocean warming on local coral reefs (Fine et al., 2013), little is known about the thermal tolerance of local H. stipulacea. In the GoA, where the tropical H. stipulacea is often found in meadows neighboring coral reefs (Mejia et al., 2016; Winters et al., 2017), loss of H. stipulacea could lead to increased sediment re-suspension, nutrient, and particulate loading into the water column (Orth et al., 2006) with adverse consequences to neighboring coral reefs due to algal overgrowth (Genin et al., 1995). Since tropical seagrass meadows also modify seawater carbonate chemistry by taking up CO₂, potentially mediating the effects of ocean acidification (Unsworth et al., 2012), loss of H. stipulacea meadows in the GoA would also enhance the decalcification process of neighboring corals (Hoegh-Guldberg et al., 2007). The fact that many reef fish species tend to hide in the structurally complex coral reefs to avoid predators, but forage in nearby structurally simple tropical seagrass meadows (Beck et al., 2006), entails that loss of seagrass meadows in the GoA will have cascading effects on local biodiversity.

Increased seawater temperatures can alter growth rates and other physiological and biochemical functions of seagrass (Short and Neckles, 1999; Beca-Carretero et al., 2018). With photosynthesis being among the first processes to be affected by elevated temperatures (Sharkey, 2005), several studies on the effects of thermal stress on seagrasses have focused on changes in photophysiological parameters (Ralph, 1998; Campbell et al., 2006; Winters et al., 2011). While studying gene expression is emerging as an important tool in seagrass thermal stress studies (Procaccini et al., 2012; Davey et al., 2016), so far, this was only applied in temperate seagrass species (Bergmann et al., 2010; Winters et al., 2011; Davey et al., 2016; Marín-Guirao et al., 2016), with no molecular studies ever performed in tropical seagrass species.

In both the eastern Mediterranean and the northern GoA, understanding the thermal tolerance of local H. stipulacea populations is crucial for seagrass conservation and management. The objective of the present study was to test whether H. stipulacea plants from native and invasive populations differ in performance and ecophysiology under thermal conditions that simulate current and future climate changes. We first analyzed temperature records and warming trends in the northern GoA (Eilat, Israel) and the eastern Mediterranean Sea (Limassol, Cyprus). We then collected plants from both native (GoA) and invasive (Limassol) H. stipulacea populations and compared their physiological, biochemical and molecular responses to current (26°C) and predicted thermal maxima (29 and 32°C) in a microcosm setup. Results suggest a rapid adaptation of the invasive population to the ongoing warming of the Mediterranean Sea. We predict that sea warming will cause vast reductions in H. stipulacea meadows growing in the GoA while facilitating their spread within the Mediterranean Sea. These results have important implications for biodiversity management and conservation in both of these regions.

Materials and Methods

Plant Materials and Habitat Descriptions

Halophila stipulacea plants (with shoots, attached roots, and rhizomes) were collected at 4 m depth from both an invasive (Limassol, eastern Mediterranean Sea, Cyprus; 34°42′20′′ N, 33°07′24′′ E; Figures 1A,B,D) and from a native population (Eilat, northern GoA, Israel; 29°32′47′′ N, 34°57′51′′ E; Figures 1A,C,E) on September 29^th 2017 and October 3^rd, respectively (n = 200 plants collected at each site; Figure 1A). The population growing in Limassol (“dream café” site; Figures 1B,D), came from a shallow area 100 m from shore that was protected by detached breakwaters. In this area we found three seagrass species: Cymodocea nodosa, residual patches of Posidonia oceanica and the invasive of H. stipulacea mixed in-between them (Figure 1D). Sediment at this site was sandy and the slope was very gentle. Similarly, the population growing in Eilat, came from a shallow meadow, 100 m from shore that is the northern-most tip of the GoA. In this area (“North Beach” site), H. stipulacea grows in extensive mono specific meadows that start at 2 m and extend to deeper waters (30 m, Winters et al., 2017). Sediment at this site was sandy and the slope was very gentle (Mejia et al., 2016). In both sites, the H. stipulacea plants that were collected, were at least 2 m away from each other in order to avoid collecting samples from the same clone. During collection, water temperatures at both sites were around 25–26°C.

Environmental Conditions at the Collection Sites

While environmental conditions at the collection site in Eilat (e.g., seawater temperatures, salinity, sediment type, Kd and nutrient levels) have been described before (Mejia et al., 2016; Winters et al., 2017) and are well-monitored in nearby sites for nearly 20 years by the National Monitoring Program at the Gulf of Eilat (NMP-Israel¹), no long-term dataset exists for anywhere near the Cyprus site. Thus, comparisons between the long-term conditions at the two sites were based on two data sets. To compare seawater salinities between the two sites, the Bio-ORACLE global dataset (Tyberghein et al., 2012; Assis et al., 2018) was used to compare minimum, maximum, and average values of water salinities over a relatively long-term data set (2000–2014, based on monthly averages; Supplementary materials S1). To compare environmental temperatures and rates of warming between study locations in Eilat and Limassol, we obtained daily sea surface temperatures (SSTs) for the period between 1982 and 2017 from the NOAA dOISST.v2 dataset at www.ngdc.noaa.gov (Smith, 2016). We used Advanced Very High-Resolution Radiometer (AVHRR) only data, due to its longer temporal span and because it has been shown to out-perform other datasets in coastal areas (Lima and Wethey, 2012). For each of the two locations, daily climatology parameters were produced by averaging corresponding daily values over the entire time span. Average warming rates were computed as the slope of the linear regression of seasonally detrended SST versus time (Lima and Wethey, 2012). Seasonal variability was removed by subtracting from each day its corresponding climatology, i.e., the average temperature for that same day computed over the entire period under analysis (Chatfield, 2016). To compensate for temporal autocorrelation that could cause an underestimation of the standard errors of the warming rates, we followed Foster and Rahmstorf (2011) and reduced the data per effective degrees of freedom, using a conservative approach where autocorrelation structure was treated as an autoregressive moving average process by applying the ARMA (1,1) model (Foster and Rahmstorf, 2011). Warming rates for February and August were computed based on averaged daily values for those months. All calculations were done in R Core Team (2018).

Mesocosm Facility

Collected plants were transferred to our seagrass dedicated microcosm facility in zip lock bags filled with seawater (Figure 2A; recently described in Oscar et al., 2018). Here plants were planted in 15 aquaria (40 cm width, 33 cm height, 45 L of seawater in each aquarium) layered with 20 L (∼7 cm high) of natural sediment (sediment collected from the shores of the GoA, Eilat; Figure 2A) and filled with artificial seawater (Red Sea Salt, Israel). Each aquarium contained 20 individuals (i.e., 10 invasive and 10 native plants) with invasive and native populations separated by plastic dividers (insert in Figure 2A). A total of 5 aquaria were assigned for each treatment and considered as experimental replicates (n = 5 in each treatment). Salinity was maintained at 40 PSU (partial salinity units, equivalent to the average salinity of the GoA’s seawater; The National Monitoring Program at the Gulf of Eilat, 2014¹) by adding purified water to compensate evaporation. In order to retain the quality of experimental seawater, 20% of the total volume of water in each aquarium was replaced by newly prepared seawater every week. This microcosm facility is fully controlled for temperature and light (ProfiLux 3.1T eX Controller, GHL; GHL, Germany). Irradiance was provided by T5 fluorescent lamps (Osram Lumilux HO 865/54W cool daylight with the color temperature of 6500 degrees Kelvin) that were set to a 10:14 h light: dark cycle, with intensities of 120 μmol photons s^–1 m^–2 measured at the bottom of the aquaria.

Water temperatures were measured automatically every 1 h with PT1000 thermal sensors (GHL, Germany) connected to two ProfiLux 3.1T eX GHL controller systems (GHL, Germany). In each one of the three treatments (detailed below), water temperature was followed over time using two sensors, situated in two different aquaria belonging to the same treatment (Figure 2A). Salinity and temperature were also checked manually every day using the WTW Multi 340i Portable Meter (Xylem Analytics, Weilheim in Oberbayern, Germany). Each one of the 15 aquaria, was fit with a submersible pump to provide gentle water movement within the aquarium.

Experimental Treatments

Following 3 weeks of the acclimation period, heating treatments were applied. These included maintaining plants at 26°C (just under the current maximal summer water temperatures in Eilat; Winters et al., 2006), and two predicted climate change scenarios, including 29 and 32°C. The experiment was started by 1 week of gradual increases of the temperatures in the control tanks (from 25 to 26°C) and in the two thermal stress treatments (from 25 to 29°C, and from 25 to 32°C, at heating rates of 1 and 0.5°CC day^–1, respectively, Figure 2B). Plants were kept for 2 weeks at these high temperature treatments (29 and 32°C), after which, water temperatures were gradually ramped down (at cooling rates similar to those mentioned above) to 26°C at the end of the experiment, with measurements at T4 taken 24 h after reaching control temperatures (Figure 2B). Plants exposed to the control treatment were kept at 26°C throughout the entire experiment (Figure 2B).

Plant Performance

Plant performance was estimated by measuring a total number of shoots (shoots plant^–1), horizontal growth rates (mm week^–1), and leaf surface area (cm²) before, during, and after the thermal stress treatment (Figure 2B). For this, three plants per population in each aquarium were randomly selected for measuring the total number of shoots every week, which were averaged into a single biological replicate (i.e., within each treatment, each aquarium = one biological replication). In addition, three plants were randomly selected per aquarium and used to assess horizontal growth rates once per week. For this, the distance from the basal leaf sheath to the tip of the three youngest leaves was measured and these three measurements were averaged into a single biological replicate. Finally, mature intact leaves (i.e., fully grown, often the third youngest shoot with leaf tip; broken leaf tips were not used) were detached from one plant per population in each aquarium once a week, scanned on a digital scanner (CanoScan LiDE 120, Canon U.S.A., Inc.) and leaf surface area was estimated by using ImageJ (Abràmoff et al., 2004).

Biochemical Measurements

Scanned leaves that were used for leaf surface area measurements (mentioned above) were subsequently used for measuring changes in chlorophyll a and b contents. Chlorophylls were extracted in 100% methanol according to Wellburn (1994) and adapted to a 96 well microplate reader (Epoch 2 Microplate Spectrophotometer, Biotek) as described in Warren (2008). After scanning, leaves were immediately weighed and 50 mg (fresh weight) of tissue was homogenized into small pieces before transferring into 1.5 mL tubes. 1 mL of 100% methanol was added into the tube, covered with aluminum foil, and kept at 4°C overnight. 200 μL of chlorophyll extract was then transferred into a 96-well plate (CELLSTAR^®, Greiner Bio-One) and the absorbance read at 470, 652, 665, and 750 nm. Chlorophyll a and b contents were calculated according to Tran et al. (2018).

Effects of thermal stress on biochemical parameters were studied over time (Figure 2B) by measuring changes in the total protein content (Elavarthi and Martin, 2010). For these measurements, 50 mg of the fourth youngest leaves were ground thoroughly to a fine powder using stainless steel beads (SSB32, Next Advance) in a 1.5 mL Eppendorf Safe-Lock Tube in a Bullet Blender tissue homogenizer (Next Advance). A volume of 0.5 mL of 0.2M potassium phosphate buffer (pH 7.8; Carlo Erba reagents) containing 0.1 mM Ethylenediaminetetraacetic acid (EDTA, Sigma-Aldrich), 0.5% (w/v) Polyvinylpyrrolidone (PVPP, Sigma-Aldrich), and 1% Protease inhibitor cocktail (Sigma-Aldrich) to neutralize phenol effect (Zhang et al., 2007) were added to the sample. The tube then was centrifuged at 13,500 g at 4°C for 20 min, and an approximately 800 μL of the supernatant was transferred to a new 1.5 mL reaction tube (Greiner Bio-One) and stored at 4°C until being used for measuring protein content. Total protein contents were determined using the Bradford assay (Bradford, 1976) adapted for microplate readers using bovine serum albumin (BSA, Sigma-Aldrich) as the protein standard according to Tran et al. (2018).

Photosynthetic Functionality

The photophysiological effects of the heat wave treatment were assessed by measuring the effective quantum yield of photosystem II (PSII) measured as ΔF/$F_m^{\prime }$ = ($F_m^{\prime }$−F/F_m) using pulse amplitude–modulated (PAM) chlorophyll fluorometer (Diving-PAM, Walz, Germany) as commonly applied in many previous studies (Ralph and Gademann, 2005; Saroussi and Beer, 2007; Winters et al., 2011). All PAM measurements were made midday, and taken underwater using the leaf distance clip on the base of the 2^nd youngest leaf (n = 5 randomly chosen shoots at each time point from each population and treatment) and after removal of epiphytes. Irradiance levels at this position ranged between 130 and 150 μmol photons m^–2 s^–1 (measured with the micro quantum sensor, Apogee, United States).

Sampling for RNA, RNA Preparation, Primer Design, and Quantitative Real-Time PCR (qRT-PCR) Assays

Fully mature leaves (2^nd youngest shoot; approximately 50–100 mg fresh weight) from each biological replicate were randomly collected on day 21, after 2 weeks of exposure to the heat-stress (T3), and day 28, at the end of the experiment (T4), 24 h after temperatures had returned to the control temperature (26°C) (Figure 2B). Samples were quickly cleaned from epiphytes and blotted dry using tissue paper, inserted into 1.5 mL Eppendorf Safe-Lock Tubes, and snap frozen immediately in liquid nitrogen before being stored in −80°C until extraction of RNA.

Total RNA was isolated from young fully mature leaves of H. stipulacea control (26°C) and heat treated (29 and 32°C) plant samples using the TRIzol method (2M guanidine thiocyanate, 4M ammonium thiocyanate, 3M sodium acetate, 38% phenol, 5% glycerol) (Kant et al., 2006). The total RNA quantity and quality were assessed with a NanoDrop spectrophotometer (ND-1000; NanoDrop Technologies) and its integrity was checked by running on a 1% agarose gel electrophoresis. Extracted total RNA was then treated with RNase-free DNase I (EN0525, Thermo Fisher Scientific) to eliminate genomic DNA. cDNA was synthesized from 1 μg of DNase I treated total RNA using the High-Capacity cDNA Reverse Transcription Kit (Catalog no. 4368813, Applied Biosystems) according to the manufacturer’s instructions. The resulting cDNA was diluted 1:10 prior to quantitative real-time PCR (qRT-PCR) assays.

The expression of heat shock protein 70 (Hsp70), Superoxidase dismutase [Mn] (SOD), Photosystem II protein D1 (PsbA), and Photosystem II protein D2 (PsbD) genes were analyzed by qRT-PCR. The gene-specific primer pairs for Hsp70 and UBQ (Ubiquitin – Reference gene) were designed with Primer Express v2.0 (Applied Biosystems, United States) using the recently sequenced and annotated H. stipulacea transcriptome (Oscar et al., unpublished data). Ubiquitin has commonly been used in seagrass studies as a reference gene and has been shown to be stable in a several seagrass species exposed to different salinity, pH, light, and thermal conditions (Serra et al., 2012; Marín-Guirao et al., 2016). Primer sequences for SOD were taken from Winters et al. (2011), while primer sequences for PsbA and PsbD were taken from Marín-Guirao et al. (2016) (see Table 1 for primer sequences).

Real-time PCR was performed with the ABI PRISM 7500 Sequence Detection System (SDS, Applied Biosystems) using SYBR Green to monitor double-stranded DNA (dsDNA) synthesis. Each reaction contained 5 μl PerfeCTa SYBR Green FastMix (Quanta Biosciences), 4 μl (20 ng) cDNA and 300 nM of gene-specific primer in a final volume of 10 μl. PCR amplifications were performed using the following conditions: 95°C for 30 s, 40 cycles of 95°C for 5 s (denaturation) and 60°C for 35 s (annealing/extension). Data were analyzed using the SDS 1.3.1 software (Applied Biosystems). To check the specificity of annealing of the primers, dissociation kinetics was performed at the end of each PCR run. The experiment was performed with three biological replicates and each reaction was performed in triplicates. The relative quantification values for each target gene were calculated by the 2^–ΔΔCT method (Livak and Schmittgen, 2001) using the Ubiquitin gene as an internal reference (Serra et al., 2012; Marín-Guirao et al., 2016). To ensure the validity of the 2^–ΔΔCT method, twofold serial dilutions of cDNA were used to create standard curves. The amplification efficiency (E) values of the target and reference genes were shown to be approximately equal, with all E > 0.91, and all linear fit R² > 0.95 (Livak and Schmittgen, 2001).

Statistical Analysis

For plant performance, biochemical and photosynthetic functionality measurements, results from within each tank (biological replicates) were averaged and considered as one experimental replicate. In total, 5 experimental replicates (n = 5 aquaria) from each of these measurements were used for statistical analyses. For gene expression results, three experimental replicates (n = 3) were used to statistically determine the effects of treatments within each gene at each time-point.

GLM (general linear model)-repeated measures analyses were conducted on plant performance, biochemical and photosynthetic functionality results from each population separately using the statistical package IBM SPSS Statistics (v. 20) to test the overall effects of different temperatures (control 26°C, treatment 29°C, and treatment 32°C). Normality of data was tested using the Shapiro–Wilk test and homogeneity of variances were verified using Levene’s (1960). When parametric assumptions were not met, data were transformed using the Johnson transformation (Johnson, 1949) by using the statistical analysis package SPC for Excel (BPI Consulting, LLC, v.5). Whenever a significant difference was detected, Tukey HSD’s multiple comparison post hoc test was applied.

Effects of different experimental temperatures on each gene of interest (GOI) were analyzed separately for each population at each time-point using Pair-wise comparison test with Permutational Multivariate Analyses of Variance (PERMANOVAs) with Primer 6 v.6.1.16 and PERMANOVA + v.1.0.6 software package (PRIMER-E, Ltd.) (Anderson et al., 2008).

Results

Environmental Conditions at the Two Sites

Comparisons between the long-term environmental conditions at the two sites that were based the Bio-ORACLE global dataset (Supplementary materials S1) demonstrated that salinity was very similar between the sites, with an annual average of 39.2 and 38.8 [Practical Salinity Scale (PSS)] for Cyprus and GoA sites, respectively (results not shown). Differences in minimal and maximal water salinity between the sites were even smaller (data not shown).

Long-term (1982–2017) analysis of daily SSTs from both sampled locations showed similar patterns in summer maxima. In both sites, highest temperatures were observed during August, with average water temperatures in August reaching 26.92 ± 0.02°C SE in Cyprus and 27.41 ± 0.02°C SE in the northern GoA (Figures 3A,B). During February, the coldest month in both locations, differences were much larger, with average water temperatures in Cyprus reaching 16.61 ± 0.01°C SE, while in the subtropical GoA, February temperatures did not go below 21.27 ± 0.01°C SE (Figures 3A,B).

In both sites, average warming rates over the multi-decade were significantly positive (regression analysis, p < 0.05), indicating that sea surface warming is taking place, but with some clear differences between the sites (Figures 3C,D). While in Cyprus both winter (February) and summer (August) warming trends were very similar (0.31 ± 0.17°C SE/decade and 0.32 ± 0.11°C SE/decade, p < 0.05, respectively; Figure 3D), in GoA August warming was more than double of that observed in February (0.46 ± 0.18°C SE/decade and 0.21 ± 0.15°C SE/decade, p < 0.05, respectively; Figure 3C). Seasonally detrended SST regression analysis continued to show that both sites have been warming significantly, although this data shows that warming has been faster in Cyprus (0.36 ± 0.06°C SE/decade, p < 0.05) compared with the northern GoA)0.26 ± 0.06°C SE/decade, p < 0.05) (Figures 3E,F).

Plant Performance

Results from the number of shoots throughout the experiment (Figures 4A,B) demonstrate the negative effects of increased water temperature on both native and invasive plants. Native plants suffered massive reductions in shoot number (Figure 4A) resulting in a significant difference between control plants and heated plants (GLM-repeated measures, p < 0.001, Table 2). Results from Tukey HSD post hoc multiple comparisons detected significant differences among the three experimental conditions, indicating that exposure to 29°C was already warm enough to stress the native plants significantly (C26 ≠ T29 ≠ T32). In contrast, while the number of shoots of the invasive plants decreased slightly over time when exposed to thermal stress (Figure 4B), no significant difference was detected with GLM-repeated measures or the Tukey HSD post hoc multiple comparisons, suggesting a higher thermal tolerance level of the invasive plants.

Similarly, increased water temperature seemed to have little or no effect on the growth rates of invasive plants as heated plants were able to maintain their growth rates similar or even slightly higher than growth rates of control plants (Figure 4D), compared with the exposure of the native plants to 32°C which dramatically slowed their growth rates (GLM-repeated measures, p < 0.001, Figure 4C and Table 2). Indeed, even during the recovery phase, native plants that were previously exposed to 32°C almost stopped growing (Figure 4C). Native plants that were exposed to 29°C heated plants, significantly reduced their number of shoots, but this reduction was much less than in the 32°C treatment (Figure 4C).

Results from the measurements of leaf surface area show big differences in the response to the increased water temperature. The fact that native leaves start bigger than their invasive counterparts corresponds to them being larger year-round also in the field (Nguyen et al., in review). In native plants, no significant difference among treatments were found (Table 2). However, in the invasive plants, extreme increased water temperature of 32°C, surprisingly caused increases to their leaf sizes (Figure 4F). Indeed, significantly larger leaves were observed in 32°C heated invasive plants (GLM-repeated measures, p = 0.022 and Tukey HSD post hoc Multiple comparisons: C26 = T29, C26 = T32, T29 ≠ T32, Table 2). This phenomenon is actually in correspondence with the growth rate results (Figure 4D) strengthening the hypothesis about a higher thermal tolerance of the invasive plants in comparison to their native counterparts.

Biochemistry

Chlorophyll a and b contents were similar in the two populations at the beginning of the experiment (T0, approximately 0.28 mg g FW^–1; Figures 5A,B). Results (Figure 5) demonstrate that the 29°C heated plants (native and invasive) maintained their biochemical contents (i.e., chlorophyll a, chlorophyll b and proteins) similar to their control plants. Indeed, in both populations, the GLM-repeated measures and Tukey HSD post hoc multiple comparisons (Table 2) showed no differences between the chlorophyll a, chlorophyll b and proteins contents of control vs. 29°C treatment (C26 = T29, Table 2).

On the other hand, the responses of native and invasive plants to the 32°C warming scenario showed big differences from T2 onwards (Figures 5A,B). While under stress of 32°C chlorophyll a and b contents degraded significantly in the native plants (from 0.28 to 0.20 mg g FW^–1, Figure 5A; GLM-repeated measures p = 0.002, Tukey HSD post hoc multiple comparisons C26 = T29 ≠ T32, Table 2), in the invasive plants, chlorophyll a and b content increased (from 0.23 to 0.50 mg g FW^–1, Figure 5B), but this increase was found to be non-significant (GLM-repeated measures p = 0.283, Tukey HSD post hoc multiple comparisons C26 = T29 = T32, Table 2).

A similar result was evident from the protein contents (Figures 5C,D). In the native plants, there was a strong trend of reduction in protein content of plants exposed to the 32°C warming scenario (GLM-repeated measures, p = 0.267, Tukey HSD post hoc multiple comparisons, C26 = T29 = T32, Table 2), indicating a possible thermal threshold response.

In contrast, compared to 26 and 29°C treatments, plants from the Mediterranean actually increased significantly their protein contents over time when exposed to the 32°C warming scenario (GLM-repeated measures, p < 0.001, Tukey HSD post hoc multiple comparisons, T26 = T29 ≠ T32, Figure 5D and Table 2).

Photosynthetic Functionality

In terms of photosynthetic functionality, results demonstrate differences between the response of native and invasive populations of H. stipulacea to the extreme thermal stress. At first, as observed in both populations, increased water temperature to 29°C did not result in any significant difference between heated plants in comparison to control ones (Figures 6A,B). Interestingly, while the extreme thermal condition (i.e., 32°C) lowered significantly the effective quantum yield in native plants (Figure 6A, GLM-repeated measures, p < 0.001, Tukey HSD post hoc multiple comparisons, T26 = T29 ≠ T32, Table 2), oppositely it had no negative effect on the invasive plants (Figure 6B). This contrasting phenomenon suggests a clear difference in thermal tolerance between the two populations.

Gene Expression

Gene expression results from this study provide the very first insight into some of the molecular responses of H. stipulacea to thermal stress. Results include the expression of three different gene categories including a molecular chaperone that is expressed in response to stress (HSP70), photosynthesis related-genes (psbD and psbA) and an antioxidant related-gene (SOD; Figure 7).

All studied genes were significantly influenced by population, treatment or their interactions (Table 3). Overall, expression patterns were lower under moderate stress (29°C) as compared to the more severe thermal stress of 32°C (Figure 7). Native plants expressed significantly more HSP70 during the stressful condition (T3) both at 29°C and significantly more at 32°C (Figure 7A), compared with the invasive plants where exposure to both temperatures did not cause a significant upregulation of HSP70 (Figure 7C). Returning plants from both populations to control temperatures at T4 entailed a shift to down-regulation of HSP70, with plants previously exposed to 32°C having significantly stronger down-regulation compared with control plants (Figures 7A,C; Pair-wise comparison, p < 0.001, Table 3). Down-regulation of plants previously exposed to 32°C was stronger in native compared with invasive plants (Figures 7B,D).

During the heat stress condition (T3), both native and invasive plants showed up-regulation of psbD and down-regulation of psbA (Figures 7A,C). The up-regulation of psbD was significant only for the invasive plants (Figures 7A,B, Pair-wise comparison, p < 0.001, Table 3). When water temperatures returned to the control level (Figures 7B,D), the 32°C invasive heated-plants shifted their expression from a very weak down-regulation during the heat stress itself, to a significant upregulation of psbA (Figure 7D, Pair-wise comparison, p < 0.001, Table 3). No change was found in the expression pattern of the native plants (Figure 7B).

Lastly, thermal-induced expression of the antioxidant-related gene (SOD) concurs with other genes showing the invasive plants were more responsive to increased water temperature. Both native and invasive plants significantly up-regulated SOD to compensate with the increased damage caused by thermal stress (Figures 7A,C), with a much stronger up-regulation found for the invasive population. Returning plants to control water temperatures (T4) resulted in reversing the upregulation of SOD found during the heat wave (T3, Figures 7A,C) to a significant down-regulation during of SOD in T4 (Figures 7B,D). In both populations, there was a stronger down-regulation of SOD during T4 in plants that were previously exposed to the higher temperatures (32°C). The fact that this downregulation of SOD was only significant in the invasive plants (Figure 7D) indicates that the invasive plants might be more thermally plastic than the native ones.

Discussion

Invasiveness of plants may be accompanied by rapid evolution due to selection in the colonized habitat (examples from land plants – Maron et al., 2004; Novy et al., 2013; Sultan et al., 2013). Here we demonstrate that plants of invasive seagrass population may have undergone an additional adaptation to the changes in the colonized habitat. Plants of H. stipulacea from the invasive population performed overall better than plants from the native population in a common garden simulated thermal stress experiment. We show that while native plants were negatively affected in photo-physiological and growth responses by thermal stress, invasive plants did not suffer and might have even benefited from it. We also show that the putative mechanism of adaptation to thermal stress could be the differential gene expression. Taken together, these results hint toward rapid evolution of H. stipulacea plants following their Lessepsian migration through the Suez Canal. This process could possibly work through strong selection on the colonizing genotypes during (e.g., the warm the ballast water of ships) or after their transport (post-introduction evolution – through exposure to many heat waves). Such a process would help to maximize the performance and establishment of the invasive population in their new environment (Molina-Montenegro et al., 2018).

Our results support the hypothesis that climate change amplifies the spread of invasive plants (Bradley et al., 2010; Zelikova et al., 2013). The long-term analysis of daily SST concurs with previous studies in the GoA (Fine et al., 2013) and the eastern Mediterranean Sea (Rilov, 2016; Ozer et al., 2017). In both of our studied sites, average SST warming during the last 35 years was higher than the global mean coastal SST trend of 0.17 ± 0.11°C/decade (Liao et al., 2015). In the GoA, the temperature has been warming at a rate 1.63 times faster than this global coastal average, but in Cyprus, warming has been occurring at even faster rates, 2.25 times faster than this global coastal average SST trend.

Under assault for at least a century, seagrasses may be nearing a crisis with respect to global sustainability (Waycott et al., 2009; Hyndes et al., 2016). In addition to increasing mean temperatures (Jordà et al., 2012), seagrass ecosystems will also face increasing climate variability that will affect the frequency, duration and intensity of summer heat waves (Schär et al., 2004; Diez et al., 2012; Oliver et al., 2018), with a projected longer duration of warm periods (Jordà et al., 2012; Oliver et al., 2018). Furthermore, these global climatic changes have been predicted to exacerbate invasions (Buckley and Csergő, 2017).

The results from our microcosm experiment suggest that native and invasive populations of H. stipulacea might differ in their evolutionary trajectories under future projections of climate changes. Plants from the native population of GoA showed signs of stress during and after the exposure to high temperatures (i.e., 32°C) in most of the measured parameters: they lost high numbers of shoots per plants (Figure 4A), reduced their growth rates (Figure 4C) and leaf sizes (Figure 4E). The plants from the invasive population of Cyprus, on the other hand, kept the number of shoots constantly (Figure 4B), their growth rates were similar to control plants (Figure 4D), and surprisingly doubled their leaf sizes under the same high temperature of 32°C; these results concur with those of Georgiou et al. (2016) that also showed increased growth rates of H. stipulacea from Cyprus exposed to increased water temperatures. Biochemical assays further suggested that GoA plants suffered at simulated high temperature and could not be able to maintain their chlorophylls as well as their protein contents (Figures 5A,C). In contrast, exposure to thermal stress caused Cyprus plants to increase significantly their photosynthetic pigments and protein (Figures 5B,D and Table 2). Application of PAM fluorometry strengthens our conclusions as significant reduction in ΔF/Fm′ was only found for native plants (Figure 6).

The results from this study suggest that native and invasive H. stipulacea plants might underwent different evolutionary pathways, and potentially will further diverge in their evolutionary trajectory given that the current trends of increased sea temperatures continue. While the plants from the northern GoA suffered greatly from simulated increases in seawater temperatures, the plants from the invasive population from the eastern Mediterranean did not suffer physiologically and seemed to take advantage of simulated seawater warming.

We speculate that if these plants are representative of larger populations, future increase in seawater temperature will threaten the native populations and benefit the invasive ones. Further exploration of this prediction should use multiple populations from both seas, in order to assess the width of the phenomenon found here in representing plants from one site in each sea.

Although It has been suggested that the range of expansion of H. stipulacea in the Mediterranean will be limited by the 15°C sea surface isotherm experienced during winter seasons (reviewed by Gambi et al., 2009), our results are in agreement with the those of Georgiou et al. (2016) that predicted that the rising winter and summer temperatures in the Mediterranean Sea (Ozer et al., 2017) will not be a limiting factor of the expansion of H. stipulacea in the Mediterranean Sea.

The work presented here, alongside that of Georgiou et al. (2016) represent the only two studies that experimentally investigated the effects of temperature on H. stipulacea, shedding light on the gap of knowledge in this ecologically important species compared with the numerous published studies on the effects of thermal stress on both temperate and tropical species (Bergmann et al., 2010; Winters et al., 2011; Tutar et al., 2017; George et al., 2018).

Our gene expression results provide some very first molecular knowledge regarding the effects of temperature on H. stipulacea. Changes in the expression level of heat shock protein (HSP70) in our study shed light on some of the molecular mechanisms behind the physiological differences in our experimental plants. Increases in expression of HSP70 as a result of thermal stress were also shown in other seagrass species including Zostera marina, Posidonia oceanica, and Cymodocea nodosa (Bergmann et al., 2010; Marín-Guirao et al., 2016). The increase in expression of HSP70 in H. stipulacea plants from the invasive population under simulated increased temperature, suggests that invasive plants might have mechanisms in place to cope with sea warming. The significantly higher expression levels of SOD in the invasive population could give these plants the protection from reactive oxygen species (ROS) needed to continue normal photosynthesis under thermal stress (Lesser, 1997; Sharkey, 2005). The fact that psbA (the chloroplast gene that codes for the photosystem II protein D1) tends in both populations to be only weakly down-regulated during the thermal stress, might indicate the possibility that photosystem II was not the main site of damage of the heated plants.

While we did not include genotyping data in this current study, a recent study on reconstructing the invasion history of the H. stipulacea applied 2bRAD sequencing to several populations, including the two populations used here (Chiquillo et al., unpublished data). Results from that study clearly demonstrate that genetic diversity in the Mediterranean Sea is actually lower than that in the GoA. The fact that the invasive Cyprus H. stipulacea population, with lower genetic diversity, did better under thermal stress than its native GoA higher diversity counterpart (Chiquillo et al., unpublished data), does not concur with those of Ehlers et al. (2008) that showed that Zostera marina populations with higher genotypic diversity survived thermal stress better than populations with reduced genotypic diversity. One possible explanation for this could be a high pre-existing genetic diversity in the native GoA range (as shown by Chiquillo et al., unpublished data) and the later filtering of genotypes that could be locally pre-adapted in the new invasive range to climate change (Molina-Montenegro et al., 2018).

With the recent doubling of Suez Canal (Galil et al., 2015, 2017) and the ongoing tropicalization and warming of the Mediterranean (Bianchi and Morri, 2003; Borghini et al., 2014) a process that is happening even faster in the eastern Mediterranean (Ozer et al., 2017), H. stipulacea could potentially become more prevalent in these waters in the coming years. This prediction is quite probable considering that conditions in the Mediterranean Sea are becoming less favorable for its temperate seagrass species (Jordà et al., 2012; Chefaoui et al., 2018; Marín-Guirao et al., 2018; Savva et al., 2018) and more welcoming to tropical species (Sghaier et al., 2014; Georgiou et al., 2016; Gerakaris et al., 2019).

In the Caribbean Sea, H. stipulacea has been demonstrating its negative characteristic by invading the entire Caribbean (Willette and Ambrose, 2009, 2012; Maréchal et al., 2013; Vera et al., 2014; Willette et al., 2014; Steiner and Willette, 2015; van Tussenbroek et al., 2016; Scheibling et al., 2018). Studies in the Caribbean have shown that H. stipulacea is physically displacing local Caribbean seagrass species (e.g., Syringodium filiforme, Halophila decipiens, and Halodule wrightii; Willette and Ambrose, 2012), thus changing the Caribbean’s seagrass landscape (Steiner and Willette, 2015). Taking into account what is already happening in the Caribbean Sea, combined with our knowledge on the changing Mediterranean Sea (Jordà et al., 2012; Rilov, 2016; Ozer et al., 2017; discussed above), our results presented here demonstrate clearly that the spreading of H. stipulacea in the Mediterranean Sea has a high potential threat to biodiversity of the Mediterranean region.

While the spread of H. stipulacea in the Mediterranean Sea is worrying, there is also a threat to the native sub-tropical populations in the GoA. The National Monitoring Program at the Gulf of Eilat (NMP-Israel) has been collecting high-resolution data of local marine ecosystems since 2004², but most of this program’s efforts have been focused on the coral reefs in the region, while local H. stipulacea meadows have only very recently been included in this program.

Results from our study highlight the importance of long-term, coordinated monitoring of this important engineering species, in both native and invasive populations of H. stipulacea. It should be noted that we used only one population from each basin. While running simulated thermal stress experiments with multiple temperatures and populations is complicated enough (3 temperatures × 2 populations in the present study), including more populations in such experiments would be necessary in the future to account for the genetic variability (Procaccini et al., 1999) and different thermal tolerances that might exist in H. stipulacea populations in both the Mediterranean and the Red Sea. Thus, we call for further studies to include multiple populations from each basin to reconfirm the results presented here.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

HN, YS, and GW conceived and designed the experiments. HN, NY, and SB performed the experiments. HN and GW analyzed the results. HN, FL, YS, and GW wrote the manuscript. All authors reviewed the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank ICA in Israel and the Israeli Ministry for Science and Technology (MOST) for supporting the Dead-Sea & Arava Science Center (ADSSC) and for supporting the seagrass mesocosm facility. We also thank the Arava International Center for Agriculture Training (AICAT) for giving financial support to HN during his M.Sc. We also thank the Manna Center Program for Food Safety & Security, at the School of Plant Sciences and Food Security Faculty of Life Sciences Tel Aviv University (TAU) for their ongoing support for the co-op between TAU and ADSSC. We also thank the Blaustein Center for Scientific Cooperation (BCSC) and the PBC Program of Israeli Council for Higher Education for postdoctoral fellowships to NY. FL was funded by FEDER funds through COMPETE (POCI-01-0145-FEDER-031088) and by National funds through FCT (PTDC/BIA-BMA/31088/2017).

Footnotes

1. ^ http://www.meteo-tech.co.il/EilatYam_data/ey_data.asp (accessed July 3, 2017)
2. ^ http://www.iui-eilat.ac.il/Research/NMPAbout.aspx (accessed July 3, 2017)

References