**Author details**

Carolyn Mahan1 and Michael Steele2 \*

1 Department of Biology, The Pennsylvania State University-Altoona College, Altoona, USA

2 Department of Biology, Wilkes University, Wilkes Barre, USA

\*Address all correspondence to<EMAIL_ADDRESS>
© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


**1. Introduction**

The oil, gas, and maritime industry are significantly impacted by biofouling, either through the costs related to the increased fuel consumption, hull cleaning, material deterioration, repainting, and corrosion [1], or even by the management intervention time and incorrect measurements in the submerged and moored sensors [2]. Ecologically, fouling events in marine environments promote species invasion and the establishment of exotic biofouling species in ports [3]. Moreover, biofouling can be a major concern in health-related problems since contamination of aquaculture facilities, such as fish cages, can occur by toxin accumulation, and air pollution may be increased through greenhouse gas emissions [4].

The non-toxic marine antifouling approaches available in the market are often expensive and not as effective as conventional biocides, which can accumulate in the marine

**Citation:** Romeu, M.J.; Gomes, L.C.; Sousa-Cardoso, F.; Morais, J.; Vasconcelos, V.; Whitehead, K.A.; Pereira, M.F.R.; Soares, O.S.G.P.; Mergulhão, F.J. How do Graphene Composite Surfaces Affect the Development and Structure of Marine Cyanobacterial Biofilms? *Coatings* **2022**, *12*, 1775. https://doi.org/10.3390/ coatings12111775

Academic Editor: Mohor Mihelˇciˇc

Received: 20 October 2022 Accepted: 17 November 2022 Published: 20 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

environment and affect non-target aquatic organisms. Inorganic ingredient-based coatings have been used for biofouling prevention, including silver nanoparticles [5], carbon nanotubes (CNTs) [6], graphene [7], and metal oxides semiconductors such as zinc oxide (ZnO) [8] and titanium dioxide (TiO2) [9]. Overall, antifouling marine paints containing nanomaterials have been reported to offer superhydrophobicity, microbial resistance, high durability, water repellency, anti-sticking, and anti-corrosive properties [10], being novel solutions for the sustainable development of the maritime industry. Advancements in materials technology have introduced carbon nanomaterials such as CNTs and graphene as a powerful approach for various applications in the marine and shipping industries [11]. Graphene consists of a single layer of carbon atoms arranged in a sp2-bonded hexagonal pattern. It is considered one of the strongest and thinnest materials available, which shows high specific surface area, electrical conductivity, and thermal stability, making it appealing for different applications [12,13]. Moreover, due to its high strength level [14], this material is a breakthrough alternative in the naval industry. All these features make the application of graphene in technical processes of maritime industries attractive, such as in water management systems, desalination, removing toxic pollutants and filtering gasses, and as a coating material [15–17]. As a coating on marine structures, besides the anti-corrosive properties, graphene can also be used in de-icing surfaces for ship operations in extremely low-temperature regions, such as the Arctic and Antarctica, due to its electrical conductivity [15].

Nanotechnology-based technologies can be of great interest in creating novel lowtoxic antifouling coatings [18]. However, analysis of the literature indicates that there is little information about the use of pristine graphene. Indeed, most in situ studies on graphene-based surfaces were performed with functionalized graphene and graphene oxide (GO) coatings [19–23]. Diatom adhesion was completely inhibited after 10 days by surfaces containing 0.36 wt% GO [24]. In turn, GO-silver nanoparticle coatings improved antibacterial and anti-algal properties [25], and showed more than 80% *Halomonas pacifica* biofilm inhibition [26]. Since the current trend is to study the potential of modified and functionalized graphene, the antibiofilm performance of graphene alone is poorly understood. Moreover, most in vitro studies are usually performed for short periods and under hydrodynamic conditions that do not mimic the real marine environment [27]. In fact, some of the in vitro studies have been performed until 24 h [6,23,25,28,29], most of them between days and weeks [7,19,20,24,30–33], but only a study performed by Fazli-Shokouhi et al. [34] extended the assay period for 3 months to evaluate the antifouling potential of graphene-based coatings. Moreover, these studies focus on organisms other than cyanobacteria, namely diatoms, algae, and macrofoulers. Therefore, the main goal of this work was to evaluate the potential of a graphene composite surface to prevent and control the development of biofilms by marine microfoulers over a long-term assay and using an in vitro platform that mimics the hydrodynamic conditions found in real marine scenarios. Cyanobacterial biofilm architecture was evaluated by three different imaging techniques: Optical Coherence Tomography (OCT), Confocal Laser Scanning Microscopy (CLSM), and Scanning Electron Microscopy (SEM). Since the epoxy resin is a commercially available coating generally used to coat the hulls of small recreational vessels [35,36] due to its exceptional physical, chemical, and mechanical properties, no safety issues, and low cost [37], pristine graphene nanoplatelets (GNP) were incorporated into this polymer matrix. Furthermore, epoxy composites exhibited high durability and resistance to fatigue and UV irradiation [38]. Surface characterization was also performed by water contact angle measurements, Optical Profilometry, and SEM.

#### **2. Materials and Methods**

Figure 1 presents the flowchart of the experimental work fully described in the upcoming sections.

**Figure 1.** Scheme of the experimental steps of the present work.


FIGURE 4 | Vertical profiles of alongshore- and time-averaged cross-shore velocities of the hydrodynamic model cases (upper panels) without wind and (middle panels) velocities averaged during the day (24–36 h), and (lower panels) night (36–48 h) in the diurnal wind cycle. Dashed line is Eulerian velocity (uE) and solid line is Lagrangian velocity (uL ). The first three panels to the left are average velocities in the rip channels (blue) and on the shoals (red), and the others are averages over the total alongshore ranges (black). Positive and negative u values correspond to offshore and onshore currents, respectively. 75 m is within the surf zone and 100 m is at the breaker line.

turbulence-driven sinking and DVM were also carried to shore (**Figure 5E**).

When phytoplankton were released in the surf zone, all cases showed cell concentrations higher in the surf zone than offshore (**Figure 6**). Sinking with sand grains enhanced retention of cells in the surf zone. Surf zone concentration of phytoplankton released within the surf zone was highest in the case of the cells with DVM and depth-varying growth (**Figure 6H**); however, cell retention rate in the surf zone in this case was not high compared to the other cases as substantial quantities of phytoplankton seemed to be carried offshore (**Figure 6**). Phytoplankton with turbulence-driven growth and sinking tended to be retained and subsequently increased in the surf zone (**Figure 6F**). Similar to the offshore-released case, cells without any traits also tended to stay in the surf zone (**Figure 6E**).

Cell concentrations in the surf zone would become much higher than offshore over time because onshore cell flux (**Figure 5**) was higher than offshore cell flux (**Figure 6**). All (except for **Figures 5A,B**) cases showed the common distribution trend of phytoplankton (per unit area), that is, highest cell concentrations were found in the rip channels, followed by

that over the shoals, and the smallest concentrations offshore (**Figures 5**, **6**).

### Comparison Between No-Wind and Diurnal Wind Case

Phytoplankton without sinking, growth, and DVM were dispersed vertically in the no-wind condition (**Figure 7A**), whereas those in the diurnal wind condition tended to be concentrated near the bottom and surface (**Figures 7C,E**). Phytoplankton concentration at the surface extended offshore in the no-wind case (**Figure 7A**). A similar extension was formed at night in the diurnal wind case (**Figure 7C**), but not pronounced like in the no-wind case. Vertical velocities at X > 150 m during the day in the diurnal wind cycle directed downward (**Figure 7F**). Offshore vertical velocities (X > 200 m) at night tended to be slightly upward (**Figure 7D**), and neutral with a little variability in the no-wind regime (**Figure 7B**).

# DISCUSSION

This modeling study suggests that turbulence-driven sinking may be an important mechanism for transport of phytoplankton into the surf zone, just like the sinking behavior of competent larvae in turbulence (Fujimura et al., 2014). Once phytoplankton enter the surf zone, their growth may be enhanced by turbulence, resulting in high cell concentration in the surf zone.

Our model showed that the diurnal wind cycle increased onshore transport (**Figures 5A,D** compared to **Figures 5E,F**) and reduced flushing rates of phytoplankton cells (**Figures 6A,D** compared to **Figures 6E,F**). These results were probably due to onshore wind, the major difference between no-wind and diurnal wind regimes. The modeled phytoplankton cells had a relatively small fall velocity (wp = −1.2 × 10−5 m s−1 ), which is 1–2 orders of magnitude smaller than the modeled vertical water velocities (**Figure 7**); thus, the vertical displacement of cells was affected by surrounding currents. In all cases, phytoplankton cells were brought to the surface by the upward flow in and near the surf zone (**Figure 7**). When wind was not applied, many of the cells were kept near the surface throughout the domain, and some cells were sinking or floating in the water column offshore (**Figure 7A**). In the diurnal wind cycle, cells floating near the surface were also carried offshore at night (**Figure 7C**), but upward currents kept them off from the water column (**Figures 7C,D**); furthermore, these cells were pushed back to the shore when the condition switched to onshore wind condition (**Figure 7E**). During the day, downward currents in the water column enhanced downward velocities of phytoplankton (**Figure 7F**), and transported via the onshore bottom currents although these were weakened due to balancing with onshore surface currents (**Figure 4**). Therefore, wind-driven onshore currents were likely to enhance phytoplankton transport toward the shore, while the constant wave-driven flow regime without wind caused higher phytoplankton flushing rates. The results were consistent with a study by Hendrickson and MacMahan (2009) who indicated that diurnal sea breeze is important for cross-shore transport.

Other physical forcing can play a role as well. Stokes drift represented by difference between uL and uE (**Figure 4**) added onshore forcing. Without Stokes drift, onshore transport of materials at the surface may be slower as shown by Fujimura et al. (2014). Bottom onshore currents and benthic streaming owing to wave forcing are also important especially when plankton are in the wave boundary layer. Alongshore variability (rip channels and shoals) can be another factor that increases onshore phytoplankton transport. Bathymetrically-controlled rip currents enhance shoaling to balance cross-shore currents. This bathymetric feature often forms a rip circulation with an eddy at its center where plankton may be trapped (Talbot and Bate, 1987; Fujimura et al., 2014; Shanks et al., 2017). Furthermore, flatter beach slopes (i.e., more dissipative beach) should be more conducive to onshore transport of plankton (Shanks et al., 2010, 2017, 2018; Fujimura et al., 2013; Morgan et al., 2016, 2017).

