Hermaphroditism in fish: incidence, distribution and associations with abiotic environmental factors

The distribution of hermaphroditism in fishes has traditionally been mainly explained by its dependence on biotic factors. However, correlates with major abiotic factors have not been investigated on a quantitative basis and at a global scale. Here, we determined the incidence of hermaphroditism in fish at the family and species level, tested the hypothesis that evolutionary relationships account for the poor presence of hermaphroditism in freshwater species, and tested the association of sexual systems with latitude, habitat type and depth. Functional hermaphroditism is reported in 8 orders, 34 families and 370 species of fishes, all teleosts. Sequential hermaphroditism predominates over simultaneous hermaphroditism at a ratio ~ 5:1 and protogyny (female-to-male sex change) predominates ~ 6:1 over protandry (male-to-female). We found 12 hermaphroditic species that can live in freshwater. However, seven of these species are from four primarily marine families while there are only five species from two mostly freshwater families. Protogynous and bi-directional sex changers have a tighter association with reef-associated tropical and subtropical habitats when compared to protandrous species, which tend to be more plastic in terms of distribution requirements. Finally, simultaneous hermaphrodite species live both in the deep sea and shallow waters in similar proportions. This study can be the basis for further research in specific groups for different purposes, including ecological and evolutionary issues as well as conservation and management of exploited species. Understanding the environmental correlates can help to forecast changes in the distribution or phenology of hermaphrodites in a global change scenario.


