A phylogenetic comparative analysis on the evolution of sequential hermaphroditism in seabreams (Teleostei: Sparidae)

The Sparids are an ideal group of fishes in which to study the evolution of sexual systems since they exhibit a great sexual diversity, from gonochorism (separate sexes) to protandrous (male-first) and protogynous (female-first) sequential hermaphroditism (sex change). According to the size-advantage model (SAM), selection should favour sex change when the second sex achieves greater reproductive success at a larger body size than the first sex. Using phylogenetic comparative methods and a sample of 68 sparid species, we show that protogyny and protandry evolve from gonochorism but evolutionary transitions between these two forms of sequential hermaphroditism are unlikely to happen. Using male gonadosomatic index (GSI) as a measure of investment in gametes and proxy for sperm competition, we find that, while gonochoristic and protogynous species support the predictions of SAM, protandrous species do not, as they exhibit higher GSI values than expected even after considering mating systems and spawning modes. We suggest that small males of protandrous species have to invest disproportionally more in sperm production than predicted not only when spawning in aggregations, with high levels of sperm competition, but also when spawning in pairs due to the need to fertilize highly fecund females, much larger than themselves. We propose that this compensatory mechanism, together with Bateman’s principles in sequential hermaphrodites, should be formally incorporated in the SAM.

mode and allometry. Although GSI, the percentage of gonad tissues on body mass, may provide an estimate of investment in the gonads irrespective of size, it is possible that at least some of the observed differences in male GSI across species are in part determined by differences in male size, if investment in the gonads is easier at larger or smaller sizes. In addition, the reported high GSI values in this study were derived from only eight protandrous species 14 , three of which were in fact those of individuals reproducing as female rather than male, or, unlike for other species in the dataset, do not correspond to the maximum recorded male GSI value for the species and so are not comparable. Finally, recent phylogenetic trees, which are essential for all comparative studies are better resolved and more comprehensive than those employed by earlier studies (e.g. 18 ), while modern phylogenetic comparative methods allow overcoming the limitations of approaches used in previous studies, such as parsimony reconstruction and phylogenetic independent contrasts, that can lead to incorrect conclusions 60,61 .
Here we have compiled the largest dataset to date of sexual systems, spawning modes and body size in the family Sparidae. Using modern phylogenetic comparative approaches, we have investigated the evolutionary history of sexual systems and tested the predictions of the SAM, considering spawning mode and mating systems, that protogynous and protandrous species should exhibit lower GSI values due to expected lower sperm competition (Table 1) than gonochoristic species, while accounting for body size. Our study therefore combines sexual systems, mating systems, sperm competition and the principles of the SAM in several sparid species, while accounting for the potential confounding effects of allometry and phylogeny.