Our data showed that the concentrations of phytoplankton were highest in the rip channels, followed by that over the shoals, and lowest offshore. This pattern is the same as the larval distribution at the same beach (Fujimura et al., 2014; Morgan et al., 2016, 2017). High concentrations of plankton in the rip channels are possibly a result of the convergence of feeder currents (**Figure 2**). Our results partially agreed with Shanks et al. (2018) that phytoplankton concentrations in the surf zone were much higher than offshore. However, contrary to our study, their data showed that cell concentrations over the shoals were lower than offshore. This discrepancy might be due to the sampling scheme. They collected samples at 1 m depth in the surf zone and throughout the water column offshore, while we averaged phytoplankton cells in the water column everywhere. Hence, we recalculated the time-averaged number of cells over the shoals at 1 m ± 0.5 m. Nonetheless, modeled phytoplankton concentrations were still higher over the shoals than offshore in all cases although shoal to offshore ratios became smaller (e.g., for DVM with depth-varying growth case, 12:1 with original calculation, 2:1 with recalculation). More plausible reason for the discrepancy between Shanks et al. (2018) and our result was a lack of significant environmental parameters and/or phytoplankton traits and further investigation will be necessary for the model to reproduce such a distribution pattern of phytoplankton.

### Limitations and Generality

Phytoplankton cells are suspended occasionally (e.g., in the ocean surface mixed layer) which has been supported by some studies showing that downward velocities of phytoplankton slowed in turbulent waters (e.g., Ruiz et al., 1996; Deleersnijder et al., 2006); nevertheless, this phenomenon can be seen only when turbulence level varies vertically (Ross, 2006; Macías et al., 2013). As the surf zone in our model was a more or less homogeneous layer with the turbulence energy dissipation rate ε > 10−5 m2 s −3 , the value used as the threshold obtained from the empirical data (Ruiz et al., 2004), the reduction of vertical velocity was ignored here. Vertical motion in turbulence is species-specific (Ruiz et al., 2004), while most studies on vertical velocities in turbulence were not conducted with real phytoplankton cells, so more observational and laboratory data will be valuable for more accurate models.

Turbulence may enhance cell growth by increasing nutrient uptake and/or exposure to light. Sullivan et al. (2003) reported that growth rates of dinoflagellate species increased with turbulence levels up to ε ≈ 10−4 m2 s −3 , but decreased at ε ≈ 10−3 m2 s −3 although they still had fairly high growth rates. If we consider their result, the model may be improved by including a growth function that varies with cross-shore varying turbulence energy dissipation rates in the range of 100–10−5 m2 s −3 (shoreline—surf zone edge). Also, the effect of turbulence on growth rate is species-specific as turbulence influences the growth rate of some species negatively (Peters and Marrasé, 2000), while others grow faster in turbulence (Davis et al., 1953; Savidge, 1981; Hondzo and Wüest, 2009). Moreover, the surf zone often receives nutrient-rich freshwater from land through the sand of the beach, and there may be spatial nutrient variability (i.e., the closer to shore, the higher in nutrients). In addition, the surf zone is a relatively shallow region, so generally phytoplankton are exposed to high irradiance, resulting in fast growth rates. Simulations of phytoplankton concentration can be improved by developing cell growth model with corresponding environmental variability.

Our DVM with depth-varying growth model was based on Talbot et al. (1990). This model case showed the most successful onshore phytoplankton migration in our simulations; however, it is not known whether this mechanism applies to any species. As mentioned previously, the ability to float by attaching to air bubbles has been observed in surf zone diatom species (Lewin and Schaefer, 1983), but offshore phytoplankton can also attach to bubbles. In fact, Shanks et al. (2018) collected foam from the sea surface in the rip current, and found that phytoplankton are often highly concentrated in the foam. Those phytoplankton consist of a lot of offshore diatoms and dinoflagellates, and a tiny percentage of surf zone species. It would have been interesting if they had determined the vertical distribution of phytoplankton and sampled at night like Talbot and Bate (1988a) performed.

The model can be improved in many aspects. It is necessary to collect species-specific information including habitat, growth rate, sinking and floating rates, cell size, and morphology. These parameters often interact with each other and the modeling is not trivial even for single species. If we consider that sand grains are sinking agent of phytoplankton, their size variability in the surf zone from observed data (e.g., Gallagher et al., 2011) and their dynamics from numerical model (e.g., Reniers et al., 2013) would help to improve the phytoplankton biophysical Lagrangian model. Likewise, if the target species interact with air bubbles, dynamics and size distributions of air bubbles produced by breaking waves in the surf zone from empirical data (e.g., Deane, 1997; Deane and Stokes, 1999, 2002) or numerical simulations (e.g., Ma et al., 2011) could be incorporated into the biophysical model. Air bubbles may be entrained by "obliquely descending eddies" (Nadaoka et al., 1989), which our model did not include. Another missing surf zone process in our model is a breaking wave roller (Feddersen, 2007; Reniers et al., 2013) that may entrain phytoplankton cells and transport them toward the shore. Moreover, wave conditions also influence onshore transport of plankton. Concentrations of phytoplankton in the surf zone at Sand City beach were higher between 0.5 and 1.0 m wave heights (Shanks et al., 2018). Most zooplankton concentrations in the surf zone tended to be negatively correlated with wave heights at Sand City beach (Morgan et al., 2017) and at a more reflective beach (Shanks et al., 2015; Fujimura et al., 2017; Morgan et al., 2017). Wave spectra may need to be applied for more realistic situations. Random wave groups generate infragravity waves and surf zone eddies that may add extra forcing to onshore transport of plankton and material within the surf zone (MacMahan et al., 2004; Reniers et al., 2010; Fujimura et al., 2014).

# CONCLUSIONS

We found that both biological and physical parameters controlled onshore delivery of phytoplankton and the concentration of the cells in the surf zone. Both wave and wind forcing play an important role in phytoplankton transport. More modeled offshore phytoplankton cells tended to reach the surf zone in the diurnal wind cycle regime than that in the no-wind condition. The cells were highly concentrated in the rip channels by current convergence and rip circulations. Turbulence-induced growth and sinking enhanced the number of phytoplankton cells in the surf zone. DVM also facilitated onshore phytoplankton transport, and depth-varying growth increased the cell concentrations in the surf zone.

The model parameters could be improved with additional empirical data. Further investigation on dynamics of waves, currents, air bubbles, and sand grain will help to develop the model. Species-specific morphological, physiological, and behavioral characteristics of phytoplankton need to be studied as well. Once the model with necessary parameters is validated, it may be coupled with other ecological model such as NPZD and

#### REFERENCES

nested in a larger scale for more practical applications, including HAB predictions and fisheries management.

#### AUTHOR CONTRIBUTIONS

Designed the investigations: AF, AR, AS, JM, SM. Developed the model: AF, AR, CP. Performed the simulations and data analysis: AF. Prepared manuscript: All authors.

#### ACKNOWLEDGMENTS

This research was supported by National Science Foundation (OCE-092735) Collaborative Research: Does coupling between the inner shelf and surf zone regulate larval supply to intertidal populations? AF is funded by National Science Foundation (1457769) Guam Ecosystems Collaboratorium.

**Conflict of Interest Statement:** 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.

Copyright © 2018 Fujimura, Reniers, Paris, Shanks, MacMahan and Morgan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Effects of Anthropogenic Stressors on Reproduction and Recruitment of Corals and Reef Organisms

#### Robert H. Richmond\*, Kaho H. Tisthammer and Narrissa P. Spies

*Kewalo Marine Laboratory, University of Hawaii at Manoa, Honolulu, HI, United States*

The persistence of populations of marine organisms depends on the success of the dual processes of reproduction and recruitment. The production of offspring alone is inconsequential unless larvae and propagules can recruit, which often entails a period of development and distribution in the water column and subsequent selection of appropriate habitats. For fish, this may mean drifting in currents before responding to particular habitat cues. For corals and other benthic invertebrates, larvae must undergo site selection, settlement and metamorphosis into the juvenile form, and survivorship is directly linked to site choice and environmental conditions. Both biotic and abiotic factors affect population replenishment success, and hence, anthropogenic influences such as pollution, sedimentation and climate change can negatively affect critical processes such as reproductive synchronization in spawning species, successful embryological development, appropriate site selection, settlement, metamorphosis and in the case of reef building corals, acquisition of the required zooxanthellae partner. Effective management practices are essential for ensuring the persistence of populations of coral reef organisms of economic, cultural and ecological value.

Keywords: corals, planulae, larval development, larval recruitment, pollution and global change, climate change impacts, metamorphic induction, anthropogenic stressors

#### Edited by: *Eric Wolanski,*

*James Cook University, Australia*

#### Reviewed by:

*Andrew Baird, James Cook University, Australia Hajime Kayanne, The University of Tokyo, Japan*

> \*Correspondence: *Robert H. Richmond<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Marine Science*

Received: *28 March 2018* Accepted: *12 June 2018* Published: *02 July 2018*

#### Citation:

*Richmond RH, Tisthammer KH and Spies NP (2018) The Effects of Anthropogenic Stressors on Reproduction and Recruitment of Corals and Reef Organisms. Front. Mar. Sci. 5:226. doi: 10.3389/fmars.2018.00226* REPRODUCTION AND RECRUITMENT

Reproduction is the process by which new individuals are formed from parental stock, and can occur through asexual or sexual means. In corals, the most common type of asexual reproduction is via fragmentation, and often involves physical separation of a portion of the tissue-covered skeleton. This can occur as a result of fish feeding on associated crabs, shrimp, bivalves or barnacles, and the driver of recruitment "behavior" is gravity, that is, the fragment sinks, settles, attaches and hopefully survives. Other means of asexual reproduction can occur via tissue sloughing or polyp bail out, and the resulting ciliated tissue bits can disperse over greater distances and exhibit selective settlement behavior (Sammarco, 1982). Larvae can also be produced via parthenogenesis, and while genetically identical to the parent colony, exhibit classic larval behaviors. Sexual reproduction involves the fusion of sperm and egg (gametes) produced by the male and female of the species respectively. Sexual reproduction results in two opportunities for increasing genetic variability: one through the contributions of two different parents to the offspring and another through the "crossing over" among chromosomes that occurs during the prophase one stage of meiosis, the reduction division that makes haploid (1n) gametes that

