Long-Term Effects of Elevated Co 2 on the Population Dynamics of The Seagrass Cymodocea Nodosa : Evidence from Volcanic Seeps

We used population reconstruction techniques to assess for the first time the population dynamics of a seagrass, Cymodocea nodosa , exposed to long-term elevated CO 2 near three volcanic seeps and compare them with reference sites away from the seeps. Under high CO 2 , the density of shoots and of individuals (apical shoots), and the vertical and horizontal elongation and production rates, were higher. Nitrogen effects on rhizome elongation and production rates and on biomass, were stronger than CO 2 as these were highest at the location where the availability of nitrogen was highest. At the seep where the availability of CO 2 was highest and nitrogen lowest, density of shoots and individuals were highest, probably due to CO 2 effects on shoot differentiation and induced reproductive output, respectively. In all three seeps there was higher short- and long-term recruitment and growth rates around zero, indicating that elevated CO 2 increases the turnover of C. nodosa shoots.


Introduction:
Seagrass photosynthesis may be limited by dissolved inorganic carbon concentration of the water column, which has led part of the research community to suggest that, seagrass could benefit from the global increase in oceanic CO2 levels (Koch et al., 2013;Russel et al., 2013;Borum et al., 2018). However, several CO2 enrichments studies (ex-situ and in-situ) on seagrasses reveal that the response is species-specific and more complex than originally thought. Relatively long experiments, such as that carried out by Alexandre et al., (2012), which lasted for five months reported positive effects on the photosynthetic production of Cymodocea nodosa under experimental CO2 enrichment conditions, but not effects on growth. Similarly, experiments with elevated CO2 enrichment conditions on Zostera marina for over one year showed no effects on specific growth rate, size or sugar content of leaves, but led to significantly higher reproductive output, below ground biomass and vegetative proliferation of new shoots (Palacios and Zimmerman, 2007). Similarly, short-term experiments have also shown increased photosynthetic rate and shoot productivity for Z. marina (Pajusalu et al., 2016;Zimmerman et al., 1997), Thalassia hemprichii (Jian et al., 2010) and Zostera muelleri (Collier et al., 2018) and increased community production for Zostera noltei (Mishra et al., 2018). Contrastingly, experiments on C. serrulata have shown no enhancement in productivity at higher CO2 (Schwarz et al., 2000;Collier et al., 2018). However, recent studies on three tropical (C. serrulata, T. hemprichii and Z, muelleri) and nine temperate (including Z. polychlamys) seagrass species showed significant increase in net productivity with increased CO2 (Ow et al., 2015;Borum et al., 2016;Collier et al., 2018).
Studies on natural CO2 seeps provide a unique opportunity to assess long-term effects of CO2, which are not possible under laboratory or in-field-controlled conditions. Studies using seeps, suggest that seagrass species can adapt to survive and enhance metabolic processes under elevated CO2 conditions (Hall-Spencer et al., 2008;Fabricius et al., 2011;Russel et al., 2013 andTakahashi et al., 2015). The photosynthetic activity of C. nodosa was stimulated near natural CO2 seeps of Vulcano, with significant increase in leaf chlorophyll-a content, maximum electron transport rate and compensation irradiance (Apostolaki et al., 2014;Olivé et al., 2017). However, this increase in efficiency of leaf metabolic process was not translated to generate higher biomass and production of the plant, probably due to nutrient limitation and grazing (Apostolaki et al., 2014) or toxic effects of seeps (Olivé et al., 2017). On the other hand, the biomass and net primary production of Halophila ovalis and C. serrulata were higher near the CO2 seeps, whereas the abundance of species increased only for C. serrulata, suggesting a species-specific response to elevated CO2 (Russell et al., 2013;Borum et al., 2016). Even though CO2 seep sites have been widely used to assess the long-term effects of elevated CO2 on benthic marine ecosystems and respective underlying mechanisms (Hall-Spencer et al., 2008;Fabricius et al., 2011;Enochs et al., 2015), care must be taken on the interpretation of results as other cofounding factors may be present, such as the emissions of heavy metals (Dando et al., 1999;Vizzini et al., 2013;Mishra et al., 2020) and sulphide (Dando et al., 1999;Capaccioni et al., 2001) which may influence the plants and population responses to elevated CO2. These environmental abiotic parameters are largely variable among seeps, whose emissions are largely dominated by CO2 (Dekov and Savelli, 2004;Varnavas and Cronan, 2005) Much of the natural CO2 seeps of Europe are concentrated in the Mediterranean Sea, in particular in the shallow waters of the active volcanic arcs in the Aegean Sea (Dando et al., 1999). These seeps provide future oceanic conditions (Hall-Spencer et al., 2008;Hall-Spencer and Rodolfo-Metalpa, 2009) with abundant seagrass at some seeps (Hall-Spencer et al., 2008;Russel et al., 2013) and sparse at others (Vizzini et al., 2013). Although these vents have attracted the relative research community over the recent past, it has never been investigated so far how these conditions will affect higher levels of seagrass biological organization other than photosynthesis and production, such as population dynamics.
Seagrass are rhizomatous plants that grow by reiteration of a limited set of modules (Borum et al., 2004). Their past growth history can therefore be reconstructed from the distinct scars left by abscised leaves on the long-lived rhizomes or the seasonal signals imprinted in the frequency and size of their modules (Short and Duarte, 2001), therefore allowing the detection and age determination of past growth (Duarte et al., 1994).Assessing population dynamics therefore allows for an in depth understanding of the response of the seagrass meadow to environmental perturbations over time.
The objective of this work is to assess for the first time the long-term responses of the meadow structure, plant growth and production, population age structure and dynamics of the seagrass C. nodosa under elevated CO2 levels near natural seeps. The seagrass C. nodosa is an endemic species that supports highly complex and biodiverse climax communities in Mediterranean Sea (Mazzella et al., 1986;Borum et al., 2004;OSPAR, 2010). We hypothesized that increased CO2 enhances plant growth and production and higher shoot turnover based on higher rates of shoot recruitment and mortality. To test this hypothesis, C. nodosa population in the vicinity and away from the influence of volcanic seeps were compared. To handle possible confounding factors of the effects of CO2 on the population dynamics of C. nodosa, we replicated the sampling effort in three seeps, two at the island of Milos (Adamas and Paleochori bay) in Greece and one at Vulcano island (Levante bay) in Italy and considered the shared responses to the three seeps as effects of elevated CO2.

