On the Spread of Microbes That Manipulate Reproduction in Marine Invertebrates

Bacterial symbionts are functionally integral to animal reproduction and development, some of which have evolved additional mechanisms to override these host programs. One habitat that is increasingly recognized to contain phylogenetically related lineages of reproductive manipulators is the ocean. The reproduction of marine invertebrates often occurs by free spawning instead of by the physical contact of copulation in terrestrial systems. We developed an integrated model to understand whether and when microbes that manipulate host reproduction by cytoplasmic incompatibility, feminization, and male killing spread within populations of free-spawning marine invertebrates. Our model supports three primary findings. First, sex ratio distortion leads to suboptimal fertilization and zygote production in planktotrophs (feeding larvae) but enhance these processes in lecithotrophs (nonfeeding larvae). Second, feminization and a combination of male killing plus enhanced growth are effective at spreading reproductive manipulators while also inducing a female-biased sex ratio. Third, the majority of free-spawning marine invertebrates could be infected across a range of life history combinations, with infections harming species with smaller eggs and longer pelagic durations while benefiting species with larger eggs and shorter pelagic durations. Together, this supports the general premise that microbes may manipulate the reproduction of free-spawning marine invertebrates (e.g., by inducing changes in developmental life history) and that these types of manipulations overlap considerably with terrestrial systems.


Introduction
Bacterial symbionts are functionally integral to animal biology, including being central to cell and tissue differentiation, reproduction, and development (McFall-Ngai 2002;Bates et al. 2006;Bright and Bulgheresi 2010;Fraune and Bosch 2010;McFall-Ngai et al. 2013;Gilbert et al. 2015). Microbes that are functionally intertwined with animal development are often beneficial to the host (McFall-Ngai 2002;Gilbert et al. 2015;Bosch and McFall-Ngai 2021). Some of these bacteria, however, have evolved additional mechanisms to override components of the host's developmental program Engelstädter and Hurst 2009a;Perlmutter and Bordenstein 2020). These microbes are often inherited through the cytoplasm of the egg and can manipulate host reproduction to favor their transmission through cytoplasmic incompatibility (where microbe-containing males are incompatible with aposymbiotic females), feminization (the conversion of genetic males to females), and male killing (microbe-induced male mortality during early or late development; Engelstädter and Hurst 2009a;Perlmutter and Bordenstein 2020;Kaur et al. 2021).
Hosts affected by microbes that manipulate reproduction are diverse and are particularly common in terrestrial nematodes and arthropods, where 160% of species are thought to be infected (Stevens et al. 2001;Duron et al. 2008;Hilgenboecker et al. 2008). The most widely recognized reproductive manipulator is Wolbachia, but bacteria within Arsenophonus, Cardinium, Rickettsia, and Spiroplasma can also interfere with host reproduction Engelstädter and Hurst 2009a;Perlmutter and Bordenstein 2020). These manipulators live extracellularly, freely inside host cells, or within vacuoles, where they can use the host's cellular machinery to induce their respective phenotypes (Skinner 1985;Callaini et al. 1994;Hurst et al. 1996;Ammar and Hogenhout 2006;Kitajima et al. 2007;Werren et al. 2008;Engelstädter and Hurst 2009a).
While much is known about these symbioses in terrestrial systems, the ocean is increasingly recognized to be home to microbes that are phylogenetically related to terrestrial reproductive manipulators (Morris et al. 2002;Perlman et al. 2006;Sunagawa et al. 2015;Harumoto and Lemaitre 2018). Dozens of marine invertebrate lineages have been found to associate with Rickettsiales, and some are believed to interfere with host reproduction (Klinges et al. 2019;Carrier et al. 2021). Marine invertebrates often reproduce by free spawning their gametes into the water column, where fertilization, embryogenesis, and larval development take place (Thorson 1950;Mileikovsky 1971;Levitan 1995). This is in contrast to terrestrial arthropods that reproduce by the direct physical contact of copulation (Engelmann 2015). Despite these vastly different life history tendencies (Strathmann 1990), the fitness of marine reproductive manipulators remains dependent on their transmission across generations. We therefore hypothesize that population-level consequences of microbe-mediated reproductive manipulations should be convergent between terrestrial and marine animal hosts (Engelstädter and Hurst 2009a;Burgess et al. 2015).
One such marine animal host is Heliocidaris, a genus of sea urchin that includes one of the most comprehensively studied transitions between the major developmental strategies among marine invertebrates (Thorson 1950;Strathmann 1985;Wray and Raff 1989;Wray and Raff 1991b;Byrne et al. 1999;Israel et al. 2016). A speciation event ∼5 million years ago resulted in sister species with alternative life history strategies: H. tuberculata is planktotrophic, while H. erythrogramma is lecithotrophic (Wray and Raff 1991a;Zigler et al. 2003). Typical of planktotrophs, H. tuberculata develops from small eggs (!100 mm in diameter) into feeding larvae that disperse for several weeks, while H. erythrogramma develops from eggs ∼53 to 86 times the volume of H. tuberculata eggs into nonfeeding larvae that lack the morphological structures for feeding and remain in the water column for 5 days (Emlet and Hoegh-Guldberg 1997;Wray and Raff 1991a).
The lecithotrophic H. erythrogramma acquired a cytoplasmically inherited Rickettsiales-Echinorickettsia raffiiin this life history switch that is most closely related to Wolbachia pipientis wMel and wAlbB (Carrier et al. 2021). The genome of E. raffii suggests that this bacterium is a nutritional mutualist and a reproductive manipulator (Carrier et al. 2021). In the former, the host is hypothesized to benefit from being provisioned essential amino acids to enhance growth, as observed in terrestrial systems (Douglas 1998;Hosokawa et al. 2010). In the latter, this bacterium is hypothesized to use effector proteins to influence the rate of male mortality and modulate fertilization (Carrier et al. 2021). Supporting the hypothesis that E. raffii manipulates host reproduction, populations of H. erythrogramma in Sydney Harbor (67.2%) and Tasmania (56.9%) are disproportionately female (Dix 1977;Carrier et al. 2021). Moreover, the youngest reproductive individuals in Sydney Harbor do not deviate from a 1∶1 female-to-male ratio, while the largest individuals have a 4∶1 female-to-male ratio (Carrier et al. 2021). There is currently no evidence that E. raffii associates with the planktotrophic H. tuberculata, indicating potential differences in whether and how reproductive manipulators may spread in these developmental life history strategies.
Understanding whether microbes that manipulate reproduction can spread and induce shifts in the population structure of marine invertebrates requires the integration of two model types. The first concerns the fertilization dynamics and zygote production of free-spawning marine invertebrates (e.g., Levitan 1991Levitan , 1993Levitan , 2004Styan 1998), and the second involves the spread and population-level consequences of cytoplasmic incompatibility, feminization, and male killing (e.g., Jansen et al. 2008;Engelstädter and Hurst 2009b;Hancock et al. 2011). Accounting for the difference in reproductive strategy between terrestrial arthropods and free-spawning marine invertebrates allows us to estimate how sex ratio distortion influences reproductive success and the likely type(s) of manipulations that occur in the sea. We first apply this general integrated model to Heliocidaris. We then develop a size-structured population model to test how a reproductive manipulator could induce a female-biased sex ratio, as is observed in H. erythrogramma (Carrier et al. 2021). Last, we develop a general model to assess when infections occur based on key components of marine invertebrate life history, limitations to such infections, and which phyla and developmental modes are particularly infectable.