fuse to become a diploid (2n) organism. No matter how the progeny are formed, they do not contribute to the population unless they successfully recruit, and that is where behavior comes in.

The persistence of populations depends on the success of both reproduction and recruitment processes. With a focus on the importance of larval behavior for successful recruitment, the following discussion will center on sexual reproduction as the primary means of producing this critical life history stage for the perpetuation of coral reefs and related ecosystems. There are six chemically mediated steps involved in the replenishment of reef populations: gamete development within mature coral colonies, synchronization of gamete release, successful fertilization of eggs by sperm, complete embryological development to the competent larval stage, settlement and metamorphic induction, and acquisition of symbiotic zooxanthellae in recruits that do not have them vertically transmitted by the parent colony (Richmond, 2014; **Figure 1**). These key processes are affected by both density dependent and density independent factors.

As with all sessile and benthic organisms, for reproductive events to be successful in the production of viable larvae, a critical density of parental stock is required. As the distance between individuals and colonies of the same species increases, the incidence of appropriate egg-sperm interactions decreases, with reduced levels of fertilization and larval production. This is called the Allee Effect, and is well known within marine organisms (Gascoigne and Lipcius, 2004). Population density is a function of environmental conditions, stressors and disturbance.

Water and bottom quality are essential parameters that affect the health of coral reefs and other marine ecosystems. Reduced water quality can stress adult populations, inhibiting gametogenesis and reducing chemical signaling between conspecifics that helps with the synchronization of spawning, the viability of gametes, fertilization success and the ability of embryos to develop into viable planula larvae. As gametes (egg and sperm) are among the most sensitive life history stages to environmental parameters, terrigenous runoff is a

FIGURE 1 | Recent coral recruit exhibiting colony growth through extra-tentacular budding and the acquisition of symbiotic zooxanthellae from the ambient environment.

serious problem that impacts reproduction of marine organisms (Richmond, 1993). Bottom quality affects recruitment success in corals and many other benthic invertebrates. Sediment cover prevents larvae from sensing key settlement cues, both chemical and textural, from preferred substrata such as crustose coralline algae. Larval metamorphic induction may be prevented by sediment and algal cover through substrate instability, the alteration of the benthic microbiome, and the alteration of bottom chemistry. In addition to the local stressors mentioned above, increased temperatures and ocean acidification tied to global climate change are factors that also affect reproduction and recruitment in a broad range of marine organisms. These anthropogenic issues will be discussed below.

#### LARVAL COMPETENCY

The competency period of larvae, the time during which they can successfully disperse, settle and metamorphose, is an important determinant of connectivity among sites and populations. The longer the competency period, the greater the distributional range of a species, and broader distributions translate into enhanced survival of species exposed to local and global level stressors with the associated threat of local extinction. The competency period of larvae is tied to nutritional mode and energy availability. Planktotrophic larvae actively feed, and often have specific structures that aid in food capture. They have a delay before they can recruit, the pre-competency period, during which time they undergo additional development and often, morphological changes. Lecithotrophic larvae survive primarily on stored energy reserves, such as lipids provided in the eggs. The term mixotrophic refers to larvae that use a combination of stored energy reserves and planktonic feeding, which can include absorption of dissolved nutrients, and in the case of corals, translocated metabolites from symbiotic zooxanthallae.

While in the plankton, larvae are not limited to drifting at the mercy of currents, but can exert a degree of control through vertical and horizontal swimming behaviors. Such movement can change the water parcels in which the larvae are contained and hence, the direction and speed of travel, ambient temperature that affects larval development and metabolic rates, and light levels including photosynthetically active radiation (PAR) and UV exposure. An extreme example is for larvae that get entrained in water parcels affected by the Kelvin Wave associated with the onset of El Nino conditions. Such larvae can be transported over hundreds to thousands of km over a period of weeks to months (Richmond, 1990). Downward vertical movement of larvae can also affect entrainment of larvae in bottom currents and eddies that enhances recruitment on coral reefs and other benthic habitats. Recent research has demonstrated that larval swimming behavior can be influenced by a number of parameters including light, gradients of salinity and temperature, as well as by auditory and chemical cues. Both fish and invertebrate larvae have been documented to respond to reef sounds that aid in the recruitment of larvae to reefs, and chemical cues can influence larval swimming and settlement behavior (Vermeij et al., 2010). While both auditory and chemical cues can attract larvae, the distance over which some, including coral planulae, can swim is limited (Figueiredo et al., 2014; Hata et al., 2017).

Fish larvae usually have an initial distribution pattern tied to currents, but after a degree of development and spinal hardening, can begin to exhibit directional swimming. Coral reef fish, as opposed to pelagic species, need to find specific habitats in order to recruit, grow and eventually reproduce, hence, cueing behavior is of paramount importance. Specific cues include sounds and smells. Waves impacting reefs and other barriers send sounds over relatively long distances. Interestingly, these reflections and refractions of wave energy are also used by traditional navigators in the Pacific Islands to identify the location of landmasses and submerged features. Over shorter distances, the sounds of crabs, snapping shrimp and fish can help larvae navigate toward healthy reefs. Dead reefs and reefs in decline may provide fewer cues to which larvae can respond (Dixson et al., 2014). The natural chemistry of reefs, tied to everything from coral mucus, algae and bacteria, also provides signals to larvae, which possess chemoreceptors. Toxicants that block these chemical signals or the larval receptors, act to impede recruitment (Peters et al., 1997).

Coral planula larvae do not randomly settle and metamorphose, but rather, select their recruitment sites (Golbuu and Richmond, 2007). Bacteria have been specifically identified as important cues for coral larval recruitment to specific substrata. Coral planulae have been found to also respond to specific crustose coralline algae (CCA), likely due to associated polysacchides and glycoglycerolipids (Tebben et al., 2015; **Figure 2**).

Coral larvae from brooding species and some spawning species (e.g., Montipora) possess symbiotic zooxanthellae that are vertically transmitted from the parent colony. These algal cells provide additional energetic input to the larvae and can extend the competency period (Richmond, 1987). While the period during which these larvae can successfully settle and metamorphose can extend from weeks to months, translating into the potential for long-distance dispersal, studies have shown that most larvae travel over short distances and many reefs are self-seeding (Sammarco et al., 1989; Jones et al., 1999; Figueiredo et al., 2013). This is due to a combination of abiotic factors (currents and associated eddies) and larval behavior, including cueing to the biotic factors of natal reefs, identified above.

Water quality remains a key concern during the planktonic phase of larvae throughout the competency period. Substratum quality is also important for settlement and metamorphosis to occur. The implications are that efforts to protect bottom and water quality in and around coral reef ecosystems are of critical importance to maintaining reef populations by supporting successful reproduction and larval recruitment.

### THE EFFECTS OF STRESSORS ON CORAL REPRODUCTION AND RECRUITMENT

#### Water Quality

The world's nearshore marine habitats are being increasingly exposed to a variety of anthropogenic stressors, such as sedimentation, eutrophication, pollution, and overfishing (Morgan et al., 2016), all of which have a profound impact on reproduction, recruitment and survival of larvae. As mentioned above, early life stages of many marine organisms are more sensitive to pollutants than adult stages (Reichelt-Brushett and Harrison, 1999), due in part to the important role that chemical cues play in substrate selection and metamorphosis of larvae in many marine organisms (Hadfield and Paul, 2001; Freckelton et al., 2017). For example, many species of coral larvae use chemical cues to settle on a substratum where specific coralline algae and/or their associated bacteria are present (Golbuu and Richmond, 2007; Vermeij et al., 2008; Tebben et al., 2015). Recently, a single bacterial metabolite, tetrabromopyrrole, was shown to induce metamorphosis in certain coral species, occurring even without attachment to a substratum (Tebben et al., 2011; Sneed et al., 2014). Chemical pollutants, therefore, affect key life stages and behaviors of corals and other marine organisms by interfering with their chemosensory functions.

FIGURE 2 | Coral planulae undergoing settlement and metamorphosis (panel 1) following exploration of the substrata (panels 2 and 3) and encountering the proper chemical/biotic cues.

FIGURE 3 | Sewage discharges along a reef on Guam, resulting in decreased recruitment of coral larvae. The freshwater effluent rises 100% of the time and bathes the adjacent shallow water reef area when there is an onshore wind. The chemical content includes personal care products, heavy metals, pesticides, pharmaceuticals, hydrocarbons and even estrogenic compounds from un-metabolized birth control pills.

A variety of chemical agents are introduced into our coastal waters via coastal development, terrestrial runoff, sewage outfalls and toxic antifouling paints from ships (**Figure 3**). Many studies have shown the negative effects of these agents on coral larvae and their behaviors (Woods et al., 2016). For example, increased heavy metal concentrations reduce fertilization rates, larval survival, settlement success, metamorphosis, and growth rates of reef-building corals. Copper (Cu) is shown to affect all of the above stages in coral larvae at relatively low concentrations (Reichelt-Brushett and Harrison, 1999, 2000; Victor and Richmond, 2005; Negri and Hoogenboom, 2011). In fact, Cu inhibits the metamorphosis of coral larvae at the lowest concentration of any metals tested so far (Reichelt-Brushett and Harrison, 1999). Exposure to iron (Fe) induced oxidative stress responses and DNA damage in coral larvae (Vijayavel et al., 2012), affecting their metabolic functions and behaviors. Exposure of coral larvae to mercury, cadmium and zinc showed neutral or mixed results (Heyward, 1988; Reichelt-Brushett and Harrison, 1999, 2005; Farina et al., 2008), suggesting interspecific differences in sensitivity to these metals.

The effects of organic compounds such as anti-fouling agents, hydrocarbons, and pesticides have been increasingly studied in marine organisms, including their impacts on reproductive and larval recruitment success. Tributylin (TBT) and Irgarol 1051 are active compounds of anti-fouling paints, and are reported to inhibit larval settlement and/or metamorphosis (Negri et al., 2002; Inoue, 2004; Negri and Marshall, 2009; Knutson et al., 2011). Since the banning of applying paints with TBT in 2003 by the International Maritime Organization, the use of the algaecide Irgarol in marine paints has increased dramatically. However, even relatively small concentrations of Irgarol can inhibit coral larval settlement (Knutson et al., 2011), which suggests that the Irgarol use may need to be reevaluated. Hydrocarbons, such as crude oil, dispersed oil, oil dispersant, and benzo[a]pyrene, are all reported to negatively affect fertilization, metamorphosis, and settlement of coral larvae (Epstein et al., 2000; Lane and Harrison, 2000; Negri and Heyward, 2000; Goodbody-Gringley et al., 2013). Oil dispersants and dispersed oil appear to have more negative effects than crude oil itself (i.e., water soluble fractions); morphological deformations and loss of normal swimming behavior of coral larvae were observed with exposure to dispersants and dispersed oil (Epstein et al., 2000).