Methods
Data collection. We used FishBase (www.fishbase.org 62 ) to gather information on the sexual systems of Sparidae, ranging from gonochorism (G), protandrous hermaphroditism (PA) and protogynous hermaphroditism (PG). We verified, and if necessary corrected, the sexual system reported in FishBase for each individual species used in this study against the primary literature 14,25,54,63,64 . We carefully revised previous assignments of sexual systems in four species, which were mainly based on the gonadal morphology of individuals collected at single or different ages 65 . Sometimes this approach cannot distinguish functional (active) hermaphrodites from non-functional hermaphrodites (i.e., individuals that despite presenting both male and female tissues reproduce only as one sex throughout their life). The assignment of the correct sex is further complicated in non-reproducing juveniles, which can present a bisexual gonadal stage. Altogether, these peculiarities can make diagnosis of sexual system extremely challenging 14,54,64 . We resolved any discrepancies from previous studies using newly published data in which care was taken not to incur in the above problems (Table S1).
Out of a total of 148 recognized sparid species, we could assign the sexual system in 68 species (Table S1). For these species, we extracted data from FishBase on male maximum body weight (in g; n = 37 spp.), male maximum total body length (in cm; n = 47 spp.) and male total body length at maturity (in cm; n = 36 spp.), to account for possible allometric effects on GSI. Finally, we extracted data from the primary literature on male GSI (n = 49) and spawning mode (pair or group spawning, based on the presence of one or more males, respectively; n = 10 spp.). When several GSI values were reported for a given species (e.g., monthly means along the year), we consistently used the highest value.
Ancestral state reconstruction. We used two molecular phylogenetic trees with time-calibrated branch lengths, an essential step for robust analyses in a phylogenetic comparative framework 66 . Specifically, we used a phylogeny of Actinopterygians 67 based on a 27-genes (6 mitochondrial and 21 nuclear genes) and a phylogeny for the family Sparidae 68 , based on three mitochondrial and two nuclear genes. These trees included 58 and 55 species, respectively, out of the 68 species with sexual system information in our dataset.  Table 1. Simplified characterization of mating systems (monogamous, polygynous, promiscuous; see text for full definitions), spawning mode (pair spawning or group spawning), sperm competition (classified as low or high) and gonadosomatic index, GSI (classified as low or high) for each sexual system (G = Gonochorism; PA = Protandry; PG = Protogyny). In bold: most common type. The prediction column summarizes the overall expectation on sperm competition (and thus GSI) for the sexual system, based on the most common mating and spawning mode, according to the SAM 23,24 . *Protogynous group spawners should be favoured as large males can produce more sperm 28 . **Promiscuous large groups are not common, but still present in protandrous species (see text for examples).
The ancestral state reconstruction infers the evolutionary history of a trait along a phylogeny given the character states of species in the tree and provides estimates of the probable character state of each node in the phylogeny. This approach is based on a Markov model of evolution for discrete traits 69 . We reconstructed the ancestral character states of sexual system using maximum likelihood (ML), setting all transition rates between G, PA and PG free to vary, i.e. they are not constrained to be of equal magnitude. We ran these analyses on both phylogenetic trees using the R package ape 70 . However, we can only report the results of the ancestral state reconstruction using the phylogenetic tree of Rabosky et al. 67 since the analysis using the Santini et al. 68 tree did not converge to a maximum likelihood solution.
Testing SAM predictions. We used phylogenetic generalized least square models (PGLS 61,71,72 ) to test the predictions of the SAM using the R package caper 73 and ML estimation, with both phylogenetic trees 67,68 . By incorporating the phylogeny, PGLS models can quantify the strength of phylogenetic signal in the data through the parameter lambda (λ), which can vary between zero (no phylogenetic signal) and one (high phylogenetic signal, whereby the species exhibit phenotypic similarity in direct proportion to their common evolutionary time) 61,72 . GSI was entered as the dependent variable, while sexual system was the independent discrete variable with three possible states (0 = G, 1 = PA, 2 = PG), and body size (as maximum length, weight or length at maturity) the independent continuous covariate. Possible allometric effects on the GSI were thus accounted for using either body length or weight as additional independent variable. Continuous variables were log 10 -transformed to meet assumptions of normality, with the exception of GSI values. Results for GSI were qualitatively similar whether this variable was log 10 -transformed, transformed with logit function or left untransformed, thus we report the results of GSI in percentage, not transformed. All model residuals were normally distributed in all analyses.
We also tested the SAM prediction within the genus Diplodus. This was the only genus that provided limited but sufficient data to consider at least two different sexual systems (G vs. PA) for statistical analysis of the relationship between sexual system and male GSI values in very closely related species. Importantly, Diplodus species have a narrow range of body sizes and thus there is less variability and potential confounding effects of allometry. However, the analysis could not be carried out in a phylogenetic context in this genus because too many species were missing from both phylogenies. Therefore, we used Student's t-test for independent samples to test differences in total male body length and GSI between the two sexual systems, and general linear model to test for the relationship between maximum male length at reproduction and GSI. In all analyses, performed in R 74 , differences were considered statistically significant when p < 0.05.