Vulcano, Italy
The third location was Levante Bay (38.40° N, 15.00° E) of Vulcano island, Italy (Fig.  1B), a well characterized area to study natural processes associated to increased CO2 levels (Boatta et al., 2013). The main underwater gas seeps are located along southern and western shores of the bay at <1 m depth (Boatta et al., 2013). The underwater gas emissions are 97-98% CO2 with 2.2% H2S close to the seeps, decreasing to less than 0.005% H2S towards the north-eastern part of the bay, where most ocean acidification research has been located (Capaccioni et al., 2001;Milazzo et al., 2014). There is a step gradient in carbonate chemistry with pH 5.65 at the main gas seeps increasing to pH 8.1, which is typical for present day Mediterranean surface seawater (Boatta et al., 2013;Howritz et al., 2015). Subtidal C. nodosa meadows dominate the northern section of the bay (Giaccone, 1969;Boatta et al., 2013). At all three locations, two sites were established a high CO2 site near the seeps (called 'seep' hereafter) and a reference site (called 'reference' hereafter) away from the influence of the vent. All sites were established at similar depths (ca. <5m).

Water sampling:
Replicated water samples (n=5) were collected at the CO2 seeps and reference sites in 100-mL Winkler bottles, in situ fixed with 20 µL mercury chloride, stored in the dark cool boxes and transported to the laboratory for total alkalinity (AT) analysis. The pHTS and temperature of the water samples were measured in the field immediately after collection using a pH meter (pH meter, Titrino Methron). In the laboratory, AT was analysed in an 80-mL sample aliquot by titration (Lab Titrino analyser, Dickson et al. 2007). We used sterilized sea water as reference materials (CRM Batch 129, accuracy-98.7%, Dickson, 2013) for AT analysis. Temperature, pHTS and AT data were used to calculate the CO2 partial pressure (pCO2) and total dissolved inorganic carbon (DIC) using the CO2SYS program developed by Pierrot and Wallace, (2006). Dissociation constants (K1 and K2) developed by Mehrbach et al., (1973) and refitted by Dickson and Miller (1987) and dissociation constant for boric acid (KB) developed by Dickson et al. (2007b) were used.