Fertilization Dynamics in Free-Spawning Marine Invertebrates
Fertilization success of free-spawning marine invertebrates is highly variable and is based on many ecological factors, resulting in the possibility for eggs within a clutch to experience both sperm limitation and polyspermy (Levitan 1995;Styan 1998;Evans and Lymbery 2020). Various gamete characteristics (e.g., egg size and sperm velocity) as well as ambient egg and sperm concentrations can be used to determine the likelihood that sperm limitation or polyspermy occur. To understand how spawning density and sex ratio may influence fertilization success and zygote production, we modified the Styan (1998) polyspermy kinetic model by converting total egg (E T ; egg/mL) and sperm (S T ; sperm/mL) concentrations to be explicitly determined by female and male density. The average number of sperm that may potentially fertilize an egg during a free-spawning event, x (eq. [5] from Styan 1998), is a function of fertilization efficiency (F e ), total sperm concentration (S T ; sperm/mL), total egg concentration (E T ; egg/mL), the biomolecular collision constant (b 0 ; mm 3 /s), and the sperm half-life (t; s): Moreover, the biomolecular collision constant (b 0 ) is determined by the cross-sectional area of the egg (j; mm 2 ) and sperm velocity (n; mm/s; eq. [2] from Styan 1998): The biomolecular collision constant and the time to elicit a block to polyspermy (t b ; s) may then be used to calculate the average number of extra sperm that may contact an already fertilized egg, which then results in polyspermy (b; eq. [13] from Styan 1998): The proportion of eggs that will be fertilized by a single sperm and successfully develop may be estimated by the following (eq. [16] from Styan 1998): Zygote density (N z,t ; zygotes/mL) may then be estimated on the basis of the probability of a monospermic zygote (eq. [4]) and total egg density (E T ): A summary of the variables and parameters used in this fertilization model are presented in table 1. Equations (4) and (5) are available in the web application at https:// mck8dg.shinyapps.io/SI_sex_ratio_distortion_marine_inverts/, where the various reproductive parameters may be altered to explore how reproductive output is influenced.

Spread of Reproductive Manipulators: No Population Structure
We developed a population model without size structure that integrated the fertilization model described above to simulate how a microbe that manipulates reproduction via either cytoplasmic incompatibility, feminization, or male killing may spread within a population of free-spawning marine invertebrates. Each scenario had four distinct groups: uninfected females (N u,f ; individuals/m 2 ), infected females (N i,f ; individuals/m 2 ), uninfected males (N u,m ; individuals/ m 2 ), and infected males (N i,m,t ; individuals/m 2 ). We used a discrete time population model with overlapping generations-where reproduction, larval development, and adult mortality are assumed to occur in that order-because ma-rine invertebrates tend to spawn seasonally (e.g., Giese 1959;Mileikovsky 1971;Strathmann 1987). Larval mortality (M z ) is assumed to be dependent on the length of the larval period and the instantaneous mortality rate (m z ), which we assumed to be constant between generations and density independent (Morgan 1995;Rumrill 1990): Adult mortality (M a ), on the other hand, is assumed to be density dependent and based on the maximum adult mortality rate (m a ), carrying capacity (k; individuals/m 2 ), total density of the population at time t (N t ; individuals/m 2 ), and the degree to which density influences mortality (i.e., the shape parameter; d). Moreover, both the current adult density and density of surviving zygotes that settled were assumed to influence density dependence, where N 0 z,t is the density of zygotes (zygotes/mL) that survived after microbe-induced mortality (i.e., during cytoplasmic incompatibility or male killing) has occurred and w (mL/m 2 ) is the settlement constant that describes the proportion of surviving zygotes that settled: The density of uninfected females in the next generation (N u,f ,t11 ; individuals/m 2 ) may then be determined by the density of surviving uninfected female zygotes (N z,u,f,t ; zygotes/mL) that settled plus the density of surviving uninfected female adults at generation t (N u,f,t ; individuals/m 2 ): The density of uninfected female zygotes (N z,u,f,t ) is directly influenced by cytoplasmic incompatibility-but not male killing or feminization-and is calculated by multiplying the number of zygotes (N z,t ) by the proportion of eggs from uninfected adult females, where c is the fecundity cost to infected females. We assumed that uninfected females produced an equal number of daughters and sons, and therefore the density of zygotes (N z,t ) is divided by two. If cytoplasmic incompatibility occurs, then the density of surviving zygotes is also dependent on the proportion of sperm from infected males and the cytoplasmic incompatibility mortality rate (m c ). Note that both the slope of eggs (E; eggs/mL/individual/m 2 ) and the slope of sperm density (S; sperm/mL/individual/m 2 ) cancel out in this equation: The density of uninfected males in the next generation (N u,m,t11 ; individuals/m 2 ) has the same general recursion equation (i.e., eq. [8]), and the density of uninfected male zygotes (N z,u,m,t ; zygotes/mL) is the same as the density of uninfected female zygotes (N z,u,f,t ) in equation (9).
Since all groups in this model have the same general recursion equations, the density of infected females in the next generation (N i,f ,t11 ; individuals/m 2 ) has the same basic structure as equation (8). The density of infected female zygotes (N z,i,f,t ; zygotes/mL), however, is unaffected by either cytoplasmic incompatibility or male killing and is only directly affected by the feminization rate. The density of infected female zygotes is calculated by multiplying the density of zygotes (N z,t ) by the proportion of infected adult females and the sex ratio biasing constant (r), which ranges from 0.5 (no sex ratio distortion) to 1.0 (complete feminization): The density of infected males in the next generation (N i,m,t11 ; individuals/m 2 ) also has the same basic structure as equation (8). The density of infected male zygotes (N z,i,m,t ; zygotes/mL) is influenced by both feminization and male killing, but not cytoplasmic incompatibility. The density of infected male zygotes has the same basic structure as equation (10), except that the complement of the sex ratio biasing constant (r) is used while also accounting for mortality due to male killing (m k ): A summary of the variables and parameters are provided in table 2, and a diagram of this model is outlined in figure S1.

