Isotopic and elemental compositions reveal density‐dependent nutrition pathways in a population of mixotrophic jellyfish

Mixotrophic organisms are increasingly recognized as important components of ecosystems, but the factors controlling their nutrition pathways (in particular their autotrophy–heterotrophy balance) are little known. Both autotrophy and heterotrophy are expected to respond to density-dependent mechanisms but not necessarily in the same direction and/or strength. We hypothesize that the autotrophy– heterotrophy balance of mixotrophic organisms might therefore be a function of population densities. To investigate this relationship, we sampled mixotrophic jellyfish holobionts (host,Mastigias papua etpisoni; symbiont, Cladocopium sp.) in a marine lake (Palau, Micronesia) on six occasions (from 2010 to 2018). Over this period, population densities varied ~100 fold. We characterized the nutrition of the holobionts using the δC and δN signatures as well as C:N ratios. δC values increased and δN values decreased with increasing population densities (respectively, R = 0.86 and 0.70, P < 0.05). Although less distinct, C:N ratios increased with increasing population densities (R = 0.59, 0.1 > P > 0.05). This indicates that the autotrophy–heterotrophy balance tends toward autotrophy when population densities increase. We propose that the availability of zooplanktonic prey is the main driver of this pattern. These results demonstrate that the autotrophy– heterotrophy balance of mixotrophic jellyfishes can be tightly regulated by density-dependent mechanisms.