Cymodocea nodosa traits:
C. nodosa shoots (n=10) samples were collected from each site using a core (20 cm diameter in Levante bay and 15 cm in Adamas and Paleochori bay) into a depth of about 30 cm in May 2013 (Levante bay) and May 2014 (Adamas and Paleochori bay). The sediment was carefully rinsed off in the field with seawater to prevent shoot dislodgement and to keep the rhizome mat intact as required for the reconstruction techniques (Duarte et al., 1994;Fourqurean et al., 2003). In the laboratory, the seagrass samples were washed again with distilled water to remove any debris and leaf epiphytes were scrapped with a plastic razor. Total and apical shoot densities (shoot m -2 ) were determined by counting the number of shoots per sample. The biomass of leaves, rhizomes (vertical and horizontal) and roots were obtained by weighing the separated tissues after oven-dried for 48 h at 60º C (dry weight, dw), and later used for calculation of production and biomass (g dw m -2 ). Above-ground (AB; leaves) to below-ground (BG; rhizomes +roots) biomass ratios (AB:BG) were also calculated. Dried plant material (leaves, rhizomes and roots) was then ground and analysed for total C and N content (% dw) in a CHN analyser (EA 1110 Model, Elemental Microanalysis Ltd, Oakhampton, Devon, UK) and C:N ratio was calculated on dw basis.
Reconstruction techniques an indirect measurement of plant growth was used to estimate the rhizome growth and production, population age structure and derived population dynamics (Duarte et al., 1994). The shoot age in plastochrone intervals (PI, the time required between the production of two consecutive new leaves) was estimated by counting the numbers of leaf scars in vertical rhizomes, plus the number of standing leaves in that shoot. The PIs are very species-specific and remarkably constant for some species as it is in the case of C. nodosa, where it varies from 26 to 33 days (Duarte et al., 1994;Short and Duarte,2001). Here we considered a PI of 29.3 days, which was estimated by Cancemi et al., (2002) for C. nodosa Mediterranean population of Ischia, Naples bay, Italy. Shoot ages were thus obtained by multiplying the number of PIs of shoots by 29.3 days. Frequency distribution of shoot ages for each population was constructed to reveal age structure and maximum longevities, and to calculate the population dynamics parameters (see below).
The vertical rhizome elongation rate (cm yr -1 ) for each site was obtained as the slope of the linear regression between the length of vertical rhizomes (from the insertion onto the horizontal rhizome to the apical meristem) and their age. Horizontal rhizome elongation rate was estimated based on the linear regression of the number of horizontal internodes in between two consecutive shoots against the age difference of those shoots. The slope (internodes yr -1 ) was then multiplied by the mean internode length (cm internode -1 ) to obtain the horizontal rhizome elongation in units cm yr -1 . The regression analyses were based on the totality of the rhizomes collected in each site with the cores plus few hand-picked horizontal rhizomes. This was not possible at Levante bay where the horizontal rhizome elongation rate was estimated based on the frequency distribution of all the measurements of the number of horizontal internodes produced per unit of time (age difference between two consecutive shoots). The horizontal growth rate is then estimated as the average of the nodes of the frequency distribution obtained for each site. With this method, it was not possible to estimate the error for the horizontal rhizome elongation rate for the Levante bay sites.
The rhizome production (g C m -2 yr -1 ) at each site was calculated as the product of the rhizome elongation rate (cm yr -1 ), the number of growing rhizomes apices per unit area (rhizome m -2 , given by non-apical and apical shoot densities for vertical and horizontal productions, respectively), the specific dry weight of rhizomes (g dw cm -1 ) and the rhizome C content (g C g -1 dw).
The long-term average shoot recruitment (R) was estimated from the shoot age structure using a negative decay model: Nx = N0 e -Rx , where Nx is the number of shoots in age class x, N0 is the number of shoots recruited into the population. This model assumes that mortality and recruitment have remained constant over the lifespan of the oldest shoots with year to year random variation around some mean value of mortality and recruitment (Fourqurean et al. 2003, Cunha andDuarte, 2005). The recruitment for the current year of sampling, i.e. present shoot recruitment rate (Ro, in units yr -1 ) was estimated as the difference of the total number of shoots and the number of shoots older than a year in the shoot population (Duarte et al., 1994). The population growth rate was given by the net shoot recruitment rate (r), estimated as: r = R0 -M, where M is the long-term mortality rate, which equals the long-term recruitment rate (R) under the assumptions of near steady state (Fourqurean et al., 2003). Population was considered growing if r is positive (R0 > R), shrinking if r is negative (R0 < R), or with the same trajectory pattern if R0 is not significantly different from R (Fourqurean et al., 2003).