A Case Study Simulation: Heliocidaris without Population Structure
We apply our model to Heliocidaris for three main reasons. First, populations of H. erythrogramma associate with an endosymbiont that is hypothesized to distort host sex ratio, and this is the only marine invertebrate known to have this type of partnership (Carrier et al. 2021). Second, these two Heliocidaris species are representative of the primary  Longo et al. 1986. developmental strategies of marine invertebrates and are very well studied, allowing us to apply species-specific values for most parameters of our model (Thorson 1950;Strathmann 1985). When species-specific parameters were unavailable, we parametrized with values from related echinoids (table 1). Finally, these species differ vastly in egg and clutch size, allowing for a broader understanding of how host reproductive biology influences the spread of microbes that manipulates the reproduction of free-spawning marine invertebrates.
Using the model described above, we simulated whether a reproductive manipulator could invade a population with each developmental life history when introduced (at 1%) using the three modes of reproductive manipulation described above in isolation. First, to model cytoplasmic incompatibility in isolation, m c was set to 0.7, m k was set to 0, and r was set to 0.5. Second, for feminization, m c and m k were set to 0 and r was set to 0.7. Last, for male killing, m c was set to 0, m k was set to 0.7, and r was set to 0.5. Pre-liminary analyses varying these values resulted in qualitatively similar results (web application at https://mck8dg .shinyapps.io/SI_sex_ratio_distortion_marine_inverts/).

A Case Study Simulation: Heliocidaris with Population Structure
We then developed a stage-based population model to determine whether and how characteristics of a reproductive manipulator (i.e., as a nutritional mutualist and reproductive manipulator) may influence the population biology and structure of H. erythrogramma. In this, we simulated the infection and spread of a reproductive manipulator in three size classes (adult diameters of 20-39 mm, 40-59 mm, and ≥60 mm), both sexes, and uninfected and infected individuals (table 3). These size classes were chosen because they represent the naturally occurring size distribution of reproductive H. erythrogramma (Williams and  1975;Dix 1977;Laegdsgaard et al. 1991;Ling et al. 2019;Carrier et al. 2021). We also performed this simulation in H. tuberculata as a proof of concept. This stage-based population model had three phases: birth, death, and transitions to the next size class. Birth followed the same fertilization dynamics described by equation (5), where the total egg (E t ; eggs/mL) and sperm (S t ; sperm/mL) concentration were determined by summing the number of each group that contributed to a successful fertilization. We assumed that the largest individuals (160 mm) produced the same number of gametes as the species-specific parameters used for the model without population structure (table 1). Individuals in the 20-39-mm and 40-59-mm size classes produced a fraction of the gametes of the largest individuals (d 1 and d 2 , respectively) that is based on the size-fecundity relationship described in equation (2) in Levitan (1991): log volume of eggs ðmLÞ p 3:56 log size ðmmÞ 2 6:59: We used median diameters for each size class (i.e., 30 mm, 50 mm, and 60 mm) in this equation, and all other reproductive parameters were assumed to not differ between size.
The density of eggs of a given infection status (E g ; eggs/ mL) is determined by the density of females (individuals/ m 2 ) in the different size classes for a given infection status (indicated by the subscript g) weighted by their respective size class contributions (d 1 and d 2 ) multiplied by the egg concentration slope E (eggs/mL/individual/m 2 ): Moreover, the density of sperm of a given infection status (S g ; sperm/mL) has the same basic structure as equation (13) except that it uses the sperm concentration slope S (sperm/ mL/individual/m 2 ): The density of uninfected female zygotes (N z,u,f,t ; zygotes/ mL) is directly influenced by cytoplasmic incompatibilitybut not male killing-and is calculated by multiplying the number of zygotes (N z,t ; zygotes/mL) by the proportion of uninfected female eggs. We assumed that uninfected females will produce equal daughters and sons, and therefore the number of zygotes (N z,t ) is divided by two. If cytoplasmic incompatibility occurs, then the density of surviving zygotes is also dependent on the cytoplasmic incompatibility mortality rate (m c ) and the proportion of infected adult male sperm: Furthermore, the density of uninfected male zygotes (N z,u,m,t ; zygotes/mL) has the same basic structure as the density of uninfected female zygotes (N z,u,f,t ).  The density of infected female zygotes (N z,i,f,t ; zygotes/ mL) is calculated by multiplying the number of zygotes (N z,t ) by the proportion of infected female eggs. We also assumed equal sex ratios at birth: The density of infected male zygotes (N z,i,m,t ; zygotes/ mL) is directly influenced by male killing and has the same basic structure as equation (16) except that mortality due to male killing (m k ) is accounted for: Surviving zygotes that settle enter the smallest size class, given their respective sex and infection status, indicated by the subscript group index, g. The density of 20-39-mm individuals in a given group (sex and infection status) in the next time step (N g,1,t11 ; individuals/m 2 ) is determined by the density of surviving adults of a given group (N g,1,t ; individuals/m 2 ) that did not transition plus the density of surviving zygotes of a given group (N z,g,t ; zygotes/mL) that settled, where M z is the larval mortality (eq. [6]), w is the settlement constant (mL/m 2 ), M a is the adult mortality (eq. [7]), and g is the transition rate to the next size class: The density of 40-59-mm individuals in a given group (sex and infection status) in the next time step (N g,2,t11 ; individuals/m 2 ) is determined by the density of surviving 20-39-mm stage class adults of a given group (N g,1,t ; individuals/m 2 ) that transition and the density of surviving 40-59-mm adults of a given group (N g,2,t ; individuals/m 2 ) that do not transition to the next size class: The density of ≥60-mm individuals in a given group (sex and infection status) in the next time step (N g,3,t11 ; individuals/m 2 ) is determined by the density of surviving 39-40-mm stage class adults of a given group (N g,2,t ; individuals/m 2 ) that transition and the density of surviving ≥60-mm adults of a given group (N g,3,t ; individuals/m 2 ): To model enhanced growth rate in each size class of infected females, g is replaced with g i such that g i 1 g. We numerically simulated this stage-based model using species-specific parameters of both Heliocidaris species for three different scenarios: (i) enhanced growth for infected females only (g i p 0:25, g p 0:1, m k p 0, m c p 0), (ii) enhanced growth for infected females plus male killing (g i p 0:25, g p 0:1, m k p 0:7, m c p 0), and (iii) enhanced growth for infected females plus male killing and cytoplasmic incompatibility (g i p 0:25, g p 0:1, m k p 0:7, m c p 0:7). As a sensitivity analysis on how enhanced growth and male killing parameters influenced model outcomes, we simulated the model at all possible combinations of male-killing rate from 0 to 0.99 at increments of 0.01 and enhanced growth from 0.11 to 0.35 at increments of 0.01.
Summary of the variables and parameters used in this model are in tables 1-3, and a diagram of this model is outlined in figure S2. All simulations were performed using the deSolve package in R (Soetaert et al. 2010), and this may be simulated independently using the web application at https://mck8dg.shinyapps.io/SI_sex_ratio_distortion _marine_inverts/.

Numerical Analyses of Reproductive Manipulators and Life History Traits
The Heliocidaris species used in this case study are broadly representative of the two major developmental strategies of free-spawning marine invertebrates. They do not, however, fully encapsulate the enormous variety of life history combinations found in the sea (Levin and Bridges 1995;Byrne 2006;Allen and Pernet 2007). We complemented our Heliocidaris case study (with no population structure) by performing numerical simulations that vary egg size (j), pelagic larval duration (t d ), and the cost of infection (c) to estimate which life history combinations may be infected by microbes that manipulate reproduction. We focus this analysis on feminization and male killing plus enhanced growth because our Heliocidaris case study indicated that these are the most likely types of reproductive manipulation for free-spawning marine invertebrates. We used the model without population structure because we were concerned with the assumptions needed to generalize the size classes of the population structure model across all of the given life history combinations. To model enhanced growth in the model without population structure, we made the cost of infection (c) negative for females associated with the microbe (i.e., infection increases fecundity).
To determine which life history characters influence host infectability, we simulated population runs at different combinations of egg diameter (from 40 to 1,500 mm with 10-mm increments) and pelagic larval duration (from 0 to 80 days with 1-day increments). The range of these parameters spans much of the diversity known in freespawning marine invertebrates (Marshall et al. 2012;Álvarez-Noriega et al. 2020). We initialized populations at an equal sex ratio with a total density of two individuals per square meter and an infection rate of 1%. We simulated populations until equilibrium was reached and set a maximum of time steps to 100,000 to account for the cycles that could happen in discrete time logistic population models (May 1974). At equilibrium we recorded the infection rate, density, sex ratio, and the time it took to reach that equilibrium.
To determine which life history combinations would result in a host benefiting or being harmed by a reproductive manipulator, we compared equilibrium density in infected populations to that in stable uninfected populations.
We overlaid these numerical analyses with the life history data for species where both egg size and the pelagic larval duration are known (Marshall et al. 2012;Álvarez-Noriega et al. 2020). We then performed separate logistic regressions (generalized linear model with a binomial family and log link) to test whether infectability was associated with host phylum (Annelida, Echinodermata, and Mollusca), developmental mode (planktotrophy and lecithotrophy), and their interaction. We excluded Cnidaria (n p 4) and Chordata (n p 1) from our logistic regressions because of their small sample size. We then performed a likelihood ratio test using a type II analysis of variance, with P values less than .05 considered to indicate statistical significance.
To test how developmental mode influenced the effect of infection on population density, we first calculated the log 2 fold change of equilibrium density from uninfected to infected and removed extreme outliers (absolute value log 2 fold change 125) or simulations that resulted in population cycles. We then performed a permutation test with 10,000 permutations (due to nonnormality of data) on the difference in mean log 2 fold change of equilibrium density between planktotrophic and lecithotrophic echinoderms. We only tested for developmental mode within Echinodermata because of much lower sample sizes in Annelida and Mollusca.
Infection by a microbial manipulator often comes at a fitness cost (Engelstädter and Hurst 2009a). To determine how this cost caused by a feminizing microbe influences infectability, we repeated simulations for a subset of egg  (4) and (5) were adopted from the Styan (1998) polyspermy model and were used to generate these graphs, with species-specific parameters from the literature (table 1). These reproductive parameters can be altered to generate additional species-specific curves using the web application at https://mck8dg.shinyapps.io/SI_sex_ratio_distortion_marine_inverts/.
sizes and pelagic larval durations while varying the feminization rate and cost. Specifically, we simulated different feminization rates (r) from 0.501 to 1.000 at increments of 0.001 and costs (c) ranging from 0.000 to 0.600 at increments of 0.001. The subset of egg sizes that we used was from 450 to 800 mm at increments of 50 mm, and the subset of pelagic larval durations was 1, 5, 9, 13, 18, and 30 days. These combinations of egg size and pelagic larval duration were usually "infectable" in a preliminary feminization rate-cost analysis using a broader range of egg sizes and pelagic larval durations from above ( fig. S3). Last, we tested which feminization rate and cost combinations resulted in populations with ≥99.9% infection.

Fertilization Dynamics in Free-Spawning Marine Invertebrates
Consistent with previous reports, the fertilization success and zygote production of free-spawning marine invertebrates is dependent on adult density (Pennington 1985;Levitan 1991) and egg size (Farley and Levitan 2001;Po-dolsky 2002). Fertilization success and zygote production are also influenced by sex ratio for both the planktotrophic Heliocidaris tuberculata and the lecithotrophic H. erythrogramma, but how they are influenced depends on the developmental strategy (figs.1, S4, S5). In H. tuberculata, fertilization success is sperm limited at low adult densities (!0.01 individuals/m 2 ) for all sex ratios until a sufficient density caused the fertilization success to reach a relative maximum. The adult density at which fertilization success reached a relative maximum depended on the sex ratio, such that the density is slightly lower in female-dominant populations (0.0177 individuals/m 2 , for a sex ratio of 0.1) and slightly higher in male-dominant populations (0.0181 individuals/m 2 , for a sex ratio of 0.9; figs. 1, S4, S5). The magnitude of this relative maximum in fertilization also depended on the sex ratio, whereby fertilization increased approximately linearly from 0% in a fully female population to 95% in a population that is 75% males and gradually plateaued toward a 100% fertilization success as the population became predominately, but not entirely, male ( fig. S5). Within the range of realistic densities (i.e., !100 individuals/m 2 ), polyspermy did not result in zero fertilization success until densities were sufficiently high (120 individuals/m 2 ) and the sex ratio was heavily male biased (i.e., ≥0.84; figs. 1, S4, S5). For most sex ratios, however, fertilization success declined slightly as density increased ( fig. 1). The interaction between density and sex ratio resulted in a greater zygote production as density increased for sex ratios below ∼0.84, even if maximum fertilization was not achieved (figs. 1, S4, S5). As the population transitioned toward female dominance, the maximum zygote production increased until 0.41, after which maximum zygote production decreased as a result of sperm limitation (figs. 1, S4, S5). The dynamics of fertilization success and zygote production for the lecithotrophic H. erythrogramma were quite different. While the density-dependent fertilization success was similar in heavily male-biased populations for both Heliocidaris species (i.e., at 0.9), this model predicted a shift toward high adult densities-instead of a reduction-in the fertilization profile as the population structure became female biased (figs. 1, S4, S5). Maximum fertilization success stayed at 100% until the population consisted solely of females (figs. 1, S5). In the near absence of a reduc-tion in fertilization success due to a shift in sex ratio, H. erythrogramma experienced an exponential increase in zygote production as the population became female dominant (figs. 1, S4, S5). Moreover, peak zygote production also shifted with the population structure and with adult density, resulting in a synergistic effect on zygote production for this developmental strategy (figs. 1, S4).

Spread of Reproductive Manipulators:
A Heliocidaris Case Study  feminization was the only mechanism that resulted in both the infection spreading and a female-biased sex ratio.
Male killing often evolves when the female offspring is provided a benefit from the sacrificed male offspring (e.g., by cannibalization; Hurst et al. 1996;Engelstädter and Hurst 2009b). This was not present in our integrated model because the dilute lifestyle of sibling larvae for free-spawning marine invertebrates makes cannibalism highly unlikely. Figure 4: Equilibrium infection rate (A) and log 2 fold change in density (B) by a microbe that induces feminization (left) or male killing plus enhanced growth (right) at different combinations of egg diameter (from 40 to 1,500 mm with 10-mm increments) and pelagic larval duration (from 0 to 80 days with 1-day increments). Numerical simulations to determine the equilibrium were performed using equations (6)-(11). For infection rate (A), light through dark blue represent 0% to 100% infected, and gray boxes indicate parameter combinations that resulted in numerical errors because of too high of a growth rate. For the log 2 fold change in equilibrium density (individuals/m 2 ; B), green through purple represent a net benefit to net cost to an infection. Moreover, gray boxes indicate parameter combinations that had less than 25% infection, had population cycles, or resulted in numerical errors because of too high of a growth rate. For this model, egg density was inversely proportional to egg diameter (Smith and Fretwell 1974). For the feminization numerical simulation (left), the cost of infection was 0 and the feminization rate was 0.7. For the male killing plus enhanced growth numerical simulations (right), the cost of infection was 21.25-representing a benefit-and the male killing rate was 0.7. Other parameters of the model were held to be the same as Heliocidaris erythrogramma (see table 1). Numerical simulations repeating this analysis but varying feminization and cost are shown in figure S3. The time it took to reach equilibrium, total population density at equilibrium, and the location of marine invertebrate species based on these two life history parameters can be explored in the web application at https://mck8dg.shinyapps.io/SI_sex_ratio_distortion_marine_inverts/. FC p fold change.
Therefore, we provided infected females with a fitness compensation by expediting growth as an alternative mechanism for endosymbiont spread (Carrier et al. 2021). Our structured population model for both Heliocidaris species found that reproductive manipulators spread when enhanced growth is relatively small and higher growth rates alone resulted in male-dominant populations, with the smaller size classes being more male biased than the largest size class (figs. 3, S7-S9). When combined with a modest degree of male killing (i.e., 0.7), the population structure for each size class became female dominant, with the magnitude of this sex ratio distortion increasing with adult size (figs. 3, S7-S9). Furthermore, enhanced growth and male killing interact to affect equilibrium sex ratios, where enhanced growth influenced smaller size classes and lead to a male-biased population while male killing became more influential in the larger size classes to sway the population toward females (S8, S9). This resulted in a nearly identical population structure to that of H. erythrogramma in Sydney Harbor but was considerably more female dominant than is observed for H. tuberculata (Carrier et al. 2021). These processes were expedited to equilibrium when combined with cytoplasmic incompatibility (figs. 3, S7, S9).

Numerical Analyses of Reproductive Manipulators and Life History Traits
Our simulations across a range of egg sizes and pelagic larval durations found that combinations of egg size and larval duration served as indicators of a potential infection, whereby species of free-spawning marine invertebrates were either fully infected or not infected, with little in between ( fig. 4). Of these two life history traits, the spread of a manipulating microbe was more limited by the pelagic larval duration than by egg size. Specifically, infection did not occur for larvae that spent more than ∼40 days in the water column but continued beyond egg diameters of 1,500 mm ( fig. 4A). The upper-bound limit of infectability was similar for a microbe that feminized the host as well as provided a fecundity benefit (i.e., enhanced growth) to the host plus male killing ( fig. 4A). Furthermore, combinations with an egg size smaller than 30 mm and a pelagic duration less than 20 days also were not infectable by microbes that feminize the host. This lower limit was not present for a microbe that elicits male killing plus enhances growth, and thus this microbial strategy could infect a wider range of life history combinations than a feminizing microbe. Figure 5: The combination of fitness cost and feminization rate that are predicted to result in populations with a ≥99.9% infection rate. These numerical simulations were run at feminization rates from 0.501 to 1.000 (with increments of 0.001) and costs from 0.000 to 0.600 (with increments of 0.001) for egg sizes from 450 to 800 mm (with increments of 50 mm) and pelagic larval durations of 1, 5, 9, 13, 18, and 30 days. Additional parameters of this simulation were held to those of Heliocidaris erythrogramma (see table 1; fig. 4).
For the numerical analysis of a feminizing microbe, we found that 68.1% of these species fell within this "infection zone" ( fig. 4A; web application at https://mck8dg.shiny apps.io/SI_sex_ratio_distortion_marine_inverts/). Phylum, but not developmental mode, was a significant predictor of infectability (phylum: x 2 2 p 8:184, P p :017; developmental mode: x 2 1 p 1:754, P p :185). Specifically, 36.8% of annelids (n p 19), 70.0% of echinoderms (n p 140), and 72.6% of molluscs (n p 73) as well as 71.4% of planktotrophs (n p 154) and 61.5% of lecithotrophs (n p 78) were within this zone. Furthermore, there was a significant interaction between phylum and developmental mode (x 2 2 p 27:6167, P ! :001), such that lecithotrophic annelids, planktotrophic echinoderms, and lecithotrophic molluscs were less likely to fall within this infection zone than the alternative developmental strategy for each phylum ( fig. S10A). Infection tended to have a negative influence on the population density of planktotrophs and a positive influence on that of lecithotrophs ( fig. 4B). Moreover, within echinoderms, infection had a more positive influence on population density for lecithotrophs than for planktotrophs (mean difference in log 2 fold change p 1:223; P p :002; fig. S11A).
For the numerical analysis of a microbe that elicits male killing plus enhances growth, we found that 83.6% of these species fell within this infection zone ( fig. 4A). Phylum was the only statistically significant predictor of infectability (phylum: x 2 2 p 14:605, P ! :001; developmental mode: x 2 1 p 2:997, P p :083; interaction: x 2 2 p 0:562, P p :756). Specifically, 100% of annelids (n p 19), 77.9% of echinoderms (n p 140), and 90.4% of molluscs (n p 73) as well as 81.8% of planktotrophs (n p 154) and 87.2% of lecithotrophs (n p 78) were within this zone ( fig. S10B). Infections tended to have a negative influence on the population density of planktotrophs and a positive influence on that of lecithotrophs ( fig. 4B), and this was particularly the case for echinoderms (mean difference in log 2 fold change p 2:308, P ! :001; fig. S11B).
The shape of this infection zone was constant across all combinations of feminization rate and infection cost ( fig. S3). After a certain cost, however, the potential for an infection across all life history combinations fell to zero, and thus this infection zone disappeared completely ( fig. S3). There was a clear boundary for the maximum cost that a reproductive manipulator can incur on its host given a fixed feminization rate, resulting in a hyperbolic relationship between these factors ( fig. 5). In other words, there is a fine balance between the feminization rate, fitness cost to the host, and whether an infection may occur.

Discussion
We integrated two types of models-the first on the fertilization dynamics of free spawning and the second on the consequences of cytoplasmic incompatibility, feminization, and male killing-to understand whether microbes that manipulate host reproduction spread within populations of free-spawning marine invertebrates and distort their sex ratio. Our model and the Heliocidaris case study support three primary findings. First, sex ratio distortion leads to suboptimal fertilization and zygote production in planktotrophs but can enhance these processes in lecithotrophs. Second, feminization and a combination of male killing plus enhanced growth are effective at spreading reproductive manipulators while also inducing a femalebiased sex ratio. Third, the majority of free-spawning marine invertebrates could be infected across a range of fitness costs and manipulation rates, with infections favoring and benefitting species with larger eggs and shorter pelagic durations. These results suggest that the life history evolution of free-spawning marine invertebrates is an underlying factor in whether microbes that manipulate reproduction may infect and spread these partnerships.

Life History Evolution and Microbes That Manipulate Reproduction
The reproductive success of free-spawning marine invertebrates is often limited by sperm availability, which may be partially offset by increasing egg size, increasing spawning density, and spawning in synchrony (Pennington 1985;Levitan et al. 1992;Oliver and Babcock 1992;Levitan and Petersen 1995;Yund 2000;Levitan 2006). An additional mechanism that may limit reproductive success is the distortion of host sex ratio toward female dominance. This population structure is common for some groups of freespawning marine invertebrates (Carrier et al. 2021), implying that it may provide some benefits. Our fertilization model is in agreement with this statement. The planktotroph in our case study experienced a lower maximum fertilization success as the sex ratio became more female biased, while maximum zygote production increased as the sex ratio became slightly female biased but decreased thereafter. Fertilization curves of the lecithotroph, on the other hand, were constant across sex ratios and shifted slightly toward higher adult densities. This, in turn, resulted in a maximum fertilization success of ∼100% for most sex ratios and an exponential increase in the maximum zygote production as sex ratio became more female biased. An increase in zygote production would benefit the host as well as the reproductive manipulator (Engelstädter and Hurst 2009a), but this fitness advantage is inconsistent across the life history strategy of the host. If a reproductive manipulator that distorts host sex ratio were to infect a planktotrophic species, then the high fecundity of this life history strategy alone should promote the rapid spread of the microbe. Transitioning from a male-dominant to a female-dominant population structure would provide an initial fitness benefit, but if the selective pressures on host reproduction are too high, then the manipulator would succumb to its own influence. For the planktotrophic life history strategy, host sex ratio appears to parallel the mutualism-parasite continuum, such that a mild reproductive manipulator could be beneficial while a highly effective one would be parasitic, with limitations to reproductive success being at fertilization. This continuum is widely observed in terrestrial arthropods (Hatcher et al. 1999;Charlat et al. 2007;Drew et al. 2021) and appears possible in the sea.
Infecting planktotrophs should also be an ineffective transmission route for microbial manipulators. Natural mortality of marine invertebrate larvae is generally dependent on the time spent in the water column. Gamete wastage (mortality) often exceeds 90%-95% because planktotrophic larvae spend weeks in the water column, and thus this would compromise the fitness of a reproductive manipulator at a rate equal to host mortality (Thorson 1950;Young and Chia 1987;Rumrill 1990;Morgan 1995;Engelstädter and Hurst 2009a). For example, if a planktotrophic larval host spends a mere 21 days in the water column, then only ∼0.1% of the reproductive manipulator cells would be transmitted to the next generation (Strathmann 1985;Rumrill 1990). This limitation is based on the pelagic duration of the larva and is reflected in our numerical analyses that predicts the infectable life history combinations (Shanks 2009). Specifically, while ∼68% (feminization) to ∼84% (enhanced growth plus male killing) of species could be infectable, those with shorter pelagic larval durations and larger maternal investment (egg size) were favored and should benefit from an infection. Maternal investment also happens to be a strong proxy for the developmental life history strategy of free-spawning marine invertebrates (Thorson 1950;Vance 1973;Emlet et al. 1987).
Lecithotrophic species produce fewer, but substantially larger, eggs than planktotrophs, and the resulting larvae spend only days in the water column (Thorson 1950;Emlet et al. 1987;Shanks 2009). If a reproductive manipulator that distorts host sex ratio were to infect a lecithotroph, then the symbiont could benefit from the eggs physical size by substantially increasing its initial titer and minimize a fitness loss because this developmental strategy spends a brief period in the water column. Distorting the sex ratio of this life history strategy could favor an exponential increase in zygote production that would be advantageous for both the host and the reproductive manipulator. Opposing that of planktotrophs, associations between lecithotrophs and reproductive manipulators could be perceived as a mutualism that would likely be more stable and increasingly beneficial as the host sex ratio is distorted toward a female-dominant population structure. Provided the differences in infection limitations and benefits to these two life history strategies for a reproductive manipulator, it may be hypothesized that these microbes could serve as selective agents to induce, mediate, and/or reinforce transitions between the developmental strategies of free-spawning marine invertebrates. Life history transitions between these developmental strategies have occurred in several major animal lineages, with rapid evolutionary shifts from planktotrophy to lecithotrophy being well documented. It is thought that an increase in the energetic content (size) of the egg relaxes the selective pressure maintaining the feeding structures and that development to metamorphosis is accelerated once these are lost (Strathmann 1978(Strathmann , 1985Wray and Raff 1991a;Raff 1992;Jaeckle 1995;Wray 1996;McEdward and Morgan 2001;Sewell and Manahan 2001;Raff and Byrne 2006;Moran et al. 2013). An increase in egg size and an accelerated time to metamorphosis would both favor the likelihood that a reproductive manipulator is transmitted to the subsequent generation. If such microbes also distort the host's sex ratio, then there would be multiple incentives to induce, mediate, and/or reinforce transitions in the developmental life history of free-spawning marine invertebrates.

Spread and Influence of Marine Reproductive Manipulators
Microbes that manipulate reproduction in terrestrial systems spread by cytoplasmic incompatibility, feminization, and male killing, but only feminization and male killing may distort host sex ratio (Engelstädter and Hurst 2009a). Our model suggests that the reproduction of free-spawning marine invertebrates may be influenced by each of these; that host sex ratio may be distorted by feminization and male killing, with the latter occurring only when the host receives some enhanced growth benefit (i.e., fecundity benefit); and that cytoplasmic incompatibility enhances the spread of these microbes. In all cases, reproductive manipulators also spread faster in planktotrophs than in lecithotrophs. Therefore, we suspect some convergence between microbe-mediated reproductive manipulation on land and in the sea. The prevailing question is, how may these occur in marine systems? Feminization involves the conversion of males to females by microbes that interact with the host's sex determination system. This is similar to sequential hermaphroditism, except that feminization is mediated by the microbial manipulator. It is common for free-spawning marine invertebrates to be sequential hermaphrodites (Sewell 1994;Collin 2006). The patterns and timing of sex change in echinoderms, for example, varies considerably and may occur in either direction, with protandry being more prevalent (see table 1 in Sewell 1994). One expectation of marine microbial manipulators is that they would interfere with the sex determination system of a sequential hermaphrodite to enable an earlier transition from male to female (Warner 1975;Munday et al. 2006). This type of reproductive manipulation has the potential to be widespread in marine invertebrates and likely depends on the capacity of these microbes to interfere with the sex determination systems for freespawning marine invertebrates.
Male killing, on the other hand, involves microbes that cause an elevated mortality during male development and can primarily help the spread of infection when females can benefit directly from the sacrificed siblings (Hurst and Jiggins 2000). This type of reproductive manipulation alone appears unlikely for free-spawning marine invertebrates for two reasons. First, embryonic mortality naturally exceeds 90%-95% and mortality due to male killing would be an additional 5%-90% (Young and Chia 1987;Rumrill 1990;Morgan 1995;Pechenik 1999;Hurst and Jiggins 2000;Engelstädter and Hurst 2009a). The combination of these mortality sources would be overwhelming and, unlike terrestrial systems, would limit microbial fitness in addition to the host. Second, life in the plankton for the offspring of free-spawning marine invertebrates is an independent endeavor, and this would not allow for a fitness compensation. We suspect that male killing may evolve in the benthic (adult) life stage, which would be a type of male killing that is not observed in terrestrial systems (Engelstädter and Hurst 2009a). Alternatively, a potential ecological "hot spot" where male killing may evolve is when offspring are brooded, thereby being confined and where sibling cannibalism is observed (e.g., echinoderms [Byrne 1996] and gastropods [Strathmann and Strathmann 2006).
An alternative mechanism to compensate for the fitness cost related to male killing is by the microbe supplementing host nutrition. Under this circumstance in our model, male killing appeared to be stable for the lecithotroph and able to elicit a shift in sex ratio toward female dominance, and such a microbe was able to infect a larger range of life history combinations than a feminizing microbe. Moreover, when modeled with enhanced growth in a sizestructured population, this combination resulted in a population structure highly similar to that observed in the field for the lecithotroph (Carrier et al. 2021). If only provided a growth benefit, then the population became male dominant, which was accelerated with cytoplasmic incompatibility. Therefore, male killing coupled with some growth (i.e., nutritional) benefit that should increase host fitness and promote the spread of the symbiont. The free-floating life of the developmental stages of marine invertebrates may thus be influenced by such reproductive manipulators.
An additional mechanism of microbial manipulation is cytoplasmic incompatibility (Engelstädter and Telschow 2009). This common type of reproductive manipulation enhances symbiont spread by eliciting embryonic mortality when infected males mate with uninfected females but avoids this when both sexes carry the compatible manipulator. Furthermore, cytoplasmic incompatibility can be either unidirectional or bidirectional, with lethality being avoided in the latter when both sexes harbor the same strain of the symbiont. This results in reproductive isolation between hosts with different infection strains and can sufficiently restrict gene flow and induce speciation events (Engelstädter and Hurst 2009a;Shropshire et al. 2020;Kaur et al. 2021). Our model suggests that such a mechanism would enable the spread of reproductive manipulators in the sea.
Although not currently interpreted as symbiontinduced cytoplasmic incompatibly, there are a diverse array of gamete incompatibilities observed in the sea. For example, there are strong male-by-female interactions and maternal effects on fertilization efficiency and zygote success in ascidians (David et al. 2016), echinoderms (Evans and Marshall 2005), and mollusks (Rawson et al. 2003). The hypothesized causal mechanisms for these examples stem from a host point of view, but these processes may also be caused by a microbial manipulator. Our model demonstrates that cytoplasmic incompatibility is a viable mechanism for a symbiont to spread in the sea, and we propose this alternative view should also be considered. We believe this is especially true when a unidirectional or bidirectional incompatibility is suspected (Shropshire et al. 2020). In the case of incompatibilities, the microbial community should then be profiled, and if a Rickettsiales is a dominant member of the gonad and egg-associated microbiome, then a reproductive manipulator may be present.

Was the Life History Evolution of Heliocidaris a Reproductive Manipulation?
If the life history evolution of Heliocidaris has been manipulated, then our model and numerical simulations would suggest that H. erythrogramma (i) has a microbe that could feminize or elicit male killing if a growth benefit is provided, (ii) has an enhanced reproductive output, and (iii) has had subsequent speciation events that led to species with larger eggs and/or shorter pelagic durations. We believe that each of these expectations have been met. First, the Echinorickettsia raffii genome has an ankyrin protein with a homologous domain to a male killing gene and a complete shikimate pathway to biosynthesize essential amino acids that may supplement host nutrition (Carrier et al. 2021). Second, adult density in female-dominant populations in Western Australia is near the modeled maximum for zygote production and this far exceeds that of the balanced (1∶1) Eastern Australia population (Dix 1977;Evans and Marshall 2005;Binks et al. 2012;Carrier et al. 2021).
Third, two subsequent speciation events occurred as H. erythrogramma spread. One of these species (H. bajulus) has large eggs that are brooded and undergo a lecithotrophic development in a thick gelatinous coat adhered to the parent (Dartnall 1972;McMillan et al. 1992;Zigler et al. 2003;Hart et al. 2011Hart et al. , 2012. The natural history and life history evolution of Heliocidaris is therefore consistent with expectations of a reproductive manipulation in the sea. Male killing would seem to be the likely type of reproductive manipulation that influenced the life history evolution of Heliocidaris. Our model, however, showed that cytoplasmic incompatibility could be viable in combination with this mechanism, and we suspect it occurs as well. First, there is unidirectional hybrid incompatibility, where H. tuberculata sperm bind to but cannot fertilize H. erythrogramma eggs unless the jelly coat is removed, allowing for a fertilization efficiency (∼85%) equal to the reciprocal cross Zigler et al. 2003). Second, while both sexes contribute to the high degree of fertilization success within populations, there is a substantial maternal effect in this system (Evans and Marshall 2005;Evans et al. 2007). Third, since spreading to Western Australia, H. erythrogramma has subspeciated to H. e. erythrogramma and H. e. armigera, and these can, but rarely do, hybridize because of asynchrony in the spawning seasons (McMillan et al. 1992;Binks 2011;Binks et al. 2012). Therefore, we hypothesize that the reproductive and population biology of H. erythrogramma is also influenced by cytoplasmic incompatibility and that marine reproductive manipulators should have cytoplasmic incompatibility factor genes (LePage et al. 2017; Kaur et al. 2021), but these have yet to be identified.

Model Limitations and a Way Forward
Modeling the fertilization and life history dynamics of freespawning marine invertebrates as well as reproductive manipulation are simplified here, and this may have influenced our results. There are three assumptions that we feel will need to be relaxed in future modeling efforts to more broadly understand how microbes that manipulate reproduction influence free-spawning marine invertebrates. First, our model was deterministic and, by definition, lacked stochasticity. Without this, the proportion of infected individuals can remain extremely small and eventually rise to fixation. Second, our model assumes a closed population, where an infection will still spread even if a population is not optimally productive if there is an advantage to the infected females. In an open system, a more productive neighboring population that was uninfected could swamp the population with a nascent infection and prevent the infection from spreading. Such a dynamic would likely be more important for planktotrophs than lecithotrophs be-cause sex ratio distortion easily moves these species off their optimum and migration rates and gene flow are likely high (McMillan et al. 1992). Third, our model does not address the possibility for variation among individuals and, thus, for the evolution of reproductive traits that influence symbiont spread. It is likely that changes in sex ratio or, more specifically, the egg-to-sperm ratio influence the selection on reproductive traits (e.g., egg size and sperm velocity). Changes in host reproductive traits may enhance or prevent the spread of a reproductive manipulator. Explicitly modeling the potential coevolution between host and symbiont is needed to test the hypothesis that symbionts may mediate transitions in the developmental life history of free-spawning marine invertebrates.

Conclusions
Our model supports the general premise that microbes may manipulate the reproduction of free-spawning marine invertebrates and that the types of manipulation overlap considerably with terrestrial systems (Hatcher et al. 1999;Charlat et al. 2007). Much of the life history landscape of free-spawning marine invertebrates may be manipulated by these microbes. The reality of a manipulation likely differs considerably for a host with different life history strategies, such that stable and beneficial shifts in population structure may be induced for lecithotrophs while there are limitations to reproductive success for planktotrophs. This begs the question of whether reproductive manipulators are omnipresent in the sea as they are on land, whether there are geographical and life history hot spots for these microbes, and whether the aquatic environment features unique ways for microbes to manipulate host reproduction.