Corals are holobionts, harboring endosymbiotic algae, zooxanthellae (Symbiodinium spp.), in their tissues. These algal symbionts provide the majority of the coral's energy needs through photosynthesis. Increased detection of herbicides and insecticides in marine environments has raised serious concerns regarding their impacts on coastal coral reef ecosystems (Markey et al., 2007; van Dam et al., 2011), since some herbicides used in agriculture and antifouling paints, such as diuron and atrazine, are photosystem II (PSII) inhibitors that directly affect zooxanthellae's functions (Owen et al., 2002; Jones et al., 2003; Jones, 2004, 2005). Our understanding about how these herbicides and insecticides affect coral hosts is still limited. The effects of pesticides found in coastal waters, on coral reproduction and physiology, are of high concern to resource managers. Our results from laboratory exposure experiments showed that the tested herbicides and insecticides had negative effects on at least one or more of coral life stages. For example, Roundup R (Monsanto Company) significantly reduced fertilization and settlement rates in corals at very low concentrations (0.0069 ppm), and caused damage to eggs at a relatively low concentration (>0.35 ppm) (Diu, 2016). Another emerging pollutant gaining attention in the marine environment is the active ingredient in some sunscreens, oxybenzone. Increased mortality, deformation, and DNA damage, as well as reductions in photosynthetic efficiency were observed in coral larvae exposed to oxybenzone (Downs et al., 2016). These observations raise concerns about the environmental safety of other personal care products that enter into coastal waters and the open ocean.

#### Bottom Quality

Water and substratum quality are critical components affecting both reproduction and recruitment of marine organisms. In order for successful population replenishment to occur, consideration must be given to the source areas, the conduits of connectivity and the recipient sites. For example, larvae formed at one site will be physiologically compromised or killed if transported across an area affected by polluted and sediment-laden watershed discharges (Victor et al., 2006). Even a drop in salinity will affect survivorship of many types of larvae, including corals and echinoderms (Richmond and Jokiel, 1984). Contaminated recipient sites and those with accumulated sediment will also be inhospitable to benthic invertebrate and fish larvae from clean source sites, as such degradation will affect critical chemically cueing.

Sedimentation affects the success of reproduction and recruitment several different ways. During spawning events, sediment, particularly the finer clay-size particles, can adhere to eggs, causing them to sink rather than float to the ocean's surface, where fertilization often occurs. Sediment is a major carrier of toxicants from land to the sea, where the change in chemistry can cause contaminant release (**Figure 4**). Such pollutants include petroleum products (polycyclic aromatic hydrocarbons–PAH's), pesticides (both insecticides and herbicides), personal care products and heavy metals. Sediment particles often serve as nuclei that aggregate organic material, including mucus and phytoplankton, which harbor bacteria, and when this material sinks to the bottom, it increases biological oxygen demand (BOD), dropping both oxygen levels and pH at the sedimentwater interface (Wolanski et al., 2003). This effect can negatively impact the survivorship and growth of larval recruits. Finally, sediment blocks key metamorphic inducers needed by many types of larvae, and most benthic epifaunal larvae cannot recruit to such unstable substrata.

Sediments impact reproductive, larval behavior and recruitment success. In additional to the chemical effects listed above, sediment has physical effects that in turn, affect the physiology of organisms. For organisms dependent upon photosynthesis for energy, metabolism, growth and reproductive allocation, the shading effects of sediment block photosynthetically active radiation (PAR) that is essential to running the light-mediated chemical reactions. Turbidity, expressed as a weight per unit volume of water (e.g., mg/l), is especially problematic for photosynthesis. Sedimentation deposition rate, a weight per unit area over a given time period (e.g., g/m2 /day) impacts both photosynthesis and energetics, due to the cost of shedding particles. Sediment can be re-suspended numerous times, so even without the addition of new materials from land-based sources, waves can result in additional impacts through reintroducing particles into the water column.

### Temperature

Coral larvae remain viable and competent within the ambient temperature range of the parent colonies (Woolsey et al., 2015). As with adult colonies, deviations of a few degrees above or below the norm can have deleterious effects, such as reducing competency or inducing mortality. The expected increase in tropical and subtropical seawater temperatures of 2–4◦C due to climate change, has major implications for coral larval survival, development, dispersal, and settlement (IPCC, 2007). Temperature is tied to larval development, and alterations

FIGURE 4 | Sediment covering the reef surface blocks metamorphic inducers and prevents larval recruitment.

in the ambient rearing temperature can have drastic effects. Abnormalities in embryogenisis of the spawning coral Acropora millepora were seen at 32◦C (4◦C above ambient temperature where gametes where collected), as well as decreased fertilization rates, with no fertilization occurring at 34◦C (Negri et al., 2007). A temperature increase reduced pre-competency periods by as much as an entire day in some coral species, which can affect the dispersal range for coral species, as well as impede the ability for devastated coral reefs to recruit new individuals and recover following bleaching events (Heyward and Negri, 2010). Model predictions suggest that increasing temperatures will serve to increase local retention of coral planulae due to reduced competency periods, and reduce connectivity among reefs separated by greater distances (Figueiredo et al., 2014). The implications are positive for locally targeted conservation efforts but raise additional concerns for reefs and populations dependent on external larval supplies for recovery and enhancing genotypic diversity.

Elevated temperature effects on coral planulae can have regional implications for population replenishment as well as global ones, as higher latitude larvae tend to develop more slowly than their tropical counterparts. An increase of even a single degree (◦C) above the normal range can increase the metabolic expenditure of stored energy, increasing premature metamorphosis and reducing recruitment success in addition to causing outright mortality (Edmunds et al., 2001; Graham et al., 2017). In aposymbiotic larvae from Diploria strigosa, a 2◦C rise in temperature from an ambient 28◦C seawater temp caused 50% mortality, and an increase to 31◦C caused 70% mortality (Bassim and Sammarco, 2003).

As mentioned, coral larvae respond to different chemical and biological cues, and must have an appropriate settlement substrate for successful settlement and metamorphosis to occur. Temperature experiments elevating seawater to 32◦C for a period of a week have been able to bleach CCA, also resulting in a shift in its microbial community (Webster et al., 2010). Acropora

millepora larval settlement and metamorphosis was reduced by 50% when these planulae were exposed to bleached CCA compared to healthy CCA from ambient seawater conditions. While a microbial shift can occur by altering the chemical compounds released by CCA, it's also important to note that the color change from pink to white when CCA bleaches can also alter larvae settlement due to a resulting lack of an appropriate spectral cue. In an experiment examining synthetic surfaces of either white or red hues, the coral larvae preferred a red surface 85% more often when compared to a white synthetic surface when both were also exposed to a chemical cue (Foster and Gilmour, 2016). The larvae themselves exhibit spectral variance by producing fluorescent proteins that change in response to environmental conditions. Fluorescent proteins are present in both larval and adult corals. While their function is not fully understood, they appear to be altered when the animal is exposed to stressors (Palmer et al., 2009). There is a demonstrated change in the expression of green fluorescent protein (GFP) in response to increased light, as well as red fluorescent protein (RFP) in response to heat(Rodriguez-Lanetty et al., 2009; Kenkel et al., 2011). The combined effect of increased water temperature and light intensity can also affect larval settlement, with a 16% decrease in larval settlement when they are exposed to both high temperatures and higher than normal light conditions (Kenkel et al., 2011).

When considering larval survival and physiology under increased temperatures, the symbiotic relationship between the animal and algal symbionts must be considered. Aposymbiotic larvae and symbiotic larvae both respond to stress in different ways, yet survivorship between temperatures of 28–32◦C is not significantly different (Baird et al., 2006). When the aposymbiotic larvae from A. millepora were exposed to temperatures above their normal range for a period of 10 h, heat shock protein gene expression was elevated but oxidative stress pathways remained unchanged (Rodriguez-Lanetty et al., 2009 ). The infection of symbionts within aposymbiotic larvae is a complex process that is just beginning to be understood under normal seawater temperatures, and may be altered under increased seawater temperatures (Weis et al., 2001). A study examining physiological responses to 4 clades of symbionts, as well as coral larvae infection, found that there were different optimum temperatures among the symbiont types where they had the ability to infect aposymbiotic coral larvae (Baird et al., 2008). Certain symbionts demonstrated lower rates of survival even at 32◦C, giving other symbiodinium strains a competitive edge. This advantage may come with tradeoffs however, as one of the most thermally tolerant strains had the slowest rate of population increase after infecting the coral larval host. This could have long-term impacts on the ability for the coral recruit to establish itself and grow, when fewer symbionts are available to provide nutrients needed for growth.

The means by which corals transmit their symbionts can affect stress in and survivorship of larvae. Hartmann et al. (2017) suggested that vertical transmission of zooxanthellae to coral eggs may exact a cost on gamete and larval survival due to the effect of increased buoyancy from lipids that result in the exposure to increased temperature and UV light at the sea surface. As many corals spawn during the rainy season in parts of the world, the sea surface can be a stressful environment with reduced salinity and increased toxicant concentrations from freshwater runoff. Hartmann et al. propose that the acquisition of symbionts after metamorphosis improves fertilization and out-crossing following spawning. The adaptive advantage of horizontal transmission of zooxanthellae is reflected in the fact that the majority of spawning corals (71%) acquire their symbionts that way.

Under conditions of elevated temperature, algal symbionts can contribute to stress and mortality in coral larvae. Yakovleva et al. (2009) found coral planulae with zooxanthellae had higher mortality rates than those lacking algal symbionts at high temperatures. Mortality was attributed to oxidative damage tied to photosynthesis and was a function of both the magnitude of increased temperature and the duration of exposure. As with exposure to other physical and chemical stressors, zooxanthellae shift from being a mutualistic partner to becoming a burden at higher temperatures.

# Ocean Acidifcation