Results
Our dataset on the sexual system of 68 sparid species across 28 genera (out of 37) includes 27 gonochoristic, 22 protogynous and 19 protandrous species (Table 2, Tables S1 and S2), with the three sexual systems roughly present in similar proportions (range~28-40%; Table 2). This represents a substantial increase in the number of species previously investigated 14 .
Reconstruction of the ancestral character state in a phylogenetic context showed that gonochorism was only slightly more likely to be the ancestral sexual system in the Sparidae family (likelihood at the root 37.4%) than protandry (33.5%) and protogyny (29.1%; Fig. 1). While both forms of sequential hermaphroditism, especially protogyny, evolved rapidly to gonochorism (PA to G: 0.035 ± 0.013; PG to G: 0.064 ± 0.033), gonochorism reverted as quickly back to protogyny (G to PG: 0.054 ± 0.039) and much less so to protandry (G to PA: 0.010 ± 0.024). Finally, the transitions between the two forms of sequential hermaphroditism were both estimated to be zero (PA to PG: 0.000 ± 0.012; PG to PA: 0.000 ± 0.027; Fig. 2).
We did not find any significant difference in total male body length (Fig. 3a), male length at maturity or maximum male body weight between sparid species with different sexual systems (Table 3). For the 46 sparid species in the tree where male GSI values were available (Table S3), GSI values of protandrous and protogynous species were higher and lower respectively than that of gonochoristic species, with statistically significant differences between protogyny and the other two sexual systems ( Fig. 3b; Table 3). Results were qualitatively similar regardless the phylogenetic tree used (Table S4) and were not influenced by allometric effects, tested using either length or weight as a covariate in the model (Table S5). Unfortunately, limited data were available on spawning mode, for both species that spawn in groups (G, n = 2; PA, n = 4; PG, n = 0) or in pairs (G: n = 1; PA: n = 2; PG: n = 1). These low numbers did not allow testing predictions for mating systems formally. However, the data appear to suggest that protandrous sparids have higher GSI values than gonochoristic species regardless of whether they spawn in groups or pairs (Fig. 4).  www.nature.com/scientificreports www.nature.com/scientificreports/ We then repeated the analyses within the genus Diplodus, which contains only gonochoristic and protandrous species with a narrower range of lengths (~25-60 cm) when compared to that of sparids as a whole (~20-200 cm). We found no significant differences in weight (t 2.98 = −0.54; p = 0.62) between gonochoristic and protandric species. However, the latter had a significantly shorter length than the former (t 4.61 = 2.71; p = 0.04; Fig. 5a). Importantly, despite being smaller in size, protandric Diplodus species had a significantly higher GSI  www.nature.com/scientificreports www.nature.com/scientificreports/ (t 4.82 = −3.19; p = 0.02) than the gonochoristic species of the same genus (Fig. 5b), confirming the results found across all sparid species. Moreover, GSI was unrelated to male total body length (F = 11.65; df=2,6; R 2 = 0.79; p = 0.246), suggesting that GSI differences between sexual systems were not determined by allometric effects.