Statistical analysis:
Data are shown as mean ± 1 standard error of the mean. Significant differences in seagrass traits, shoot density, biomass and rhizome production rates among sites ( factor with 2 levels: CO2 seeps and reference) and locations ( factors with three levels: Adamas, Paleochori bay and Levante bay) were investigated using two-way ANOVA (results given by the p-value associated to each fixed factor), after testing for homogeneity of variances (Fligner test) and normality of distribution (Shaphiro-Wilks test). The Tukey's multiple comparison test was applied to determine significant differences among factor levels. When ANOVA assumptions were not verified even after root-squared or Log10 transformation, comparison between sites and locations were performed by non-parametric Kruskal-Wallis rank sum tests (results given be the p-value associated to each fixed factor) and the post hoc Dunn's test. Differences in the vertical and horizontal rhizome elongation rates between reference and CO2 sites at each location were tested by ANCOVA analysis, in which the rhizome trait (rhizome length for vertical rhizome rate and number of internodes for horizontal one) was fit to rhizome age using site and its interaction as fixed factor and their slopes compared. The confidence limits of the exponent coefficient from exponential decay regression model used to estimate the long-term average recruitment rate (R ) were used to test if it was different from the present recruitment rate (Ro; Fourqurean et al., 2003). Significance levels was considered at p < 0.05.

Seawater carbonate chemistry
Seawater carbonate chemistry varied between CO2 seeps and reference sites in a similar way for all locations, being pCO2 1.7 to 6-fold higher and pH 0.2 to 0.7 units lower at the seep sites than at reference sites (Table 1). DIC presented similar values at all three seep locations. Much higher CO2 concentration in the water column was observed at Adamas than at the other two locations (Table 1).

Cymodocea nodosa traits
The C content (35.1 ± 0.2 % vs 32.1 ± 0.7 %) at the seeps were higher than the reference site of Adamas and significantly different (Table 2), as expected due to the higher CO2 availability (Table 1). At Paleochori bay, the rhizome N content near the seep was 2.3fold lower than at the reference sites and consequently the C:N ratio was 2-fold higher at the seeps (Table 2). In general, the C and N contents of C. nodosa varied significantly among locations. The C content of C. nodosa roots were higher at Adamas than of the other locations, whereas the C and N contents of leaves (and rhizome N) were higher at the seeps of Levante bay.
C. nodosa total and apical shoot density were higher at the CO2 seeps than at the reference sites (Fig.2). Whereas total shoot density did not differ among locations, apical shoot density did, with higher densities at Adamas (Fig. 2B). C. nodosa biomass was not affected by seeps (Table 3). Meadows at Levante bay exhibited 4-fold higher biomass than those at the other locations. AB:BG differed among locations, being higher at Paleochori bay (Table 3).
The vertical rhizome elongation rate of C. nodosa was significantly higher at the CO2 seep than at the reference site in all locations and ranged from 0.58 ± 0.02 cm yr -1 (reference site at Paleochori bay) to 1.85± 0.10 cm yr -1 (CO2 seep at Levante bay; Fig. 3A). The horizonal rhizome elongation rate was only higher at the CO2 seep site than at the reference site at Paleochori bay (Fig. 3B). C. nodosa population at Adamas presented the lowest horizontal elongation rates (<5cm yr -1 ), whereas the highest were found at Levante bay (> 10 cm yr -1 ; Fig.  3B). The vertical rhizome production rate was 2.1 to 2.7-fold higher and horizontal rhizome production was 1.5 to 2.3-fold higher at CO2 seeps than at reference sites ( Fig.3C and D).
C. nodosa shoot longevity was higher at Levante bay (7-8 years), whereas at Adamas and Paleochori bay maximum shoot ages reached about 3 years. No differences between CO2 seep and reference sites were patent (Table 4). In contrast, the present and long-term recruitment rates were generally higher in populations growing at the CO2 seeps than those at the reference sites (Table 4). CO2 seeps also affected the frequency distribution of C. nodosa shoots, which consistently showed a higher number of younger plants (1 and 2 years old) at the seeps, whereas the number of older plants (more than 3 years) were higher at the reference sites in all locations (Fig.4). The net population growth rates were generally around zero, indicating stable population except at Paleochori bay reference site, where the population was actively growing (growth rate=0.3, Table 4).