Introduction
Hermaphroditism, defined as the presence of the male and female function (i.e., sperm and egg production, respectively) in the same individual, either sequentially or simultaneously, is present in the major taxonomic divisions of plants and is common in several Metazoans (Lewis 1942;Leonard 2013). Most invertebrate hermaphrodites are simultaneous and a rough estimate puts the number of hermaphrodite species * 65,000, or about 5-6% of animal species, a Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/ s11160-021- 09681-9. figure that increases to * 30% if insects are excluded (Jarne and Auld 2006).
The best model to explain when sequential hermaphroditism is favored is the size-advantage model, which predicts that sex change will occur when the reproductive success of one sex increases more rapidly with size (or age) than the reproductive success of the opposite sex (Ghiselin 1969;Warner 1975;Charnov 1982). On the other hand, the bestknown and supported theoretical framework to explain simultaneous hermaphroditism is the low-density model, which suggests that this sexual system is associated with the low probability of finding a partner (Tomlinson 1966;Ghiselin 1969;Charnov 1979;Charnov 1982).
Regarding the benefits of hermaphroditism, it has been accepted that sexual systems should be very sensitive to ecological parameters such as mate availability and reproductive assurance (Jarne and Auld 2006), mate search-efficiency (Eppley and Jesson 2008) and low encounter probability (Leonard 2010). However, the above models do not explain well the distribution of sexual systems in metazoans (Williams 1975). Instead, phylogeny reflects the large-scale distribution of sexual systems better than ecology in plants and metazoans (Leonard 2013(Leonard 2018. The difference between theory predictions and actual distribution is known as ' 'Williams' paradox'' (Leonard 199020132018. Fish are the only vertebrate group where hermaphroditism is present, with two main types. The first main type is sequential hermaphroditism, where individuals can: a) first mature as females, change sex and function as males for the rest of their lives (protogynous species), b) first mature as males, change sex and function as females for the rest of their lives (protandrous species), or c) can revert to the previous sex after change in low density situations (bi-directional sequential hermaphroditism; Sunobe and Nakazono 1993;Munday et al. 2010;Manabe et al. 2013). The second main type is simultaneous hermaphroditism, where both male and female gametes are produced by the same individual, either at the same time or within a very short period of time (Atz 1964;Yamamoto 1969).
Over the years, several studies have attempted to determine the incidence of hermaphroditism in fish. It was estimated to occur in 2% of 25,000 fish species (Pauly 2004a). Later, it was concluded to be distributed in 20 families in 9 orders (Mank et al. 2006). A plot has been presented displaying about half of the fish species (the ones available in The Tree of Sex) as hermaphroditic (Bachtrog et al. 2014). Discrepancies among estimations may be due to the fact that the diagnosis of hermaphroditism is not straightforward (Sadovy and Shapiro 1987). Criteria for diagnosing functional hermaphroditism in fish should be based on a detailed histological analysis of the gonads in various stages of sexual transition, simultaneous occurrence of mature gonadal tissues and observations of functional sex change (Sadovy and Shapiro 1987;Sadovy and Domeier 2005;Sadovy de Mitcheson and Liu 2008). Following these criteria, and after an extensive review of the literature, functional hermaphroditism was confirmed in 7 orders, 27 families (6% of all fish families) and 94 genera, but without estimation of actual number of species for each type of hermaphroditism (Sadovy de Mitcheson and Liu 2008). More recent studies have investigated the incidence of hermaphroditism in specific groups, e.g., gobies (Cole 2010;Manabe et al. 2013;Sunobe et al. 2017), serranids (Erisman and Hastings 2011), polynemids (Shihab et al. 2017;Butler et al. 2018) and sparids (Pla et al. 2020) to name some examples, but a global picture, based on current phylogenetic relationships, and with number of species for each major type of hermaphroditism taking into account only confirmed cases is missing. These last points, i.e., based on currently accepted phylogeny and only on confirmed species are very important. In this regard, recently Kuwamura et al. (2020) presented an annotated list of hermaphroditism in fish, with 485 species, but in several of these species functional hermaphroditism had not been confirmed. Further, the phylogenetic relationships were based on the classification of Nelson et al. (2016), which nowadays is considered outdated in many aspects. Nevertheless, the work of Kuwamura et al. (2020) is relevant because the analysis of the different forms of hermaphroditism with the mating system provides support for the size-advantage model. Based on the above, the first aim of this study was to provide an update of the different types of sexual systems in fish, focusing on the incidence of hermaphroditism at the species level, based on currently accepted phylogenetic relationships, and considering only confirmed species.
As explained above, at a small-scale, in fish the most important ecological drivers of hermaphroditism are social structure and mating system (Warner 1988;Munday et al. 2006;Godwin 2009). A comprehensive study in various fish lineages confirmed that social and mating systems, not evolutionary history, were the main drivers (Erisman et al. 2013;Kuwamura et al. 2020). However, one of the conclusions of the Sadovy de Mitcheson and Liu (2008) review was the need to also pay attention to the contribution of abiotic factors to better understand the distribution of hermaphroditism in fish because these, when compared to the influence of biotic factors, have received less attention (Sadovy de Mitcheson and Liu 2008). In fact, several studies show that fine-scale variation in the physical environment can influence community dynamics (Menke and Holway 2006) and this could affect reproduction. In the case of hermaphrodites, increased temperature due to global warming could affect the synthesis of sex steroids and influence the temporal dynamics of sex change (Pankhurst and Munday 2011). In this regard, in aquatic ecology four important abiotic factors are salinity, latitude (temperature), habitat type and depth (Margalef 1998).
Functional hermaphroditism has traditionally been considered either absent or very rare in freshwater habitats (Sadovy de Mitcheson and Liu 2008;Pavlov et al. 2009;Wootton and Smith 2014;Kuwamura et al. 2020). The number of hermaphroditic species in freshwater was calculated to represent only * 3% of the total number of hermaphroditic species, in clear contrast with the fact that freshwater fish species account for about half the number of all fish species (Kuwamura et al. 2020). The underlying cause is unknown but explanations for this phenomenon have relied on morphological aspects such as egg size. Thus, compared to marine fish, freshwater fish tend to spawn fewer and larger demersal eggs and larvae (Freedman and Noakes 2002;Pavlov et al. 2009). It has been argued that anatomical differences in the reproductive ducts may accentuate the physical differences in male and female anatomy, forming a sort of barrier for the development of hermaphroditism (Warner 1978;Sadovy de Mitcheson and Liu 2008). Also, parental care, common in freshwater fish, may reduce differences in reproductive success among sexes, as proposed by the size-advantage model, thus not favoring hermaphroditism (Sadovy de Mitcheson and Liu 2008). However, to the best of our knowledge, these hypotheses have not formally been tested, perhaps because more data needs to be collected. In any case, the low preponderance of hermaphroditism in freshwater fishes remains to be explained since it constitutes a gap in knowledge. Recently, the origin of vertebrates has been placed in shallow intertidal and subtidal environments ). However, trait reconstructions suggest that all extant marine Actinopterygians (ray-finned) fishes, the most speciesrich clade of marine vertebrates, containing 96% of fish species with roughly similar number of marine and freshwater species, were derived from a freshwater ancestor (Vega and Wiens 2012). Both marine and freshwater habitats are dominated by two clades, the Ostariophysi and the Percomorpha, although, interestingly, the most species-rich families of freshwater fishes (e.g., Cyprinidae, Cichlidae) do not contain hermaphrodite species (Sadovy de Mitcheson and Liu 2008). Knowledge on the association of the different sexual systems according to water salinity can provide opportunities to test hypotheses of possible adaptive significance of variations in sexual system (Sadovy de Mitcheson and Liu 2008). However, to the best of our knowledge, the distribution of the different sexual systems at species level in different water salinities has never been addressed on a numerical basis and on an evolutionary perspective. Thus, the second aim of this study was to test the hypothesis that evolutionary relationships account for the low presence of hermaphroditism in freshwater.
Latitude and temperature are closely linked. Temperature influences many aspects of aquatic ectotherms like fish such as survival, distribution, growth and reproductive output (Conover 1992;Heibo et al. 2005;Ruttenberg et al. 2005;Trip et al. 2008). It has been speculated that sex change at high latitudes is physiologically limited (Trip et al. 2011) because of the high associated energetic costs (Sadovy de Mitcheson and Liu 2008), despite that some protogynous species are found at high latitudes, e.g., Odax pullus (Trip et al. 2011) and Labrus bergylta, (Muncaster et al. 2010). The study of the distribution of the different sexual systems between each latitude region can help in future comparative analyses of ecological and geographical traits (Sadovy de Mitcheson and Liu 2008). In this regard, Pauly (2004b) showed a preponderance of hermaphrodites close to tropical areas. Thus, sequential hermaphroditism, especially protogyny, is associated with shallow tropical marine reef habitats (Warner 1984;Sadovy de Mitcheson and Liu 2008). This has been explained by ecological stability for the abundance of food (Barlow 1975) and predictable mating opportunities that these habitats offer, but whether this richness is a consequence of the higher abundance of species at low latitudes is not known. Also, and to the best of our knowledge, a quantitative assessment of the distribution of the different forms of hermaphroditism with different types of habitats has never been made. On the other hand, simultaneous hermaphroditism has been associated with the deep-sea (Ghiselin 1969;Baldwin and Johnson 1996;Davis and Fielitz 2010) in accordance with the low-density model (Ghiselin 1969). The classical example concerns several families of Aulopiformes (Davis and Fielitz 2010). However, simultaneous hermaphrodites are also present in shallow waters mostly represented by members of the family Serranidae but the relative distribution of simultaneous hermaphrodite species has not been addressed in a quantitative manner either. Thus, based on the above, the third aim of this study was to quantitatively test the association of the different sexual systems according to latitude, habitat and depth, taking into account phylogenetic relationships.

Creation of an ad hoc comprehensive database
We collected information on the sexual system of fishes from the primary literature, taking also into account previous general reviews (Sadovy de Mitcheson and Liu 2008), reviews on specific groups or families, e.g., Sparidae (Pla et al. 2020); Serranidae (Erisman and Hastings 2011). In total we searched 379 papers (references in the Supplementary information). We also obtained information on the following environmental variables, collected from FishBase (www.fishbase.com; Froese and Pauly 2018): salinity, latitude, habitat, and depth (see below for further details). Importantly, in this study we only consider as being functional hermaphrodites those species where evidence has been based on the criteria of Sadovy and Shapiro (1987) as in Sadovy de Mitcheson and Liu (2008), but not in those cases where evidence has been inferred and is not direct, or when functional hermaphroditism has not been fully confirmed. For taxonomic classification level higher than genus and species, as well as the most updated counts of valid species, we used Eschmeyer's Catalog of Fishes database (Fricke et al. 2019): http://researcharchive.calacademy. org/research/ichthyology/catalog/fishcatmain.asp. For evolutionary relationships, we used a comprehensive and updated phylogenetic tree of Actinopterygian fishes (Rabosky et al. 2018). The information provided here represents an expanded, curated and validated version of the dataset of Pla (2019) with only cases that could be confirmed.

Statistical analysis
In this study, the sexual system was always the dependent variable. Sexual system was treated as nominal categorical variable that, for a given species, could take one of the following values: gonochorism, protogyny, protandry, bi-directional sex change or simultaneous hermaphroditism. Because of the handful number of androdioecious species (species consisting of males and simultaneous hermaphrodites) in fish (Costa 2016;Petersen and Fischer 1986;Petersen 1990), this sexual system was not contemplated here. On the other hand, the independent variables or abiotic environmental factors considered were: salinity, latitude, habitat, and depth.
Regarding salinity, the presence or absence of a given species in each one of the three major types of water salinities (freshwater, brackish water and saltwater) in FishBase was coded as a binary variable with the values of 1 or 0, respectively. This allowed to account for the possibility that one species could be present in more than one type of water salinity. We then tested the association of the different sexual systems with the different water salinities types employing a multinomial generalized linear model (GLM) using the library nnet in R (Ripley et al. 2016).
Latitude was treated as a nominal categorical variable with the following categories: tropical, subtropical, temperate, and boreal. Each species was assigned to only one of the latitude categories, the most frequent one in case of being present in more than one. We tested the association of the different types of sexual systems with latitude categories with a multinomial GLM.
Habitat was also treated as a nominal categorical variable with the following possible categories: bathydemersal, bathypelagic, benthopelagic, demersal and reef-associated. The pelagic habitat was not considered because no hermaphrodites were present in this habitat. Each species was assigned to only one of the habitat categories, the most frequent one in case of being present in more than one. We tested the association of different types of hermaphroditism with these different marine habitats and to do so we used a multinomial GLM, as described above for salinity.
As for depth, we took into account depth data (all in meters) provided by FishBase: depth shallow, the shallowest depth at which a given species is found; depth common shallow, the mean upper depth at which a given species is commonly found; depth common deep, the mean bottom depth at which a given species is commonly found; and depth deep, the maximum depth reported for a given species. We tested the distribution of simultaneous hermaphroditism as a function of depth. To that end, a multinomial GLM on log depth deep was used. All data for the hermaphroditic species considered in this study can be found in the online Supplementary information. Importantly, because of the breakdown of the different types of hermaphroditism (numbers of species in the range of * 20 to * 245 depending on type; see the ''Summary tables'' in the database of the Supplementary information) combined with the low number of records for some of the abiotic factors, it was not possible, despite our attempts, to properly correct for phylogenetic signal, understood as the tendency for related species to resemble each other more for a given trait, in our case sexual system, than they resemble species drawn at random from the tree (Orzack and Sober 2001;Blomberg and Garland 2002). We thus were aware that this is a limitation that could potentially affect the results of the second and third objectives of this study, but not the first one, and therefore have taken every precaution we could in the interpretation of the results.
To check the sensitivity of the GLM results to the fact that certain categories are much more abundant than others (e.g., fishes are particularly abundant at tropical latitudes, in coral reefs and demersal habitats), we randomized the original data set by keeping the number of species per environmental factor category fixed and resampling the sexual system (50,000 times, without replacement). A Chi-square statistic was recalculated at each iteration to examine whether the contingency table observed could be within any of the randomly generated contingency tables by the resampling algorithm.

Incidence and distribution of the different forms of hermaphroditism in fishes
A summary of the incidence of the different types of hermaphroditism in fish and their classification in orders and families is shown in Table 1. A detailed list of all the confirmed functional hermaphrodite fish species considered in this study and their supporting primary literature can be found in the Supplementary information.

Distribution of hermaphroditism according to water salinity
We compiled information on salinity preference for the 370 hermaphroditic species (details in the Supplementary information). We did not include Squalius alburnoides, Henicorhynchus lobatus and Crenicara punctulatum because lack of enough evidence to consider that hermaphroditism is their usual form of Of the 244 protogynous species, 216 species (88.5%) were exclusively found in saltwater, 23 species (9.4%) both in saltwater and brackish water, none exclusively in brackish water, 1 species (0.4%), Lophogobius cyprinoides, in saltwater, brackish water and freshwater, 3 species (1.2%), Synbranchus marmoratus, Ophisternon bengalense and Monopterus albus, both in brackish water and freshwater, and 1 species (0.4%), Monopterus boueti, exclusively in freshwater. Regarding the 39 protandrous species, 20 species (74.4%) were exclusively found in saltwater, 12 species (30.7%) both in saltwater and brackish water, none exclusively in brackish water, 7 species (18%) in saltwater, brackish water and freshwater, none both in brackish water and freshwater, and none exclusively in freshwater. Regarding the 22 bi-directional species, 21 species (95.5%) were exclusively found in saltwater, 1 species (4.5%), Priolepis cincta, both in saltwater and brackish water, and none in any of the remaining combinations. Finally, of the 65 simultaneous hermaphrodite species, 62 species (95.4%) were exclusively found in saltwater, 3 species (4.6%), Chlorophthalmus agassizi, Diplectrum radiale and Serranus psittacinus, both in saltwater and brackish water, and none in any of the remaining combinations (Fig. 2a). Recall we are not considering here the androdioecious species (e.g., from the genus Kryptolebias). Thus, we found a total of twelve species of hermaphrodites present in freshwater (Supplementary information) of which seven species were from four mostly seawater families (Gobiidae, Centropomidae, Polynemidae and Sparidae) that can be found also in brackish water and fresh water, and five species from two families that mainly live in freshwater: the catadromous protandrous Lates calcarifer (Latidae), and the protogynous Monopterus boueti, M. albus, Ophisternon bengalense and Synbranchus marmoratus (Synbranchidae).
Since a given species could be present in more than one type of water salinity, Table 2 shows the number of species of each type of hermaphroditism according Bi-directional (7.2%) Fig. 1 Graphical representation of the relative abundance of the major sexual systems among teleost fishes. Sample sizes: Gonochorism, n = 5,268 (*) please note, this takes into account only the gonochoristic species from families that also contain hermaphrodite species, not the total number of extant species that are gonochores; Hermaphroditism, n = 370; Sequential hermaphroditism, n = 305; Simultaneous hermaphroditism, n = 65; Protogyny, n = 244; Protandry, n = 39; Bi-directional hermaphroditism, n = 22. Parthenogenesis (i.e., unisexuality, n = 2) and androdioecy (n = 3) were not included due to the very low number of species and difficulty of visualization

(c)
to water salinity. The ratio of species occurrence in saltwater, brackish water and fresh water was 48:5.4:1 for protogynous species and * 6:3:1 for protandrous species. In accordance, the results of the multinomial GLM showed hermaphrodites were underrepresented in freshwater (significantly so for protogynous and simultaneous hermaphrodites; P \ 0.001). Conversely, protandric hermaphrodites were relatively more frequent in brackish waters when compared to the other forms of hermaphroditism (P \ 0.001). Saltwater showed proportionally higher representation of bidirectional, protandric and simultaneous hermaphrodites (P \ 0.001; Table 3).
From an evolutionary point of view, basal osteichthyans do not contain hermaphroditic species. Hermaphroditism is first found in basal teleosts of the families Muraenidae (two simultaneous hermaphrodite species) and Gonostomatidae (one protandrous species), all marine (Fig. 2b). Most hermaphrodite species belong to families of derived teleosts, including the twelve hermaphrodite species found in freshwater (Fig. 2c).

Hermaphroditism and latitude
We obtained data on the most common latitude distribution for all the 370 hermaphroditic species. Eighteen species, all of them simultaneous hermaphrodites, had a circumglobal distribution. Of the 352 remaining species, 225 (* 64%) inhabited tropical regions. The number of hermaphrodite species was inversely related with latitude. Protogynous hermaphrodites predominated in all but the boreal regions, where simultaneous hermaphrodites were the most abundant type (Fig. 3a). Analysis of the latitudinal distribution as a function of the sexual system showed significant differences. Thus, GLM (tropical as reference) showed that sequential hermaphrodites, specifically protogynous species, were significantly more frequent in the tropical and subtropical latitudes (P \ 0.001). Conversely, the latitudinal distribution of simultaneous hermaphrodites was significantly different from that of the rest of hermaphrodites (P \ 0.001) while protandrous hermaphrodites were underrepresented (Table 4).  Fig. 2 Distribution of the different forms of hermaphroditism according to water salinity. a Numerical distribution of the different types of hermaphroditism in freshwater, brackish water and saltwater. Numbers in italics outside the circles represent the total number of species of a given sexual system found in freshwater (left, in green), brackish water (top, in purple) and saltwater (right, in blue). Numbers inside the colored circles represent the actual number of species in each type of salinity. To confirm that higher incidence of hermaphroditism is not due to the high number of species in the tropics, we randomized the original data set by keeping the number of species per latitudinal category fixed and resampling the type of reproduction (50,000 times, without replacement). The chi-square statistic was recalculated at each iteration. The expected distribution of the chi square statistic is approximately 0. As shown in Supplementary Fig. 1, the chi-square statistic estimated from the original data set does not fall within the values of the 50,000 randomized data sets (P \ 0.001).

The association of hermaphroditism with different marine habitats
Most hermaphrodites, particularly sequential hermaphrodites, are associated with reefs and, to a much lower degree, with demersal and benthopelagic habitats. Their presence in bathypelagic and bathydemersal habitats is much lower (Fig. 3b). In particular, protogynous hermaphrodites had a high probability to be found in reefs and low probability to be found in the benthopelagic (P \ 0.001). However, protandry had a similar presence in reefs and demersal (* 40% each; n = 16 and n = 15, respectively), with fewer species present in benthopelagic (16%, n = 6) and bathypelagic (3%, n = 1) habitats (no statistical difference).  Bi-directional hermaphrodites were mostly associated with reefs (86%, n = 19) and demersal (14%, n = 3) habitats, and significantly less abundant in the benthopelagic and bathydemersal. On the other hand, simultaneous hermaphroditism was found to be significantly reef-associated (n = 20, 32.2%; P \ 0.01) and the predominant form in bathypelagic (n = 22; 35.5%, P \ 0.001) and bathydemersal (n = 7; 11.3%; P \ 0.01) habitats (Table 5; Fig. 3b).

The association of hermaphroditism with depth
Of the 65 species with verified simultaneous hermaphroditism, only members of the order Aulopiformes inhabit the deep-sea. In contrast, we found that the maximum depth attained by simultaneous hermaphrodite species from other orders, i.e., Perciformes (26 species of the genus Serranus, Hypoplectrus and Diplectrum), Gobiiformes (two Priolepis sp.) and Anguilliformes (two Gymnothorax sp.), are not found beyond the first * 500 m (Fig. 4a). GLM results showed that simultaneous hermaphroditism is relatively more frequent with increasing depth (P \ 0.001). (Table 6). However, binning depth values in ranges of 200 m, showed that 25 hermaphroditic species live within the first 200 m and 31 hermaphroditic species (i.e., * 50% of the total) within the first 400 m of the water column (Fig. 4b).
The mean depth values of these species were 74.6 m for the 0-200 m range (n = 25 species); 337.6 m for the 201-400 m range (n = 6 species) and 2,392 m for [ 400 m (n = 31 species Fig. 4c). These results indicate that although there is an increase with depth, the number of simultaneous hermaphrodite species The Chi-square statistic refers to the result of the randomization test whose maximum depth is restricted to shallow waters (B 200 m) is similar to those whose maximum depth is [ 400 m.

Incidence and distribution of the different forms of hermaphroditism in fishes
Of the eight orders with hermaphroditism present, six (Anguilliformes, Aulopiformes, Perciformes, Scorpaeniformes, Stomiiformes, Synbranchiformes) coincide with six of the seven orders confirmed by Sadovy de Mitcheson and Liu (2008) and seven of the sixteen orders reported by Kuwamura et al. (2020). However,  for the Cyprinodontiformes, androdioecy rather than hermaphroditism, is the sexual system accepted and this distinction is here taken into account (see below). Two orders containing hermaphrodite species (Centrarchiformes, Gobiiformes) were previously considered as families by Sadovy de Mitcheson and Liu (2008) and two orders (Labriformes and Spariformes) by Kuwamura et al. (2020) are considered families in this study. Thus, there is no net increase in the number of orders in this study with respect to that of Sadovy de Mitcheson and Liu (2008). The presence of hermaphroditism in the order Siluriformes (Avise and Mank 2009) was not included by Sadovy de Mitcheson and Liu (2008) nor in the present study because of lack of direct evidence. The same applies to the orders Cichliformes, Clupeiformes, Cypriniformes, Moroniformes, Tetraodontiformes and Trachiniformes from the study by Kuwamura et al. 2002. Overall, the 370 hermaphrodite species identified so far represent 7% of the total of the 5,268 extant species in the teleost families containing at least one hermaphrodite species, but only 1.1% of the 33,786 extant species of fish if we include the rest of the families. These figures are in line with previous estimations. Thus, based on information on 25,000 species drawn from FishBase the number of hermaphrodites fish species was estimated * 2% (Pauly 2004a, b). In their review, Sadovy de Mitcheson and Liu (2008) did not provide a number of hermaphrodite species but estimated that hermaphroditism was present in at least 6% of teleost families. We find that the percent of hermaphrodite species varies greatly (\ 1% to 100%) among families but these figures must be taken with caution because of the low number of species for which there is information on sexual system in some families. Of note: regardless of whether one takes into account the five families with the higher number of hermaphrodite species or the families where hermaphroditism is present that contain the higher number (C 100) of species, the proportion of hermaphrodite species within a family is at most 24.1% (Sparidae) but usually much lower.
Recently, Kuwamura et al. (2020) reported a total of 462 hermaphrodite species in fish. Of these, 152 species differ from our database and at least a third (36%) show discrepancies with our list for various reasons. First, the evidence used by Kuwamura et al. (2020) was based for some species on the same old bibliographic references lacking direct or functional evidence that led Sadovy de Mitcheson and Liu (2008) and ourselves not to consider those species as hermaphrodites. Some of these species include Rhinomuraena quaesita or several Gonostomatidae species such as Cyclothone atraria, C. microdon, Gonostoma elongatum and Sigmops bathyphilum. Analysis of the gonadal histology or type of gametes present in the gonads did not provide sufficient evidence to confidently classify them as protandric hermaphrodites (Badcock and Merrett 1976;Sshen et al. 1979;Fisher 1983;Nemoto 1985 1987;Badcock 1986). Regarding the number of families, Kuwamura et al. (2020) report 41 families containing hermaphrodite species, of which seven are different from ours. When checking the sexual systems of the species (1-3 species per family at most) belonging to the families Balistidae, Clupeidae, Cobitidae, Malacanthidae, Moronidae, Poeciliidae and Terapontidae, it was found that the references provided were the same as those used in the Sadovy de Mitcheson and Liu (2008) study, which are proposed but not confirmed as hermaphrodites. Finally, the Creediidae family was not present in the Sadovy de Mitcheson and Liu (2008) paper but, again, there is no reliable evidence (Langston 2004;Shitamitsu and Sunobe 2017). Therefore, this increase in the number of families lacks proper support.
On the other hand, at least a dozen species taken as hermaphrodites in Kuwamura et al. (2020) are confirmed gonochoristic. These species include, for example, some serranids such as Paralabrax maculatofasciatus (Sadovy and Domeier 2005) or epinephelids such as Epinephelus striatus (Sadovy and Colin 1995); sparids such as Diplodus vulgaris (Gonçalves and Erzini 2000), Pagrus major (Matsuyama et al. 1988;Booth and Buxton 1997) or labrids such as Symphodus tinca (Warner and Lejeune 1985) and Bodianus eclancheri (Hoffman 1980). These assignations, coupled with the fact that the total number of species per family used by Kuwamura et al. (2020) was based on Nelson et al. (2016), and that both these numbers and the phylogenetic classification of the teleosts have changed, explains some of the discrepancies in their work. For example, they report 84.2% of the members of the family Sparidae being hermaphrodites. This number is a gross overestimation because is based on a total number of 38 sparid species, when their actual figure is * 150 species (Pla et al. 2020). Therefore, the study of Kuwamura et al. (2020) must be taken with these considerations in mind.
Focusing only on the two major types of hermaphroditism, sequential hermaphroditism (82.4%) predominates over simultaneous hermaphroditism (17.6%), is in sharp contrast with the situation in plants and invertebrates, where simultaneous hermaphroditism predominates (Leonard 2013). Even if we added the 158 species of simultaneous hermaphrodites of the order Aulopiformes mentioned elsewhere (Baldwin and Johnson 1996;Erisman et al. 2013) but for which we did not find explicit evidence or confirmation, sequential hermaphroditism would still be predominant in fish. The reason for the low numbers of simultaneous hermaphrodites in fish is not known but it could be related to the difficulty of handling sex hormones with opposite functions in the same gonad and with two functional roles at the same time (Warner 1978;Devlin and Nagahama 2002).
Androdioecy, a mixed sexual system where individuals in a population can be simultaneous hermaphrodites or males, but never females, has so far only been described in three species of the genus Kryptolebias (K. ocellatus, K. marmoratus and K. hermaphroditus) of the Rivulidae family (Costa et al. 2010;Costa 2016). In addition, there are two species of Serranus (S. baldwini and S. psittacinus, the latter often referred to as S. fasciatus in previous studies) that have also been considered androdioecious (Erisman et al. 2013). However, unlike Kryptolebias, that in addition to simultaneous hermaphrodites can also contain primary males under certain environmental conditions (Harrington 1971), populations of these serranid species consist of simultaneous hermaphrodites, where a dominant individual will eventually end up losing its female function and behave as a male in a harem of hermaphrodites (Petersen and Fischer 1986;Petersen 1990). Thus, although populations of the species of the genus Kryptolebias and Serranus referred to above can end up consisting of males and simultaneous hermaphrodites, in this study, S. baldwini and S. psittacinus were considered simultaneous hermaphrodites, not androdioecious, in line with had been done before by Petersen and Fischer (1986) and Petersen (1990).
Bi-directional sex change is thought to be related to the small size, short life cycle and elevated predation risk (Patzner et al. 2011). Bi-directional sex change has been found naturally in the wild occurring mainly in species of the family Gobiidae (21 species), although instances have been also reported in other families such as Cirrhitidae and Pseudochromidae (Kadota et al. 2012;Kuwamura et al. 2015). However, sometimes evidence has been obtained mainly in laboratory settings after hormonal or social induction as reported, for example, in several species of the family Pomacanthidae (Hioki and Suzuki 1996), Epinephelidae (Tanaka 1990), Labridae (Kuwamura et al. 2002(Kuwamura et al. 2007Ohta et al. 2003) and Pseudochromidae (Wittenrich and Munday 2005;Kuwamura et al. 2015). Many of these bi-directional sex changers (e.g., Cirrhitichthys falco) are mainly protogynous in the wild. Thus, in the present study, only the most common sexual system has been considered in these cases, i.e., protogyny. The species thus concerned have been highlighted in the Supplementary information to denote that they can also be bi-directional under certain circumstances. Adding these species, the incidence of bi-directional sequential hermaphroditism would rise from 7.2% to 12.7%, a figure similar to that of the number of protandrous species.

Hermaphroditism in freshwater fishes
Hermaphroditism in freshwater fishes has been traditionally considered to be absent (Wootton and Smith 2014) or rare (Sadovy de Mitcheson and Liu 2008;Pavlov et al. 2009). In our study we attempt to provide the first explicit assessment of the number of hermaphrodite fish species present in freshwater. There are 12 species of hermaphrodites in freshwater (the androdioecious species mentioned above were excluded). These 12 species (3.2% of all hermaphroditic species) belong to six families (Centropomidae, Gobiidae, Polynemidae, Sparidae, Latidae, and Synbranchidae) representing * 18% of all families containing hermaphrodites. This might suggest diversity in hermaphroditic freshwater fishes. However, seven of these species belong to families that are typically marine or marine-brackish water (the first four families named above). These species are euryhaline and can live also in freshwater. Thus, there is only five protogynous species that belong to two typically freshwater families: one is the catadromous protandrous Lates calcarifer (Moore 1979;Guiguen et al. 1994), family Latidae, and four protogynous species of the family Synbranchidae: Monopterus boueti, M. albus, Ophisternon bengalense and Synbranchus marmoratus (Liem 1963;Liem 1968). However, only M. boueti is exclusively freshwater. There are no simultaneous hermaphrodites that exclusively live in freshwater, again if members of the Rivulidae family with androdioecy are not considered. These low numbers do not allow proper correction for phylogenetic signal but given the variety of the families to which the species of hermaphrodites found in freshwater belong it is unlikely that results would be affected.
It has been argued that morphological (e.g., the large eggs typically produced by freshwater fishes, which could accentuate sex-related differences in the anatomy of the reproductive system) and behavioral aspects (e.g., nest building and parental care) could represent a barrier to the development of hermaphrodites (Warner 1978;Sadovy de Mitcheson and Liu 2008). This would be so despite the fact that teleosts have higher potential, compared to other vertebrates, for developing hermaphroditism (Adolfi et al. 2018) due to their gonadal plasticity (Oldfield 2005;Piferrer 2021). In this regard, Lates calcarifer, produces eggs of 0.7-0.8 mm diameter, i.e., like those of many pelagic marine fishes and in Monopterus boueti, males build a nest and guard the eggs, which does not seem to represent any barrier for the presence of hermaphroditism. Thus, we propose that evolutionary relationships primarily account for the presence or absence of hermaphroditism in freshwater fishes although the influence of morphological or behavioral aspects as mentioned above cannot be ruled out. ''Primary'' freshwater fish, according to accepted evolutionary classifications (Witten and Huysseune 2009), and exemplified by families such as Cyprinidae and Salmonidae, do not contain hermaphrodite species. All the six families that contain hermaphroditic species that are found in freshwater belong to what can be called the more derived teleosts, where sequential hermaphroditism evolved. Together, these results suggest that hermaphroditism in freshwater species is ultimately influenced by evolutionary history. Comparative studies on the anatomical and physiological aspects that negatively might influence the development of hermaphroditism in freshwater are needed and should examine the families, particularly those comprising species living in different types of salinities, and all the species described in the present study, given their low number. These studies could shed more light on why hermaphroditism is not common in freshwater when compared to marine habitats.
The association of the different sexual systems according to latitude, habitat and depth Our results confirm the preponderance of sequential hermaphroditism, mainly protogynous hermaphrodites, in rich and complex habitats such as coral reefs. Our permutation tests show that this is not due to the higher number of species present in some latitudes or habitats. Fifteen out of the sixteen families containing protogynous hermaphrodites have species that live in reefs. The exception is the Synbranchidae, one of the two freshwater families with hermaphrodite species. The 45 protogynous species (* 20% of all protogynous) not associated with reefs live in demersal and benthopelagic habitats. Although the social system of most of these species is not as well-known as those associated with reefs, the haremic system is present in some (Cole 1983;Kline et al. 2011), suggesting that adaptation to these habitats is possible as long as they maintain their social system, the main driver for protogyny (Warner 1988;Munday et al. 2006).
Regarding latitudinal distribution, protogyny is mainly associated with tropical reefs because of habitat stability and potential for polygamy of adults, among other possible environmental factors (DeMartini and Sikkel 2006) that allows the territoriality and sedentary life required by this system (Warner 1984). However, around 5% protogynous species have adapted to temperate environments, which are regarded as more unpredictable environments when compared to reefs (DeMartini and Sikkel 2006). But even so, few protogynous species such as Centropyge interruptus and Labrus bergylta inhabit temperate environments and exhibit haremic system (Moyer and Nakazono 1978;Muncaster et al. 2013), again suggesting that adaptation to these latitudes is also possible as long as they maintain their social system. We detected a clear association of protogyny with relatively shallow waters in reef-associated habitat and their absence in bathypelagic and bathydemersal habitats (Fig. 3b).
The distribution of protandrous species, although also most abundant in tropical latitudes (64% of species) and in reef habitats, is not as skewed in favor of those habitats as that of protogynous species since they have a similar proportion between reefs and demersal habitats (* 40% of species in each) with a few species in benthopelagic habitats. In addition, protandrous species have a comparatively more association with brackish water environments as discussed above.
Sequential bi-directional species show, broadly speaking, the same distribution as protogynous species, perhaps even more associated with relatively shallow waters (to * 60 m depth) and tropical reefs. This tight association may be due to their general small size, short life cycle and elevated risk of predation (Patzner et al. 2011), which does not favor dispersion. In summary, among sequential hermaphrodites, protogynous and bi-directional sex changers have a somehow tighter association with specific habitats when compared to protandrous species, which tend to be more plastic in terms of habitat/distribution requirements. This difference is probably due to the social and mating system characteristic of protogynous species, which normally require complex habitats such as coral reefs that favor the maintenance also of social structure such as harems where a dominant male secures access to a group of females.
Simultaneous hermaphroditism evolved independently in at least four orders (Aulopiformes, Anguilliformes, Gobiiformes, and Perciformes). Thus, it first appeared in basal teleosts, making it the oldest lineage of hermaphrodites (Davis and Fielitz 2010). Simultaneous hermaphroditism, exemplified by members of several families of Aulopiformes, has been associated with the deep-sea (Ghiselin 1969;Baldwin and Johnson 1996;Davis and Fielitz 2010) in accordance with the low-density model (Ghiselin 1969), but simultaneous hermaphrodites of the families Serranidae and Gobiidae inhabit shallower waters. We, therefore, tested the distribution of simultaneous hermaphrodites with depth. We did not consider sequential hermaphrodites because the most abundant type of them, protogyny, are reef-associated, essentially absent or totally absent form bathypelagic or bathydemersal habitats, as our data shows (Fig. 3b). The absence in certain depth categories made us consider that there was no need to test depth with sequential hermaphroditism.
Our analysis showed a similar number of simultaneous hermaphrodites in shallow (B 200 m) waters and depths [ 400 m (ratio 1:1.2). Even assuming that all of the 126 species of Aulopiformes mentioned in Baldwin and Johnson (1996) for which we could not confirm their sexual system lived exclusively in deep waters, the distribution of simultaneous hermaphroditism in fish is not only associated with the deep sea. In fish, simultaneous hermaphroditism has a significant probability of being found in both reefs and bathypelagic habitats, in the latter being more likely than other sexual systems. Thus, while it is true that simultaneous hermaphroditism is the only type of hermaphroditism present in the deep sea, the opposite is not true, i.e., simultaneous hermaphrodites live both in the deep sea and in shallow waters in roughly similar proportions based on the number of confirmed species at present.
In summary, we have provided the first attempt to quantify, utilizing only species which are confirmed hermaphrodites, the incidence and distribution of the different forms of hermaphroditism in fish focusing on the influence on abiotic factors, an area underexplored when compared to the influence of the social environment mating behavior and spawning mode. Of note, our data lacks control for phylogenetic signal, which could introduce some bias in the analysis of the role of environmental correlates we show. We could not correct for phylogeny due to the limited number of species present in some forms of hermaphroditism. However, close inspection of the data where in principle correction for phylogeny could potentially be useful, e.g., data presented in Fig. 4a indicates an even distribution of species among different depths, suggesting that correcting for phylogenetic signal would not significantly alter the results. The randomization tests showed that the associations we found between the different sexual systems and latitude, habitat and depth are indeed real and not influenced by the higher abundance of species in some latitudes, habitats or depths. This however, does not diminish the importance of biotic factors in explaining the distribution of hermaphroditism. In fact, both biotic and abiotic factors must act together for making the presence of hermaphroditism possible.
To conclude, the major insights that our study bring to the knowledge of hermaphroditism in fishes can be summarized as follows: 1) we have attempted, for the first time, to quantify at the species level the phylogenetic distribution of hermaphroditism in fish, considering only species in which it has been confirmed, 2) we show the higher presence of the most abundant form of hermaphroditism in some latitudes and habitats is not merely due to the fact that these habitats are the most rich in species abundance, 3) the scarce presence of hermaphroditism in freshwater fishes seems better explained by evolutionary relationships rather than by morphological or anatomical constraints, although these cannot be excluded, and 4) simultaneous hermaphroditism appears as common in shallow waters and complex environments as it is in deep-water ecosystems with very low encounter probabilities. This study can serve as the basis for further research in specific groups for different purposes, including ecological and evolutionary studies. Theories that account for environmental heterogeneity predicts that biotic interactions determine distribution limits under benign climatic conditions while abiotic limitations set distribution limits under harsh climatic conditions (Gaston 2009;Louthan et al. 2015). Therefore, our study can also be useful to better understand possible distribution shifts in response to global change, mainly due to increases in temperatures that could affect the reproductive dynamics of hermaphrodites and dispersion of larvae, hence risking mismatches with other environmental factors. This is particularly relevant taking into account that many hermaphrodite fish species sustain important fisheries.
Author contributions FP conceived the study; FP and SP designed the study; SP collected the data; SP and FM analyzed the data; SP and FP wrote the manuscript. All authors revised and approved the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Data availability statement All the data used in the present study is contained in the online Supplementary information.