Discussion
With a larger dataset of sexual systems in the family Sparidae than previously used, this study reveals that protandry and protogyny can evolve from gonochorism, although the ancestral state remains still uncertain in this family. Importantly, we find that transitions between the two forms of sequential hermaphroditism are unlikely, if ever, to occur. We find strong support for the SAM predictions, incorporating mating system and sperm competition, that protogynous species should exhibit lower levels of sperm competition relative to gonochoristic species (Table 1), as quantified by their low GSI values, consistent with their mating systems that allow large males to monopolize access to fertile females. However, we find no evidence in support of similar predictions in protandrous sparids, i.e., low levels of sperm competition because they are expected to spawn in pairs or small groups (Table 1). Unexpectedly, protandrous species have similar GSI values to those of gonochoristic species and higher GSI than protogynous species, regardless of their mating system. Below we propose how a compensatory mechanism, together with mating system and spawning mode, may explain this unexpected finding.
Using twice as many species relative to an earlier study, recent molecular, time-calibrated phylogenies and modern phylogenetic comparative approaches, our study shows that gonochorism is only marginally more likely to be the ancestral state in this family. We find that gonochorism can evolve into both protogyny and protandry. . Asterisks indicate statistically significant differences with the following equivalence: *p < 0.05.   However, sequential hermaphroditism is an evolutionary unstable state as it reverts quickly back to gonochorism, suggesting that both types of sequential hermaphroditism are costlier to sustain than gonochorism. These results may explain why, despite hermaphroditism being anatomically and physiologically possible in fish in contrast to other vertebrates 75 , gonochorism predominates among fishes 64 . Female has been often considered the 'default sex' in fishes because, even in some gonochoristic species, all individuals, exhibit ovarian differentiation at the early stages of development, a process that is subsequently halted in the individuals that become males 63,76 . Thus, it is perhaps not surprising that our analysis reveals that the evolutionary transition rate from gonochorism to protandry is very low and that transitions from protogyny to protandry and vice versa are unlikely to occur; once canalized towards initial development as a female, it may be too costly to switch the developmental pathway to male-first sex changer (and vice versa).  www.nature.com/scientificreports www.nature.com/scientificreports/ We tested whether the sparids conform to SAM predictions when incorporating sperm competition, mating system and spawning mode as previously done in protogynous and gonochoristic epinephelids 35 . Specifically, we tested whether gonochoristic species, which often spawn in large groups or aggregations and are typically characterized by intense sperm competition, have a higher GSI than species with either types of sequential hermaphroditism. Sperm competition is indeed less intense in haremic species (often protogynous), where the presence of few dominant large males drastically reduce the interaction between sperm of different males 24,35,40 . It is also generally accepted that protandrous hermaphrodites normally reproduce in small, random mating groups (no size-assortative) or in strictly monogamous pairs 23,24,77 as, for example, in most clownfish (family Pomacentridae) of the genus Amphiprion, such as A. melanopus and A. percula [78][79][80] . In both mating systems, this would imply a low degree of sperm competition and thus low values of GSI, as predicted by the SAM [22][23][24] . As expected, we found that the average male GSI in gonochoristic sparids was >3%, in agreement with previous reports 14,35,40 , and significantly higher than the mean GSI of protogynous sparids (≤2%). These results, therefore, support the SAM prediction, incorporating sperm competition, in protogynous species (Table 1). However, our study also reveals that protandrous sparids have, on average, the highest male average GSI value (3.6%), regardless of their mating system. Thus, even when mating in pairs or small groups, protandrous males invest heavily in the gonads, indeed even more than gonochoristic species that mate in large aggregations with intense sperm competition. Further, we demonstrate that these results are not determined by differences in body size, as we find no allometric effects on GSI values, neither across all species, nor within one genus (Diplodus) containing closely related species of similar size with different mating systems. Altogether, our study unambiguously demonstrates that, while gonochoristic and protogynous sparids conform to a broader version of the SAM including sperm competition, spawning mode and mating system, protandrous sparids do not. We suggest that these results may be explained by both spawning mode and a compensatory mechanism determined by high sexual size dimorphism (SSD) in protandrous species. In protandrous species, sexual size dimorphism has been confirmed in several species such as Diplodus annularis and Lithognathus mormyrus 81,82 . However, when we checked male GSI in relation to SSD in the few protandrous fishes for which there is information we found a weak (r 2 = 0.3) positive relationship that nevertheless was not significant probably due to low sample size (data not shown). Therefore, the high GSI in protandrous species is not, in principle, explained by the magnitude of SSD. However, more data to would be needed to confirm this observation.