Discussion:
Shallow water CO2 seeps may act as natural analogues of future oceanic scenarios. The presence of seagrass at the CO2 seeps of Adamas and Paleochori bay of Milos island, Greece and of Levante bay at Vulcano island, Italy, provided an opportunity to test the effects of longterm CO2 enrichment on seagrass population dynamics. Confirming our initial hypothesis, we observed a positive response of C. nodosa population to high CO2 at the three seep sites, with higher total and apical shoot density, and higher vertical and horizontal rhizome elongation and production rates than at the reference sites.
Higher shoot densities near seeps suggests that elevated CO2 promotes the differentiation of vertical rhizomes from the horizontal rhizome meristem. It is not clear how CO2 may influence the differentiation of vertical rhizomes, but the stimulation of meristem activity and tissue-specific responses under elevated CO2 has been well described in terrestrial grasses (Kinsman et al., 1997;Masle, 2000). The positive influence of high CO2 on the density of C. nodosa individuals is probably related to increased reproductive output as observed by Palacios and Zimmerman (2007) in Z. marina.
At the Levante bay, where more N is available as indicated by the general higher concentration of N in the plant tissues, both vertical and horizontal rhizome grew faster. C. nodosa at this location showed also higher biomass, lower AB:BG and higher shoot longevity, indicating more favourable growth conditions. At Adamas, where the N content in leaves, rhizomes and roots were lower, the highest CO2 availability was not reflected at the growth level. The N content of C. nodosa leaves at Adamas and Paleochori bay were very low, below the threshold limit of 1.8%, that indicates nitrogen limitation in seagrasses (Duarte, 1990). Nutrient limitation of seagrass growth under high CO2 has been often invoked (Alexandre et al., 2012;Apostolaki et al., 2014) even though Campbell and Fourqurean (2018) failed to demonstrate the influence of nutrient enrichment on seagrass response to elevated CO2. The response of seagrasses to the synergistic effects of CO2 and nutrients are complex and often mediated by epiphyte growth (Martí nez-Crego et al., 2014). Negative effects of high N content due to epiphyte shading has been observed for P. oceanica meadows at Ischia CO2 seeps of Italy (Ravaglioli et al., 2017).
The density of C. nodosa increased under the influence of high CO2 at all three seep locations. However, contrasting evidence of no significant differences in density of C. nodosa at the CO2 seeps of Levante bay, Vulcano has been reported by Apostolaki et al., (2014) and Vizzini et al., (2019). In other species, significant positive effects on density near the seeps were found, such as for P. oceanica at Ischia (Hall-Spencer et al., 2008), C. rotundata at Papua New Guinea (Fabricius et al., 2011;Takahashi et al., 2015) and C. serrulata at Papua New Guinea (Russel et al., 2013). Even though we did not find any significant site effects on C. nodosa biomass, as well as Apostolaki et al. (2014), positive effects of seeps have been documented on other Cymodocea spp. (Fabricius et al., 2011;Russel et al., 2013;Takahashi et al., 2015). The effects of seeps are less clear in the biomass distribution between above-and below-ground plant parts. We did not find any significant differences, contrasting with Fabricius et al. (2011) andRussel et al. (2013), who found increase of below-ground biomass of C. rotundata and C. serrulata, respectively, near CO2 seeps. This was also observed in the one-year long, experimental CO2 enrichment study of Z. marina (Palacios and Zimmerman, 2007). The vertical and horizontal elongation rates of C. nodosa varied from 0.59 to 1.85 cm yr -1 and 3.9 to 12.8 cm yr -1 respectively. The vertical elongations observed lie around the average of those reported elsewhere (1.4 cm yr -1 , Marbà and Duarte, 1998), whereas the horizontal elongations are much lower than the average value reported in Marbà and Duarte, (1998), 40 cm yr -1 . The horizontal elongation of C. nodosa varies enormously, from 7 to 204 cm yr-1, being the higher values typical of the "runners", i.e., the long horizontal rhizomes located in the border of the meadows, which are exploring new grounds. Our observations, were made in the middle of the meadows and thus are consistent with the lower values of variation.
To our best knowledge, this is the first report of the effects of increased CO2 on the population dynamics of a seagrass. Observations confirmed the initial hypothesis and showed that populations near the seeps in all three locations had higher short-and long-term recruitment and growth rates around zero, indicating that the turnover of C. nodosa shoots under high CO2 will be higher. On the other hand, the shoot longevity was more positively affected by nutrient availability than by elevated CO2. The maximum shoot longevity of 8.3 years observed at Levante bay reference site, which was higher than previously reported (7.6 yr at Ria Formosa lagoon, Cunha and Duarte, 2005;4 yr at Ria Formosa, Cabaç o et al., 2010) In conclusion, our findings show that under long-term CO2 enrichment an increase of C. nodosa densities (both shoot and apical), an increase of elongation rates and production of rhizomes (both vertical and horizontal), a faster turnover of shoots, but no changes in the population biomass are expected. However, nitrogen effects on rhizome elongation and production rates and on biomass were stronger than CO2. Our results provide a picture of putative future alterations of the population structure and dynamics of C. nodosa meadows. We emphasize the importance of using replicated CO2 seeps to capture the common responses of seagrasses to long-term CO2 enrichment and its synergistic effects with the environmental fluctuation among seeps.