Reduced oceanic pH attributed to increasing levels of atmospheric CO2 is a concern for the future of coral reefs as the structural foundation of these ecosystems results from the deposition of CaCO3 (aragonite). Chua et al. (2013), performed a study of the effects of pH on coral planulae using four relevant exposure treatments, and found no consistent effects on embryological development, larval survivorship and metamorphosis. They posit that larvae containing zooxanthellae may be more sensitive to the effects of decreased pH through metabolic effects on the symbionts. As coral planulae (generally) do not calcify while in the plankton, ocean acidification is more likely to directly affect larvae that do, such as those of sea urchins (Kurihara, 2008). However, the long-term effects of OA on coral recruitment and survivorship are expected to increase and, based on current predictions, be substantial. Reduced pH does affect the recruitment and growth of crustose coralline algae, which has been shown to be a preferred settlement substratum for many coral species (Kuffner et al., 2008). As small coral recruits can be outcompeted for space by a number of benthic organisms including algae, sponges and ascidians, reduced colony growth rates will affect survivorship as OA increases.

# SURVIVORSHIP–SELECTION/ADAPTATION VS. ACCLIMATIZATION

Intraspecific variation, or inter-individual differences, in traits arise from genetic predispositions, past experiences and physiological states (Guest et al., 2012; Nanninga and Berumen, 2014), and is increasingly recognized as an important ecological factor that affects population dynamics and evolution of species, including reproductive success (Benton et al., 2006; Bolnick et al., 2011). In reef-building corals, significant differences in larval success has been observed among different genotypes (Meyer et al., 2009; Baums et al., 2013; Miller, 2014; Miller et al., 2016), as well as due to the parental experiences of exposure to stressful conditions (Putnam and Gates, 2015). This tells us that the larval stages are subject to strong selective pressures from their environmental conditions. Since persistence of populations, and ultimately species, depends on producing better adapted genotypes, the larval stages play a key role in coral population dynamics, connectivity, and survival. The fitness of existing individuals depends on their acclimatization potential, and many larval behaviors likely contribute to increasing their fitness; examples include selecting an optimal environment to settle (e.g., Davies et al., 2014; Dixson et al., 2014; Doropoulos et al., 2016), using reverse metamorphosis to escape from adverse conditions (Richmond, 1985; Negri et al., 2005), and potentially benefitting from increased size by forming a chimeric aggregate (Amar et al., 2008). When selective pressures exceed the acclimatization capacity of an individual, the individual will be eliminated from the population (i.e., selection). However, since no individual larva will experience exactly the same conditions, one evolutionary response is to produce genetically variable offspring so that the population will have enough adaptive capacity to respond to environmental heterogeneity and unpredictability. Intraspecific variation in reproductive success may, therefore, be a result of such an evolutionary strategy. Looking forward, corals and many other marine organisms will face strong selective forces from rising sea surface temperatures and ocean acidification (Webster et al., 2013; Baums et al., 2014). The larval stages known to be vulnerable to environmental stressors (Polato et al., 2010) will especially be subject to this strong selection, resulting in adaptation playing an essential role in how coral populations respond to changing environmental conditions. Standing genetic variation fuels such short-term adaptation (Stapley et al., 2010; Voolstra et al., 2011; Shinzato et al., 2012). In this regard, loss of genetic diversity already seen in some threatened coral populations (Baums et al., 2013) poses a great concern, as it reduces the coral's future reproductive and adaptive capacity. Human activities, including those responsible for global climate change, are predicted to reduce overall genetic

#### REFERENCES

diversity in populations and species (Templeton et al., 2001; Pauls et al., 2013), which further threatens the persistence of coral reef ecosystems.

## SUMMARY AND CONCLUSIONS

In a world of mounting anthropogenic stressors affecting the world's oceans and all life contained within, the ability of populations of organisms to survive, reproduce, evolve and grow is threatened at the local, regional and global levels. An understanding of the reproductive and recruitment behavior of marine organisms is essential to addressing the challenges of resource sustainability through responsive and effective management of the human activities responsible for resource declines. New tools such as genomics, proteomics and transcriptomics allow for the study of behaviors at the cellular and molecular levels, which when combined with more traditional studies of physiology and organismal behavior, have the potential to provide key guidance. The needs are pressing if there is to be a sound legacy of ocean resources for future generations.

# AUTHOR CONTRIBUTIONS

RR organized the manuscript, coordinated sections with the coauthors, wrote several sections, provided figures and edited versions. KT wrote several sections and edited versions. NS wrote several sections, contributed figures and edited versions.

# FUNDING

Funds partially supporting this research were provided by NOAA's Protected Resources Division, NOAA Fisheries' Pacific Islands Regional Office (NOAA grant # NA17NMF4520144) and the National Fish and Wildlife Foundation (NFWF grant # 34413).

stress tolerance in the reef building coral, Acropora palmata. Coral Reefs 32, 703–717. doi: 10.1007/s00338-013-1012-6

Palmer, C. V., Modi, C. K., and Mydlarz, L. D. (2009). Coral fluorescent proteins as antioxidants. PLoS ONE 4:e7298. doi: 10.1371/journal.pone.0007298

millepora by tetrabromopyrrole isolated from a Pseudoalteromonas Bacterium. PLoS ONE 6:e19082. doi: 10.1371/journal.pone.0019082

**Conflict of Interest Statement:** 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.

The reviewer AB and handling Editor declared their shared affiliation.

Copyright © 2018 Richmond, Tisthammer and Spies. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Biophysical Simulations Support Schooling Behavior of Fish Larvae Throughout Ontogeny

Igal Berenshtein1,2,3 \*, Claire B. Paris 3 , Hezi Gildor 4 , Erick Fredj 5 , Yael Amitai 4 and Moshe Kiflawi 6,2

*1 Department of Life Sciences, Eilat Campus, Ben-Gurion University of the Negev, Beer-Sheva, Israel, 2 The Interuniversity Institute for Marine Sciences of Eilat, Eilat, Israel, 3 Department of Ocean Sciences, Rosenstiel School of Marine and Atmospheric Science, University of Miami, Miami, FL, United States, 4 Institute of Earth Sciences, The Hebrew University, Jerusalem, Israel, 5 Department of Computer Science, Jerusalem Institute of Technology, Jerusalem, Israel, 6 Department of Life-Sciences, Ben-Gurion University of the Negev, Beer-Sheva, Israel*

Schooling is very common in adult and juvenile fish, but has been rarely studied during the larval stage. Recent otolith micro-chemistry studies of coral reef fish have demonstrated that cohorts of larvae can move through similar paths and settle within a few meters one from another. However, little is known about the processes involved in the formation and maintenance of these cohorts. Here we use a biophysical modeling approach to examine whether local hydrodynamics, various individual behaviors, or larval schooling can explain cohesive patterns observed for *Neopomacentrus miryae* in the Gulf of Aqaba/Eilat (Red Sea), and whether schooling is feasible in terms of initial encounter probability and cohesiveness maintenance. We then examine the consequences of schooling behavior on larval settlement success and connectivity. Our results indicate that: (1) Schooling behavior is necessary for generating cohesive dispersal patterns, (2) Initial larval encounter of newly-hatched larvae is plausible, depending mainly on initial larval densities and patchiness, and (3) schooling behavior increases the rate of larval settlement while decreasing the percentage of realized connections. Together with mounting evidence of cohesive dispersal, this numerical study demonstrates that larval schooling throughout the pelagic phase is realistic, and has a significant effect on settlement success and connectivity patterns. Future research is needed to understand the mechanisms of fission-fusion dynamics of larval cohorts and their effect on dispersal. Our findings should be considered in future efforts of larval dispersal models, specifically in the context of marine connectivity and the planning of marine protected area networks.

Keywords: larval dispersal, schooling, cohorts, coral reef fish larvae, cohesive dispersal, settlement, connectivity, Gulf of Aqaba/Eilat

#### INTRODUCTION

Adult and juvenile fish often swim in groups. Such schooling behavior can help reduce drag and conserve energy (Weihs, 1975), improve orientation (Codling et al., 2007), reduce the risk of predation (Shaw, 1978), and enhance the detection of food and mates (Shaw, 1978). On the other hand, schooling fish may suffer higher mortality and lower growth rates due to competition within the group (Hixon and Jones, 2005). The prevalence of schooling suggests that the overall balance would seem to be positive (Shaw, 1978).

#### Edited by:

*Alberto Basset, University of Salento, Italy*

#### Reviewed by:

*Ulisses Miranda Azeiteiro, University of Aveiro, Portugal Kevin Sorochan, Dalhousie University, Canada*

\*Correspondence: *Igal Berenshtein<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Marine Science*

Received: *23 December 2017* Accepted: *03 July 2018* Published: *25 July 2018*

#### Citation:

*Berenshtein I, Paris CB, Gildor H, Fredj E, Amitai Y and Kiflawi M (2018) Biophysical Simulations Support Schooling Behavior of Fish Larvae Throughout Ontogeny. Front. Mar. Sci. 5:254. doi: 10.3389/fmars.2018.00254*

The advantages of schooling may similarly benefit the pelagic larval phase that characterizes the life-cycle of most coral reef fishes (Leis, 2006). These larvae are thought to actively group, mainly during late larval stages, once their sensory and movement capacities have matured enough to support this behavior (Shaw, 1978; Arvedlund and Kavanagh, 2009). Specifically, larval schooling may be similar to adult fish schooling, where collective behavior emerges from simple "individual-fish" rules based on the location of the nearest neighbor/s (Couzin et al., 2002).

Little is known about the frequency of schooling behavior across species and its ontogenetic dynamics. Previous studies indicated that Aldrichetta forsteri and Gobiosoma boscifish larvae school at early stages when their size is ∼5 mm and their fins are not fully developed yet (Breitburg, 1991; Kingsford and Tricklebank, 1991). However schooling was also reported at later larval and post-larval stages for different species (reviewed in Leis, 2006).

Mounting evidence, derived from the analysis of otolith micro-chemistry (Ben-Tzvi et al., 2012; Shima and Swearer, 2016) and genetic relatedness (Bernardi et al., 2012; Selwyn et al., 2016), suggests that cohorts of coastal marine fish may be sharing a substantial portion of their dispersal histories (by "cohort" we mean a group of larvae of same age). While it is tempting to relate these patterns to larval behavior, they may also result from aggregation by prevailing hydrological features and/or local cohesive Lagrangian structures (Siegel et al., 2008). The relative contribution of biological, physical, and bio-physical processes to the generation and maintenance of such cohesive dispersal pattern, i.e., shared dispersal paths among individuals within the same cohort, is largely unexplored.

This gap of knowledge is critical, given the demographic and ecological ramifications of the dispersive fate of marine larvae (Levin, 2006). Specifically, the dispersal outcome of schooling larvae is likely to differ from individual larvae as larval schools are characterized by faster swimming speeds, and higher precision of orientation (Irisson et al., 2015). Distinguishing between individual and schooling dispersal is important since dispersal and connectivity are central elements affecting population dynamics, and are often considered when designing marine protected areas networks (Cowen and Sponaugle, 2009). Additional importance of larval schooling lies in its potential effect on the ecology of pelagic ecosystems via planktonic predator-prey dynamics (McGurk, 1986; Frank et al., 1993; Bradbury et al., 2003), planktonic patchiness enhancement (Bradbury et al., 2003), species diversity (Hill, 1973).

In the current study we used a biophysical model of the Gulf of Aqaba (GoA) to test which aspect of larval behavior is necessary and sufficient to reproduce patterns of cohesive dispersal reported for Neopomacentrus miryae in the GoA (Ben-Tzvi et al., 2012). Briefly, Ben-Tzvi et al. (2012) have demonstrated that newly settled N. miryae cohorts shared a highly similar otolith micro-chemistry pattern within each cohort compared to between cohorts. In the current study, we tested whether different scenarios generate cohesive cohorts by examining passive drift, as well as by invoking larval behaviors such as simple random walk (SRW), auto-correlated directional swimming (CRW; Berenshtein et al., 2017), ontogenetic vertical migration (OVM; Paris and Cowen, 2004; Irisson et al., 2010; Huebert et al., 2011), biased correlated directional swimming (BCRW; Staaterman et al., 2012; Mouritsen et al., 2013), and schooling behavior (Irisson et al., 2015). We then examined the plausibility of larval encounter during the early pelagic phase; their capacity to overcome turbulent velocities and remain grouped; and finally, we quantified the relative changes in settlementsuccess and connectivity in schooling vs. non-schooling dispersal.

# METHODS

# Study Site

The Gulf of Aqaba is a semi-enclosed elongated extension of the Red Sea (**Figure 1A**), stretching across ∼180 km (**Figure 1B**), with maximal width and depth of ∼27 km and 1,850 m respectively. The GoA circulation is characterized by deep mixing during winter (Feb-March), and stratification during the rest of the year, with a strong density-driven exchange flow with the Red sea during April–August (Biton and Gildor, 2011). The GoA circulation is also dominated by a chain of several sub-meso-scale eddies along its main axis (Biton and Gildor, 2011).

## Study Fish

Neopomacentrus miryae is endemic to the Red Sea and is highly abundant in the GoA. The ecology of this species was described by Fishelson et al. (1974) referring to it as Abudefduf azysron (Dor and Allen, 1974). Adult N. miryae are planktivorous, found in large stationary schools along vertical walls of forereefs, with males guarding and aerating demersal eggs. N. miryae spawning season spans between November-May, and larvae are found at depths of 0-100 m in the water column, with highest abundance at the top 25 m (Kimmerling et al., 2017).

#### Biophysical Model of the GoA

The biophysical model of the GoA (Berenshtein et al., 2017) was composed of 3 parts: First, the winds at 10 m above sea level were computed using an atmospheric model- Weather Research and Forecasting model (WRF; Skamarock et al., 2001), Next, wind data was used to force a 3D high resolution oceanographic model- Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997), Finally, an offline particle tracking algorithm (modified from Fredj et al., 2016) made use of the 3D circulation data combined with larval biological traits to compute larval trajectories. The resolution of the MITgcm data was of 1 h (temporal), 300 m (horizontal), and 32 depth layers most of which were concentrated in the upper 300 m (Biton and Gildor, 2011; Berenshtein et al., 2017). The temporal resolution of the particle tracking model was of 10 min (time-step length), with a spatial resolution of centimeters due to a continuous 3D grid, and spatial interpolation of larval positions. The detailed description of the model can be found in Berenshtein et al. (2017).