Variables
Some protandrous sparid species like Acanthopagrus berda, Sarpa salpa 83 , Diplodus capensis 84 and D. annularis 81 , spawn in aggregations 59 , similarly to many gonochoristic species 14 and should therefore experience intense sperm competition, leading to a high GSI value. Indeed, we often tend to oversimplify the complexity of sequential hermaphroditism: not all protogynous species are haremic, as not all protandrous species mate in pairs. A recent study 28,85 has revealed a broader variation in effective population size in protogynous species that differs from the more limited expectations obtained when all protogynous species were considered haremic by default. Instead, some protogynous species were found to be group spawners, altering the expectations that all protogynous species should have low effective population size due to their supposed haremic system. Similarly, not all protandrous species are monogamous or mate in pairs and small groups. Furthermore, many of the protandrous sparids that engage mostly or exclusively in pair mating, including Sparus aurata 86 and Rhabdosargus sarba 58 , exhibit a surprisingly high GSI of 4.4% which, albeit lower than values of group spawning sparids (5.5%), is still higher than most gonochoristic species (3%). Therefore, while spawning in aggregations can explain the high GSI of many gonochoristic and protogynous species, mating system and spawning mode alone cannot explain the high GSI consistently found in protandrous sparids. This is not unique to the sparids. For example, the majority of damselfishes (family Pomacentridae) reproduce in pairs and, as predicted by SAM 14 , exhibit lower GSI values (max. GSI < 1% 87 ); however, other protandrous damselfishes, such as the yellowtail clownfish Amphiprion clarkii, reproduce in pairs 88 but exhibit unexpectedly high GSI (4.14%) 89 . This corroborates the suggestion that protandrous species can exhibit high GSI values, regardless of mating system.
Here we propose that it is precisely the nature of protandry what explains the high GSI in protandrous sparids with pair mating. Specifically, protandrous males (first sex) are always smaller than females. Thus, given that fecundity increases with size in females 90 , small protandrous males need to produce large amounts of sperm to effectively fertilize highly fecund females much larger than themselves, even when mating in pairs and under conditions of low levels of sperm competition. To do so successfully, they need to invest disproportionally in the gonads. For small-sized protandrous males, there might be a body size threshold below which the GSI has to increase to ensure enough sperm production to fertilize the eggs produced by the larger females. Consistently, evidence of sperm production adjustment in relation to the amount of eggs to be fertilized has been reported in a coral reef fish, Thalassoma bifasciatum 91 . Further, in many species with alternative male mating strategies, smaller "sneaker" males invest more in the gonads than larger territorial males 38 due to the competitive environment as well as to compensate for their reduced size. Taken together, this evidence supports our suggestion that protandrous males have higher GSI than expected because they need to ensure the fertilization of eggs produced by much larger females than themselves.
The potential limitation in the fertilization capacity of small males is in agreement with the current debate questioning Bateman's principles [92][93][94] . Briefly, Bateman principles state that, due to the smaller cost of producing sperm when compared to eggs: (i) male reproductive success (RS) increases with mate number whereas female RS does not; (ii) males have greater variance in RS than females, and (iii) the sex with the greater variance in RS undergoes stronger sexual selection 92,94 . Bateman's principles have been used to test the strength of sexual selection between the sexes in gonochoristic 95 and simultaneous hermaphroditic species, although whether they really apply to the latter is currently highly debated 93 . Surprisingly, although the presence of sexual selection in sequential hermaphrodites has been recognized with the assumption that it is stronger on the most abundant sex 96 , Bateman's principles have not been formally tested. We argue that it is the combination of the existence of male-skewed sex ratios and male-male competition, on the one hand, and the problem for small males of Scientific RepoRtS | (2020) 10:3606 | https://doi.org/10.1038/s41598-020-60376-w www.nature.com/scientificreports www.nature.com/scientificreports/ producing enough sperm to fertilize large females, on the other hand, that together explain the high GSI in males of protandrous sparids, depending on the social/mating system. In fact, there is evidence that sperm depletion is a problem for many males across several taxa (reviewed in 94 ). For example, in the simultaneous hermaphrodite polychaete Ophryotrocha diadema, small protandrous males can have difficulties fertilizing a full clutch of eggs 97 . Sperm depletion has been also documented in fish 98,99 . Moreover, in external spawners as most fishes are, mating rates (how many females can be fertilized) and ejaculation rates (how many sperm should be released in the water) should also be considered potential causes of higher GSI 100 . Thus, a higher than expected GSI may be related not only to sperm competition, which would relate to strong sexual selection in males in protandrous species spawning in groups, but also to a physiological compensatory mechanism that allows males to fertilize the many eggs that large females produce.
To test these ideas further, field studies are needed to corroborate whether the protandrous sparids with the highest GSI values are those that spawn in large groups or aggregations rather than in random matings or pairs. Furthermore, laboratory experiments aimed at determining the actual fertilization capacity of small-sized protandrous sparids, specifically fertilization rates with different amounts of sperm, will be key to determine whether there is indeed a size threshold below which the GSI needs to increase in order to ensure fertilization of the eggs released by the larger females. This evidence would advance substantially our understanding of the relationship between sexual systems and mating systems in this diverse family in particular and in teleosts in general.
To conclude, this study provides the most updated analysis of the distribution and incidence of different sexual systems in the family Sparidae, an ideal model taxon in which to study the evolution of sexual systems and SAM predictions. We have found that both protandry and protogyny can evolve from gonochorism and back, but evolutionary transitions between the two types of sequential hermaphroditism are unlikely if ever to occur. We show that the predictions of the SAM incorporating mating system, spawning mode and sperm competition, hold well for protogynous and gonochoristic species. In contrast, protandrous species do not conform to theoretical expectations. The high GSI values of some protandrous sparids suggest that males compete to fertilize the eggs of females while others mate in pairs in the absence of male-male competition but still invest greatly in gonad tissue. We suggest that this is due to a compensatory mechanism that is intrinsic to protandry: boosting male investment in the gonads to ensure successful fertilization of the considerable number of eggs released by highly fecund females that are much larger than protandrous males.

Data availability
The datasets supporting this article have been uploaded as part of the supplementary material.