Acknowledgement:
This work was funded by Erasmus Mundus MARES program coordinated by Ghent University (FPA 2011-0016). This study received Portuguese national funds from the FCT-Foundation for Science and Technology through project UID/04326/2020 granted to CCMAR, and the fellowship SFRH/BPD/119344/2016 granted to C.B.d.1.S.We are grateful for the support of HCMR, Greece, particularly to Thanos Dailianis and Julius Glampedakis in field sampling and Veronica Santinelli in laboratory work. We are thankful to Andrew Tonkin and Robert Clough at Plymouth University, UK for helping in laboratory analysis.

Figure 2.
Cymodocea nodosa total and apical shoot density (shoots m -2 , mean ± standard error) at the reference and CO2 seep sites at Adamas, Paleochori bay and Levante bay locations. Summary of the 2-way ANOVA, showing the level of statistical significance (***p < 0.001, **p < 0.01, *p < 0.05) of the location, site and their interactions as fixed factors. Uppercase letters on bars represent post hoc Tukey pairwise groupings that indicate differences when the effect of location factor was significant.

Figure 3.
Cymodocea nodosa vertical and horizontal rhizome elongation rate (cm yr -1 ) and production rates (g C m -2 yr -1 ; mean ± standard error) at the reference and CO2 seep sites at locations Adamas, Paleochori bay and Levante bay. It was not possible to estimate the error for the horizontal rhizome elongation rate (see Methods). Asterisks indicate the level of statistical significance (***p < 0.001, **p < 0.01, *p < 0.05) of the site factor in the ANCOVA analysis for the elongation rates, and of the site, location and site x location interaction factor in the 2-way ANOVA for the production rates.