# Backward and Forward-in-Time Passive Simulations

For our simulations we used 5 real recruitment events which occurred at 3 sites in the northwestern coast of the GoA: North Beach (NB), Inter-University Institute (IUI), and Princess Beach (PB) (**Figure 1C**; Ben-Tzvi et al., 2012). In Ben-Tzvi et al. (2012), newly-recruited fish were collected with hand-nets, and their age and day of recruitment was determined using otolith analysis. Newly-recruited fish that shared the same site, age and recruitment day were considered of the same cohort (Ben-Tzvi et al., 2012).

To examine whether the local hydrodynamics (advection and turbulence) could generate patterns of cohort cohesiveness similar to those observed in N. miryae (Ben-Tzvi et al., 2012), we performed backward and forward-in-time simulations using the dates and locations of the collected cohorts and their estimated PLDs (**Figure 2C**; Ben-Tzvi et al., 2012).

The rational of this approach was that if currents were to generate a cohesive dispersal pattern, they would have to be sufficiently strong and persistent to overcome the system's stochasticity. Such pattern would result in similar trajectories in backward and forward-in-time simulations, connecting between the settlement sites, and the hatching locations on the reef.

In the backtracking simulation (Sim. #1, **Table 1**) 5,000 passive virtual larvae were released from each of the recruitment sites at their corresponding dates. In the forward simulation (Sim. #2, **Table 1**), 20 virtual larvae were released from each of the locations of virtual larvae which ended up on the reef in the backtracking simulation, releasing a total of 4,480 virtual larvae (Sim. #2, **Table 1**).

We then visually examined the simulated trajectories for indication of cohort-unique dispersal trajectories, imposed solely by local hydrodynamics.

# Pairwise Distances and Probability of Adjacent Settlement

To determine if various larval behaviors can generate the observed cohesive dispersal pattern we released pairs of virtual larvae from the entire fringing reef of the GoA (41 reef polygons in **Figure 1B**) and tracked them forward in time (Sim. #3–8, **Table 1**). For each simulation, 1,000 pairs of larvae, characterized by identical hatching times and locations, were released daily at midnight across 9 release days (15th−23th of March 2004). Larval trajectories were computed in the biophysical model according to advection, turbulence, and larval behavior (Berenshtein et al., 2017). Mortality was not applied in the biophysical simulations.

To quantitatively assess the occurrence of cohesive dispersal we computed: (1) larval pairwise distances, i.e., the distance separating paired virtual larvae across the PLD, and (2) the probability of sampling 2 settled nearest neighbors, which hatched at the same location and time. While the former corresponds to the observed similarity in larval paths deducted from the trace-element of otolith rings, the latter corresponds to the similarity in the site of origin deducted from the traceelement of otolith cores (Ben-Tzvi et al., 2012).

Since very little is known about the behavior of N. miryae larvae, we chose to apply the commonly observed and modeled larval behaviors for other species within the Pomacentridae family (Leis and Carson-Ewart, 2003; Irisson et al., 2009, 2015; Mouritsen et al., 2013; Berenshtein et al., 2014). Behaviors, which include movement toward a common direction (e.g., OVM and BCRW), are expected to produce convergence of larval paths, which would be expressed in reduced pairwise distances compared to passive simulations. The following behavioral scenarios were implemented (Sim. #3–8, **Table 1**).

TABLE 1 | Biophysical simulations scenarios- model configuration and behavior of virtual larvae.

*SRW, simple random walk; CRW, auto-correlated random walk; BCRW, biased correlated random walk; OVM, ontogenetic vertical migration. In all simulations turbulence was applied as a simple random walk (SRW) component (Berenshtein et al., 2017).* \**For the schooling scenario (#8), the pairs of larvae were moving together by definition. Therefore the effective number of simulated trajectories is half of the released number of larvae.*

#### Simple Random Walk (SRW)

For this scenario, ontogenetic swimming speeds were applied according to Fisher (2005):

$$u''\_{\
u} = u''\_{\
u} + 10^{\frac{\log(t)}{\log(PLD)} \cdot \log\left(u''\_s - u''\_h\right)}\tag{1}$$

Age (t) and PLD are expressed in hours.

Hatching speeds were parameterized with (u ′′ h) 2 cm s−1 , reaching maximum swimming speeds (u ′′ s) of 20 cm s−1 , which are typical values for Pomacentridae (Leis and Carson-Ewart, 1997; Fisher, 2005). Here, the movement directions were drawn from a random distribution, therefore un-correlated and traditionally termed "Simple Random Walk" (SRW). Note, that in all the following behaviors, larval swimming speeds are modeled similarly (Equation 1), with the exception of the passive scenario, were larvae are immobile.

#### Auto-Correlated Directional Swimming (CRW)

Auto-correlated directional swimming for individual larvae was simulated as a CRW process, such that larval heading at a given time step is correlated to the previous one.

$$
\theta\_{\mathbf{i}} = \theta\_{\mathbf{i}-1} + \delta\theta \tag{2}
$$

Where θi and θi−1 are the current and previous headings, and δθ is the turning angle, which is drawn from a von-Mises distribution.

$$f(\delta\theta | kappa) = \frac{1}{2\pi I\_0(\kappa)} e^{kappa} \tag{3}$$

Where I0 is a modified Bessel function of order 0, and Kappa is the concentration parameter of the von-Mises distribution (Codling et al., 2004; Staaterman et al., 2012; Berenshtein et al., 2017). Briefly, the von-Mises Kappa controls the modeled precision of orientation, such that high Kappa values would result in directional trajectories, compared to low Kappa values. For example, Kappa = 0, 2, and 4 will result in a total decay of directional correlation after 10, 100, and 240 min respectively (Berenshtein et al., 2017). Kappa values for individual larvae were estimated from the experimental data of Irisson et al. (2015) for the con-familial Chromis atripectoralis using maximum likelihood estimation for von-Misses distribution (R core team, 2013) for each headings sequence (i.e., each individual Drifting In-Situ Chamber experiment). Following this analysis Kappa was set to 2.5.

#### Ontogenetic Vertical Migration (OVM)

The occurrence of OVM was not documented for N. miryae but was demonstrated in other pomacentrids, with larvae shown to migrate deeper with age (Paris and Cowen, 2004; Irisson et al., 2009; Huebert et al., 2011). This behavior was implemented in our model (Sim. #6, **Table 1**), such that larvae at the age of 0-5 days post-hatch (DPH) maintained a depth range of 0-25 m, and larvae at the age of 6–27 DPH maintained a depth range of 25- 100 m. This is achieved by a vertical movement speed of 2 cm s−1 . Larvae older than 27 DPH were within their competence window, and settled to the reef upon encounter.

#### Biased-Correlated Directional Swimming (BCRW)

Biased-correlated directional swimming for individual larvae was simulated as a biased correlated random walk process (BCRW) such that larval headings are drawn from a von-Mises distribution centered around a pre-defined azimuth; we used 17◦ , which is the axial direction of the GoA (i.e., larvae that swim toward this direction will end up at the northern tip of the GoA). To date, there is no empirical information about a common swimming directions of larvae in the GoA (in contrast to other locations, e.g., Bottesch et al., 2016), but using the axial direction of the GoA makes biological sense since it minimizes the coastal encounters of the larvae, and therefore minimizes their risk of predation. Kappa values for BCRW individual larvae were estimated from the experimental data of Irisson et al. (2015) (see explanation in CRW section above) resulting in kappa = 1.

#### Schooling Behavior

For schooling behavior, we simulated the most simplified scenario in which pairs of larvae were moving together keeping the exact same positions and paths. Therefore, the effective number of virtual larvae released is half from that of the other strategies (**Table 1**). We have computed the Kappa parameter of the von-Mises distribution of turning angles from Irisson et al. (2015) data for the con-familial C. atripectoralis. Mean Kappa value for CRW computed from that data equaled to 4.4 (see explanation in CRW section above). The mean terminal swimming speed of schooling larvae was set to 21.4 cm s−1 which is 7% faster than the individual terminal swimming speed (Irisson et al., 2015).

#### Passive

In addition the behaviors mentioned above, we have included a passive scenario to provide a baseline for the effect of the various behaviors on the pairwise distances.

#### Encounter Probability Estimation

For schooling behavior to occur, larvae need to encounter one another. We therefore sought to examine the probability of larvae encountering each other in the pelagic environment soon after hatching (within 1 day).

To address this question, we used two approaches. The first is spatiotemporally explicit, using 3D currents and turbulence from our biophysical model. In the second approach, we use mathematical models modified from the classical planktonic predator-prey encounter rate and probability models (Gerritsen and Strickler, 1977; Kiørboe and MacKenzie, 1995).

### Encounter Probability Using the 3D Spatiotemporally Explicit Model

The first approach raised several computational challenges. First, it required a high temporal resolution (<1 s), whereas the temporal resolution of our model is of 10 min. Second, it required realistic densities of larvae, which can be orders of magnitude higher than could be tracked in biophysical models across the entire GoA. Third, to compute all possible pairwise distances across time requires immense computational power.

We therefore created a sub-sample of 100 × 100 × 100 m (length × width × depth) of the GoA, centered around 29◦ 29.4′ N 34◦ 57.4′ E. In this space we followed either 2,600 or 26,000 passive virtual larvae, which correspond to the lower and midlevel values of measured densities of the con-familial Stegastes partitus larvae (<6 DPH, Paris and Cowen, 2004). The larvae, were released on the 15th of March 2004 at 12:00 p.m., and tracked for 10 min at a temporal resolution of 0.1 s. Initial positions of the passive larvae were generated at levels of patchiness [Lloyd's index of patchiness (L); (Lloyd, 1967)] that ranged from random distribution (L = 1) to moderately clumped (L = 40) distribution (Bradbury et al., 2003).

Random distributions were generated by selecting random XYZ coordinates within the cubic domain. Clumped distributions were generated using an iterative algorithm which produced the initial positions of virtual larvae around predefined number of centroids (Np1 = 20), randomly drawn from a uniform distribution within a pre-defined distance (Rp1 = 20 m) from a given patch centroid. For each iteration, to compute Lloyd's index of patchiness, 100 samples (7 × 7 × 7 m) were randomly placed in the domain. Then Npi and the Rpi were iteratively modified until reaching the desired range of patchiness. The sample's volume (343 m3 ) is within the range of typical ichthyoplankton samples (Kimmerling et al., 2017). Visualization examples of the spatial patterns are presented in Supplementary Figure S1.

We assumed a detection distance of 0.1 m, which corresponds to the visual acuity in fish larvae (Renee Lara, 2001). Encounter probability was computed as the proportion of larvae, out of the total number released, that came within 0.1 m of at least one other larvae throughout the 10 min simulation.

## Encounter Probability Using the Classic Encounter Model

For the classic approach we used the encounter rate model from Kiørboe and MacKenzie (1995) between planktonic predator and prey. Here the encounter rate (Er) was a function of larval density (D; using 0.026 and 0.0026 ind m−3 ), detection distance (R; using 0.1 m), and swimming speed (v; using values of 0 and 0.02 m s−1 ; Fisher, 2005).

$$Er = \pi DR^2 \sqrt{(2\nu^2 + 2\nu^2)}\tag{4}$$

Turbulent velocity (w) was computed following Kiørboe and MacKenzie (1995), using turbulent dissipation rates (ε) of 10−7 and 10−4 m2 s −3 which represent the range of turbulent dissipation rates for the shelf and coastal zones (Kiorboe and Saiz, 1995).

$$
\omega = 1.9(\epsilon \text{R})^{1/3} \tag{5}
$$

Encounter probability was computed following Gerritsen and Strickler (1977), where encounters were assumed to be randomly occurring events (Poisson distribution). Therefore the probability encountering at least one larvae (PE) was formulated as a function of the time interval (t; using 10 min and 24 h).

$$P\_E \left(\Delta t\right) = 1 - e^{-Er\Delta t} \frac{\left(Er\Delta t\right)^0}{0!} = 1 - e^{-Er\Delta t} \tag{6}$$

The term e −Ert in Equation (6) represents the probability encountering 0 larvae.

### Assessment of the Capacity of Newly Hatched Larvae to Overcome Turbulence and Remain Grouped

To assess the capacity of adjacent larvae to overcome horizontal turbulence and stay grouped, we computed turbulent velocity using Kiørboe and MacKenzie (1995) approach (Equation 5), and compared it to the combined velocity of two larvae swimming one toward the other.

# Assessment of Grouping Tendency From Kimmerling et al. (2017) Data

The grouping tendency of N. miryae and Chromis viridis was assessed using the sampling data from Kimmerling et al. (2017). We focused on these two species as they were the species examined in Ben-Tzvi et al. (2012), with N. miryae demonstrating cohesive dispersal pattern, while C. viridis did not. The frequency and number of individuals occurring in the ichthyoplankton samples was recorded, and their Lloyds Index of Patchiness was computed.

# Settlement Success and Connectivity Network Density

Settlement-success was computed as the percent of virtual larvae that have reached a reef polygon within their competence window (last 48 h of their PLD; Berenshtein et al., 2017). Networkdensity was computed as the percent of realized connections out of all possible connections between the different polygons in the GoA domain. Network-density considered only whether or not connections occurred, without considering the connections strength.

Settlement success and network density were computed per release event (N = 9; Sim #5 and #8, **Table 1**). Paired t-tests were used to examine the significance of differences in terms of settlement-success and network-density between schooling vs. individual CRW larvae. Shapiro-Wilk test was used to examine the normality of the paired differences. All simulations and statistical analyses were run on MATLAB R2014b (The MathWorks Inc., Natick, MA, USA) and R (R core team, 2013).

# RESULTS

# Backward and Forward-in-Time Passive Simulations

Backtracking of passive virtual larvae from the sites and times of settlement by N. miryae (**Figure 2A**) resulted in a wide spread of trajectories around the northern two-thirds of the GoA; with a high degree of overlap between cohorts (**Figure 2A**). Similarly, forward-in-time trajectories, released from reef sources from the backtracking simulations, did not show any channeling patterns toward the observed settlement sites, but instead were widespread across most of the GoA area (**Figure 2B**). Both backward and forward in-time trajectories showed no evidence for prevailing currents that could generate cohort cohesion (Ben-Tzvi et al., 2012).

#### Pairwise Distances and Probability of Adjacent Settlement

The analysis of pairwise-distances indicated that all but the schooling scenarios resulted in a high degree of dispersion. Mean pairwise distances at the end of the PLD ranged between 13 and 38 km, for pairs of larvae that hatched at the same location and time (**Figure 3**). Behaviors which include movement toward similar horizontal or vertical directions among larvae (BCRW and OVM) resulted in path convergence in comparison to the passive scenario, while dispersive behaviors which lacked "common movement direction" (SRW and CRW) resulted in greater path divergence in comparison to the passive scenario. Yet, even for the most convergent scenario (BCRW), the pairwise distances were far from generating patterns of cohesive dispersal, i.e., cohort members that settled within a few meters one from another (Ben-Tzvi et al., 2012). By definition, pairwise-distances under the schooling behavior scenario equaled zero. All behaviors are characterized by a nearly linear increase of pairwise distances with slopes slightly decreasing toward the end of the PLD due to the bounded domain of the GoA. An exception is BCRW, which shows a considerably earlier and steeper decline in the pairwise distances due to convergence of larvae toward the northern tip of the GoA.

Similarly, the probability of sampling a nearest neighbor pair which hatched at the same location and time for all but the "schooling behavior" strategy was extremely low (<0.00004). Naturally, in the schooling behavior scenario, as the larval pairs hatched, moved, and settled together, the probability for pairwise settlement is maximal.

#### Larval Encounter Probabilities

The probability of larval encounter per 10 min, from the 3D biophysical model, ranged between 0.0046 and 0.49 at low larval densities, and 0.046-0.78 at high larval densities, depending the degree of initial larval patchiness (**Figure 4**; for spatial distribution see Supplementary Figure S1).

Encounter probabilities computed with the classic encounter model (Kiørboe and MacKenzie, 1995) ranged from 0.001 to 0.99

Magenta asterisks in (A,B) represent particle release locations. (C) Pelagic Larval Durations (PLD) of the recruiting cohorts from April-March 2004, which were deducted from otolith analysis in Ben-Tzvi et al. (2012). Maps in (A,B) were generated using MATLAB R2014b (The MathWorks Inc., Natick, MA, USA), and are rotated views of the GoA (73◦ Clockwise).

depending on larval density, time span, turbulent velocity, and swimming speeds (**Table 2**). Across 24 h, for high larval densities, the encounter probability of swimming larvae is 0.2 and 0.99 for low and high turbulent dissipation rates respectively (**Table 2**, Supplementary Figure S2). Increase in larval density, time span, turbulent velocity, and swimming speeds resulted in an increase in encounter probabilities.

# Assessment of the Capacity of Newly Hatched Larvae to Overcome Turbulence and Remain Grouped

Computation of turbulent velocity (Equation 5) for dissipation rate of 10−4 m2 s −3 resulted in a turbulent velocity of w = 0.041 m s−1 . Given maximal larval speeds during the first 24 h after hatching is 0.02–0.03 m s−1 , larvae that swim toward

FIGURE 3 | Distances between pairs of virtual larvae that hatched at the same locations and times throughout the pelagic larval duration (PLD) for the following strategies: passive, simple random walk (SRW), auto-correlated directional swimming (CRW), ontogenetic vertical migration (OVM), biased correlated random walk (BCRW), and schooling behavior (simulation #3–8, Table 1). Shaded areas indicate 99% confidence intervals (*N* = 32,046 and 28,414 for the passive and other strategies respectively).

each other would be able to overcome turbulent velocity and remain grouped. Dissipation rate of 10−4 m2 s −3 represents the upper limit of coastal and shelf dissipation rates, suggesting that such larvae, which can swim since hatching, will be able to overcome turbulent velocity and remain grouped at nearly all

29◦29.4′ N 34◦57.4′ E during the 15th of March 2004 at 12:00–12:10 p.m.

# oceanic conditions, from early stages of their PLD.

# Assessment of Grouping Tendency From Kimmerling et al. (2017) Data

For N. miryae, 90% of the samples contained two or more larvae, compared to only 29% in C. viridis. Lloyd's index of patchiness for N. miryae was 64 compared to 27 for C. viridis.

# Settlement Success and Connectivity Network Density

Settlement success of schooling larvae was significantly higher than that of non-schooling larvae (pairwise t-test, t = 26.3, df = 8, p < 0.01, **Figure 5**). Connectivity network density in schooling larvae was significantly lower than that achieved in the non-schooling scenario (pairwise t-test, t = −35.2, df = 8, p < 0.01, **Figure 5**).

# DISCUSSION

Of the simulation scenarios considered in the current study, patterns of cohesive dispersal could only be reproduced by schooling behavior (**Figure 3**). This would seem to suggest that newly hatched N. miryae may actively form long-lived schools.

To school, larvae must encounter one another, yet schooling does not seem to occur immediately after hatching for the following reasons: (1) There is no indication for synchronous hatching from Pomacenridae demersal eggs, (2) There is no genetic structure between larval cohorts (Ben-Tzvi et al., 2012), and (3) Larval swimming velocities at hatching are not necessarily sufficient to overcome turbulent velocity. Therefore, N. miryae encounter their conspecifics shortly—but not immediately—after hatching.

The larval encounter analyses computed in this study demonstrated feasible probabilities for pelagic larval encounter within 24 h (**Table 2**, **Figure 4**). However, both encounter computation methods encompass some limitations. The mathematical model assumes larval distributions are random rather than clumped. The latter has been documented for fish larvae (e.g., Bradbury et al., 2003, L>3; Frank et al., 1993, L>2), and for later-stage larvae of N. miryae in the GoA (Kimmerling

TABLE 2 | Sensitivity analysis of larval encounter rates and probabilities as a function of larval density, time span, larval swimming speeds, and turbulent velocities.

*Encounter rates and turbulent velocities were computed according to Kiørboe and MacKenzie (1995), and encounter probabilities were computed according to Gerritsen and Strickler (1977), for larval encounter with at least one other larva. Rows are sorted according to increasing encounter probabilities values. For graphical representation of this data see Supplementary Figure S2.*

et al., 2017, L = 64). The 3D—spatio-temporally explicit model is limited as well, since advection, and turbulence transport the larvae away from their initially confined cubic domain to waters which are empty of larvae. This alters the shape and the volume of the waters encompassing the larval patch, resulting in a modification of the effective densities across time. Despite these limitations, both methods indicate that larval encounter within the first 24 h post-hatching can be feasible.

Even when larvae encounter each other they still need to overcome turbulent velocity to keep together. The fact that larval swimming velocities soon after hatching were similar to the upper range of turbulent velocities (Dissipation rate of 10−4 m3 s −2 ; **Table 2**) suggests that fish larvae can maintain schooling in most oceanic conditions.

To keep together, larvae may utilize their vision and lateral line senses similarly to adult fish (Partridge and Pitcher, 1980), and their olfactory sense similarly to settlement stage larvae (Lecchini et al., 2014). They may also use their auditory sense, as fish larvae can produce sounds during the night (Staaterman et al., 2014).

The fact that larvae within cohorts stayed together through their PLD (Ben-Tzvi et al., 2012) suggests a high cohort fidelity and a lack of (or little) mixing with other cohorts or sporadic individuals (i.e., cohesive dispersal). This, in turn, suggests that N. miryae larvae would form scarce and clumped patches in the pelagic environment.

The spatial scale of cohesive dispersal is unknown due to the limited information about the spatial heterogeneity of traceelements used in Ben-Tzvi et al. (2012). However, using the available information, we can make a simplistic assessment of the spatial extent of newly formed cohorts, which results in length scales of tens to hundreds of meters (Supplementary Data S1). At these small scales, convergent oceanographic features, which could potentially aggregate larvae but are too small to be resolved by the MITgcm model, could occur (e.g., Gildor et al., 2009). Yet, such features are rare and persistent across no more than a few days, and therefore cannot produce the observed cohesive pattern across the PLD (Ben-Tzvi et al., 2012).

Schooling behavior, along with other biological and physical factors (e.g., Fronts, pycnoclines, and heterogeneous distribution of larva predator and prey) can contribute to larval patchiness (e.g., McGurk, 1986; Frank et al., 1993). Ichthyoplankton samples from Kimmerling et al. (2017) suggest a stronger tendency of N. miryae to be found in groups compared to C. viridis, whose otolith analysis lacked any evidence for cohesive dispersal (Ben-Tzvi et al., 2012). Note however, that the larvae of these two species are characterized by different spatio-temporal dynamics which may also effect their tendency to be found in groups (Kimmerling et al., 2017). Naturally, Interplay between biological and physical factors can effect schooling dynamics, as more patchy populations have a higher probability of encounter, and therefore facilitate schooling.

Previous studies demonstrated that the extent of patchiness often follows a "U" shaped curve across ontogeny (for example in Capelin and Sandlance larvae/eggs; Bradbury et al., 2003). Initial high patchiness is attributed to the coherence of larvae/eggs patches from mass-spawning/hatching at specific areas; the following decrease is attributed mainly to advection and mortality; and the increase toward the end of the PLD is attributed to active larval behavior (Bradbury et al., 2003). While Frank et al. (1993) attributed that increase to the occurrence of prey (Sagitta elegans), schooling behavior may have also played a role in these dynamics.

In our current biophysical model, mortality was not applied. In other modeling efforts, mortality is mostly applied as a constant rate in time and space due to limited empirical information. Such application of mortality is not expected to effect the results and the conclusions of this study (e.g., Berenshtein et al., 2017).

Increased schooling capacity can affect fitness both positively and negatively. Positive effects include increased swimming speeds and orientation capacities (Irisson et al., 2015), which are manifested in increased settlement success (**Figure 5**; Berenshtein et al., 2017); while negative effects may include higher probability of predator encounter due to increased mobility (swimming velocity and directional persistence; Visser and Kiørboe, 2006), increased competition for food, and the reduction of variation in effective dispersal paths (i.e., connectivity; **Figure 5**).

The increased settlement success of schooling larvae is attributed to higher swimming speeds and a more precise orientation (Irisson et al., 2015), which lead to an increase in displacement associated with larval behavior (Berenshtein et al., 2017). This, in turn, allows the larvae to depart from entraining currents (e.g., sub-mesoscale eddies) and increases their probability of reaching the coast and settling (Berenshtein et al., 2017). The decrease in connectivity network density (**Figure 5**) is attributed to the reduction of the variation in dispersal trajectories. In other words, a population of larvae which disperses in pairs for example, can produce half of the dispersal paths generated by individually dispersing larvae.

These changes in settlement success and connectivity are important as these factors are central in population dynamics, affecting the biogeography and the genetic structure of populations in multiple marine species (Cowen and Sponaugle, 2009). Moreover, larval supply from external (non-local) sources is essential for sustaining populations, i.e., sites which receive larvae from multiple sources would be less susceptible for a population collapse compared to sites which are dependent on fewer sources (Cowen et al., 2006).As schooling in fish larvae might be very common (reviewed in Leis, 2006), in-situ and laboratory studies are needed to examine ontogenetic larval collective behavior, and specifically their fission-fusion dynamics across their PLD. This, in conjunction with ichthyoplankton analyses of species abundances throughout ontogeny, could set the basis for a realistic incorporation of schooling behavior in biophysical models.

In conclusion, this study supports the occurrence of schooling of N. miryae in the GoA across the PLD for realistic reconstruction of observed settlement patterns. In addition, our study demonstrates the effects of schooling behavior on larval settlement success and connectivity. The dynamics and characteristics of larval schooling should be further studied, and its application should be considered in larval dispersal models when managing fisheries and determining marine protected areas.

# AUTHOR CONTRIBUTIONS

IB, HG, YA, EF, and MK participated in data analysis, EF adjusted and run the atmospheric model, and helped with the particle tracking algorithm, HG and YA adjusted and run the oceanographic model, IB, CP, and MK participated in conceiving, designing, and drafting the study.

## FUNDING

This study was supported by the United States-Israel Bi-national Science Foundation (BSF grant 2008/144 to MK 1210 and CP).

# ACKNOWLEDGMENTS

We thank Dr. Ofer Ben-Tzvi for providing details about the field sampling. We thank Dr. Robin Faillettaz for support with graphics. We thank the IUI staff and the Paris lab for their

#### REFERENCES

contribution in various aspects of this study. We greatly thank two reviewers who substantially improved the quality of our manuscript. This study was supported by the United States–Israel Bi-national Science Foundation (BSF grant 2008/144 to MK and CP).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2018.00254/full#supplementary-material

**Conflict of Interest Statement:** 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.

Copyright © 2018 Berenshtein, Paris, Gildor, Fredj, Amitai and Kiflawi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.