INTRODUCTION
Density-dependent processes are key mechanisms that have profound impacts on the biology and ecology of organisms (e.g., May et al. 1974). One central aspect of organisms that can be affected by density dependence is their nutrition. This is true for both heterotrophs (e.g., increased competition for prey, Svanbäck and Bolnick 2007;disruption of feeding, Brown et al. 1994, Kelly et al. 1996 and autotrophs (e.g., increased shading or increased competition for nutrients, Stevenson et al. 1991, Blundell and Peart 2004, Sjøtun et al. 2006. Such density-dependent impacts on organisms' nutrition can then have important ecological consequences for population (e.g., through impact on recruitment, Sjøtun et al. 2006) or food-web dynamics (e.g., use of alternative prey; Svanbäck and Bolnick 2007).
Density-dependent effects on nutrition have been well documented in heterotrophs and autotrophs. However, there is substantial lack of knowledge about density-dependent dynamics of mixotrophic organisms. Mixotrophy, can be defined as a nutrition mode where an organism, or an holobiont (i.e., symbionts + host), produces organic matter through autotrophy (photosynthesis or chemosynthesis) and acquire it directly from its environment (e.g., osmotrophy, phagotrophy). This nutrition strategy is ubiquitous in both aquatic and terrestrial ecosystems (Selosse et al. 2017). Well known are coral reefs where the mixotrophy of the coral holobiont allows for recycling of nutrients, which supports rich and complex ecosystems in otherwise nutrient-poor, oligotrophic, environments (Muscatine and Porter 1977). More generally, mixotrophy has critical roles in all ecosystems from terrestrial plant nutrition (Schmidt et al. 2013) to planktonic mixotrophic protists (Stoecker et al. 2017).
The contribution of autotrophic and heterotrophic nutrition to total nutrition (i.e., for carbon, nitrogen, or other compounds) in mixotrophic organisms is variable, representing a full spectrum of nutrition strategies from dominantly autotrophic to dominantly heterotrophic. A given mixotrophic species might move along this spectrum as a function of environmental conditions and/or biotic interactions (Wilken et al. 2013, Fox et al. 2019. Then, the position of a mixotrophic organism along the autotrophyto-heterotrophy spectrum has important implications such as its position in the food web, its contribution to primary productivity, and thereby for biogeochemical cycling (see Stoecker et al. 2017). Therefore, to understand the ecological roles of mixotrophs in ecosystem functioning, it is of crucial importance to identify and characterize the biological and ecological factors that determine their nutrition pathways (Selosse et al. 2017).
Evidence for the importance of density-dependent processes for mixotrophic organisms is known from studies on corals, where such effects are related to the availability of suitable space and the settlement of recruits (e.g., Vermeij andSandin 2008, Doropoulos et al. 2017). However, less is known on how density-dependent processes impact the access to nutritive resources of mixotrophic organisms (but see Tanner 2002). It is important to stress here that predicting the effect of density dependence on the nutrition of mixotrophic organisms is not trivial a priori. Indeed, autotrophic and heterotrophic nutrition pathways can be affected in a different direction or strength by density-dependent processes. Therefore, density-dependent processes could ultimately lead to a shift in the autotrophyheterotrophy balance of mixotrophic organisms.
Here, we study the nutrition of the mixotrophic, zooxanthellate, jellyfish Mastigias papua etpisoni holobiont (host, Mastigias papua etpisoni, Mastigiidae, Rhizostomeae, Scyphozoa, see Dawson 2005, Souza andDawson 2018; symbiont, Cladocopium sp., Symbiodiniaceae, Suessiales, Dinophyceae, see LaJeunesse in Vega de Luna et al. 2019) in Ongeim'l Tketau a marine lake from Palau (Micronesia, Western Pacific Ocean; see Hamner et al. 1982). These marine lakes, especially Ongeim'l Tketau with its jellyfish population, offer an opportunity as model systems for the investigation of densitydependent dynamics impacting the nutrition of mixotrophic organisms. Zooxanthellate jellyfishes in general, including Mastigias, can show important variability in their nutrition (reviewed in Djeghri et al. 2019). The population of medusae in Ongeim'l Tketau can be large (in the order of tens of individuals per m 2 , Hamner et al. 1982, Cimino et al. 2018 and are a key driver of community dynamics (Wollrab et al. 2019). The medusae are normally present year-round and recruit more or less continuously (Hamner et al. 1982). However, the population fluctuates dramatically in relation to El Niño events (Dawson et al. 2001, Martin et al. 2006. Mastigias papua medusae in Palau's marine lakes are known to be able to exert strong predatory pressure on copepods (Wollrab et al. 2019) and the simple planktonic food web of Ongeim'l Tketau (only two small, dominantly herbivorous, species of copepods: Bestiolina similis and Oithona dissimilis, Hamner et al. 1982, Saitoh et al. 2011, see Benedetti et al. 2016 for copepod sizes and trophic regimes) limits the possibilities of diet switching. Finally, due to its enclosed nature, Ongeim'l Tketau is less influenced by dispersal through advection than more open systems (e.g., coastal areas), making density estimates much more accurate.
To investigate the potential link between population density and mixotrophy, we sampled Mastigias papua etpisoni medusae from Ongeim'l Tketau over six different years (between 2010 and 2018) including large variations in population densities. We characterized medusae nutrition through the use of their stable isotope signatures (δ 13 C and δ 15 N) and C:N ratios. In photosymbiotic organisms such as Mastigias papua etpisoni holobionts, δ 13 C values tend to increase with the relative contribution of autotrophy and decrease with the relative contribution of heterotrophy, whereas the reverse is true of δ 15 N (reviewed in Ferrier-Pagès and Leal 2018, see Djeghri et al. 2020 for an application to zooxanthellate jellyfishes). Conversely, C:N ratios tend to increase with the relative contribution of autotrophy and decrease with the relative contribution of heterotrophy (Djeghri et al. 2020 and references therein). We therefore hypothesize that density-dependent process would affect the mixotrophic nutrition of the Mastigias papua etpisoni holobiont leading to changes in its autotrophy-heterotrophy balance. We expect to detect these changes by observing coherent variations of the δ 13 C, δ 15 N, and C:N ratios of the medusae as a function of their population densities.

Elemental and isotopic composition of medusae and copepods
Mastigias papua etpisoni medusae were sampled by hand in Ongeim'l Tketau during six sampling missions in August-September 2010, July-August 2011, August 2013, April 2015, April 2016, and August-September 2018. The medusae were brought back alive in the laboratory and left for 2-3 h in clean seawater to allow for the removal of prey from their oral arms. Then, a piece of oral arm (including symbionts, therefore representative of the holobiont) was dissected and oven-dried for 24 h at 60°C. The dried samples were then ground to a powder and inserted in tin capsules for elemental and isotopic analyses. The number of individuals sampled varied from one year to another (n = 12, 5, 5, 10, 6, and 18 for, respectively, 2010, 2011, 2013, 2015, 2016, and 2018). The carbon and nitrogen mean masses of the medusae samples analyzed were of, respectively, 147.0 AE 83.5 µg and 28.1 AE 17.3 µg (mean AE SD). Additionally, zooplankton samples were collected during each sampling using different plankton nets as a function of material availability and years (mesh sizes comprised between 63 and 200 µm). Copepods were then isolated using forceps under a dissecting microscope and inserted in tin capsules. Importantly, this isolation allowed for the correction of potential sampling bias induced by the use of different plankton nets. The tin capsules containing copepods were then oven-dried for 24 h at 60°C. The number of samples analyzed varied from one year to another (n = 2 for 2010, 2011, and 2016, n = 1 for 2013 and n = 3 for 2015 and 2018). The carbon and nitrogen mean masses of the Elemental and isotopic analyses of the medusae samples from 2018 were analyzed using an elemental analyzer (Thermo Scientific EA Flash 2000) coupled to a gas isotope ratio mass spectrometer (Thermo Scientific DELTA V Plus) at the Stable Isotopes Facility of the Pôle Spectrométrie Océan (PSO-IUEM, Plouzané, France). All other samples were analyzed by a high-sensitivity elemental analyzer (Hansen and Sommer 2007) connected to a gas isotope ratio mass spectrometer (DELTAplus Advantage; Thermo Fisher Scientific, Waltham, Massachusetts, USA) at GEO-MAR, Helmholtz Centre for Ocean Research, Kiel, Germany. The high sensitivity of the GEO-MAR mass spectrometer allowed for the accurate estimation of the isotopic signatures of low-mass samples (most notably, the copepod samples; see Hansen and Sommer 2007). Stable isotope values are expressed as permil (‰) using the δ notation (normalized to Vienna Pee Dee Belemnite and atmospheric N 2 for, respectively, carbon and nitrogen): where X is the element measured, X H is the amount of heavy isotope, and X L is the amount of light isotope from the samples (X sample ) and standard (X std ). An internal standard-casein (Analysentechnik, Germany) in the Pôle Spectrométrie Océan, and acetanilide (Merck, Germany) in GEOMAR-was measured every sixth sample within each sample run. The overall standard deviations measured in the Pôle Spectrométrie Océan were of AE0.12‰ and AE0.24‰ for, respectively, nitrogen (N) and carbon (C). At GEOMAR, the overall standard deviation for the low measurement range 2.5-8.0 µg N and 5.0-80 µg C was AE0.25‰ and AE0.2‰, respectively. The overall standard deviation for the higher measurement range 3-15 µg N and 10.0-140 µg C was AE0.2‰ and AE0.15‰, respectively. To make the different years comparable, the δ 13 C and δ 15 N values of copepods were subtracted from the δ 13 C and δ 15 N values of medusae of the corresponding year, thus correcting for the potential shifts in the isotopic baseline (hereafter referred to as "corrected δ 13 C" and "corrected δ 15 N"). C:N values were expressed as mass ratios.

Monitoring of Ongeim'l Tketau and of the medusae population
The Mastigias papua etpisoni population density in Ongeim'l Tketau is evaluated monthly by the Coral Reef Research Foundation (CRRF) through net sampling (15 sampling points, replicated three times during the day, using a 50 cm diameter and 1-mm mesh size zooplankton net, Martin et al. 2006, Cimino et al. 2018). In addition to population density, the medusae are measured (umbrella diameter to the nearest cm) resolving the size structure of the population. The population data used here are the evaluation the closest in time to the sampling date for isotopic and elemental composition.
Additionally, the temperature in Ongeim'l Tketau is monitored monthly at two sites (Eastern and Western basins) at 1-m depth intervals between 0 and 12 m (plus one measure at 0.5 m since April 2016). These temperatures were averaged representing the mean mixolimnion temperature (following Martin et al. 2006).

Statistics
Previous experimental and observational studies (Ferrier-Pagès and Leal 2018, Djeghri et al. 2020 and references therein) suggest that if δ 13 C, δ 15 N, and C:N ratios of photosymbiotic organisms change in response to a shift in the autotrophy-heterotrophy balance some specific correlations patterns should be seen: δ 13 C and C:N ratios should be positively correlated with each other but negatively correlated with δ 15 N. This was tested by computing the pairwise Pearson's correlation coefficients between the corrected δ 13 C, corrected δ 15 N, and C:N ratios of the medusae.
Then, the relationship between Mastigias papua etpisoni population density and their corrected δ 13 C, corrected δ 15 N, and C:N ratios was tested by performing linear regressions using the mean isotopic or elemental value of each year as the variable and log 10 -transformed population density as the covariate. As sampling effort varied from year to year, each year was weighted by the number of medusae sampled for isotopic and elemental measures to improve confidence estimates. Additionally, a possible relationship v www.esajournals.org between temperature and the corrected δ 13 C, corrected δ 15 N, and C:N ratios of the medusae was investigated using the same protocol.

RESULTS AND DISCUSSION
In this work, we hypothesize that densitydependent processes could affect the nutrition of the present model organism, the Mastigias papua etpisoni medusae holobiont. Particularly, we expect that the strength and direction of the effect of density-dependent processes on autotrophy and heterotrophy would not necessarily be similar leading to changes in the autotrophyheterotrophy balance of the medusae.
To characterize the nutrition of Mastigias papua etpisoni holobiont in Ongeim'l Tketau, we used δ 13 C and δ 15 N signatures along with C:N ratios. These three indicators can provide information on the relative importance of autotrophy and heterotrophy in photosymbiotic organisms (see, e.g., Ferrier-Pagès and Leal 2018, Djeghri et al. 2020) but can also be affected by other events such as a switch in preferred prey or nutrient. Such events may act as confounding factors possibly impeding the interpretation of our indicators as representatives of the autotrophyheterotrophy balance of medusae. However, if the δ 13 C, δ 15 N, and C:N ratios represent well the autotrophy-heterotrophy balance of a photosymbiotic organism, a precise correlation pattern is expected: The δ 13 C and C:N ratios should be positively correlated with each other while being negatively correlated with the δ 15 N (see Ferrier-Pagès and Leal 2018, Djeghri et al. 2020 and references therein). This is what is observed here ( Table 2). It is unlikely that this correlation pattern would have been observed if the isotopic and elemental compositions of the medusae were dominated by confounding factors. For instance, changes in preferred prey may have changed the δ 13 C (e.g., prey from a different origin) and the δ 15 N (e.g., prey from a different trophic level). However, such shift must not necessarily result in a negative correlation between δ 13 C and δ 15 N signatures. It is also unclear how such a change would have impacted the C:N ratios of the medusae. In contrast, a change in the autotrophy-heterotrophy balance of the medusae provides a simple and well-proofed explanation to this correlation pattern. Although this does not exclude that additional confounding factors such as prey switching may have played a role, this suggests that their role is probably minor.
When plotted against population density data of Mastigias papua etpisoni in Ongeim'l Tketau (Fig. 1), the different indicators used are congruent. Corrected δ 13 C values increased significantly with log 10 -transformed population density (R 2 = 0.86, P < 0.01; Fig. 1a), whereas corrected δ 15 N values decreased significantly with log 10 -transformed population density (R 2 = 0.70, P < 0.05; Fig. 1b). C:N ratios tended to increase with log 10 -transformed population density although this trend was not significant at α = 0.05 (i.e., 0.05 < P < 0.1; Fig. 1c). Therefore, the indicators used all suggest that the denser the Mastigias papua etpisoni population, the more the autotrophic its nutrition and vice versa. Variations in lake water temperature could potentially be a confounding factor. Indeed, mixotrophic organisms can tend toward a dominantly heterotrophic nutrition when temperatures increase (Wilken et al. 2013). For Mastigias papua etpisoni specifically, high temperatures can result in bleaching (expulsion of the zooxanthellae, Dawson et al. 2001). Additionally, the medusae population density in Ongeim'l Tketau is negatively correlated with temperatures (high temperatures correlate with lower population densities; Martin et al. 2006). Therefore, it could be expected that the variations seen in Mastigias papua etpisoni nutrition are due to temperature variations rather than density-dependent processes. However, we found no significant correlations between the corrected δ 13 C, corrected δ 15 N, and C:N ratios of Mastigias papua etpisoni and the variations of the mean mixolimnion temperature in Ongeim'l Tketau (R 2 < 0.05 and P > 0.7 in all cases; Appendix S1: Fig. S1). This suggests that temperature is not a strong driver of Mastigias papua etpisoni nutrition in this study as compared to density-dependent mechanisms that can indeed alter the autotrophy-heterotrophy balance of Mastigias papua etpisoni (Fig. 1). Two main mechanisms, not mutually exclusive, might shape the observed density-dependent patterns. On one hand, high population densities might favor autotrophy. On the other hand, low population densities might favor heterotrophy. An increase in the contribution of autotrophy with population density could be related to a better protection from light through self-shading potentially diminishing photoinhibition. Photoinhibition is indeed known to alter autotrophy in other mixotrophic organisms such as corals (e.g., Winters et al. 2003). However, Mastigias papua medusae are known to have complex behavior and can adjust for depth as a function of light intensity, suggesting that they can mitigate photoinhibition (Dawson and Hamner 2003). Moreover, and more generally, high population densities tend to have negative effects on photosynthesis through self-shading and/or competition for inorganic nutrients (e.g., Wang et al. 2005, Sjøtun et al. 2006). Therefore, a strong positive impact of low population densities on heterotrophy appears more likely as the driving factor of the observed patterns. The ultimate reason could be an increased prey availability at low medusae population densities. Mixotrophic medusae could decrease zooplankton densities through either bottom-up or top-down dynamics (see Pitt et al. 2009, West et al. 2009). Bottom-up processes result from mixotrophic medusae competing with phytoplankton for nutrients, thereby reducing available prey for zooplankton. Topdown processes result from mixotrophic medusae directly controlling zooplankton populations through predation. These two processes are not mutually exclusive and most probably work in parallel. In the specific case of Palau's marine lakes, however, top-down mechanisms, and cascading effects on food-web structure, appear to be the main driver (Wollrab et al. 2019). Furthermore, the body plan of medusae makes them very efficient zooplankton predators (Acuña et al. 2011). When population densities are high, medusae are able to reduce dramatically zooplanktonic prey densities sometimes resulting in density-dependent limitation of their own growth (e.g., Lucas 2001, Goldstein andRiisgård 2016). For Ongeim'l Tketau, it can be estimated that, when numerous (e.g., in 2013 or in 2015 with~100 ind./m 2 ), the medusae would clear~600% of the lake water column of its copepods each day (estimation derived from the size-based clearance rates given by García and Durbin 1993 for Phyllorhiza punctata medusae, Mastigiidae, and assuming a 14 m deep mixolimnion). This estimation represents a maximum possible clearance effect, assuming a homogeneous and overlapping distribution of medusae and copepods, which might not always be the  12, 5, 5, 10, 6, and 18 for, respectively, 2010, 2011, 2013, 2015, 2016, and 2018). Empty symbols are the individual data from each year; full symbols are the means of each sampling event. Dotted lines are 95% CI around regressions. case (Hamner et al. 1982, Dawson andHamner 2003). By contrast, when the medusae population is low (e.g., in 2016 or 2018, <10 ind./m 2 ) less than 10% of the water column would be cleared of its copepods each day. It is thus likely that more prey per individual would be available when the medusae population densities, and hence the predation pressures, are low.
It is also interesting to note that the mixotrophic nutrition of Mastigias papua etpisoni-and of mixotrophic organisms more generallymight help mitigate the negative effects of density-dependent processes. If either autotrophic or heterotrophic nutrition is particularly negatively impacted by density-dependent effects, the other nutrition mode might buffer or compensate such negative effects. For instance, scleractinian corals may compensate for the loss of autotrophic nutritive input due to bleaching by an increase in heterotrophic nutrition (Hughes and Grottoli 2013). Evidence for the opposite mechanism-that is, compensation of decreased heterotrophy through increased autotrophy-is, however, more ambiguous (see, e.g., Hoogenboom et al. 2010). Even without compensation, the existence of a second source of nutrition, unaffected, or less affected, by density-dependent processes, could act as a buffer against negative effects such as-as hypothesized here-increased competition for prey. In the specific case of Mastigias papua in Palau's marine lakes, this could help to explain the very large, and rather stable over years, population densities reached (in the order of tens of individuals per m 2 , Hamner et al. 1982, Cimino et al. 2018. In contrast, heterotrophic medusae may suffer more strongly from density-dependent processes reducing their body size (either by limiting growth or by shrinking) as a consequence of food limitation (Lucas 2001, Pitt et al. 2014. This is in line with the idea that the so-called R* (the amount of resources at which birth rate equals mortality, Tilman 1990) of mixotrophic organisms for prey would be substantially lower than the ones of comparable heterotrophic organisms (see, e.g., Rothhaupt 1996, Ptacnik 2016. To conclude, our results show that the contribution of autotrophy and heterotrophy to nutrition is significantly linked to population density in a directed and predictable way in Mastigias papua etpisoni medusae holobionts. There is no reason why similar dynamics should not operate in populations of other mixotrophic organisms, although this remains to be demonstrated. Such shifts in nutrition could in turn affect the role played by mixotrophic organisms in ecosystem functioning (see, e.g., Stoecker et al. 2017 for planktonic food webs). More generally, we argue that such generally overlooked density-dependent processes must be taken into account to increase our understanding of the roles and consequences of mixotrophic nutrition in aquatic and terrestrial ecosystems.

ACKNOWLEDGMENTS
We thank Pat Colin, Lori J Bell, Michael N Dawson, Mathew Mesubed, and Emilio Basilius for assistance in logistics and fieldwork. We particularly remember and are forever grateful for the invaluable contribution of Emilio Basilius  to research conducted in Palau's marine lakes-working in the marine lakes since 1995, and Emilio was always leading the trail with energy and enthusiasm. The Stable Isotope Facility of the "Pôle Spectrométrie Océan" (PSO, Plouzané, France) is supported by funding from the University of Brest (UBO), CNRS, IFREMER, and IRD. This work was supported by the "Laboratoire d'Excellence" LabexMER (ANR-10-LABX-19) and co-funded by a grant from the French government under the program "Investissements d'Avenir." We are grateful to two anonymous reviewers for their valuable comments who helped raising the quality of this manuscript. Finally, we thank Palau's National Bureau of Marine Resources and Koror State Government for permitting our research in Palau. v www.esajournals.org